35 	 * We have xu <= 2^31-1 and f <= 2^31-1.  in br_i31_montymul()
36 	 * Thus, if d <= 2*m-1 on input, then:  in br_i31_montymul()
37 	 *   2*m-1 + 2*(2^31-1)*m <= (2^32)*m-1  in br_i31_montymul()
41 	 * which can thus be only 0 or 1.  in br_i31_montymul()
43 	size_t len, len4, u, v;  in br_i31_montymul()  local
53 		 *   d[v+1] <= 2^31-1  in br_i31_montymul()
54 		 *   xu*y[v+1] <= (2^31-1)*(2^31-1)  in br_i31_montymul()
55 		 *   f*m[v+1] <= (2^31-1)*(2^31-1)  in br_i31_montymul()
56 		 *   r <= 2^32-1  in br_i31_montymul()
57 		 *   (2^31-1) + 2*(2^31-1)*(2^31-1) + (2^32-1) = 2^63 - 2^31  in br_i31_montymul()
58 		 * After division by 2^31, the new r is then at most 2^32-1  in br_i31_montymul()
72 		xu = x[u + 1];  in br_i31_montymul()
73 		f = MUL31_lo((d[1] + MUL31_lo(x[u + 1], y[1])), m0i);  in br_i31_montymul()
76 		for (v = 0; v < len4; v += 4) {  in br_i31_montymul()
79 			z = (uint64_t)d[v + 1] + MUL31(xu, y[v + 1])  in br_i31_montymul()
80 				+ MUL31(f, m[v + 1]) + r;  in br_i31_montymul()
82 			d[v + 0] = (uint32_t)z & 0x7FFFFFFF;  in br_i31_montymul()
83 			z = (uint64_t)d[v + 2] + MUL31(xu, y[v + 2])  in br_i31_montymul()
84 				+ MUL31(f, m[v + 2]) + r;  in br_i31_montymul()
86 			d[v + 1] = (uint32_t)z & 0x7FFFFFFF;  in br_i31_montymul()
87 			z = (uint64_t)d[v + 3] + MUL31(xu, y[v + 3])  in br_i31_montymul()
88 				+ MUL31(f, m[v + 3]) + r;  in br_i31_montymul()
90 			d[v + 2] = (uint32_t)z & 0x7FFFFFFF;  in br_i31_montymul()
91 			z = (uint64_t)d[v + 4] + MUL31(xu, y[v + 4])  in br_i31_montymul()
92 				+ MUL31(f, m[v + 4]) + r;  in br_i31_montymul()
94 			d[v + 3] = (uint32_t)z & 0x7FFFFFFF;  in br_i31_montymul()
96 		for (; v < len; v ++) {  in br_i31_montymul()
99 			z = (uint64_t)d[v + 1] + MUL31(xu, y[v + 1])  in br_i31_montymul()
100 				+ MUL31(f, m[v + 1]) + r;  in br_i31_montymul()
102 			d[v] = (uint32_t)z & 0x7FFFFFFF;  in br_i31_montymul()
106 		 * Since the new dh can only be 0 or 1, the addition of  in br_i31_montymul()