* 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.
* %sccs.include.redist.c%
#if defined(LIBC_SCCS) && !defined(lint)
static char sccsid
[] = "@(#)muldi3.c 8.1 (Berkeley) %G%";
#endif /* LIBC_SCCS and not lint */
* 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)
* 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.
* uv = (2^2n + 2^n) (u1 v1) +
* (2^n) (u1 v0 - u1 v1 + u0 v1 - 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
* 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)
* (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
);
union uu u
, v
, low
, prod
;
register u_long high
, mid
, udiff
, vdiff
;
register int negall
, negmid
;
* Get u and v such that u, v >= 0. When this is finished,
* u1, u0, v1, and v0 will be directly accessible through the
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
);
* 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.
negmid
= 0, udiff
= u1
- u0
;
negmid
= 1, udiff
= u0
- u1
;
vdiff
= v1
- v0
, negmid
^= 1;
* Assemble the final product.
prod
.ul
[H
] = high
+ (negmid
? -mid
: mid
) + low
.ul
[L
] +
return (negall
? -prod
.q
: prod
.q
);
* 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
* 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
* splits into high and low longs as HHALF(l) and LHUP(l) respectively.
__lmulq(u_long u
, u_long v
)
u_long u1
, u0
, v1
, v0
, udiff
, vdiff
, high
, mid
, low
;
u_long prodh
, prodl
, was
;
/* This is the same small-number optimization as before. */
udiff
= u1
- u0
, neg
= 0;
udiff
= u0
- u1
, neg
= 1;
vdiff
= v1
- v0
, neg
^= 1;
/* prod = (high << 2N) + (high << N); */
prodh
= high
+ HHALF(high
);
/* if (neg) prod -= mid << N; else prod += mid << N; */
prodh
-= HHALF(mid
) + (prodl
> was
);
prodh
+= HHALF(mid
) + (prodl
< was
);
prodh
+= HHALF(low
) + (prodl
< was
);
if ((prodl
+= low
) < low
)
/* return 4N-bit product */