faf6fd8c81932dfe58580ca23abe2af059080a4f
1 /*
2  * Copyright (C) 2007 Michael Brown <mbrown@fensystems.co.uk>.
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License as
6  * published by the Free Software Foundation; either version 2 of the
7  * License, or any later version.
8  *
9  * This program is distributed in the hope that it will be useful, but
10  * WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
12  * General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software
16  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
17  */
19 /** @file
20  *
21  * 64-bit division
22  *
23  * The x86 CPU (386 upwards) has a divl instruction which will perform
24  * unsigned division of a 64-bit dividend by a 32-bit divisor.  If the
25  * resulting quotient does not fit in 32 bits, then a CPU exception
26  * will occur.
27  *
28  * Unsigned integer division is expressed as solving
29  *
30  *   x = d.q + r                        0 <= q, 0 <= r < d
31  *
32  * given the dividend (x) and divisor (d), to find the quotient (q)
33  * and remainder (r).
34  *
35  * The x86 divl instruction will solve
36  *
37  *   x = d.q + r                        0 <= q, 0 <= r < d
38  *
39  * given x in the range 0 <= x < 2^64 and 1 <= d < 2^32, and causing a
40  * hardware exception if the resulting q >= 2^32.
41  *
42  * We can therefore use divl only if we can prove that the conditions
43  *
44  *   0 <= x < 2^64
45  *   1 <= d < 2^32
46  *        q < 2^32
47  *
48  * are satisfied.
49  *
50  *
51  * Case 1 : 1 <= d < 2^32
52  * ======================
53  *
54  * We express x as
55  *
56  *   x = xh.2^32 + xl                   0 <= xh < 2^32, 0 <= xl < 2^32  (1)
57  *
58  * i.e. split x into low and high dwords.  We then solve
59  *
60  *   xh = d.qh + r'                     0 <= qh, 0 <= r' < d            (2)
61  *
62  * which we can do using a divl instruction since
63  *
64  *   0 <= xh < 2^64                     since 0 <= xh < 2^32 from (1)   (3)
65  *
66  * and
67  *
68  *   1 <= d < 2^32                      by definition of this Case      (4)
69  *
70  * and
71  *
72  *   d.qh = xh - r'                     from (2)
73  *   d.qh <= xh                         since r' >= 0 from (2)
74  *   qh <= xh                           since d >= 1 from (2)
75  *   qh < 2^32                          since xh < 2^32 from (1)        (5)
76  *
77  * Having obtained qh and r', we then solve
78  *
79  *   ( r'.2^32 + xl ) = d.ql + r        0 <= ql, 0 <= r < d             (6)
80  *
81  * which we can do using another divl instruction since
82  *
83  *   xl <= 2^32 - 1                     from (1), so
84  *   r'.2^32 + xl <= ( r' + 1 ).2^32 - 1
85  *   r'.2^32 + xl <= d.2^32 - 1         since r' < d from (2)
86  *   r'.2^32 + xl < d.2^32                                              (7)
87  *   r'.2^32 + xl < 2^64                since d < 2^32 from (4)         (8)
88  *
89  * and
90  *
91  *   1 <= d < 2^32                      by definition of this Case      (9)
92  *
93  * and
94  *
95  *   d.ql = ( r'.2^32 + xl ) - r        from (6)
96  *   d.ql <= r'.2^32 + xl               since r >= 0 from (6)
97  *   d.ql < d.2^32                      from (7)
98  *   ql < 2^32                          since d >= 1 from (2)           (10)
99  *
100  * This then gives us
101  *
102  *   x = xh.2^32 + xl                   from (1)
103  *   x = ( d.qh + r' ).2^32 + xl        from (2)
104  *   x = d.qh.2^32 + ( r'.2^32 + xl )
105  *   x = d.qh.2^32 + d.ql + r           from (3)
106  *   x = d.( qh.2^32 + ql ) + r                                         (11)
107  *
108  * Letting
109  *
110  *   q = qh.2^32 + ql                                                   (12)
111  *
112  * gives
113  *
114  *   x = d.q + r                        from (11) and (12)
115  *
116  * which is the solution.
117  *
118  *
119  * This therefore gives us a two-step algorithm:
120  *
121  *   xh = d.qh + r'                     0 <= qh, 0 <= r' < d            (2)
122  *   ( r'.2^32 + xl ) = d.ql + r        0 <= ql, 0 <= r < d             (6)
123  *
124  * which translates to
125  *
126  *   %edx:%eax = 0:xh
127  *   divl d
128  *   qh = %eax
129  *   r' = %edx
130  *
131  *   %edx:%eax = r':xl
132  *   divl d
133  *   ql = %eax
134  *   r = %edx
135  *
136  * Note that if
137  *
138  *   xh < d
139  *
140  * (which is a fast dword comparison) then the first divl instruction
141  * can be omitted, since the answer will be
142  *
143  *   qh = 0
144  *   r = xh
145  *
146  *
147  * Case 2 : 2^32 <= d < 2^64
148  * =========================
149  *
150  * We first express d as
151  *
152  *   d = dh.2^k + dl                    2^31 <= dh < 2^32,
153  *                                      0 <= dl < 2^k, 1 <= k <= 32     (1)
154  *
155  * i.e. find the highest bit set in d, subtract 32, and split d into
156  * dh and dl at that point.
157  *
158  * We then express x as
159  *
160  *   x = xh.2^k + xl                    0 <= xl < 2^k                   (2)
161  *
162  * giving
163  *
164  *   xh.2^k = x - xl                    from (2)
165  *   xh.2^k <= x                        since xl >= 0 from (1)
166  *   xh.2^k < 2^64                      since xh < 2^64 from (1)
167  *   xh < 2^(64-k)                                                      (3)
168  *
169  * We then solve the division
170  *
171  *   xh = dh.q' + r'                            0 <= r' < dh            (4)
172  *
173  * which we can do using a divl instruction since
174  *
175  *   0 <= xh < 2^64                     since x < 2^64 and xh < x
176  *
177  * and
178  *
179  *   1 <= dh < 2^32                     from (1)
180  *
181  * and
182  *
183  *   dh.q' = xh - r'                    from (4)
184  *   dh.q' <= xh                        since r' >= 0 from (4)
185  *   dh.q' < 2^(64-k)                   from (3)                        (5)
186  *   q'.2^31 <= dh.q'                   since dh >= 2^31 from (1)       (6)
187  *   q'.2^31 < 2^(64-k)                 from (5) and (6)
188  *   q' < 2^(33-k)
189  *   q' < 2^32                          since k >= 1 from (1)           (7)
190  *
191  * This gives us
192  *
193  *   xh.2^k = dh.q'.2^k + r'.2^k        from (4)
194  *   x - xl = ( d - dl ).q' + r'.2^k    from (1) and (2)
195  *   x = d.q' + ( r'.2^k + xl ) - dl.q'                                 (8)
196  *
197  * Now
198  *
199  *  r'.2^k + xl < r'.2^k + 2^k          since xl < 2^k from (2)
200  *  r'.2^k + xl < ( r' + 1 ).2^k
201  *  r'.2^k + xl < dh.2^k                since r' < dh from (4)
202  *  r'.2^k + xl < ( d - dl )            from (1)                        (9)
203  *
204  *
205  * (missing)
206  *
207  *
208  * This gives us two cases to consider:
209  *
210  * case (a):
211  *
212  *   dl.q' <= ( r'.2^k + xl )                                           (15a)
213  *
214  * in which case
215  *
216  *   x = d.q' + ( r'.2^k + xl - dl.q' )
217  *
218  * is a direct solution to the division, since
219  *
220  *   r'.2^k + xl < d                    from (9)
221  *   ( r'.2^k + xl - dl.q' ) < d        since dl >= 0 and q' >= 0
222  *
223  * and
224  *
225  *   0 <= ( r'.2^k + xl - dl.q' )       from (15a)
226  *
227  * case (b):
228  *
229  *   dl.q' > ( r'.2^k + xl )                                            (15b)
230  *
231  * Express
232  *
233  *  x = d.(q'-1) + ( r'.2^k + xl ) + ( d - dl.q' )
234  *
235  *
236  * (missing)
237  *
238  *
239  * special case: k = 32 cannot be handled with shifts
240  *
241  * (missing)
242  *
243  */
245 #include <stdint.h>
246 #include <assert.h>
248 typedef uint64_t UDItype;
250 struct uint64 {
251         uint32_t l;
252         uint32_t h;
253 };
255 static inline void udivmod64_lo ( const struct uint64 *x,
256                                   const struct uint64 *d,
257                                   struct uint64 *q,
258                                   struct uint64 *r ) {
259         uint32_t r_dash;
261         q->h = 0;
262         r->h = 0;
263         r_dash = x->h;
265         if ( x->h >= d->l ) {
266                 __asm__ ( "divl %2"
267                           : "=&a" ( q->h ), "=&d" ( r_dash )
268                           : "g" ( d->l ), "0" ( x->h ), "1" ( 0 ) );
269         }
271         __asm__ ( "divl %2"
272                   : "=&a" ( q->l ), "=&d" ( r->l )
273                   : "g" ( d->l ), "0" ( x->l ), "1" ( r_dash ) );
274 }
276 static void udivmod64 ( const struct uint64 *x,
277                         const struct uint64 *d,
278                         struct uint64 *q,
279                         struct uint64 *r ) {
281         if ( d->h == 0 ) {
282                 udivmod64_lo ( x, d, q, r );
283         } else {
284                 assert ( 0 );
285                 while ( 1 ) {};
286         }
287 }
289 /**
290  * 64-bit division with remainder
291  *
292  * @v x                 Dividend
293  * @v d                 Divisor
294  * @ret r               Remainder
295  * @ret q               Quotient
296  */
297 UDItype __udivmoddi4 ( UDItype x, UDItype d, UDItype *r ) {
298         UDItype q;
299         UDItype *_x = &x;
300         UDItype *_d = &d;
301         UDItype *_q = &q;
302         UDItype *_r = r;
304         udivmod64 ( ( struct uint64 * ) _x, ( struct uint64 * ) _d,
305                     ( struct uint64 * ) _q, ( struct uint64 * ) _r );
307         assert ( ( x == ( ( d * q ) + (*r) ) ) );
308         assert ( (*r) < d );
310         return q;
311 }
313 /**
314  * 64-bit division
315  *
316  * @v x                 Dividend
317  * @v d                 Divisor
318  * @ret q               Quotient
319  */
320 UDItype __udivdi3 ( UDItype x, UDItype d ) {
321         UDItype r;
322         return __udivmoddi4 ( x, d, &r );
323 }
325 /**
326  * 64-bit modulus
327  *
328  * @v x                 Dividend
329  * @v d                 Divisor
330  * @ret q               Quotient
331  */
332 UDItype __umoddi3 ( UDItype x, UDItype d ) {
333         UDItype r;
334         __udivmoddi4 ( x, d, &r );
335         return r;
336 }