3 /*-
4  * SPDX-License-Identifier: BSD-3-Clause
10  * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
11  * contributed to Berkeley.
22  *    may be used to endorse or promote products derived from this software
26  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
30  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
57  * Now add 2^n u1 v1 to the first term and subtract it from the middle,
58  * and add 2^n u0 v0 to the last term and subtract it from the middle.
62  *	         (2^n)    (u1 v0 - u1 v1 + u0 v1 - u0 v0)  +
67  *	uv = (2^2n + 2^n) (u1 v1)  +			[u1v1 = high]
68  *		 (2^n)    (u1 - u0) (v0 - v1)  +	[(u1-u0)... = mid]
69  *	       (2^n + 1)  (u0 v0)			[u0v0 = low]
71  * The terms (u1 v1), (u1 - u0) (v0 - v1), and (u0 v0) can all be done
73  * of (u1 - u0) or (v0 - v1) may be negative.)
83  * bits, but we are going to discard the upper 64.  In other words,
87  *	(2^n)(high) + (2^n)(mid) + (2^n + 1)(low)
91  *	(2^n)(high + mid + low) + low
93  * Furthermore, `high' and `mid' can be computed mod 2^n, as any factor
94  * of 2^n in either one will also vanish.  Only `low' need be computed
103 	union uu u, v, low, prod;  in __muldi3()  local
104 	u_int high, mid, udiff, vdiff;  in __muldi3()  local
119 		u.q = -a, negall = 1;  in __muldi3()
123 		v.q = -b, negall ^= 1;  in __muldi3()
136 		 * any upper bits in high and mid, so we can use native  in __muldi3()
139 		low.q = __lmulq(u0, v0);  in __muldi3()
142 			negmid = 0, udiff = u1 - u0;  in __muldi3()
144 			negmid = 1, udiff = u0 - u1;  in __muldi3()
146 			vdiff = v0 - v1;  in __muldi3()
148 			vdiff = v1 - v0, negmid ^= 1;  in __muldi3()
151 		high = u1 * v1;  in __muldi3()
156 		prod.ul[H] = high + (negmid ? -mid : mid) + low.ul[L] +  in __muldi3()
157 		    low.ul[H];  in __muldi3()
158 		prod.ul[L] = low.ul[L];  in __muldi3()
160 	return (negall ? -prod.q : prod.q);  in __muldi3()
168  * Multiply two 2N-bit ints to produce a 4N-bit quad, where N is half
169  * the number of bits in an int (whatever that is---the code below
170  * does not care as long as quad.h does its part of the bargain---but
175  * we can get away with native multiplication---none of our input terms
178  * Note that, for u_int l, the quad-precision result
182  * splits into high and low ints as HHALF(l) and LHUP(l) respectively.
187 	u_int u1, u0, v1, v0, udiff, vdiff, high, mid, low;  in __lmulq()  local
197 	low = u0 * v0;  in __lmulq()
199 	/* This is the same small-number optimization as before. */  in __lmulq()
201 		return (low);  in __lmulq()
204 		udiff = u1 - u0, neg = 0;  in __lmulq()
206 		udiff = u0 - u1, neg = 1;  in __lmulq()
208 		vdiff = v0 - v1;  in __lmulq()
210 		vdiff = v1 - v0, neg ^= 1;  in __lmulq()
213 	high = u1 * v1;  in __lmulq()
215 	/* prod = (high << 2N) + (high << N); */  in __lmulq()
216 	prodh = high + HHALF(high);  in __lmulq()
217 	prodl = LHUP(high);  in __lmulq()
219 	/* if (neg) prod -= mid << N; else prod += mid << N; */  in __lmulq()
222 		prodl -= LHUP(mid);  in __lmulq()
223 		prodh -= HHALF(mid) + (prodl > was);  in __lmulq()
230 	/* prod += low << N */  in __lmulq()
232 	prodl += LHUP(low);  in __lmulq()
233 	prodh += HHALF(low) + (prodl < was);  in __lmulq()
234 	/* ... + low; */  in __lmulq()
235 	if ((prodl += low) < low)  in __lmulq()
238 	/* return 4N-bit product */  in __lmulq()