Inside Bitcoin

501*η=-δ*. In this representation, negating *δ* corresponds to negating *η*, and incrementing
502*δ* corresponds to decrementing *η*. This allows us to remove the negation in the *c1*
503computation:
504
505```python
506    # Compute a mask c1 for eta < 0, and compute the conditional negation x of f:
507    c1 = eta >> 63
508    x = (f ^ c1) - c1
509    # Compute a mask c2 for odd g, and conditionally add x to g:
510    c2 = -(g & 1)
511    g += x & c2
512    # Compute a mask c for (eta < 0) and odd (input) g, and use it to conditionally negate eta,
513    # and add g to f:
514    c3 = c1 & c2
515    eta = (eta ^ c3) - c3
516    f += g & c3
517    # Incrementing delta corresponds to decrementing eta.
518    eta -= 1
519    g >>= 1
520```
521
522A variant of divsteps with better worst-case performance can be used instead: starting *&delta;* at
523*1/2* instead of *1*. This reduces the worst case number of iterations to *590* for *256*-bit inputs
524(which can be shown using convex hull analysis). In this case, the substitution *&zeta;=-(&delta;+1/2)*
525is used instead to keep the variable integral. Incrementing *&delta;* by *1* still translates to
526decrementing *&zeta;* by *1*, but negating *&delta;* now corresponds to going from *&zeta;* to *-(&zeta;+1)*, or
527*~&zeta;*. Doing that conditionally based on *c3* is simply:
528
529```python
530    ...
531    c3 = c1 & c2
532    zeta ^= c3
533    ...
534```
535
536By replacing the loop in `divsteps_n_matrix` with a variant of the divstep code above (extended to
537also apply all *f* operations to *u*, *v* and all *g* operations to *q*, *r*), a constant-time version of
538`divsteps_n_matrix` is obtained. The full code will be in section 7.
539
540These bit fiddling tricks can also be used to make the conditional negations and additions in
541`update_de` and `normalize` constant-time.
542
543
544## 6. Variable-time optimizations
545
546In section 5, we modified the `divsteps_n_matrix` function (and a few others) to be constant time.
547Constant time operations are only necessary when computing modular inverses of secret data. In
548other cases, it slows down calculations unnecessarily. In this section, we will construct a
549faster non-constant time `divsteps_n_matrix` function.
550
551To do so, first consider yet another way of writing the inner loop of divstep operations in
552`gcd` from section 1. This decomposition is also explained in the paper in section 8.2. We use
553the original version with initial *&delta;=1* and *&eta;=-&delta;* here.
554
555```python
556for _ in range(N):
557    if g & 1 and eta < 0:
558        eta, f, g = -eta, g, -f
559    if g & 1:
560        g += f
561    eta -= 1
562    g >>= 1
563```
564
565Whenever *g* is even, the loop only shifts *g* down and decreases *&eta;*. When *g* ends in multiple zero
566bits, these iterations can be consolidated into one step. This requires counting the bottom zero
567bits efficiently, which is possible on most platforms; it is abstracted here as the function
568`count_trailing_zeros`.
569
570```python
571def count_trailing_zeros(v):
572    """
573    When v is zero, consider all N zero bits as "trailing".
574    For a non-zero value v, find z such that v=(d<<z) for some odd d.
575    """
576    if v == 0:
577        return N
578    else:
579        return (v & -v).bit_length() - 1
580
581i = N # divsteps left to do
582while True:
583    # Get rid of all bottom zeros at once. In the first iteration, g may be odd and the following
584    # lines have no effect (until "if eta < 0").
585    zeros = min(i, count_trailing_zeros(g))
586    eta -= zeros
587    g >>= zeros
588    i -= zeros
589    if i == 0:
590        break
591    # We know g is odd now
592    if eta < 0:
593        eta, f, g = -eta, g, -f
594    g += f
595    # g is even now, and the eta decrement and g shift will happen in the next loop.
596```
597
598We can now remove multiple bottom *0* bits from *g* at once, but still need a full iteration whenever
599there is a bottom *1* bit. In what follows, we will get rid of multiple *1* bits simultaneously as
600well.
601
602Observe that as long as *&eta; &geq; 0*, the loop does not modify *f*. Instead, it cancels out bottom
603bits of *g* and shifts them out, and decreases *&eta;* and *i* accordingly - interrupting only when *&eta;*
604becomes negative, or when *i* reaches *0*. Combined, this is equivalent to adding a multiple of *f* to
605*g* to cancel out multiple bottom bits, and then shifting them out.
606
607It is easy to find what that multiple is: we want a number *w* such that *g+w&thinsp;f* has a few bottom
608zero bits. If that number of bits is *L*, we want *g+w&thinsp;f mod 2<sup>L</sup> = 0*, or *w = -g/f mod 2<sup>L</sup>*. Since *f*
609is odd, such a *w* exists for any *L*. *L* cannot be more than *i* steps (as we'd finish the loop before
610doing more) or more than *&eta;+1* steps (as we'd run `eta, f, g = -eta, g, -f` at that point), but
611apart from that, we're only limited by the complexity of computing *w*.
612
613This code demonstrates how to cancel up to 4 bits per step:
614
615```python
616NEGINV16 = [15, 5, 3, 9, 7, 13, 11, 1] # NEGINV16[n//2] = (-n)^-1 mod 16, for odd n
617i = N
618while True:
619    zeros = min(i, count_trailing_zeros(g))
620    eta -= zeros
621    g >>= zeros
622    i -= zeros
623    if i == 0:
624        break
625    # We know g is odd now
626    if eta < 0:
627        eta, f, g = -eta, g, -f
628    # Compute limit on number of bits to cancel
629    limit = min(min(eta + 1, i), 4)
630    # Compute w = -g/f mod 2**limit, using the table value for -1/f mod 2**4. Note that f is
631    # always odd, so its inverse modulo a power of two always exists.
632    w = (g * NEGINV16[(f & 15) // 2]) % (2**limit)
633    # As w = -g/f mod (2**limit), g+w*f mod 2**limit = 0 mod 2**limit.
634    g += w * f
635    assert g % (2**limit) == 0
636    # The next iteration will now shift out at least limit bottom zero bits from g.
637```
638
639By using a bigger table more bits can be cancelled at once. The table can also be implemented
640as a formula. Several formulas are known for computing modular inverses modulo powers of two;
641some can be found in Hacker's Delight second edition by Henry S. Warren, Jr. pages 245-247.
642Here we need the negated modular inverse, which is a simple transformation of those:
643
644- Instead of a 3-bit table:
645  - *-f* or *f ^ 6*
646- Instead of a 4-bit table:
647  - *1 - f(f + 1)*
648  - *-(f + (((f + 1) & 4) << 1))*
649- For larger tables the following technique can be used: if *w=-1/f mod 2<sup>L</sup>*, then *w(w&thinsp;f+2)* is
650  *-1/f mod 2<sup>2L</sup>*. This allows extending the previous formulas (or tables). In particular we
651  have this 6-bit function (based on the 3-bit function above):
652  - *f(f<sup>2</sup> - 2)*
653
654This loop, again extended to also handle *u*, *v*, *q*, and *r* alongside *f* and *g*, placed in
655`divsteps_n_matrix`, gives a significantly faster, but non-constant time version.
656
657
658## 7. Final Python version
659
660All together we need the following functions:
661
662- A way to compute the transition matrix in constant time, using the `divsteps_n_matrix` function
663  from section 2, but with its loop replaced by a variant of the constant-time divstep from
664  section 5, extended to handle *u*, *v*, *q*, *r*:
665
666```python
667def divsteps_n_matrix(zeta, f, g):
668    """Compute zeta and transition matrix t after N divsteps (multiplied by 2^N)."""
669    u, v, q, r = 1, 0, 0, 1 # start with identity matrix
670    for _ in range(N):
671        c1 = zeta >> 63
672        # Compute x, y, z as conditionally-negated versions of f, u, v.
673        x, y, z = (f ^ c1) - c1, (u ^ c1) - c1, (v ^ c1) - c1
674        c2 = -(g & 1)
675        # Conditionally add x, y, z to g, q, r.
676        g, q, r = g + (x & c2), q + (y & c2), r + (z & c2)
677        c1 &= c2                     # reusing c1 here for the earlier c3 variable
678        zeta = (zeta ^ c1) - 1       # inlining the unconditional zeta decrement here
679        # Conditionally add g, q, r to f, u, v.
680        f, u, v = f + (g & c1), u + (q & c1), v + (r & c1)
681        # When shifting g down, don't shift q, r, as we construct a transition matrix multiplied
682        # by 2^N. Instead, shift f's coefficients u and v up.
683        g, u, v = g >> 1, u << 1, v << 1
684    return zeta, (u, v, q, r)
685```
686
687- The functions to update *f* and *g*, and *d* and *e*, from section 2 and section 4, with the constant-time
688  changes to `update_de` from section 5:
689
690```python
691def update_fg(f, g, t):
692    """Multiply matrix t/2^N with [f, g]."""
693    u, v, q, r = t
694    cf, cg = u*f + v*g, q*f + r*g
695    return cf >> N, cg >> N
696
697def update_de(d, e, t, M, Mi):
698    """Multiply matrix t/2^N with [d, e], modulo M."""
699    u, v, q, r = t
700    d_sign, e_sign = d >> 257, e >> 257
701    md, me = (u & d_sign) + (v & e_sign), (q & d_sign) + (r & e_sign)
702    cd, ce = (u*d + v*e) % 2**N, (q*d + r*e) % 2**N
703    md -= (Mi*cd + md) % 2**N
704    me -= (Mi*ce + me) % 2**N
705    cd, ce = u*d + v*e + M*md, q*d + r*e + M*me
706    return cd >> N, ce >> N
707```
708
709- The `normalize` function from section 4, made constant time as well:
710
711```python
712def normalize(sign, v, M):
713    """Compute sign*v mod M, where v in (-2*M,M); output in [0,M)."""
714    v_sign = v >> 257
715    # Conditionally add M to v.
716    v += M & v_sign
717    c = (sign - 1) >> 1
718    # Conditionally negate v.
719    v = (v ^ c) - c
720    v_sign = v >> 257
721    # Conditionally add M to v again.
722    v += M & v_sign
723    return v
724```
725
726- And finally the `modinv` function too, adapted to use *&zeta;* instead of *&delta;*, and using the fixed
727  iteration count from section 5:
728
729```python
730def modinv(M, Mi, x):
731    """Compute the modular inverse of x mod M, given Mi=1/M mod 2^N."""
732    zeta, f, g, d, e = -1, M, x, 0, 1
733    for _ in range((590 + N - 1) // N):
734        zeta, t = divsteps_n_matrix(zeta, f % 2**N, g % 2**N)
735        f, g = update_fg(f, g, t)
736        d, e = update_de(d, e, t, M, Mi)
737    return normalize(f, d, M)
738```
739
740- To get a variable time version, replace the `divsteps_n_matrix` function with one that uses the
741  divsteps loop from section 5, and a `modinv` version that calls it without the fixed iteration
742  count:
743
744```python
745NEGINV16 = [15, 5, 3, 9, 7, 13, 11, 1] # NEGINV16[n//2] = (-n)^-1 mod 16, for odd n
746def divsteps_n_matrix_var(eta, f, g):
747    """Compute eta and transition matrix t after N divsteps (multiplied by 2^N)."""
748    u, v, q, r = 1, 0, 0, 1
749    i = N
750    while True: