bitcoinsecp256k1_modinv/modinv32.rs
1/*!
2 | This file implements modular inversion
3 | based on the paper "Fast constant-time
4 | gcd computation and modular inversion"
5 | by Daniel J. Bernstein and Bo-Yin Yang.
6 |
7 | For an explanation of the algorithm,
8 | see doc/safegcd_implementation.md.
9 | This file contains an implementation
10 | for N=30, using 30-bit signed limbs
11 | represented as int32_t.
12 |
13 */
14
15crate::ix!();
16
17
18
19//-------------------------------------------[.cpp/bitcoin/src/secp256k1/src/modinv32.h]
20
21/**
22 | A signed 30-bit limb representation
23 | of integers.
24 |
25 | Its value is sum(v[i] * 2^(30*i), i=0..8).
26 |
27 */
28pub struct ModInv32Signed30 {
29 v: [i32; 9],
30}
31
32pub struct ModInv32ModInfo {
33
34 /**
35 | The modulus in signed30 notation, must
36 | be odd and in [3, 2^256].
37 |
38 */
39 modulus: ModInv32Signed30,
40
41 /**
42 | modulus^{-1} mod 2^30
43 |
44 */
45 modulus_inv30: u32,
46}
47
48//-------------------------------------------[.cpp/bitcoin/src/secp256k1/src/modinv32_impl.h]
49
50#[cfg(VERIFY)]
51pub const SIGNED30_ONE: ModInv32Signed30 = {{1}};
52
53/**
54 | Compute a*factor and put it in r. All
55 | but the top limb in r will be in range [0,2^30).
56 |
57 */
58#[cfg(VERIFY)]
59pub fn modinv32_mul_30(
60 r: *mut ModInv32Signed30,
61 a: *const ModInv32Signed30,
62 alen: i32,
63 factor: i32) {
64
65 todo!();
66 /*
67 const int32_t M30 = (int32_t)(UINT32_MAX >> 2);
68 int64_t c = 0;
69 int i;
70 for (i = 0; i < 8; ++i) {
71 if (i < alen) c += (int64_t)a->v[i] * factor;
72 r->v[i] = (int32_t)c & M30; c >>= 30;
73 }
74 if (8 < alen) c += (int64_t)a->v[8] * factor;
75 VERIFY_CHECK(c == (int32_t)c);
76 r->v[8] = (int32_t)c;
77 */
78}
79
80/**
81 | Return -1 for a<b*factor, 0 for a==b*factor,
82 | 1 for a>b*factor. A consists of alen
83 | limbs; b has 9.
84 |
85 */
86#[cfg(VERIFY)]
87pub fn modinv32_mul_cmp_30(
88 a: *const ModInv32Signed30,
89 alen: i32,
90 b: *const ModInv32Signed30,
91 factor: i32) -> i32 {
92
93 todo!();
94 /*
95 int i;
96 modinv32_signed30 am, bm;
97 modinv32_mul_30(&am, a, alen, 1); /* Normalize all but the top limb of a. */
98 modinv32_mul_30(&bm, b, 9, factor);
99 for (i = 0; i < 8; ++i) {
100 /* Verify that all but the top limb of a and b are normalized. */
101 VERIFY_CHECK(am.v[i] >> 30 == 0);
102 VERIFY_CHECK(bm.v[i] >> 30 == 0);
103 }
104 for (i = 8; i >= 0; --i) {
105 if (am.v[i] < bm.v[i]) return -1;
106 if (am.v[i] > bm.v[i]) return 1;
107 }
108 return 0;
109 */
110}
111
112/**
113 | Take as input a signed30 number in range
114 | (-2*modulus,modulus), and add a multiple of the
115 | modulus to it to bring it to range [0,modulus).
116 |
117 | If sign < 0, the input will also be negated in
118 | the process.
119 |
120 | The input must have limbs in range
121 | (-2^30,2^30). The output will have limbs in
122 | range [0,2^30).
123 */
124pub fn modinv32_normalize_30(
125 r: *mut ModInv32Signed30,
126 sign: i32,
127 modinfo: *const ModInv32ModInfo) {
128
129 todo!();
130 /*
131 const int32_t M30 = (int32_t)(UINT32_MAX >> 2);
132 int32_t r0 = r->v[0], r1 = r->v[1], r2 = r->v[2], r3 = r->v[3], r4 = r->v[4],
133 r5 = r->v[5], r6 = r->v[6], r7 = r->v[7], r8 = r->v[8];
134 int32_t cond_add, cond_negate;
135
136 #ifdef VERIFY
137 /* Verify that all limbs are in range (-2^30,2^30). */
138 int i;
139 for (i = 0; i < 9; ++i) {
140 VERIFY_CHECK(r->v[i] >= -M30);
141 VERIFY_CHECK(r->v[i] <= M30);
142 }
143 VERIFY_CHECK(modinv32_mul_cmp_30(r, 9, &modinfo->modulus, -2) > 0); /* r > -2*modulus */
144 VERIFY_CHECK(modinv32_mul_cmp_30(r, 9, &modinfo->modulus, 1) < 0); /* r < modulus */
145 #endif
146
147 /* In a first step, add the modulus if the input is negative, and then negate if requested.
148 * This brings r from range (-2*modulus,modulus) to range (-modulus,modulus). As all input
149 * limbs are in range (-2^30,2^30), this cannot overflow an int32_t. Note that the right
150 * shifts below are signed sign-extending shifts (see assumptions.h for tests that that is
151 * indeed the behavior of the right shift operator). */
152 cond_add = r8 >> 31;
153 r0 += modinfo->modulus.v[0] & cond_add;
154 r1 += modinfo->modulus.v[1] & cond_add;
155 r2 += modinfo->modulus.v[2] & cond_add;
156 r3 += modinfo->modulus.v[3] & cond_add;
157 r4 += modinfo->modulus.v[4] & cond_add;
158 r5 += modinfo->modulus.v[5] & cond_add;
159 r6 += modinfo->modulus.v[6] & cond_add;
160 r7 += modinfo->modulus.v[7] & cond_add;
161 r8 += modinfo->modulus.v[8] & cond_add;
162 cond_negate = sign >> 31;
163 r0 = (r0 ^ cond_negate) - cond_negate;
164 r1 = (r1 ^ cond_negate) - cond_negate;
165 r2 = (r2 ^ cond_negate) - cond_negate;
166 r3 = (r3 ^ cond_negate) - cond_negate;
167 r4 = (r4 ^ cond_negate) - cond_negate;
168 r5 = (r5 ^ cond_negate) - cond_negate;
169 r6 = (r6 ^ cond_negate) - cond_negate;
170 r7 = (r7 ^ cond_negate) - cond_negate;
171 r8 = (r8 ^ cond_negate) - cond_negate;
172 /* Propagate the top bits, to bring limbs back to range (-2^30,2^30). */
173 r1 += r0 >> 30; r0 &= M30;
174 r2 += r1 >> 30; r1 &= M30;
175 r3 += r2 >> 30; r2 &= M30;
176 r4 += r3 >> 30; r3 &= M30;
177 r5 += r4 >> 30; r4 &= M30;
178 r6 += r5 >> 30; r5 &= M30;
179 r7 += r6 >> 30; r6 &= M30;
180 r8 += r7 >> 30; r7 &= M30;
181
182 /* In a second step add the modulus again if the result is still negative, bringing r to range
183 * [0,modulus). */
184 cond_add = r8 >> 31;
185 r0 += modinfo->modulus.v[0] & cond_add;
186 r1 += modinfo->modulus.v[1] & cond_add;
187 r2 += modinfo->modulus.v[2] & cond_add;
188 r3 += modinfo->modulus.v[3] & cond_add;
189 r4 += modinfo->modulus.v[4] & cond_add;
190 r5 += modinfo->modulus.v[5] & cond_add;
191 r6 += modinfo->modulus.v[6] & cond_add;
192 r7 += modinfo->modulus.v[7] & cond_add;
193 r8 += modinfo->modulus.v[8] & cond_add;
194 /* And propagate again. */
195 r1 += r0 >> 30; r0 &= M30;
196 r2 += r1 >> 30; r1 &= M30;
197 r3 += r2 >> 30; r2 &= M30;
198 r4 += r3 >> 30; r3 &= M30;
199 r5 += r4 >> 30; r4 &= M30;
200 r6 += r5 >> 30; r5 &= M30;
201 r7 += r6 >> 30; r6 &= M30;
202 r8 += r7 >> 30; r7 &= M30;
203
204 r->v[0] = r0;
205 r->v[1] = r1;
206 r->v[2] = r2;
207 r->v[3] = r3;
208 r->v[4] = r4;
209 r->v[5] = r5;
210 r->v[6] = r6;
211 r->v[7] = r7;
212 r->v[8] = r8;
213
214 #ifdef VERIFY
215 VERIFY_CHECK(r0 >> 30 == 0);
216 VERIFY_CHECK(r1 >> 30 == 0);
217 VERIFY_CHECK(r2 >> 30 == 0);
218 VERIFY_CHECK(r3 >> 30 == 0);
219 VERIFY_CHECK(r4 >> 30 == 0);
220 VERIFY_CHECK(r5 >> 30 == 0);
221 VERIFY_CHECK(r6 >> 30 == 0);
222 VERIFY_CHECK(r7 >> 30 == 0);
223 VERIFY_CHECK(r8 >> 30 == 0);
224 VERIFY_CHECK(modinv32_mul_cmp_30(r, 9, &modinfo->modulus, 0) >= 0); /* r >= 0 */
225 VERIFY_CHECK(modinv32_mul_cmp_30(r, 9, &modinfo->modulus, 1) < 0); /* r < modulus */
226 #endif
227 */
228}
229
230/**
231 | Data type for transition matrices (see
232 | section 3 of explanation).
233 |
234 | t = [ u v ]
235 | [ q r ]
236 |
237 */
238pub struct ModInv32Trans2x2 {
239 u: i32,
240 v: i32,
241 q: i32,
242 r: i32,
243}
244
245/**
246 | Compute the transition matrix and zeta for 30
247 | divsteps.
248 |
249 | Input: zeta: initial zeta
250 | f0: bottom limb of initial f
251 | g0: bottom limb of initial g
252 | Output: t: transition matrix
253 | Return: final zeta
254 |
255 | Implements the divsteps_n_matrix function from
256 | the explanation.
257 */
258pub fn modinv32_divsteps_30(
259 zeta: i32,
260 f0: u32,
261 g0: u32,
262 t: *mut ModInv32Trans2x2) -> i32 {
263
264 todo!();
265 /*
266 /* u,v,q,r are the elements of the transformation matrix being built up,
267 * starting with the identity matrix. Semantically they are signed integers
268 * in range [-2^30,2^30], but here represented as unsigned mod 2^32. This
269 * permits left shifting (which is UB for negative numbers). The range
270 * being inside [-2^31,2^31) means that casting to signed works correctly.
271 */
272 uint32_t u = 1, v = 0, q = 0, r = 1;
273 uint32_t c1, c2, f = f0, g = g0, x, y, z;
274 int i;
275
276 for (i = 0; i < 30; ++i) {
277 VERIFY_CHECK((f & 1) == 1); /* f must always be odd */
278 VERIFY_CHECK((u * f0 + v * g0) == f << i);
279 VERIFY_CHECK((q * f0 + r * g0) == g << i);
280 /* Compute conditional masks for (zeta < 0) and for (g & 1). */
281 c1 = zeta >> 31;
282 c2 = -(g & 1);
283 /* Compute x,y,z, conditionally negated versions of f,u,v. */
284 x = (f ^ c1) - c1;
285 y = (u ^ c1) - c1;
286 z = (v ^ c1) - c1;
287 /* Conditionally add x,y,z to g,q,r. */
288 g += x & c2;
289 q += y & c2;
290 r += z & c2;
291 /* In what follows, c1 is a condition mask for (zeta < 0) and (g & 1). */
292 c1 &= c2;
293 /* Conditionally change zeta into -zeta-2 or zeta-1. */
294 zeta = (zeta ^ c1) - 1;
295 /* Conditionally add g,q,r to f,u,v. */
296 f += g & c1;
297 u += q & c1;
298 v += r & c1;
299 /* Shifts */
300 g >>= 1;
301 u <<= 1;
302 v <<= 1;
303 /* Bounds on zeta that follow from the bounds on iteration count (max 20*30 divsteps). */
304 VERIFY_CHECK(zeta >= -601 && zeta <= 601);
305 }
306 /* Return data in t and return value. */
307 t->u = (int32_t)u;
308 t->v = (int32_t)v;
309 t->q = (int32_t)q;
310 t->r = (int32_t)r;
311 /* The determinant of t must be a power of two. This guarantees that multiplication with t
312 * does not change the gcd of f and g, apart from adding a power-of-2 factor to it (which
313 * will be divided out again). As each divstep's individual matrix has determinant 2, the
314 * aggregate of 30 of them will have determinant 2^30. */
315 VERIFY_CHECK((int64_t)t->u * t->r - (int64_t)t->v * t->q == ((int64_t)1) << 30);
316 return zeta;
317 */
318}
319
320/**
321 | Compute the transition matrix and eta for 30
322 | divsteps (variable time).
323 |
324 | Input: eta: initial eta
325 | f0: bottom limb of initial f
326 | g0: bottom limb of initial g
327 | Output: t: transition matrix
328 | Return: final eta
329 |
330 | Implements the divsteps_n_matrix_var function
331 | from the explanation.
332 */
333pub fn modinv32_divsteps_30_var(
334 eta: i32,
335 f0: u32,
336 g0: u32,
337 t: *mut ModInv32Trans2x2) -> i32 {
338
339 todo!();
340 /*
341 /* inv256[i] = -(2*i+1)^-1 (mod 256) */
342 static const uint8_t inv256[128] = {
343 0xFF, 0x55, 0x33, 0x49, 0xC7, 0x5D, 0x3B, 0x11, 0x0F, 0xE5, 0xC3, 0x59,
344 0xD7, 0xED, 0xCB, 0x21, 0x1F, 0x75, 0x53, 0x69, 0xE7, 0x7D, 0x5B, 0x31,
345 0x2F, 0x05, 0xE3, 0x79, 0xF7, 0x0D, 0xEB, 0x41, 0x3F, 0x95, 0x73, 0x89,
346 0x07, 0x9D, 0x7B, 0x51, 0x4F, 0x25, 0x03, 0x99, 0x17, 0x2D, 0x0B, 0x61,
347 0x5F, 0xB5, 0x93, 0xA9, 0x27, 0xBD, 0x9B, 0x71, 0x6F, 0x45, 0x23, 0xB9,
348 0x37, 0x4D, 0x2B, 0x81, 0x7F, 0xD5, 0xB3, 0xC9, 0x47, 0xDD, 0xBB, 0x91,
349 0x8F, 0x65, 0x43, 0xD9, 0x57, 0x6D, 0x4B, 0xA1, 0x9F, 0xF5, 0xD3, 0xE9,
350 0x67, 0xFD, 0xDB, 0xB1, 0xAF, 0x85, 0x63, 0xF9, 0x77, 0x8D, 0x6B, 0xC1,
351 0xBF, 0x15, 0xF3, 0x09, 0x87, 0x1D, 0xFB, 0xD1, 0xCF, 0xA5, 0x83, 0x19,
352 0x97, 0xAD, 0x8B, 0xE1, 0xDF, 0x35, 0x13, 0x29, 0xA7, 0x3D, 0x1B, 0xF1,
353 0xEF, 0xC5, 0xA3, 0x39, 0xB7, 0xCD, 0xAB, 0x01
354 };
355
356 /* Transformation matrix; see comments in modinv32_divsteps_30. */
357 uint32_t u = 1, v = 0, q = 0, r = 1;
358 uint32_t f = f0, g = g0, m;
359 uint16_t w;
360 int i = 30, limit, zeros;
361
362 for (;;) {
363 /* Use a sentinel bit to count zeros only up to i. */
364 zeros = ctz32_var(g | (UINT32_MAX << i));
365 /* Perform zeros divsteps at once; they all just divide g by two. */
366 g >>= zeros;
367 u <<= zeros;
368 v <<= zeros;
369 eta -= zeros;
370 i -= zeros;
371 /* We're done once we've done 30 divsteps. */
372 if (i == 0) break;
373 VERIFY_CHECK((f & 1) == 1);
374 VERIFY_CHECK((g & 1) == 1);
375 VERIFY_CHECK((u * f0 + v * g0) == f << (30 - i));
376 VERIFY_CHECK((q * f0 + r * g0) == g << (30 - i));
377 /* Bounds on eta that follow from the bounds on iteration count (max 25*30 divsteps). */
378 VERIFY_CHECK(eta >= -751 && eta <= 751);
379 /* If eta is negative, negate it and replace f,g with g,-f. */
380 if (eta < 0) {
381 uint32_t tmp;
382 eta = -eta;
383 tmp = f; f = g; g = -tmp;
384 tmp = u; u = q; q = -tmp;
385 tmp = v; v = r; r = -tmp;
386 }
387 /* eta is now >= 0. In what follows we're going to cancel out the bottom bits of g. No more
388 * than i can be cancelled out (as we'd be done before that point), and no more than eta+1
389 * can be done as its sign will flip once that happens. */
390 limit = ((int)eta + 1) > i ? i : ((int)eta + 1);
391 /* m is a mask for the bottom min(limit, 8) bits (our table only supports 8 bits). */
392 VERIFY_CHECK(limit > 0 && limit <= 30);
393 m = (UINT32_MAX >> (32 - limit)) & 255U;
394 /* Find what multiple of f must be added to g to cancel its bottom min(limit, 8) bits. */
395 w = (g * inv256[(f >> 1) & 127]) & m;
396 /* Do so. */
397 g += f * w;
398 q += u * w;
399 r += v * w;
400 VERIFY_CHECK((g & m) == 0);
401 }
402 /* Return data in t and return value. */
403 t->u = (int32_t)u;
404 t->v = (int32_t)v;
405 t->q = (int32_t)q;
406 t->r = (int32_t)r;
407 /* The determinant of t must be a power of two. This guarantees that multiplication with t
408 * does not change the gcd of f and g, apart from adding a power-of-2 factor to it (which
409 * will be divided out again). As each divstep's individual matrix has determinant 2, the
410 * aggregate of 30 of them will have determinant 2^30. */
411 VERIFY_CHECK((int64_t)t->u * t->r - (int64_t)t->v * t->q == ((int64_t)1) << 30);
412 return eta;
413 */
414}
415
416/**
417 | Compute (t/2^30) * [d, e] mod modulus,
418 | where t is a transition matrix for 30 divsteps.
419 |
420 | On input and output, d and e are in range
421 | (-2*modulus,modulus). All output limbs will be
422 | in range
423 |
424 | (-2^30,2^30).
425 |
426 | This implements the update_de function from the
427 | explanation.
428 */
429pub fn modinv32_update_de_30(
430 d: *mut ModInv32Signed30,
431 e: *mut ModInv32Signed30,
432 t: *const ModInv32Trans2x2,
433 modinfo: *const ModInv32ModInfo) {
434
435 todo!();
436 /*
437 const int32_t M30 = (int32_t)(UINT32_MAX >> 2);
438 const int32_t u = t->u, v = t->v, q = t->q, r = t->r;
439 int32_t di, ei, md, me, sd, se;
440 int64_t cd, ce;
441 int i;
442 #ifdef VERIFY
443 VERIFY_CHECK(modinv32_mul_cmp_30(d, 9, &modinfo->modulus, -2) > 0); /* d > -2*modulus */
444 VERIFY_CHECK(modinv32_mul_cmp_30(d, 9, &modinfo->modulus, 1) < 0); /* d < modulus */
445 VERIFY_CHECK(modinv32_mul_cmp_30(e, 9, &modinfo->modulus, -2) > 0); /* e > -2*modulus */
446 VERIFY_CHECK(modinv32_mul_cmp_30(e, 9, &modinfo->modulus, 1) < 0); /* e < modulus */
447 VERIFY_CHECK((labs(u) + labs(v)) >= 0); /* |u|+|v| doesn't overflow */
448 VERIFY_CHECK((labs(q) + labs(r)) >= 0); /* |q|+|r| doesn't overflow */
449 VERIFY_CHECK((labs(u) + labs(v)) <= M30 + 1); /* |u|+|v| <= 2^30 */
450 VERIFY_CHECK((labs(q) + labs(r)) <= M30 + 1); /* |q|+|r| <= 2^30 */
451 #endif
452 /* [md,me] start as zero; plus [u,q] if d is negative; plus [v,r] if e is negative. */
453 sd = d->v[8] >> 31;
454 se = e->v[8] >> 31;
455 md = (u & sd) + (v & se);
456 me = (q & sd) + (r & se);
457 /* Begin computing t*[d,e]. */
458 di = d->v[0];
459 ei = e->v[0];
460 cd = (int64_t)u * di + (int64_t)v * ei;
461 ce = (int64_t)q * di + (int64_t)r * ei;
462 /* Correct md,me so that t*[d,e]+modulus*[md,me] has 30 zero bottom bits. */
463 md -= (modinfo->modulus_inv30 * (uint32_t)cd + md) & M30;
464 me -= (modinfo->modulus_inv30 * (uint32_t)ce + me) & M30;
465 /* Update the beginning of computation for t*[d,e]+modulus*[md,me] now md,me are known. */
466 cd += (int64_t)modinfo->modulus.v[0] * md;
467 ce += (int64_t)modinfo->modulus.v[0] * me;
468 /* Verify that the low 30 bits of the computation are indeed zero, and then throw them away. */
469 VERIFY_CHECK(((int32_t)cd & M30) == 0); cd >>= 30;
470 VERIFY_CHECK(((int32_t)ce & M30) == 0); ce >>= 30;
471 /* Now iteratively compute limb i=1..8 of t*[d,e]+modulus*[md,me], and store them in output
472 * limb i-1 (shifting down by 30 bits). */
473 for (i = 1; i < 9; ++i) {
474 di = d->v[i];
475 ei = e->v[i];
476 cd += (int64_t)u * di + (int64_t)v * ei;
477 ce += (int64_t)q * di + (int64_t)r * ei;
478 cd += (int64_t)modinfo->modulus.v[i] * md;
479 ce += (int64_t)modinfo->modulus.v[i] * me;
480 d->v[i - 1] = (int32_t)cd & M30; cd >>= 30;
481 e->v[i - 1] = (int32_t)ce & M30; ce >>= 30;
482 }
483 /* What remains is limb 9 of t*[d,e]+modulus*[md,me]; store it as output limb 8. */
484 d->v[8] = (int32_t)cd;
485 e->v[8] = (int32_t)ce;
486 #ifdef VERIFY
487 VERIFY_CHECK(modinv32_mul_cmp_30(d, 9, &modinfo->modulus, -2) > 0); /* d > -2*modulus */
488 VERIFY_CHECK(modinv32_mul_cmp_30(d, 9, &modinfo->modulus, 1) < 0); /* d < modulus */
489 VERIFY_CHECK(modinv32_mul_cmp_30(e, 9, &modinfo->modulus, -2) > 0); /* e > -2*modulus */
490 VERIFY_CHECK(modinv32_mul_cmp_30(e, 9, &modinfo->modulus, 1) < 0); /* e < modulus */
491 #endif
492 */
493}
494
495/**
496 | Compute (t/2^30) * [f, g],
497 | where t is a transition matrix for 30 divsteps.
498 |
499 | This implements the update_fg function from the
500 | explanation.
501 */
502pub fn modinv32_update_fg_30(
503 f: *mut ModInv32Signed30,
504 g: *mut ModInv32Signed30,
505 t: *const ModInv32Trans2x2) {
506
507 todo!();
508 /*
509 const int32_t M30 = (int32_t)(UINT32_MAX >> 2);
510 const int32_t u = t->u, v = t->v, q = t->q, r = t->r;
511 int32_t fi, gi;
512 int64_t cf, cg;
513 int i;
514 /* Start computing t*[f,g]. */
515 fi = f->v[0];
516 gi = g->v[0];
517 cf = (int64_t)u * fi + (int64_t)v * gi;
518 cg = (int64_t)q * fi + (int64_t)r * gi;
519 /* Verify that the bottom 30 bits of the result are zero, and then throw them away. */
520 VERIFY_CHECK(((int32_t)cf & M30) == 0); cf >>= 30;
521 VERIFY_CHECK(((int32_t)cg & M30) == 0); cg >>= 30;
522 /* Now iteratively compute limb i=1..8 of t*[f,g], and store them in output limb i-1 (shifting
523 * down by 30 bits). */
524 for (i = 1; i < 9; ++i) {
525 fi = f->v[i];
526 gi = g->v[i];
527 cf += (int64_t)u * fi + (int64_t)v * gi;
528 cg += (int64_t)q * fi + (int64_t)r * gi;
529 f->v[i - 1] = (int32_t)cf & M30; cf >>= 30;
530 g->v[i - 1] = (int32_t)cg & M30; cg >>= 30;
531 }
532 /* What remains is limb 9 of t*[f,g]; store it as output limb 8. */
533 f->v[8] = (int32_t)cf;
534 g->v[8] = (int32_t)cg;
535 */
536}
537
538
539/**
540 | Compute (t/2^30) * [f, g],
541 | where t is a transition matrix for 30 divsteps.
542 |
543 | Version that operates on a variable number of
544 | limbs in f and g.
545 |
546 | This implements the update_fg function from the
547 | explanation in modinv64_impl.h.
548 */
549pub fn modinv32_update_fg_30_var(
550 len: i32,
551 f: *mut ModInv32Signed30,
552 g: *mut ModInv32Signed30,
553 t: *const ModInv32Trans2x2) {
554
555 todo!();
556 /*
557 const int32_t M30 = (int32_t)(UINT32_MAX >> 2);
558 const int32_t u = t->u, v = t->v, q = t->q, r = t->r;
559 int32_t fi, gi;
560 int64_t cf, cg;
561 int i;
562 VERIFY_CHECK(len > 0);
563 /* Start computing t*[f,g]. */
564 fi = f->v[0];
565 gi = g->v[0];
566 cf = (int64_t)u * fi + (int64_t)v * gi;
567 cg = (int64_t)q * fi + (int64_t)r * gi;
568 /* Verify that the bottom 62 bits of the result are zero, and then throw them away. */
569 VERIFY_CHECK(((int32_t)cf & M30) == 0); cf >>= 30;
570 VERIFY_CHECK(((int32_t)cg & M30) == 0); cg >>= 30;
571 /* Now iteratively compute limb i=1..len of t*[f,g], and store them in output limb i-1 (shifting
572 * down by 30 bits). */
573 for (i = 1; i < len; ++i) {
574 fi = f->v[i];
575 gi = g->v[i];
576 cf += (int64_t)u * fi + (int64_t)v * gi;
577 cg += (int64_t)q * fi + (int64_t)r * gi;
578 f->v[i - 1] = (int32_t)cf & M30; cf >>= 30;
579 g->v[i - 1] = (int32_t)cg & M30; cg >>= 30;
580 }
581 /* What remains is limb (len) of t*[f,g]; store it as output limb (len-1). */
582 f->v[len - 1] = (int32_t)cf;
583 g->v[len - 1] = (int32_t)cg;
584 */
585}
586
587/**
588 | Compute the inverse of x modulo modinfo->modulus,
589 | and replace x with it (constant time
590 | in x).
591 |
592 | Same as secp256k1_modinv32_var, but
593 | constant time in x (not in the modulus).
594 |
595 */
596pub fn modinv32(
597 x: *mut ModInv32Signed30,
598 modinfo: *const ModInv32ModInfo) {
599
600 todo!();
601 /*
602 /* Start with d=0, e=1, f=modulus, g=x, zeta=-1. */
603 modinv32_signed30 d = {{0}};
604 modinv32_signed30 e = {{1}};
605 modinv32_signed30 f = modinfo->modulus;
606 modinv32_signed30 g = *x;
607 int i;
608 int32_t zeta = -1; /* zeta = -(delta+1/2); delta is initially 1/2. */
609
610 /* Do 20 iterations of 30 divsteps each = 600 divsteps. 590 suffices for 256-bit inputs. */
611 for (i = 0; i < 20; ++i) {
612 /* Compute transition matrix and new zeta after 30 divsteps. */
613 modinv32_trans2x2 t;
614 zeta = modinv32_divsteps_30(zeta, f.v[0], g.v[0], &t);
615 /* Update d,e using that transition matrix. */
616 modinv32_update_de_30(&d, &e, &t, modinfo);
617 /* Update f,g using that transition matrix. */
618 #ifdef VERIFY
619 VERIFY_CHECK(modinv32_mul_cmp_30(&f, 9, &modinfo->modulus, -1) > 0); /* f > -modulus */
620 VERIFY_CHECK(modinv32_mul_cmp_30(&f, 9, &modinfo->modulus, 1) <= 0); /* f <= modulus */
621 VERIFY_CHECK(modinv32_mul_cmp_30(&g, 9, &modinfo->modulus, -1) > 0); /* g > -modulus */
622 VERIFY_CHECK(modinv32_mul_cmp_30(&g, 9, &modinfo->modulus, 1) < 0); /* g < modulus */
623 #endif
624 modinv32_update_fg_30(&f, &g, &t);
625 #ifdef VERIFY
626 VERIFY_CHECK(modinv32_mul_cmp_30(&f, 9, &modinfo->modulus, -1) > 0); /* f > -modulus */
627 VERIFY_CHECK(modinv32_mul_cmp_30(&f, 9, &modinfo->modulus, 1) <= 0); /* f <= modulus */
628 VERIFY_CHECK(modinv32_mul_cmp_30(&g, 9, &modinfo->modulus, -1) > 0); /* g > -modulus */
629 VERIFY_CHECK(modinv32_mul_cmp_30(&g, 9, &modinfo->modulus, 1) < 0); /* g < modulus */
630 #endif
631 }
632
633 /* At this point sufficient iterations have been performed that g must have reached 0
634 * and (if g was not originally 0) f must now equal +/- GCD of the initial f, g
635 * values i.e. +/- 1, and d now contains +/- the modular inverse. */
636 #ifdef VERIFY
637 /* g == 0 */
638 VERIFY_CHECK(modinv32_mul_cmp_30(&g, 9, &SIGNED30_ONE, 0) == 0);
639 /* |f| == 1, or (x == 0 and d == 0 and |f|=modulus) */
640 VERIFY_CHECK(modinv32_mul_cmp_30(&f, 9, &SIGNED30_ONE, -1) == 0 ||
641 modinv32_mul_cmp_30(&f, 9, &SIGNED30_ONE, 1) == 0 ||
642 (modinv32_mul_cmp_30(x, 9, &SIGNED30_ONE, 0) == 0 &&
643 modinv32_mul_cmp_30(&d, 9, &SIGNED30_ONE, 0) == 0 &&
644 (modinv32_mul_cmp_30(&f, 9, &modinfo->modulus, 1) == 0 ||
645 modinv32_mul_cmp_30(&f, 9, &modinfo->modulus, -1) == 0)));
646 #endif
647
648 /* Optionally negate d, normalize to [0,modulus), and return it. */
649 modinv32_normalize_30(&d, f.v[8], modinfo);
650 *x = d;
651 */
652}
653
654/**
655 | Compute the inverse of x modulo modinfo->modulus,
656 | and replace x with it (variable time).
657 |
658 -------------------------
659 | Replace x with its modular inverse mod
660 | modinfo->modulus. x must be in range
661 | [0, modulus).
662 |
663 | If x is zero, the result will be zero as
664 | well. If not, the inverse must exist
665 | (i.e., the gcd of x and modulus must be
666 | 1). These rules are automatically satisfied
667 | if the modulus is prime.
668 |
669 | On output, all of x's limbs will be in
670 | [0, 2^30).
671 |
672 */
673pub fn modinv32_var(
674 x: *mut ModInv32Signed30,
675 modinfo: *const ModInv32ModInfo) {
676
677 todo!();
678 /*
679 /* Start with d=0, e=1, f=modulus, g=x, eta=-1. */
680 modinv32_signed30 d = {{0, 0, 0, 0, 0, 0, 0, 0, 0}};
681 modinv32_signed30 e = {{1, 0, 0, 0, 0, 0, 0, 0, 0}};
682 modinv32_signed30 f = modinfo->modulus;
683 modinv32_signed30 g = *x;
684 #ifdef VERIFY
685 int i = 0;
686 #endif
687 int j, len = 9;
688 int32_t eta = -1; /* eta = -delta; delta is initially 1 (faster for the variable-time code) */
689 int32_t cond, fn, gn;
690
691 /* Do iterations of 30 divsteps each until g=0. */
692 while (1) {
693 /* Compute transition matrix and new eta after 30 divsteps. */
694 modinv32_trans2x2 t;
695 eta = modinv32_divsteps_30_var(eta, f.v[0], g.v[0], &t);
696 /* Update d,e using that transition matrix. */
697 modinv32_update_de_30(&d, &e, &t, modinfo);
698 /* Update f,g using that transition matrix. */
699 #ifdef VERIFY
700 VERIFY_CHECK(modinv32_mul_cmp_30(&f, len, &modinfo->modulus, -1) > 0); /* f > -modulus */
701 VERIFY_CHECK(modinv32_mul_cmp_30(&f, len, &modinfo->modulus, 1) <= 0); /* f <= modulus */
702 VERIFY_CHECK(modinv32_mul_cmp_30(&g, len, &modinfo->modulus, -1) > 0); /* g > -modulus */
703 VERIFY_CHECK(modinv32_mul_cmp_30(&g, len, &modinfo->modulus, 1) < 0); /* g < modulus */
704 #endif
705 modinv32_update_fg_30_var(len, &f, &g, &t);
706 /* If the bottom limb of g is 0, there is a chance g=0. */
707 if (g.v[0] == 0) {
708 cond = 0;
709 /* Check if all other limbs are also 0. */
710 for (j = 1; j < len; ++j) {
711 cond |= g.v[j];
712 }
713 /* If so, we're done. */
714 if (cond == 0) break;
715 }
716
717 /* Determine if len>1 and limb (len-1) of both f and g is 0 or -1. */
718 fn = f.v[len - 1];
719 gn = g.v[len - 1];
720 cond = ((int32_t)len - 2) >> 31;
721 cond |= fn ^ (fn >> 31);
722 cond |= gn ^ (gn >> 31);
723 /* If so, reduce length, propagating the sign of f and g's top limb into the one below. */
724 if (cond == 0) {
725 f.v[len - 2] |= (uint32_t)fn << 30;
726 g.v[len - 2] |= (uint32_t)gn << 30;
727 --len;
728 }
729 #ifdef VERIFY
730 VERIFY_CHECK(++i < 25); /* We should never need more than 25*30 = 750 divsteps */
731 VERIFY_CHECK(modinv32_mul_cmp_30(&f, len, &modinfo->modulus, -1) > 0); /* f > -modulus */
732 VERIFY_CHECK(modinv32_mul_cmp_30(&f, len, &modinfo->modulus, 1) <= 0); /* f <= modulus */
733 VERIFY_CHECK(modinv32_mul_cmp_30(&g, len, &modinfo->modulus, -1) > 0); /* g > -modulus */
734 VERIFY_CHECK(modinv32_mul_cmp_30(&g, len, &modinfo->modulus, 1) < 0); /* g < modulus */
735 #endif
736 }
737
738 /* At this point g is 0 and (if g was not originally 0) f must now equal +/- GCD of
739 * the initial f, g values i.e. +/- 1, and d now contains +/- the modular inverse. */
740 #ifdef VERIFY
741 /* g == 0 */
742 VERIFY_CHECK(modinv32_mul_cmp_30(&g, len, &SIGNED30_ONE, 0) == 0);
743 /* |f| == 1, or (x == 0 and d == 0 and |f|=modulus) */
744 VERIFY_CHECK(modinv32_mul_cmp_30(&f, len, &SIGNED30_ONE, -1) == 0 ||
745 modinv32_mul_cmp_30(&f, len, &SIGNED30_ONE, 1) == 0 ||
746 (modinv32_mul_cmp_30(x, 9, &SIGNED30_ONE, 0) == 0 &&
747 modinv32_mul_cmp_30(&d, 9, &SIGNED30_ONE, 0) == 0 &&
748 (modinv32_mul_cmp_30(&f, len, &modinfo->modulus, 1) == 0 ||
749 modinv32_mul_cmp_30(&f, len, &modinfo->modulus, -1) == 0)));
750 #endif
751
752 /* Optionally negate d, normalize to [0,modulus), and return it. */
753 modinv32_normalize_30(&d, f.v[len - 1], modinfo);
754 *x = d;
755 */
756}