Inside Bitcoin

251    limb_t u = 1, v = 0, q = 0, r = 1;
252    // The number of divsteps still left.
253    int i = SIGNED_LIMB_SIZE;
254    while (true) {
255        /* Use a sentinel bit to count zeros only up to i. */
256        int zeros = std::countr_zero(g | (MAX_LIMB << i));
257        /* Perform zeros divsteps at once; they all just divide g by two. */
258        g >>= zeros;
259        u <<= zeros;
260        v <<= zeros;
261        eta -= zeros;
262        i -= zeros;
263         /* We're done once we've performed SIGNED_LIMB_SIZE divsteps. */
264        if (i == 0) break;
265        /* If eta is negative, negate it and replace f,g with g,-f. */
266        if (eta < 0) {
267            limb_t tmp;
268            eta = -eta;
269            tmp = f; f = g; g = -tmp;
270            tmp = u; u = q; q = -tmp;
271            tmp = v; v = r; r = -tmp;
272        }
273        /* eta is now >= 0. In what follows we're going to cancel out the bottom bits of g. No more
274         * than i can be cancelled out (as we'd be done before that point), and no more than eta+1
275         * can be done as its sign will flip once that happens. */
276        int limit = ((int)eta + 1) > i ? i : ((int)eta + 1);
277        /* m is a mask for the bottom min(limit, 8) bits (our table only supports 8 bits). */
278        limb_t m = (MAX_LIMB >> (LIMB_SIZE - limit)) & 255U;
279        /* Find what multiple of f must be added to g to cancel its bottom min(limit, 8) bits. */
280        limb_t w = (g * NEGINV256[(f >> 1) & 127]) & m;
281        /* Do so. */
282        g += f * w;
283        q += u * w;
284        r += v * w;
285    }
286    out.u = (signed_limb_t)u;
287    out.v = (signed_limb_t)v;
288    out.q = (signed_limb_t)q;
289    out.r = (signed_limb_t)r;
290    return eta;
291}
292
293/** Apply matrix t/2^SIGNED_LIMB_SIZE to vector [d,e], modulo modulus.
294 *
295 * On input and output, d and e are in range 1-2*modulus..modulus-1.
296 */
297inline void UpdateDE(Num3072Signed& d, Num3072Signed& e, const SignedMatrix& t)
298{
299    const signed_limb_t u = t.u, v=t.v, q=t.q, r=t.r;
300
301    /* [md,me] start as zero; plus [u,q] if d is negative; plus [v,r] if e is negative. */
302    signed_limb_t sd = d.limbs[SIGNED_LIMBS - 1] >> (LIMB_SIZE - 1);
303    signed_limb_t se = e.limbs[SIGNED_LIMBS - 1] >> (LIMB_SIZE - 1);
304    signed_limb_t md = (u & sd) + (v & se);
305    signed_limb_t me = (q & sd) + (r & se);
306    /* Begin computing t*[d,e]. */
307    signed_limb_t di = d.limbs[0], ei = e.limbs[0];
308    signed_double_limb_t cd = (signed_double_limb_t)u * di + (signed_double_limb_t)v * ei;
309    signed_double_limb_t ce = (signed_double_limb_t)q * di + (signed_double_limb_t)r * ei;
310    /* Correct md,me so that t*[d,e]+modulus*[md,me] has SIGNED_LIMB_SIZE zero bottom bits. */
311    md -= (MODULUS_INVERSE * limb_t(cd) + md) & MAX_SIGNED_LIMB;
312    me -= (MODULUS_INVERSE * limb_t(ce) + me) & MAX_SIGNED_LIMB;
313    /* Update the beginning of computation for t*[d,e]+modulus*[md,me] now md,me are known. */
314    cd -= (signed_double_limb_t)1103717 * md;
315    ce -= (signed_double_limb_t)1103717 * me;
316    /* Verify that the low SIGNED_LIMB_SIZE bits of the computation are indeed zero, and then throw them away. */
317    Assume((cd & MAX_SIGNED_LIMB) == 0);
318    Assume((ce & MAX_SIGNED_LIMB) == 0);
319    cd >>= SIGNED_LIMB_SIZE;
320    ce >>= SIGNED_LIMB_SIZE;
321    /* Now iteratively compute limb i=1..SIGNED_LIMBS-2 of t*[d,e]+modulus*[md,me], and store them in output
322     * limb i-1 (shifting down by SIGNED_LIMB_SIZE bits). The corresponding limbs in modulus are all zero,
323     * so modulus/md/me are not actually involved here. */
324    for (int i = 1; i < SIGNED_LIMBS - 1; ++i) {
325        di = d.limbs[i];
326        ei = e.limbs[i];
327        cd += (signed_double_limb_t)u * di + (signed_double_limb_t)v * ei;
328        ce += (signed_double_limb_t)q * di + (signed_double_limb_t)r * ei;
329        d.limbs[i - 1] = (signed_limb_t)cd & MAX_SIGNED_LIMB; cd >>= SIGNED_LIMB_SIZE;
330        e.limbs[i - 1] = (signed_limb_t)ce & MAX_SIGNED_LIMB; ce >>= SIGNED_LIMB_SIZE;
331    }
332    /* Compute limb SIGNED_LIMBS-1 of t*[d,e]+modulus*[md,me], and store it in output limb SIGNED_LIMBS-2. */
333    di = d.limbs[SIGNED_LIMBS - 1];
334    ei = e.limbs[SIGNED_LIMBS - 1];
335    cd += (signed_double_limb_t)u * di + (signed_double_limb_t)v * ei;
336    ce += (signed_double_limb_t)q * di + (signed_double_limb_t)r * ei;
337    cd += (signed_double_limb_t)md << FINAL_LIMB_MODULUS_BITS;
338    ce += (signed_double_limb_t)me << FINAL_LIMB_MODULUS_BITS;
339    d.limbs[SIGNED_LIMBS - 2] = (signed_limb_t)cd & MAX_SIGNED_LIMB; cd >>= SIGNED_LIMB_SIZE;
340    e.limbs[SIGNED_LIMBS - 2] = (signed_limb_t)ce & MAX_SIGNED_LIMB; ce >>= SIGNED_LIMB_SIZE;
341    /* What remains goes into output limb SINGED_LIMBS-1 */
342    d.limbs[SIGNED_LIMBS - 1] = (signed_limb_t)cd;
343    e.limbs[SIGNED_LIMBS - 1] = (signed_limb_t)ce;
344}
345
346/** Apply matrix t/2^SIGNED_LIMB_SIZE to vector (f,g).
347 *
348 * The matrix t must be chosen such that t*(f,g) results in multiples of 2^SIGNED_LIMB_SIZE.
349 * This is the case for matrices computed by ComputeDivstepMatrix().
350 */
351inline void UpdateFG(Num3072Signed& f, Num3072Signed& g, const SignedMatrix& t, int len)
352{
353    const signed_limb_t u = t.u, v=t.v, q=t.q, r=t.r;
354
355    signed_limb_t fi, gi;
356    signed_double_limb_t cf, cg;
357    /* Start computing t*[f,g]. */
358    fi = f.limbs[0];
359    gi = g.limbs[0];
360    cf = (signed_double_limb_t)u * fi + (signed_double_limb_t)v * gi;
361    cg = (signed_double_limb_t)q * fi + (signed_double_limb_t)r * gi;
362    /* Verify that the bottom SIGNED_LIMB_BITS bits of the result are zero, and then throw them away. */
363    Assume((cf & MAX_SIGNED_LIMB) == 0);
364    Assume((cg & MAX_SIGNED_LIMB) == 0);
365    cf >>= SIGNED_LIMB_SIZE;
366    cg >>= SIGNED_LIMB_SIZE;
367    /* Now iteratively compute limb i=1..SIGNED_LIMBS-1 of t*[f,g], and store them in output limb i-1 (shifting
368     * down by SIGNED_LIMB_BITS bits). */
369    for (int i = 1; i < len; ++i) {
370        fi = f.limbs[i];
371        gi = g.limbs[i];
372        cf += (signed_double_limb_t)u * fi + (signed_double_limb_t)v * gi;
373        cg += (signed_double_limb_t)q * fi + (signed_double_limb_t)r * gi;
374        f.limbs[i - 1] = (signed_limb_t)cf & MAX_SIGNED_LIMB; cf >>= SIGNED_LIMB_SIZE;
375        g.limbs[i - 1] = (signed_limb_t)cg & MAX_SIGNED_LIMB; cg >>= SIGNED_LIMB_SIZE;
376    }
377    /* What remains is limb SIGNED_LIMBS of t*[f,g]; store it as output limb SIGNED_LIMBS-1. */
378    f.limbs[len - 1] = (signed_limb_t)cf;
379    g.limbs[len - 1] = (signed_limb_t)cg;
380
381}
382} // namespace
383
384Num3072 Num3072::GetInverse() const
385{
386    // Compute a modular inverse based on a variant of the safegcd algorithm:
387    // - Paper: https://gcd.cr.yp.to/papers.html
388    // - Inspired by this code in libsecp256k1:
389    //   https://github.com/bitcoin-core/secp256k1/blob/master/src/modinv32_impl.h
390    // - Explanation of the algorithm:
391    //   https://github.com/bitcoin-core/secp256k1/blob/master/doc/safegcd_implementation.md
392
393    // Local variables d, e, f, g:
394    // - f and g are the variables whose gcd we compute (despite knowing the answer is 1):
395    //   - f is always odd, and initialized as modulus
396    //   - g is initialized as *this (called x in what follows)
397    // - d and e are the numbers for which at every step it is the case that:
398    //   - f = d * x mod modulus; d is initialized as 0
399    //   - g = e * x mod modulus; e is initialized as 1
400    Num3072Signed d, e, f, g;
401    e.limbs[0] = 1;
402    // F is initialized as modulus, which in signed limb representation can be expressed
403    // simply as 2^3072 + -MAX_PRIME_DIFF.
404    f.limbs[0] = -MAX_PRIME_DIFF;
405    f.limbs[FINAL_LIMB_POSITION] = ((limb_t)1) << FINAL_LIMB_MODULUS_BITS;
406    g.FromNum3072(*this);
407    int len = SIGNED_LIMBS; //!< The number of significant limbs in f and g
408    signed_limb_t eta = -1; //!< State to track knowledge about ratio of f and g
409    // Perform divsteps on [f,g] until g=0 is reached, keeping (d,e) synchronized with them.
410    while (true) {
411        // Compute transformation matrix t that represents the next SIGNED_LIMB_SIZE divsteps
412        // to apply. This can be computed from just the bottom limb of f and g, and eta.
413        SignedMatrix t;
414        eta = ComputeDivstepMatrix(eta, f.limbs[0], g.limbs[0], t);
415        // Apply that transformation matrix to the full [f,g] vector.
416        UpdateFG(f, g, t, len);
417        // Apply that transformation matrix to the full [d,e] vector (mod modulus).
418        UpdateDE(d, e, t);
419
420        // Check if g is zero.
421        if (g.limbs[0] == 0) {
422            signed_limb_t cond = 0;
423            for (int j = 1; j < len; ++j) {
424                cond |= g.limbs[j];
425            }
426            // If so, we're done.
427            if (cond == 0) break;
428        }
429
430        // Check if the top limbs of both f and g are both 0 or -1.
431        signed_limb_t fn = f.limbs[len - 1], gn = g.limbs[len - 1];
432        signed_limb_t cond = ((signed_limb_t)len - 2) >> (LIMB_SIZE - 1);
433        cond |= fn ^ (fn >> (LIMB_SIZE - 1));
434        cond |= gn ^ (gn >> (LIMB_SIZE - 1));
435        if (cond == 0) {
436            // If so, drop the top limb, shrinking the size of f and g, by
437            // propagating the sign to the previous limb.
438            f.limbs[len - 2] |= (limb_t)f.limbs[len - 1] << SIGNED_LIMB_SIZE;
439            g.limbs[len - 2] |= (limb_t)g.limbs[len - 1] << SIGNED_LIMB_SIZE;
440            --len;
441        }
442    }
443    // At some point, [f,g] will have been rewritten into [f',0], such that gcd(f,g) = gcd(f',0).
444    // This is proven in the paper. As f started out being modulus, a prime number, we know that
445    // gcd is 1, and thus f' is 1 or -1.
446    Assume((f.limbs[0] & MAX_SIGNED_LIMB) == 1 || (f.limbs[0] & MAX_SIGNED_LIMB) == MAX_SIGNED_LIMB);
447    // As we've maintained the invariant that f = d * x mod modulus, we get d/f mod modulus is the
448    // modular inverse of x we're looking for. As f is 1 or -1, it is also true that d/f = d*f.
449    // Normalize d to prepare it for output, while negating it if f is negative.
450    d.Normalize(f.limbs[len - 1] >> (LIMB_SIZE  - 1));
451    Num3072 ret;
452    d.ToNum3072(ret);
453    return ret;
454}
455
456void Num3072::Multiply(const Num3072& a)
457{
458    limb_t c0 = 0, c1 = 0, c2 = 0;
459    Num3072 tmp;
460
461    /* Compute limbs 0..N-2 of this*a into tmp, including one reduction. */
462    for (int j = 0; j < LIMBS - 1; ++j) {
463        limb_t d0 = 0, d1 = 0, d2 = 0;
464        mul(d0, d1, this->limbs[1 + j], a.limbs[LIMBS + j - (1 + j)]);
465        for (int i = 2 + j; i < LIMBS; ++i) muladd3(d0, d1, d2, this->limbs[i], a.limbs[LIMBS + j - i]);
466        mulnadd3(c0, c1, c2, d0, d1, d2, MAX_PRIME_DIFF);
467        for (int i = 0; i < j + 1; ++i) muladd3(c0, c1, c2, this->limbs[i], a.limbs[j - i]);
468        extract3(c0, c1, c2, tmp.limbs[j]);
469    }
470
471    /* Compute limb N-1 of a*b into tmp. */
472    assert(c2 == 0);
473    for (int i = 0; i < LIMBS; ++i) muladd3(c0, c1, c2, this->limbs[i], a.limbs[LIMBS - 1 - i]);
474    extract3(c0, c1, c2, tmp.limbs[LIMBS - 1]);
475
476    /* Perform a second reduction. */
477    muln2(c0, c1, MAX_PRIME_DIFF);
478    for (int j = 0; j < LIMBS; ++j) {
479        addnextract2(c0, c1, tmp.limbs[j], this->limbs[j]);
480    }
481
482    assert(c1 == 0);
483    assert(c0 == 0 || c0 == 1);
484
485    /* Perform up to two more reductions if the internal state has already
486     * overflown the MAX of Num3072 or if it is larger than the modulus or
487     * if both are the case.
488     * */
489    if (this->IsOverflow()) this->FullReduce();
490    if (c0) this->FullReduce();
491}
492
493void Num3072::SetToOne()
494{
495    this->limbs[0] = 1;
496    for (int i = 1; i < LIMBS; ++i) this->limbs[i] = 0;
497}
498
499void Num3072::Divide(const Num3072& a)
500{