1 /*-
2 * Copyright (c) 1992, 1993
3 * The Regents of the University of California. All rights reserved.
4 *
5 * This software was developed by the Computer Systems Engineering group
6 * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
7 * contributed to Berkeley.
8 *
9 * Redistribution and use in source and binary forms, with or without
10 * modification, are permitted provided that the following conditions
11 * are met:
12 * 1. Redistributions of source code must retain the above copyright
13 * notice, this list of conditions and the following disclaimer.
14 * 2. Redistributions in binary form must reproduce the above copyright
15 * notice, this list of conditions and the following disclaimer in the
16 * documentation and/or other materials provided with the distribution.
17 * 4. Neither the name of the University nor the names of its contributors
18 * may be used to endorse or promote products derived from this software
19 * without specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
22 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
25 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
27 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
30 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
31 * SUCH DAMAGE.
32 */
33
34 #ifndef _LIBKERN_QUAD_H_
35 #define _LIBKERN_QUAD_H_
36
37 /*
38 * Quad arithmetic.
39 *
40 * This library makes the following assumptions:
41 *
42 * - The type long long (aka quad_t) exists.
43 *
44 * - A quad variable is exactly twice as long as `long'.
45 *
46 * - The machine's arithmetic is two's complement.
47 *
48 * This library can provide 128-bit arithmetic on a machine with 128-bit
49 * quads and 64-bit longs, for instance, or 96-bit arithmetic on machines
50 * with 48-bit longs.
51 */
52 /*
53 #include <sys/cdefs.h>
54 #include <sys/types.h>
55 #include <sys/limits.h>
56 #include <sys/syslimits.h>
57 */
58
59 #include <limits.h>
60 typedef long long quad_t;
61 typedef unsigned long long u_quad_t;
62 typedef unsigned long u_long;
63 #define CHAR_BIT __CHAR_BIT__
64
65 /*
66 * Define the order of 32-bit words in 64-bit words.
67 * For little endian only.
68 */
69 #define _QUAD_HIGHWORD 1
70 #define _QUAD_LOWWORD 0
71
72 /*
73 * Depending on the desired operation, we view a `long long' (aka quad_t) in
74 * one or more of the following formats.
75 */
76 union uu {
77 quad_t q; /* as a (signed) quad */
78 quad_t uq; /* as an unsigned quad */
79 long sl[2]; /* as two signed longs */
80 u_long ul[2]; /* as two unsigned longs */
81 };
82
83 /*
84 * Define high and low longwords.
85 */
86 #define H _QUAD_HIGHWORD
87 #define L _QUAD_LOWWORD
88
89 /*
90 * Total number of bits in a quad_t and in the pieces that make it up.
91 * These are used for shifting, and also below for halfword extraction
92 * and assembly.
93 */
94 #define QUAD_BITS (sizeof(quad_t) * CHAR_BIT)
95 #define LONG_BITS (sizeof(long) * CHAR_BIT)
96 #define HALF_BITS (sizeof(long) * CHAR_BIT / 2)
97
98 /*
99 * Extract high and low shortwords from longword, and move low shortword of
100 * longword to upper half of long, i.e., produce the upper longword of
101 * ((quad_t)(x) << (number_of_bits_in_long/2)). (`x' must actually be u_long.)
102 *
103 * These are used in the multiply code, to split a longword into upper
104 * and lower halves, and to reassemble a product as a quad_t, shifted left
105 * (sizeof(long)*CHAR_BIT/2).
106 */
107 #define HHALF(x) ((x) >> HALF_BITS)
108 #define LHALF(x) ((x) & ((1 << HALF_BITS) - 1))
109 #define LHUP(x) ((x) << HALF_BITS)
110
111 typedef unsigned int qshift_t;
112
113 quad_t __ashldi3(quad_t, qshift_t);
114 quad_t __ashrdi3(quad_t, qshift_t);
115 int __cmpdi2(quad_t a, quad_t b);
116 quad_t __divdi3(quad_t a, quad_t b);
117 quad_t __lshrdi3(quad_t, qshift_t);
118 quad_t __moddi3(quad_t a, quad_t b);
119 u_quad_t __qdivrem(u_quad_t u, u_quad_t v, u_quad_t *rem);
120 u_quad_t __udivdi3(u_quad_t a, u_quad_t b);
121 u_quad_t __umoddi3(u_quad_t a, u_quad_t b);
122 int __ucmpdi2(u_quad_t a, u_quad_t b);
123 quad_t __divmoddi4(quad_t a, quad_t b, quad_t *rem);
124
125 #endif /* !_LIBKERN_QUAD_H_ */
126
127 #if defined (_X86_) && !defined (__x86_64__)
128 /*
129 * Shift a (signed) quad value left (arithmetic shift left).
130 * This is the same as logical shift left!
131 */
132 quad_t
__ashldi3(a,shift)133 __ashldi3(a, shift)
134 quad_t a;
135 qshift_t shift;
136 {
137 union uu aa;
138
139 aa.q = a;
140 if (shift >= LONG_BITS) {
141 aa.ul[H] = shift >= QUAD_BITS ? 0 :
142 aa.ul[L] << (shift - LONG_BITS);
143 aa.ul[L] = 0;
144 } else if (shift > 0) {
145 aa.ul[H] = (aa.ul[H] << shift) |
146 (aa.ul[L] >> (LONG_BITS - shift));
147 aa.ul[L] <<= shift;
148 }
149 return (aa.q);
150 }
151
152 /*
153 * Shift a (signed) quad value right (arithmetic shift right).
154 */
155 quad_t
__ashrdi3(a,shift)156 __ashrdi3(a, shift)
157 quad_t a;
158 qshift_t shift;
159 {
160 union uu aa;
161
162 aa.q = a;
163 if (shift >= LONG_BITS) {
164 long s;
165
166 /*
167 * Smear bits rightward using the machine's right-shift
168 * method, whether that is sign extension or zero fill,
169 * to get the `sign word' s. Note that shifting by
170 * LONG_BITS is undefined, so we shift (LONG_BITS-1),
171 * then 1 more, to get our answer.
172 */
173 s = (aa.sl[H] >> (LONG_BITS - 1)) >> 1;
174 aa.ul[L] = shift >= QUAD_BITS ? s :
175 aa.sl[H] >> (shift - LONG_BITS);
176 aa.ul[H] = s;
177 } else if (shift > 0) {
178 aa.ul[L] = (aa.ul[L] >> shift) |
179 (aa.ul[H] << (LONG_BITS - shift));
180 aa.sl[H] >>= shift;
181 }
182 return (aa.q);
183 }
184
185 /*
186 * Return 0, 1, or 2 as a <, =, > b respectively.
187 * Both a and b are considered signed---which means only the high word is
188 * signed.
189 */
190 int
__cmpdi2(a,b)191 __cmpdi2(a, b)
192 quad_t a, b;
193 {
194 union uu aa, bb;
195
196 aa.q = a;
197 bb.q = b;
198 return (aa.sl[H] < bb.sl[H] ? 0 : aa.sl[H] > bb.sl[H] ? 2 :
199 aa.ul[L] < bb.ul[L] ? 0 : aa.ul[L] > bb.ul[L] ? 2 : 1);
200 }
201
202 /*
203 * Divide two signed quads.
204 * ??? if -1/2 should produce -1 on this machine, this code is wrong
205 */
206 quad_t
__divdi3(a,b)207 __divdi3(a, b)
208 quad_t a, b;
209 {
210 u_quad_t ua, ub, uq;
211 int neg;
212
213 if (a < 0)
214 ua = -(u_quad_t)a, neg = 1;
215 else
216 ua = a, neg = 0;
217 if (b < 0)
218 ub = -(u_quad_t)b, neg ^= 1;
219 else
220 ub = b;
221 uq = __qdivrem(ua, ub, (u_quad_t *)0);
222 return (neg ? -uq : uq);
223 }
224
225 /*
226 * Shift an (unsigned) quad value right (logical shift right).
227 */
228 quad_t
__lshrdi3(a,shift)229 __lshrdi3(a, shift)
230 quad_t a;
231 qshift_t shift;
232 {
233 union uu aa;
234
235 aa.q = a;
236 if (shift >= LONG_BITS) {
237 aa.ul[L] = shift >= QUAD_BITS ? 0 :
238 aa.ul[H] >> (shift - LONG_BITS);
239 aa.ul[H] = 0;
240 } else if (shift > 0) {
241 aa.ul[L] = (aa.ul[L] >> shift) |
242 (aa.ul[H] << (LONG_BITS - shift));
243 aa.ul[H] >>= shift;
244 }
245 return (aa.q);
246 }
247
248 /*
249 * Return remainder after dividing two signed quads.
250 *
251 * XXX
252 * If -1/2 should produce -1 on this machine, this code is wrong.
253 */
254 quad_t
__moddi3(a,b)255 __moddi3(a, b)
256 quad_t a, b;
257 {
258 u_quad_t ua, ub, ur;
259 int neg;
260
261 if (a < 0)
262 ua = -(u_quad_t)a, neg = 1;
263 else
264 ua = a, neg = 0;
265 if (b < 0)
266 ub = -(u_quad_t)b;
267 else
268 ub = b;
269 (void)__qdivrem(ua, ub, &ur);
270 return (neg ? -ur : ur);
271 }
272
273
274 /*
275 * Multiprecision divide. This algorithm is from Knuth vol. 2 (2nd ed),
276 * section 4.3.1, pp. 257--259.
277 */
278
279 #define B (1 << HALF_BITS) /* digit base */
280
281 /* Combine two `digits' to make a single two-digit number. */
282 #define COMBINE(a, b) (((u_long)(a) << HALF_BITS) | (b))
283
284 /* select a type for digits in base B: use unsigned short if they fit */
285 #if ULONG_MAX == 0xffffffff && USHRT_MAX >= 0xffff
286 typedef unsigned short digit;
287 #else
288 typedef u_long digit;
289 #endif
290
291 /*
292 * Shift p[0]..p[len] left `sh' bits, ignoring any bits that
293 * `fall out' the left (there never will be any such anyway).
294 * We may assume len >= 0. NOTE THAT THIS WRITES len+1 DIGITS.
295 */
296 static void
__shl(register digit * p,register int len,register int sh)297 __shl(register digit *p, register int len, register int sh)
298 {
299 register int i;
300
301 for (i = 0; i < len; i++)
302 p[i] = LHALF(p[i] << sh) | (p[i + 1] >> (HALF_BITS - sh));
303 p[i] = LHALF(p[i] << sh);
304 }
305
306 /*
307 * __qdivrem(u, v, rem) returns u/v and, optionally, sets *rem to u%v.
308 *
309 * We do this in base 2-sup-HALF_BITS, so that all intermediate products
310 * fit within u_long. As a consequence, the maximum length dividend and
311 * divisor are 4 `digits' in this base (they are shorter if they have
312 * leading zeros).
313 */
314 u_quad_t
__qdivrem(uq,vq,arq)315 __qdivrem(uq, vq, arq)
316 u_quad_t uq, vq, *arq;
317 {
318 union uu tmp;
319 digit *u, *v, *q;
320 register digit v1, v2;
321 u_long qhat, rhat, t;
322 int m, n, d, j, i;
323 digit uspace[5], vspace[5], qspace[5];
324
325 /*
326 * Take care of special cases: divide by zero, and u < v.
327 */
328 if (vq == 0) {
329 /* divide by zero. */
330 static volatile const unsigned int zero = 0;
331
332 tmp.ul[H] = tmp.ul[L] = 1 / zero;
333 if (arq)
334 *arq = uq;
335 return (tmp.q);
336 }
337 if (uq < vq) {
338 if (arq)
339 *arq = uq;
340 return (0);
341 }
342 u = &uspace[0];
343 v = &vspace[0];
344 q = &qspace[0];
345
346 /*
347 * Break dividend and divisor into digits in base B, then
348 * count leading zeros to determine m and n. When done, we
349 * will have:
350 * u = (u[1]u[2]...u[m+n]) sub B
351 * v = (v[1]v[2]...v[n]) sub B
352 * v[1] != 0
353 * 1 < n <= 4 (if n = 1, we use a different division algorithm)
354 * m >= 0 (otherwise u < v, which we already checked)
355 * m + n = 4
356 * and thus
357 * m = 4 - n <= 2
358 */
359 tmp.uq = uq;
360 u[0] = 0;
361 u[1] = HHALF(tmp.ul[H]);
362 u[2] = LHALF(tmp.ul[H]);
363 u[3] = HHALF(tmp.ul[L]);
364 u[4] = LHALF(tmp.ul[L]);
365 tmp.uq = vq;
366 v[1] = HHALF(tmp.ul[H]);
367 v[2] = LHALF(tmp.ul[H]);
368 v[3] = HHALF(tmp.ul[L]);
369 v[4] = LHALF(tmp.ul[L]);
370 for (n = 4; v[1] == 0; v++) {
371 if (--n == 1) {
372 u_long rbj; /* r*B+u[j] (not root boy jim) */
373 digit q1, q2, q3, q4;
374
375 /*
376 * Change of plan, per exercise 16.
377 * r = 0;
378 * for j = 1..4:
379 * q[j] = floor((r*B + u[j]) / v),
380 * r = (r*B + u[j]) % v;
381 * We unroll this completely here.
382 */
383 t = v[2]; /* nonzero, by definition */
384 q1 = u[1] / t;
385 rbj = COMBINE(u[1] % t, u[2]);
386 q2 = rbj / t;
387 rbj = COMBINE(rbj % t, u[3]);
388 q3 = rbj / t;
389 rbj = COMBINE(rbj % t, u[4]);
390 q4 = rbj / t;
391 if (arq)
392 *arq = rbj % t;
393 tmp.ul[H] = COMBINE(q1, q2);
394 tmp.ul[L] = COMBINE(q3, q4);
395 return (tmp.q);
396 }
397 }
398
399 /*
400 * By adjusting q once we determine m, we can guarantee that
401 * there is a complete four-digit quotient at &qspace[1] when
402 * we finally stop.
403 */
404 for (m = 4 - n; u[1] == 0; u++)
405 m--;
406 for (i = 4 - m; --i >= 0;)
407 q[i] = 0;
408 q += 4 - m;
409
410 /*
411 * Here we run Program D, translated from MIX to C and acquiring
412 * a few minor changes.
413 *
414 * D1: choose multiplier 1 << d to ensure v[1] >= B/2.
415 */
416 d = 0;
417 for (t = v[1]; t < B / 2; t <<= 1)
418 d++;
419 if (d > 0) {
420 __shl(&u[0], m + n, d); /* u <<= d */
421 __shl(&v[1], n - 1, d); /* v <<= d */
422 }
423 /*
424 * D2: j = 0.
425 */
426 j = 0;
427 v1 = v[1]; /* for D3 -- note that v[1..n] are constant */
428 v2 = v[2]; /* for D3 */
429 do {
430 register digit uj0, uj1, uj2;
431
432 /*
433 * D3: Calculate qhat (\^q, in TeX notation).
434 * Let qhat = min((u[j]*B + u[j+1])/v[1], B-1), and
435 * let rhat = (u[j]*B + u[j+1]) mod v[1].
436 * While rhat < B and v[2]*qhat > rhat*B+u[j+2],
437 * decrement qhat and increase rhat correspondingly.
438 * Note that if rhat >= B, v[2]*qhat < rhat*B.
439 */
440 uj0 = u[j + 0]; /* for D3 only -- note that u[j+...] change */
441 uj1 = u[j + 1]; /* for D3 only */
442 uj2 = u[j + 2]; /* for D3 only */
443 if (uj0 == v1) {
444 qhat = B;
445 rhat = uj1;
446 goto qhat_too_big;
447 } else {
448 u_long nn = COMBINE(uj0, uj1);
449 qhat = nn / v1;
450 rhat = nn % v1;
451 }
452 while (v2 * qhat > COMBINE(rhat, uj2)) {
453 qhat_too_big:
454 qhat--;
455 if ((rhat += v1) >= B)
456 break;
457 }
458 /*
459 * D4: Multiply and subtract.
460 * The variable `t' holds any borrows across the loop.
461 * We split this up so that we do not require v[0] = 0,
462 * and to eliminate a final special case.
463 */
464 for (t = 0, i = n; i > 0; i--) {
465 t = u[i + j] - v[i] * qhat - t;
466 u[i + j] = LHALF(t);
467 t = (B - HHALF(t)) & (B - 1);
468 }
469 t = u[j] - t;
470 u[j] = LHALF(t);
471 /*
472 * D5: test remainder.
473 * There is a borrow if and only if HHALF(t) is nonzero;
474 * in that (rare) case, qhat was too large (by exactly 1).
475 * Fix it by adding v[1..n] to u[j..j+n].
476 */
477 if (HHALF(t)) {
478 qhat--;
479 for (t = 0, i = n; i > 0; i--) { /* D6: add back. */
480 t += u[i + j] + v[i];
481 u[i + j] = LHALF(t);
482 t = HHALF(t);
483 }
484 u[j] = LHALF(u[j] + t);
485 }
486 q[j] = qhat;
487 } while (++j <= m); /* D7: loop on j. */
488
489 /*
490 * If caller wants the remainder, we have to calculate it as
491 * u[m..m+n] >> d (this is at most n digits and thus fits in
492 * u[m+1..m+n], but we may need more source digits).
493 */
494 if (arq) {
495 if (d) {
496 for (i = m + n; i > m; --i)
497 u[i] = (u[i] >> d) |
498 LHALF(u[i - 1] << (HALF_BITS - d));
499 u[i] = 0;
500 }
501 tmp.ul[H] = COMBINE(uspace[1], uspace[2]);
502 tmp.ul[L] = COMBINE(uspace[3], uspace[4]);
503 *arq = tmp.q;
504 }
505
506 tmp.ul[H] = COMBINE(qspace[1], qspace[2]);
507 tmp.ul[L] = COMBINE(qspace[3], qspace[4]);
508 return (tmp.q);
509 }
510
511 /*
512 * Return 0, 1, or 2 as a <, =, > b respectively.
513 * Neither a nor b are considered signed.
514 */
515 int
__ucmpdi2(a,b)516 __ucmpdi2(a, b)
517 u_quad_t a, b;
518 {
519 union uu aa, bb;
520
521 aa.uq = a;
522 bb.uq = b;
523 return (aa.ul[H] < bb.ul[H] ? 0 : aa.ul[H] > bb.ul[H] ? 2 :
524 aa.ul[L] < bb.ul[L] ? 0 : aa.ul[L] > bb.ul[L] ? 2 : 1);
525 }
526
527 /*
528 * Divide two unsigned quads.
529 */
530 u_quad_t
__udivdi3(a,b)531 __udivdi3(a, b)
532 u_quad_t a, b;
533 {
534
535 return (__qdivrem(a, b, (u_quad_t *)0));
536 }
537
538 /*
539 * Return remainder after dividing two unsigned quads.
540 */
541 u_quad_t
__umoddi3(a,b)542 __umoddi3(a, b)
543 u_quad_t a, b;
544 {
545 u_quad_t r;
546
547 (void)__qdivrem(a, b, &r);
548 return (r);
549 }
550
551 /*
552 * Divide two signed quads.
553 * This function is new in GCC 7.
554 */
555 quad_t
__divmoddi4(a,b,rem)556 __divmoddi4(a, b, rem)
557 quad_t a, b, *rem;
558 {
559 u_quad_t ua, ub, uq, ur;
560 int negq, negr;
561
562 if (a < 0)
563 ua = -(u_quad_t)a, negq = 1, negr = 1;
564 else
565 ua = a, negq = 0, negr = 0;
566 if (b < 0)
567 ub = -(u_quad_t)b, negq ^= 1;
568 else
569 ub = b;
570 uq = __qdivrem(ua, ub, &ur);
571 if (rem)
572 *rem = (negr ? -ur : ur);
573 return (negq ? -uq : uq);
574 }
575
576 #else
577 static int __attribute__((unused)) dummy;
578 #endif /*deined (_X86_) && !defined (__x86_64__)*/
579
580