mod.gno
14.52 Kb · 608 lines
1package uint256
2
3import (
4 "math/bits"
5)
6
7// Reciprocal computes the 320-bit reciprocal estimate used by Barrett reduction.
8//
9// The implementation is specialized for a four-word modulus with m[3] non-zero
10// (2^192 <= m < 2^256). For a modulus whose most-significant word is zero it
11// returns an all-zero estimate instead of running the refinement steps.
12//
13// Parameters:
14// - m: Four-word unsigned modulus; its most-significant word selects the supported range.
15//
16// Returns:
17// - mu: Five little-endian uint64 words containing the reciprocal estimate for m.
18func Reciprocal(m *Uint) (mu [5]uint64) {
19 if m[3] == 0 {
20 return mu
21 }
22
23 s := bits.LeadingZeros64(m[3]) // Replace with leadingZeros(m) for general case
24 p := 255 - s // floor(log_2(m)), m>0
25
26 // 0 or a power of 2?
27
28 // Check if at least one bit is set in m[2], m[1] or m[0],
29 // or at least two bits in m[3]
30
31 if m[0]|m[1]|m[2]|(m[3]&(m[3]-1)) == 0 {
32
33 mu[4] = 18446744073709551615 >> uint(p&63)
34 mu[3] = 18446744073709551615
35 mu[2] = 18446744073709551615
36 mu[1] = 18446744073709551615
37 mu[0] = 18446744073709551615
38
39 return mu
40 }
41
42 // Maximise division precision by left-aligning divisor
43
44 var (
45 y Uint // left-aligned copy of m
46 r0 uint32 // estimate of 2^31/y
47 )
48
49 y.Lsh(m, uint(s)) // 1/2 < y < 1
50
51 // Extract most significant 32 bits
52
53 yh := uint32(y[3] >> 32)
54
55 if yh == 0x80000000 { // Avoid overflow in division
56 r0 = 0xffffffff
57 } else {
58 r0, _ = bits.Div32(0x80000000, 0, yh)
59 }
60
61 // First iteration: 32 -> 64
62
63 t1 := uint64(r0) // 2^31/y
64 t1 *= t1 // 2^62/y^2
65 t1, _ = bits.Mul64(t1, y[3]) // 2^62/y^2 * 2^64/y / 2^64 = 2^62/y
66
67 r1 := uint64(r0) << 32 // 2^63/y
68 r1 -= t1 // 2^63/y - 2^62/y = 2^62/y
69 r1 *= 2 // 2^63/y
70
71 if (r1 | (y[3] << 1)) == 0 {
72 r1 = 18446744073709551615
73 }
74
75 // Second iteration: 64 -> 128
76
77 // square: 2^126/y^2
78 a2h, a2l := bits.Mul64(r1, r1)
79
80 // multiply by y: e2h:e2l:b2h = 2^126/y^2 * 2^128/y / 2^128 = 2^126/y
81 b2h, _ := bits.Mul64(a2l, y[2])
82 c2h, c2l := bits.Mul64(a2l, y[3])
83 d2h, d2l := bits.Mul64(a2h, y[2])
84 e2h, e2l := bits.Mul64(a2h, y[3])
85
86 b2h, c := bits.Add64(b2h, c2l, 0)
87 e2l, c = bits.Add64(e2l, c2h, c)
88 e2h, _ = bits.Add64(e2h, 0, c)
89
90 _, c = bits.Add64(b2h, d2l, 0)
91 e2l, c = bits.Add64(e2l, d2h, c)
92 e2h, _ = bits.Add64(e2h, 0, c)
93
94 // subtract: t2h:t2l = 2^127/y - 2^126/y = 2^126/y
95 t2l, b := bits.Sub64(0, e2l, 0)
96 t2h, _ := bits.Sub64(r1, e2h, b)
97
98 // double: r2h:r2l = 2^127/y
99 r2l, c := bits.Add64(t2l, t2l, 0)
100 r2h, _ := bits.Add64(t2h, t2h, c)
101
102 if (r2h | r2l | (y[3] << 1)) == 0 {
103 r2h = 18446744073709551615
104 r2l = 18446744073709551615
105 }
106
107 // Third iteration: 128 -> 192
108
109 // square r2 (keep 256 bits): 2^190/y^2
110 a3h, a3l := bits.Mul64(r2l, r2l)
111 b3h, b3l := bits.Mul64(r2l, r2h)
112 c3h, c3l := bits.Mul64(r2h, r2h)
113
114 a3h, c = bits.Add64(a3h, b3l, 0)
115 c3l, c = bits.Add64(c3l, b3h, c)
116 c3h, _ = bits.Add64(c3h, 0, c)
117
118 a3h, c = bits.Add64(a3h, b3l, 0)
119 c3l, c = bits.Add64(c3l, b3h, c)
120 c3h, _ = bits.Add64(c3h, 0, c)
121
122 // multiply by y: q = 2^190/y^2 * 2^192/y / 2^192 = 2^190/y
123
124 x0 := a3l
125 x1 := a3h
126 x2 := c3l
127 x3 := c3h
128
129 var q0, q1, q2, q3, q4, t0 uint64
130
131 q0, _ = bits.Mul64(x2, y[0])
132 q1, t0 = bits.Mul64(x3, y[0])
133 q0, c = bits.Add64(q0, t0, 0)
134 q1, _ = bits.Add64(q1, 0, c)
135
136 t1, _ = bits.Mul64(x1, y[1])
137 q0, c = bits.Add64(q0, t1, 0)
138 q2, t0 = bits.Mul64(x3, y[1])
139 q1, c = bits.Add64(q1, t0, c)
140 q2, _ = bits.Add64(q2, 0, c)
141
142 t1, t0 = bits.Mul64(x2, y[1])
143 q0, c = bits.Add64(q0, t0, 0)
144 q1, c = bits.Add64(q1, t1, c)
145 q2, _ = bits.Add64(q2, 0, c)
146
147 t1, t0 = bits.Mul64(x1, y[2])
148 q0, c = bits.Add64(q0, t0, 0)
149 q1, c = bits.Add64(q1, t1, c)
150 q3, t0 = bits.Mul64(x3, y[2])
151 q2, c = bits.Add64(q2, t0, c)
152 q3, _ = bits.Add64(q3, 0, c)
153
154 t1, _ = bits.Mul64(x0, y[2])
155 q0, c = bits.Add64(q0, t1, 0)
156 t1, t0 = bits.Mul64(x2, y[2])
157 q1, c = bits.Add64(q1, t0, c)
158 q2, c = bits.Add64(q2, t1, c)
159 q3, _ = bits.Add64(q3, 0, c)
160
161 t1, t0 = bits.Mul64(x1, y[3])
162 q1, c = bits.Add64(q1, t0, 0)
163 q2, c = bits.Add64(q2, t1, c)
164 q4, t0 = bits.Mul64(x3, y[3])
165 q3, c = bits.Add64(q3, t0, c)
166 q4, _ = bits.Add64(q4, 0, c)
167
168 t1, t0 = bits.Mul64(x0, y[3])
169 q0, c = bits.Add64(q0, t0, 0)
170 q1, c = bits.Add64(q1, t1, c)
171 t1, t0 = bits.Mul64(x2, y[3])
172 q2, c = bits.Add64(q2, t0, c)
173 q3, c = bits.Add64(q3, t1, c)
174 q4, _ = bits.Add64(q4, 0, c)
175
176 // subtract: t3 = 2^191/y - 2^190/y = 2^190/y
177 _, b = bits.Sub64(0, q0, 0)
178 _, b = bits.Sub64(0, q1, b)
179 t3l, b := bits.Sub64(0, q2, b)
180 t3m, b := bits.Sub64(r2l, q3, b)
181 t3h, _ := bits.Sub64(r2h, q4, b)
182
183 // double: r3 = 2^191/y
184 r3l, c := bits.Add64(t3l, t3l, 0)
185 r3m, c := bits.Add64(t3m, t3m, c)
186 r3h, _ := bits.Add64(t3h, t3h, c)
187
188 // Fourth iteration: 192 -> 320
189
190 // square r3
191
192 a4h, a4l := bits.Mul64(r3l, r3l)
193 b4h, b4l := bits.Mul64(r3l, r3m)
194 c4h, c4l := bits.Mul64(r3l, r3h)
195 d4h, d4l := bits.Mul64(r3m, r3m)
196 e4h, e4l := bits.Mul64(r3m, r3h)
197 f4h, f4l := bits.Mul64(r3h, r3h)
198
199 b4h, c = bits.Add64(b4h, c4l, 0)
200 e4l, c = bits.Add64(e4l, c4h, c)
201 e4h, _ = bits.Add64(e4h, 0, c)
202
203 a4h, c = bits.Add64(a4h, b4l, 0)
204 d4l, c = bits.Add64(d4l, b4h, c)
205 d4h, c = bits.Add64(d4h, e4l, c)
206 f4l, c = bits.Add64(f4l, e4h, c)
207 f4h, _ = bits.Add64(f4h, 0, c)
208
209 a4h, c = bits.Add64(a4h, b4l, 0)
210 d4l, c = bits.Add64(d4l, b4h, c)
211 d4h, c = bits.Add64(d4h, e4l, c)
212 f4l, c = bits.Add64(f4l, e4h, c)
213 f4h, _ = bits.Add64(f4h, 0, c)
214
215 // multiply by y
216
217 x1, x0 = bits.Mul64(d4h, y[0])
218 x3, x2 = bits.Mul64(f4h, y[0])
219 t1, t0 = bits.Mul64(f4l, y[0])
220 x1, c = bits.Add64(x1, t0, 0)
221 x2, c = bits.Add64(x2, t1, c)
222 x3, _ = bits.Add64(x3, 0, c)
223
224 t1, t0 = bits.Mul64(d4h, y[1])
225 x1, c = bits.Add64(x1, t0, 0)
226 x2, c = bits.Add64(x2, t1, c)
227 x4, t0 := bits.Mul64(f4h, y[1])
228 x3, c = bits.Add64(x3, t0, c)
229 x4, _ = bits.Add64(x4, 0, c)
230 t1, t0 = bits.Mul64(d4l, y[1])
231 x0, c = bits.Add64(x0, t0, 0)
232 x1, c = bits.Add64(x1, t1, c)
233 t1, t0 = bits.Mul64(f4l, y[1])
234 x2, c = bits.Add64(x2, t0, c)
235 x3, c = bits.Add64(x3, t1, c)
236 x4, _ = bits.Add64(x4, 0, c)
237
238 t1, t0 = bits.Mul64(a4h, y[2])
239 x0, c = bits.Add64(x0, t0, 0)
240 x1, c = bits.Add64(x1, t1, c)
241 t1, t0 = bits.Mul64(d4h, y[2])
242 x2, c = bits.Add64(x2, t0, c)
243 x3, c = bits.Add64(x3, t1, c)
244 x5, t0 := bits.Mul64(f4h, y[2])
245 x4, c = bits.Add64(x4, t0, c)
246 x5, _ = bits.Add64(x5, 0, c)
247 t1, t0 = bits.Mul64(d4l, y[2])
248 x1, c = bits.Add64(x1, t0, 0)
249 x2, c = bits.Add64(x2, t1, c)
250 t1, t0 = bits.Mul64(f4l, y[2])
251 x3, c = bits.Add64(x3, t0, c)
252 x4, c = bits.Add64(x4, t1, c)
253 x5, _ = bits.Add64(x5, 0, c)
254
255 t1, t0 = bits.Mul64(a4h, y[3])
256 x1, c = bits.Add64(x1, t0, 0)
257 x2, c = bits.Add64(x2, t1, c)
258 t1, t0 = bits.Mul64(d4h, y[3])
259 x3, c = bits.Add64(x3, t0, c)
260 x4, c = bits.Add64(x4, t1, c)
261 x6, t0 := bits.Mul64(f4h, y[3])
262 x5, c = bits.Add64(x5, t0, c)
263 x6, _ = bits.Add64(x6, 0, c)
264 t1, t0 = bits.Mul64(a4l, y[3])
265 x0, c = bits.Add64(x0, t0, 0)
266 x1, c = bits.Add64(x1, t1, c)
267 t1, t0 = bits.Mul64(d4l, y[3])
268 x2, c = bits.Add64(x2, t0, c)
269 x3, c = bits.Add64(x3, t1, c)
270 t1, t0 = bits.Mul64(f4l, y[3])
271 x4, c = bits.Add64(x4, t0, c)
272 x5, c = bits.Add64(x5, t1, c)
273 x6, _ = bits.Add64(x6, 0, c)
274
275 // subtract
276 _, b = bits.Sub64(0, x0, 0)
277 _, b = bits.Sub64(0, x1, b)
278 r4l, b := bits.Sub64(0, x2, b)
279 r4k, b := bits.Sub64(0, x3, b)
280 r4j, b := bits.Sub64(r3l, x4, b)
281 r4i, b := bits.Sub64(r3m, x5, b)
282 r4h, _ := bits.Sub64(r3h, x6, b)
283
284 // Multiply candidate for 1/4y by y, with full precision
285
286 x0 = r4l
287 x1 = r4k
288 x2 = r4j
289 x3 = r4i
290 x4 = r4h
291
292 q1, q0 = bits.Mul64(x0, y[0])
293 q3, q2 = bits.Mul64(x2, y[0])
294 q5, q4 := bits.Mul64(x4, y[0])
295
296 t1, t0 = bits.Mul64(x1, y[0])
297 q1, c = bits.Add64(q1, t0, 0)
298 q2, c = bits.Add64(q2, t1, c)
299 t1, t0 = bits.Mul64(x3, y[0])
300 q3, c = bits.Add64(q3, t0, c)
301 q4, c = bits.Add64(q4, t1, c)
302 q5, _ = bits.Add64(q5, 0, c)
303
304 t1, t0 = bits.Mul64(x0, y[1])
305 q1, c = bits.Add64(q1, t0, 0)
306 q2, c = bits.Add64(q2, t1, c)
307 t1, t0 = bits.Mul64(x2, y[1])
308 q3, c = bits.Add64(q3, t0, c)
309 q4, c = bits.Add64(q4, t1, c)
310 q6, t0 := bits.Mul64(x4, y[1])
311 q5, c = bits.Add64(q5, t0, c)
312 q6, _ = bits.Add64(q6, 0, c)
313
314 t1, t0 = bits.Mul64(x1, y[1])
315 q2, c = bits.Add64(q2, t0, 0)
316 q3, c = bits.Add64(q3, t1, c)
317 t1, t0 = bits.Mul64(x3, y[1])
318 q4, c = bits.Add64(q4, t0, c)
319 q5, c = bits.Add64(q5, t1, c)
320 q6, _ = bits.Add64(q6, 0, c)
321
322 t1, t0 = bits.Mul64(x0, y[2])
323 q2, c = bits.Add64(q2, t0, 0)
324 q3, c = bits.Add64(q3, t1, c)
325 t1, t0 = bits.Mul64(x2, y[2])
326 q4, c = bits.Add64(q4, t0, c)
327 q5, c = bits.Add64(q5, t1, c)
328 q7, t0 := bits.Mul64(x4, y[2])
329 q6, c = bits.Add64(q6, t0, c)
330 q7, _ = bits.Add64(q7, 0, c)
331
332 t1, t0 = bits.Mul64(x1, y[2])
333 q3, c = bits.Add64(q3, t0, 0)
334 q4, c = bits.Add64(q4, t1, c)
335 t1, t0 = bits.Mul64(x3, y[2])
336 q5, c = bits.Add64(q5, t0, c)
337 q6, c = bits.Add64(q6, t1, c)
338 q7, _ = bits.Add64(q7, 0, c)
339
340 t1, t0 = bits.Mul64(x0, y[3])
341 q3, c = bits.Add64(q3, t0, 0)
342 q4, c = bits.Add64(q4, t1, c)
343 t1, t0 = bits.Mul64(x2, y[3])
344 q5, c = bits.Add64(q5, t0, c)
345 q6, c = bits.Add64(q6, t1, c)
346 q8, t0 := bits.Mul64(x4, y[3])
347 q7, c = bits.Add64(q7, t0, c)
348 q8, _ = bits.Add64(q8, 0, c)
349
350 t1, t0 = bits.Mul64(x1, y[3])
351 q4, c = bits.Add64(q4, t0, 0)
352 q5, c = bits.Add64(q5, t1, c)
353 t1, t0 = bits.Mul64(x3, y[3])
354 q6, c = bits.Add64(q6, t0, c)
355 q7, c = bits.Add64(q7, t1, c)
356 q8, _ = bits.Add64(q8, 0, c)
357
358 // Final adjustment
359
360 // subtract q from 1/4
361 _, b = bits.Sub64(0, q0, 0)
362 _, b = bits.Sub64(0, q1, b)
363 _, b = bits.Sub64(0, q2, b)
364 _, b = bits.Sub64(0, q3, b)
365 _, b = bits.Sub64(0, q4, b)
366 _, b = bits.Sub64(0, q5, b)
367 _, b = bits.Sub64(0, q6, b)
368 _, b = bits.Sub64(0, q7, b)
369 _, b = bits.Sub64(uint64(1)<<62, q8, b)
370
371 // decrement the result
372 x0, t := bits.Sub64(r4l, 1, 0)
373 x1, t = bits.Sub64(r4k, 0, t)
374 x2, t = bits.Sub64(r4j, 0, t)
375 x3, t = bits.Sub64(r4i, 0, t)
376 x4, _ = bits.Sub64(r4h, 0, t)
377
378 // commit the decrement if the subtraction underflowed (reciprocal was too large)
379 if b != 0 {
380 r4h, r4i, r4j, r4k, r4l = x4, x3, x2, x1, x0
381 }
382
383 // Shift to correct bit alignment, truncating excess bits
384
385 p = (p & 63) - 1
386
387 x0, c = bits.Add64(r4l, r4l, 0)
388 x1, c = bits.Add64(r4k, r4k, c)
389 x2, c = bits.Add64(r4j, r4j, c)
390 x3, c = bits.Add64(r4i, r4i, c)
391 x4, _ = bits.Add64(r4h, r4h, c)
392
393 if p < 0 {
394 r4h, r4i, r4j, r4k, r4l = x4, x3, x2, x1, x0
395 p = 0 // avoid negative shift below
396 }
397
398 {
399 r := uint(p) // right shift
400 l := uint(64 - r) // left shift
401
402 x0 = (r4l >> r) | (r4k << l)
403 x1 = (r4k >> r) | (r4j << l)
404 x2 = (r4j >> r) | (r4i << l)
405 x3 = (r4i >> r) | (r4h << l)
406 x4 = (r4h >> r)
407 }
408
409 if p > 0 {
410 r4h, r4i, r4j, r4k, r4l = x4, x3, x2, x1, x0
411 }
412
413 mu[0] = r4l
414 mu[1] = r4k
415 mu[2] = r4j
416 mu[3] = r4i
417 mu[4] = r4h
418
419 return mu
420}
421
422// reduce4 computes the least non-negative residue of x modulo m
423//
424// requires a four-word modulus (m[3] > 1) and its inverse (mu)
425func reduce4(x [8]uint64, m *Uint, mu [5]uint64) (z Uint) {
426 // NB: Most variable names in the comments match the pseudocode for
427 // Barrett reduction in the Handbook of Applied Cryptography.
428
429 // q1 = x/2^192
430
431 x0 := x[3]
432 x1 := x[4]
433 x2 := x[5]
434 x3 := x[6]
435 x4 := x[7]
436
437 // q2 = q1 * mu; q3 = q2 / 2^320
438
439 var q0, q1, q2, q3, q4, q5, t0, t1, c uint64
440
441 q0, _ = bits.Mul64(x3, mu[0])
442 q1, t0 = bits.Mul64(x4, mu[0])
443 q0, c = bits.Add64(q0, t0, 0)
444 q1, _ = bits.Add64(q1, 0, c)
445
446 t1, _ = bits.Mul64(x2, mu[1])
447 q0, c = bits.Add64(q0, t1, 0)
448 q2, t0 = bits.Mul64(x4, mu[1])
449 q1, c = bits.Add64(q1, t0, c)
450 q2, _ = bits.Add64(q2, 0, c)
451
452 t1, t0 = bits.Mul64(x3, mu[1])
453 q0, c = bits.Add64(q0, t0, 0)
454 q1, c = bits.Add64(q1, t1, c)
455 q2, _ = bits.Add64(q2, 0, c)
456
457 t1, t0 = bits.Mul64(x2, mu[2])
458 q0, c = bits.Add64(q0, t0, 0)
459 q1, c = bits.Add64(q1, t1, c)
460 q3, t0 = bits.Mul64(x4, mu[2])
461 q2, c = bits.Add64(q2, t0, c)
462 q3, _ = bits.Add64(q3, 0, c)
463
464 t1, _ = bits.Mul64(x1, mu[2])
465 q0, c = bits.Add64(q0, t1, 0)
466 t1, t0 = bits.Mul64(x3, mu[2])
467 q1, c = bits.Add64(q1, t0, c)
468 q2, c = bits.Add64(q2, t1, c)
469 q3, _ = bits.Add64(q3, 0, c)
470
471 t1, _ = bits.Mul64(x0, mu[3])
472 q0, c = bits.Add64(q0, t1, 0)
473 t1, t0 = bits.Mul64(x2, mu[3])
474 q1, c = bits.Add64(q1, t0, c)
475 q2, c = bits.Add64(q2, t1, c)
476 q4, t0 = bits.Mul64(x4, mu[3])
477 q3, c = bits.Add64(q3, t0, c)
478 q4, _ = bits.Add64(q4, 0, c)
479
480 t1, t0 = bits.Mul64(x1, mu[3])
481 q0, c = bits.Add64(q0, t0, 0)
482 q1, c = bits.Add64(q1, t1, c)
483 t1, t0 = bits.Mul64(x3, mu[3])
484 q2, c = bits.Add64(q2, t0, c)
485 q3, c = bits.Add64(q3, t1, c)
486 q4, _ = bits.Add64(q4, 0, c)
487
488 t1, t0 = bits.Mul64(x0, mu[4])
489 _, c = bits.Add64(q0, t0, 0)
490 q1, c = bits.Add64(q1, t1, c)
491 t1, t0 = bits.Mul64(x2, mu[4])
492 q2, c = bits.Add64(q2, t0, c)
493 q3, c = bits.Add64(q3, t1, c)
494 q5, t0 = bits.Mul64(x4, mu[4])
495 q4, c = bits.Add64(q4, t0, c)
496 q5, _ = bits.Add64(q5, 0, c)
497
498 t1, t0 = bits.Mul64(x1, mu[4])
499 q1, c = bits.Add64(q1, t0, 0)
500 q2, c = bits.Add64(q2, t1, c)
501 t1, t0 = bits.Mul64(x3, mu[4])
502 q3, c = bits.Add64(q3, t0, c)
503 q4, c = bits.Add64(q4, t1, c)
504 q5, _ = bits.Add64(q5, 0, c)
505
506 // Drop the fractional part of q3
507
508 q0 = q1
509 q1 = q2
510 q2 = q3
511 q3 = q4
512 q4 = q5
513
514 // r1 = x mod 2^320
515
516 x0 = x[0]
517 x1 = x[1]
518 x2 = x[2]
519 x3 = x[3]
520 x4 = x[4]
521
522 // r2 = q3 * m mod 2^320
523
524 var r0, r1, r2, r3, r4 uint64
525
526 r4, r3 = bits.Mul64(q0, m[3])
527 _, t0 = bits.Mul64(q1, m[3])
528 r4, _ = bits.Add64(r4, t0, 0)
529
530 t1, r2 = bits.Mul64(q0, m[2])
531 r3, c = bits.Add64(r3, t1, 0)
532 _, t0 = bits.Mul64(q2, m[2])
533 r4, _ = bits.Add64(r4, t0, c)
534
535 t1, t0 = bits.Mul64(q1, m[2])
536 r3, c = bits.Add64(r3, t0, 0)
537 r4, _ = bits.Add64(r4, t1, c)
538
539 t1, r1 = bits.Mul64(q0, m[1])
540 r2, c = bits.Add64(r2, t1, 0)
541 t1, t0 = bits.Mul64(q2, m[1])
542 r3, c = bits.Add64(r3, t0, c)
543 r4, _ = bits.Add64(r4, t1, c)
544
545 t1, t0 = bits.Mul64(q1, m[1])
546 r2, c = bits.Add64(r2, t0, 0)
547 r3, c = bits.Add64(r3, t1, c)
548 _, t0 = bits.Mul64(q3, m[1])
549 r4, _ = bits.Add64(r4, t0, c)
550
551 t1, r0 = bits.Mul64(q0, m[0])
552 r1, c = bits.Add64(r1, t1, 0)
553 t1, t0 = bits.Mul64(q2, m[0])
554 r2, c = bits.Add64(r2, t0, c)
555 r3, c = bits.Add64(r3, t1, c)
556 _, t0 = bits.Mul64(q4, m[0])
557 r4, _ = bits.Add64(r4, t0, c)
558
559 t1, t0 = bits.Mul64(q1, m[0])
560 r1, c = bits.Add64(r1, t0, 0)
561 r2, c = bits.Add64(r2, t1, c)
562 t1, t0 = bits.Mul64(q3, m[0])
563 r3, c = bits.Add64(r3, t0, c)
564 r4, _ = bits.Add64(r4, t1, c)
565
566 // r = r1 - r2
567
568 var b uint64
569
570 r0, b = bits.Sub64(x0, r0, 0)
571 r1, b = bits.Sub64(x1, r1, b)
572 r2, b = bits.Sub64(x2, r2, b)
573 r3, b = bits.Sub64(x3, r3, b)
574 r4, b = bits.Sub64(x4, r4, b)
575
576 // if r<0 then r+=m
577
578 if b != 0 {
579 r0, c = bits.Add64(r0, m[0], 0)
580 r1, c = bits.Add64(r1, m[1], c)
581 r2, c = bits.Add64(r2, m[2], c)
582 r3, c = bits.Add64(r3, m[3], c)
583 r4, _ = bits.Add64(r4, 0, c)
584 }
585
586 // while (r>=m) r-=m
587
588 for {
589 // q = r - m
590 q0, b = bits.Sub64(r0, m[0], 0)
591 q1, b = bits.Sub64(r1, m[1], b)
592 q2, b = bits.Sub64(r2, m[2], b)
593 q3, b = bits.Sub64(r3, m[3], b)
594 q4, b = bits.Sub64(r4, 0, b)
595
596 // if borrow break
597 if b != 0 {
598 break
599 }
600
601 // r = q
602 r4, r3, r2, r1, r0 = q4, q3, q2, q1, q0
603 }
604
605 z[3], z[2], z[1], z[0] = r3, r2, r1, r0
606
607 return z
608}
