diff options
Diffstat (limited to 'sys/libkern/muldi3.c')
-rw-r--r-- | sys/libkern/muldi3.c | 244 |
1 files changed, 0 insertions, 244 deletions
diff --git a/sys/libkern/muldi3.c b/sys/libkern/muldi3.c deleted file mode 100644 index 922f742..0000000 --- a/sys/libkern/muldi3.c +++ /dev/null @@ -1,244 +0,0 @@ -/*- - * Copyright (c) 1992, 1993 - * The Regents of the University of California. All rights reserved. - * - * This software was developed by the Computer Systems Engineering group - * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and - * contributed to Berkeley. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * 3. All advertising materials mentioning features or use of this software - * must display the following acknowledgement: - * This product includes software developed by the University of - * California, Berkeley and its contributors. - * 4. Neither the name of the University nor the names of its contributors - * may be used to endorse or promote products derived from this software - * without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - * - * $FreeBSD$ - */ - -#include <libkern/quad.h> - -/* - * Multiply two quads. - * - * Our algorithm is based on the following. Split incoming quad values - * u and v (where u,v >= 0) into - * - * u = 2^n u1 * u0 (n = number of bits in `u_long', usu. 32) - * - * and - * - * v = 2^n v1 * v0 - * - * Then - * - * uv = 2^2n u1 v1 + 2^n u1 v0 + 2^n v1 u0 + u0 v0 - * = 2^2n u1 v1 + 2^n (u1 v0 + v1 u0) + u0 v0 - * - * Now add 2^n u1 v1 to the first term and subtract it from the middle, - * and add 2^n u0 v0 to the last term and subtract it from the middle. - * This gives: - * - * uv = (2^2n + 2^n) (u1 v1) + - * (2^n) (u1 v0 - u1 v1 + u0 v1 - u0 v0) + - * (2^n + 1) (u0 v0) - * - * Factoring the middle a bit gives us: - * - * uv = (2^2n + 2^n) (u1 v1) + [u1v1 = high] - * (2^n) (u1 - u0) (v0 - v1) + [(u1-u0)... = mid] - * (2^n + 1) (u0 v0) [u0v0 = low] - * - * The terms (u1 v1), (u1 - u0) (v0 - v1), and (u0 v0) can all be done - * in just half the precision of the original. (Note that either or both - * of (u1 - u0) or (v0 - v1) may be negative.) - * - * This algorithm is from Knuth vol. 2 (2nd ed), section 4.3.3, p. 278. - * - * Since C does not give us a `long * long = quad' operator, we split - * our input quads into two longs, then split the two longs into two - * shorts. We can then calculate `short * short = long' in native - * arithmetic. - * - * Our product should, strictly speaking, be a `long quad', with 128 - * bits, but we are going to discard the upper 64. In other words, - * we are not interested in uv, but rather in (uv mod 2^2n). This - * makes some of the terms above vanish, and we get: - * - * (2^n)(high) + (2^n)(mid) + (2^n + 1)(low) - * - * or - * - * (2^n)(high + mid + low) + low - * - * Furthermore, `high' and `mid' can be computed mod 2^n, as any factor - * of 2^n in either one will also vanish. Only `low' need be computed - * mod 2^2n, and only because of the final term above. - */ -static quad_t __lmulq(u_long, u_long); - -quad_t -__muldi3(a, b) - quad_t a, b; -{ - union uu u, v, low, prod; - register u_long high, mid, udiff, vdiff; - register int negall, negmid; -#define u1 u.ul[H] -#define u0 u.ul[L] -#define v1 v.ul[H] -#define v0 v.ul[L] - - /* - * Get u and v such that u, v >= 0. When this is finished, - * u1, u0, v1, and v0 will be directly accessible through the - * longword fields. - */ - if (a >= 0) - u.q = a, negall = 0; - else - u.q = -a, negall = 1; - if (b >= 0) - v.q = b; - else - v.q = -b, negall ^= 1; - - if (u1 == 0 && v1 == 0) { - /* - * An (I hope) important optimization occurs when u1 and v1 - * are both 0. This should be common since most numbers - * are small. Here the product is just u0*v0. - */ - prod.q = __lmulq(u0, v0); - } else { - /* - * Compute the three intermediate products, remembering - * whether the middle term is negative. We can discard - * any upper bits in high and mid, so we can use native - * u_long * u_long => u_long arithmetic. - */ - low.q = __lmulq(u0, v0); - - if (u1 >= u0) - negmid = 0, udiff = u1 - u0; - else - negmid = 1, udiff = u0 - u1; - if (v0 >= v1) - vdiff = v0 - v1; - else - vdiff = v1 - v0, negmid ^= 1; - mid = udiff * vdiff; - - high = u1 * v1; - - /* - * Assemble the final product. - */ - prod.ul[H] = high + (negmid ? -mid : mid) + low.ul[L] + - low.ul[H]; - prod.ul[L] = low.ul[L]; - } - return (negall ? -prod.q : prod.q); -#undef u1 -#undef u0 -#undef v1 -#undef v0 -} - -/* - * Multiply two 2N-bit longs to produce a 4N-bit quad, where N is half - * the number of bits in a long (whatever that is---the code below - * does not care as long as quad.h does its part of the bargain---but - * typically N==16). - * - * We use the same algorithm from Knuth, but this time the modulo refinement - * does not apply. On the other hand, since N is half the size of a long, - * we can get away with native multiplication---none of our input terms - * exceeds (ULONG_MAX >> 1). - * - * Note that, for u_long l, the quad-precision result - * - * l << N - * - * splits into high and low longs as HHALF(l) and LHUP(l) respectively. - */ -static quad_t -__lmulq(u_long u, u_long v) -{ - u_long u1, u0, v1, v0, udiff, vdiff, high, mid, low; - u_long prodh, prodl, was; - union uu prod; - int neg; - - u1 = HHALF(u); - u0 = LHALF(u); - v1 = HHALF(v); - v0 = LHALF(v); - - low = u0 * v0; - - /* This is the same small-number optimization as before. */ - if (u1 == 0 && v1 == 0) - return (low); - - if (u1 >= u0) - udiff = u1 - u0, neg = 0; - else - udiff = u0 - u1, neg = 1; - if (v0 >= v1) - vdiff = v0 - v1; - else - vdiff = v1 - v0, neg ^= 1; - mid = udiff * vdiff; - - high = u1 * v1; - - /* prod = (high << 2N) + (high << N); */ - prodh = high + HHALF(high); - prodl = LHUP(high); - - /* if (neg) prod -= mid << N; else prod += mid << N; */ - if (neg) { - was = prodl; - prodl -= LHUP(mid); - prodh -= HHALF(mid) + (prodl > was); - } else { - was = prodl; - prodl += LHUP(mid); - prodh += HHALF(mid) + (prodl < was); - } - - /* prod += low << N */ - was = prodl; - prodl += LHUP(low); - prodh += HHALF(low) + (prodl < was); - /* ... + low; */ - if ((prodl += low) < low) - prodh++; - - /* return 4N-bit product */ - prod.ul[H] = prodh; - prod.ul[L] = prodl; - return (prod.q); -} |