29  * without the usual one-word header. Value is split into 31-bit words,
30  * each stored in a 32-bit slot (top bit is zero) in little-endian
37  * Negate big integer conditionally. The value consists of 'len' words,
43 cond_negate(uint32_t *a, size_t len, uint32_t ctl)  in cond_negate()  argument
49 	xm = -ctl >> 1;  in cond_negate()
50 	for (k = 0; k < len; k ++) {  in cond_negate()
63  *   if neg = 1, then -m <= a < 0
71 finish_mod(uint32_t *a, size_t len, const uint32_t *m, uint32_t neg)  in finish_mod()  argument
83 	for (k = 0; k < len; k ++) {  in finish_mod()
88 		cc = (aw - mw - cc) >> 31;  in finish_mod()
97 	xm = -neg >> 1;  in finish_mod()
98 	ym = -(neg | (1 - cc));  in finish_mod()
100 	for (k = 0; k < len; k ++) {  in finish_mod()
105 		aw = aw - mw - cc;  in finish_mod()
113  *   a <- (a*pa+b*pb)/(2^31)
114  *   b <- (a*qa+b*qb)/(2^31)
115  * The division is assumed to be exact (i.e. the low word is dropped).
126 co_reduce(uint32_t *a, uint32_t *b, size_t len,  in co_reduce()  argument
135 	for (k = 0; k < len; k ++) {  in co_reduce()
144 		 *   0 <= wa <= 2^31 - 1  in co_reduce()
145 		 *   0 <= wb <= 2^31 - 1  in co_reduce()
146 		 *   |cca| <= 2^32 - 1  in co_reduce()
148 		 *   |za| <= (2^31-1)*(2^32) + (2^32-1) = 2^63 - 1  in co_reduce()
150 		 * Thus, the new value of cca is such that |cca| <= 2^32 - 1.  in co_reduce()
158 			a[k - 1] = za & 0x7FFFFFFF;  in co_reduce()
159 			b[k - 1] = zb & 0x7FFFFFFF;  in co_reduce()
164 		 * right-shift; since, in C, right-shifting a signed  in co_reduce()
165 		 * negative value is implementation-defined, we use a  in co_reduce()
171 		tta = (tta ^ M) - M;  in co_reduce()
172 		ttb = (ttb ^ M) - M;  in co_reduce()
177 	a[len - 1] = (uint32_t)cca;  in co_reduce()
178 	b[len - 1] = (uint32_t)ccb;  in co_reduce()
182 	cond_negate(a, len, nega);  in co_reduce()
183 	cond_negate(b, len, negb);  in co_reduce()
189  *   a <- (a*pa+b*pb)/(2^31) mod m
190  *   b <- (a*qa+b*qb)/(2^31) mod m
192  * m0i is equal to -1/m[0] mod 2^31.
199 co_reduce_mod(uint32_t *a, uint32_t *b, size_t len,  in co_reduce_mod()  argument
211 	for (k = 0; k < len; k ++) {  in co_reduce_mod()
227 			a[k - 1] = (uint32_t)za & 0x7FFFFFFF;  in co_reduce_mod()
228 			b[k - 1] = (uint32_t)zb & 0x7FFFFFFF;  in co_reduce_mod()
234 		tta = (tta ^ M) - M;  in co_reduce_mod()
235 		ttb = (ttb ^ M) - M;  in co_reduce_mod()
240 	a[len - 1] = (uint32_t)cca;  in co_reduce_mod()
241 	b[len - 1] = (uint32_t)ccb;  in co_reduce_mod()
245 	 *   -m <= a < 2*m  in co_reduce_mod()
246 	 *   -m <= b < 2*m  in co_reduce_mod()
248 	 * The top word of 'a' and 'b' may have a 32-th bit set.  in co_reduce_mod()
251 	finish_mod(a, len, m, (uint32_t)((uint64_t)cca >> 63));  in co_reduce_mod()
252 	finish_mod(b, len, m, (uint32_t)((uint64_t)ccb >> 63));  in co_reduce_mod()
276 	 *   - If a is even, then a <- a/2 and u <- u/2 mod m.  in br_i31_moddiv()
277 	 *   - Otherwise, if b is even, then b <- b/2 and v <- v/2 mod m.  in br_i31_moddiv()
278 	 *   - Otherwise, if a > b, then a <- (a-b)/2 and u <- (u-v)/2 mod m.  in br_i31_moddiv()
279 	 *   - Otherwise, b <- (b-a)/2 and v <- (v-u)/2 mod m.  in br_i31_moddiv()
287 	 * if m has bit length k bits, then 2k-2 steps are sufficient.  in br_i31_moddiv()
290 	 * Though complexity is quadratic in the size of m, the bit-by-bit  in br_i31_moddiv()
303 	 * the division being exact.  in br_i31_moddiv()
307 	 * pa with -pa, and pb with -pb. The total length of a and b is  in br_i31_moddiv()
317 	size_t len, k;  in br_i31_moddiv()  local
321 	len = (m[0] + 31) >> 5;  in br_i31_moddiv()
323 	b = a + len;  in br_i31_moddiv()
325 	v = b + len;  in br_i31_moddiv()
326 	memcpy(a, y + 1, len * sizeof *y);  in br_i31_moddiv()
327 	memcpy(b, m + 1, len * sizeof *m);  in br_i31_moddiv()
328 	memset(v, 0, len * sizeof *v);  in br_i31_moddiv()
334 	for (num = ((m[0] - (m[0] >> 5)) << 1) + 30; num >= 30; num -= 30) {  in br_i31_moddiv()
346 		 * (a[j] << 31) + a[j - 1], and (b[j] << 31) + b[j - 1].  in br_i31_moddiv()
350 		c0 = (uint32_t)-1;  in br_i31_moddiv()
351 		c1 = (uint32_t)-1;  in br_i31_moddiv()
356 		j = len;  in br_i31_moddiv()
357 		while (j -- > 0) {  in br_i31_moddiv()
367 			c0 &= (((aw | bw) + 0x7FFFFFFF) >> 31) - (uint32_t)1;  in br_i31_moddiv()
402 			 *   a <- (a-b)/2 if: a is odd, b is odd, a_hi > b_hi  in br_i31_moddiv()
403 			 *   b <- (b-a)/2 if: a is odd, b is odd, a_hi <= b_hi  in br_i31_moddiv()
404 			 *   a <- a/2 if: a is even  in br_i31_moddiv()
405 			 *   b <- b/2 if: a is odd, b is even  in br_i31_moddiv()
409 			 * non-multiplication by 2.  in br_i31_moddiv()
417 			 * so we inline a 64-bit version here.  in br_i31_moddiv()
419 			rz = b_hi - a_hi;  in br_i31_moddiv()
442 			a_lo -= b_lo & -cAB;  in br_i31_moddiv()
443 			a_hi -= b_hi & -(uint64_t)cAB;  in br_i31_moddiv()
444 			pa -= qa & -(int64_t)cAB;  in br_i31_moddiv()
445 			pb -= qb & -(int64_t)cAB;  in br_i31_moddiv()
446 			b_lo -= a_lo & -cBA;  in br_i31_moddiv()
447 			b_hi -= a_hi & -(uint64_t)cBA;  in br_i31_moddiv()
448 			qa -= pa & -(int64_t)cBA;  in br_i31_moddiv()
449 			qb -= pb & -(int64_t)cBA;  in br_i31_moddiv()
454 			a_lo += a_lo & (cA - 1);  in br_i31_moddiv()
455 			pa += pa & ((int64_t)cA - 1);  in br_i31_moddiv()
456 			pb += pb & ((int64_t)cA - 1);  in br_i31_moddiv()
457 			a_hi ^= (a_hi ^ (a_hi >> 1)) & -(uint64_t)cA;  in br_i31_moddiv()
458 			b_lo += b_lo & -cA;  in br_i31_moddiv()
459 			qa += qa & -(int64_t)cA;  in br_i31_moddiv()
460 			qb += qb & -(int64_t)cA;  in br_i31_moddiv()
461 			b_hi ^= (b_hi ^ (b_hi >> 1)) & ((uint64_t)cA - 1);  in br_i31_moddiv()
467 		r = co_reduce(a, b, len, pa, pb, qa, qb);  in br_i31_moddiv()
468 		pa -= pa * ((r & 1) << 1);  in br_i31_moddiv()
469 		pb -= pb * ((r & 1) << 1);  in br_i31_moddiv()
470 		qa -= qa * (r & 2);  in br_i31_moddiv()
471 		qb -= qb * (r & 2);  in br_i31_moddiv()
472 		co_reduce_mod(u, v, len, pa, pb, qa, qb, m + 1, m0i);  in br_i31_moddiv()
483 	for (k = 1; k < len; k ++) {  in br_i31_moddiv()