Inside Bitcoin

751        zeros = min(i, count_trailing_zeros(g))
752        eta, i = eta - zeros, i - zeros
753        g, u, v = g >> zeros, u << zeros, v << zeros
754        if i == 0:
755            break
756        if eta < 0:
757            eta, f, u, v, g, q, r = -eta, g, q, r, -f, -u, -v
758        limit = min(min(eta + 1, i), 4)
759        w = (g * NEGINV16[(f & 15) // 2]) % (2**limit)
760        g, q, r = g + w*f, q + w*u, r + w*v
761    return eta, (u, v, q, r)
762
763def modinv_var(M, Mi, x):
764    """Compute the modular inverse of x mod M, given Mi = 1/M mod 2^N."""
765    eta, f, g, d, e = -1, M, x, 0, 1
766    while g != 0:
767        eta, t = divsteps_n_matrix_var(eta, f % 2**N, g % 2**N)
768        f, g = update_fg(f, g, t)
769        d, e = update_de(d, e, t, M, Mi)
770    return normalize(f, d, Mi)
771```
772
773## 8. From GCDs to Jacobi symbol
774
775We can also use a similar approach to calculate Jacobi symbol *(x | M)* by keeping track of an
776extra variable *j*, for which at every step *(x | M) = j (g | f)*. As we update *f* and *g*, we
777make corresponding updates to *j* using
778[properties of the Jacobi symbol](https://en.wikipedia.org/wiki/Jacobi_symbol#Properties):
779* *((g/2) | f)* is either *(g | f)* or *-(g | f)*, depending on the value of *f mod 8* (negating if it's *3* or *5*).
780* *(f | g)* is either *(g | f)* or *-(g | f)*, depending on *f mod 4* and *g mod 4* (negating if both are *3*).
781
782These updates depend only on the values of *f* and *g* modulo *4* or *8*, and can thus be applied
783very quickly, as long as we keep track of a few additional bits of *f* and *g*. Overall, this
784calculation is slightly simpler than the one for the modular inverse because we no longer need to
785keep track of *d* and *e*.
786
787However, one difficulty of this approach is that the Jacobi symbol *(a | n)* is only defined for
788positive odd integers *n*, whereas in the original safegcd algorithm, *f, g* can take negative
789values. We resolve this by using the following modified steps:
790
791```python
792        # Before
793        if delta > 0 and g & 1:
794            delta, f, g = 1 - delta, g, (g - f) // 2
795
796        # After
797        if delta > 0 and g & 1:
798            delta, f, g = 1 - delta, g, (g + f) // 2
799```
800
801The algorithm is still correct, since the changed divstep, called a "posdivstep" (see section 8.4
802and E.5 in the paper) preserves *gcd(f, g)*. However, there's no proof that the modified algorithm
803will converge. The justification for posdivsteps is completely empirical: in practice, it appears
804that the vast majority of nonzero inputs converge to *f=g=gcd(f<sub>0</sub>, g<sub>0</sub>)* in a
805number of steps proportional to their logarithm.
806
807Note that:
808- We require inputs to satisfy *gcd(x, M) = 1*, as otherwise *f=1* is not reached.
809- We require inputs *x &neq; 0*, because applying posdivstep with *g=0* has no effect.
810- We need to update the termination condition from *g=0* to *f=1*.
811
812We account for the possibility of nonconvergence by only performing a bounded number of
813posdivsteps, and then falling back to square-root based Jacobi calculation if a solution has not
814yet been found.
815
816The optimizations in sections 3-7 above are described in the context of the original divsteps, but
817in the C implementation we also adapt most of them (not including "avoiding modulus operations",
818since it's not necessary to track *d, e*, and "constant-time operation", since we never calculate
819Jacobi symbols for secret data) to the posdivsteps version.
820