Inside Bitcoin

1# The safegcd implementation in libsecp256k1 explained
2
3This document explains the modular inverse and Jacobi symbol implementations in the `src/modinv*.h` files.
4It is based on the paper
5["Fast constant-time gcd computation and modular inversion"](https://gcd.cr.yp.to/papers.html#safegcd)
6by Daniel J. Bernstein and Bo-Yin Yang. The references below are for the Date: 2019.04.13 version.
7
8The actual implementation is in C of course, but for demonstration purposes Python3 is used here.
9Most implementation aspects and optimizations are explained, except those that depend on the specific
10number representation used in the C code.
11
12## 1. Computing the Greatest Common Divisor (GCD) using divsteps
13
14The algorithm from the paper (section 11), at a very high level, is this:
15
16```python
17def gcd(f, g):
18    """Compute the GCD of an odd integer f and another integer g."""
19    assert f & 1  # require f to be odd
20    delta = 1     # additional state variable
21    while g != 0:
22        assert f & 1  # f will be odd in every iteration
23        if delta > 0 and g & 1:
24            delta, f, g = 1 - delta, g, (g - f) // 2
25        elif g & 1:
26            delta, f, g = 1 + delta, f, (g + f) // 2
27        else:
28            delta, f, g = 1 + delta, f, (g    ) // 2
29    return abs(f)
30```
31
32It computes the greatest common divisor of an odd integer *f* and any integer *g*. Its inner loop
33keeps rewriting the variables *f* and *g* alongside a state variable *δ* that starts at *1*, until
34*g=0* is reached. At that point, *|f|* gives the GCD. Each of the transitions in the loop is called a
35"division step" (referred to as divstep in what follows).
36
37For example, *gcd(21, 14)* would be computed as:
38- Start with *δ=1 f=21 g=14*
39- Take the third branch: *δ=2 f=21 g=7*
40- Take the first branch: *δ=-1 f=7 g=-7*
41- Take the second branch: *δ=0 f=7 g=0*
42- The answer *|f| = 7*.
43
44Why it works:
45- Divsteps can be decomposed into two steps (see paragraph 8.2 in the paper):
46  - (a) If *g* is odd, replace *(f,g)* with *(g,g-f)* or (f,g+f), resulting in an even *g*.
47  - (b) Replace *(f,g)* with *(f,g/2)* (where *g* is guaranteed to be even).
48- Neither of those two operations change the GCD:
49  - For (a), assume *gcd(f,g)=c*, then it must be the case that *f=a c* and *g=b c* for some integers *a*
50    and *b*. As *(g,g-f)=(b c,(b-a)c)* and *(f,f+g)=(a c,(a+b)c)*, the result clearly still has
51    common factor *c*. Reasoning in the other direction shows that no common factor can be added by
52    doing so either.
53  - For (b), we know that *f* is odd, so *gcd(f,g)* clearly has no factor *2*, and we can remove
54    it from *g*.
55- The algorithm will eventually converge to *g=0*. This is proven in the paper (see theorem G.3).
56- It follows that eventually we find a final value *f'* for which *gcd(f,g) = gcd(f',0)*. As the
57  gcd of *f'* and *0* is *|f'|* by definition, that is our answer.
58
59Compared to more [traditional GCD algorithms](https://en.wikipedia.org/wiki/Euclidean_algorithm), this one has the property of only ever looking at
60the low-order bits of the variables to decide the next steps, and being easy to make
61constant-time (in more low-level languages than Python). The *δ* parameter is necessary to
62guide the algorithm towards shrinking the numbers' magnitudes without explicitly needing to look
63at high order bits.
64
65Properties that will become important later:
66- Performing more divsteps than needed is not a problem, as *f* does not change anymore after *g=0*.
67- Only even numbers are divided by *2*. This means that when reasoning about it algebraically we
68  do not need to worry about rounding.
69- At every point during the algorithm's execution the next *N* steps only depend on the bottom *N*
70  bits of *f* and *g*, and on *δ*.
71
72
73## 2. From GCDs to modular inverses
74
75We want an algorithm to compute the inverse *a* of *x* modulo *M*, i.e. the number a such that *a x=1
76mod M*. This inverse only exists if the GCD of *x* and *M* is *1*, but that is always the case if *M* is
77prime and *0 < x < M*. In what follows, assume that the modular inverse exists.
78It turns out this inverse can be computed as a side effect of computing the GCD by keeping track
79of how the internal variables can be written as linear combinations of the inputs at every step
80(see the [extended Euclidean algorithm](https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm)).
81Since the GCD is *1*, such an algorithm will compute numbers *a* and *b* such that a&thinsp;x + b&thinsp;M = 1*.
82Taking that expression *mod M* gives *a&thinsp;x mod M = 1*, and we see that *a* is the modular inverse of *x
83mod M*.
84
85A similar approach can be used to calculate modular inverses using the divsteps-based GCD
86algorithm shown above, if the modulus *M* is odd. To do so, compute *gcd(f=M,g=x)*, while keeping
87track of extra variables *d* and *e*, for which at every step *d = f/x (mod M)* and *e = g/x (mod M)*.
88*f/x* here means the number which multiplied with *x* gives *f mod M*. As *f* and *g* are initialized to *M*
89and *x* respectively, *d* and *e* just start off being *0* (*M/x mod M = 0/x mod M = 0*) and *1* (*x/x mod M
90= 1*).
91
92```python
93def div2(M, x):
94    """Helper routine to compute x/2 mod M (where M is odd)."""
95    assert M & 1
96    if x & 1: # If x is odd, make it even by adding M.
97        x += M
98    # x must be even now, so a clean division by 2 is possible.
99    return x // 2
100
101def modinv(M, x):
102    """Compute the inverse of x mod M (given that it exists, and M is odd)."""
103    assert M & 1
104    delta, f, g, d, e = 1, M, x, 0, 1
105    while g != 0:
106        # Note that while division by two for f and g is only ever done on even inputs, this is
107        # not true for d and e, so we need the div2 helper function.
108        if delta > 0 and g & 1:
109            delta, f, g, d, e = 1 - delta, g, (g - f) // 2, e, div2(M, e - d)
110        elif g & 1:
111            delta, f, g, d, e = 1 + delta, f, (g + f) // 2, d, div2(M, e + d)
112        else:
113            delta, f, g, d, e = 1 + delta, f, (g    ) // 2, d, div2(M, e    )
114        # Verify that the invariants d=f/x mod M, e=g/x mod M are maintained.
115        assert f % M == (d * x) % M
116        assert g % M == (e * x) % M
117    assert f == 1 or f == -1  # |f| is the GCD, it must be 1
118    # Because of invariant d = f/x (mod M), 1/x = d/f (mod M). As |f|=1, d/f = d*f.
119    return (d * f) % M
120```
121
122Also note that this approach to track *d* and *e* throughout the computation to determine the inverse
123is different from the paper. There (see paragraph 12.1 in the paper) a transition matrix for the
124entire computation is determined (see section 3 below) and the inverse is computed from that.
125The approach here avoids the need for 2x2 matrix multiplications of various sizes, and appears to
126be faster at the level of optimization we're able to do in C.
127
128
129## 3. Batching multiple divsteps
130
131Every divstep can be expressed as a matrix multiplication, applying a transition matrix *(1/2 t)*
132to both vectors *[f, g]* and *[d, e]* (see paragraph 8.1 in the paper):
133
134```
135  t = [ u,  v ]
136      [ q,  r ]
137
138  [ out_f ] = (1/2 * t) * [ in_f ]
139  [ out_g ] =             [ in_g ]
140
141  [ out_d ] = (1/2 * t) * [ in_d ]  (mod M)
142  [ out_e ]               [ in_e ]
143```
144
145where *(u, v, q, r)* is *(0, 2, -1, 1)*, *(2, 0, 1, 1)*, or *(2, 0, 0, 1)*, depending on which branch is
146taken. As above, the resulting *f* and *g* are always integers.
147
148Performing multiple divsteps corresponds to a multiplication with the product of all the
149individual divsteps' transition matrices. As each transition matrix consists of integers
150divided by *2*, the product of these matrices will consist of integers divided by *2<sup>N</sup>* (see also
151theorem 9.2 in the paper). These divisions are expensive when updating *d* and *e*, so we delay
152them: we compute the integer coefficients of the combined transition matrix scaled by *2<sup>N</sup>*, and
153do one division by *2<sup>N</sup>* as a final step:
154
155```python
156def divsteps_n_matrix(delta, f, g):
157    """Compute delta and transition matrix t after N divsteps (multiplied by 2^N)."""
158    u, v, q, r = 1, 0, 0, 1 # start with identity matrix
159    for _ in range(N):
160        if delta > 0 and g & 1:
161            delta, f, g, u, v, q, r = 1 - delta, g, (g - f) // 2, 2*q, 2*r, q-u, r-v
162        elif g & 1:
163            delta, f, g, u, v, q, r = 1 + delta, f, (g + f) // 2, 2*u, 2*v, q+u, r+v
164        else:
165            delta, f, g, u, v, q, r = 1 + delta, f, (g    ) // 2, 2*u, 2*v, q  , r
166    return delta, (u, v, q, r)
167```
168
169As the branches in the divsteps are completely determined by the bottom *N* bits of *f* and *g*, this
170function to compute the transition matrix only needs to see those bottom bits. Furthermore all
171intermediate results and outputs fit in *(N+1)*-bit numbers (unsigned for *f* and *g*; signed for *u*, *v*,
172*q*, and *r*) (see also paragraph 8.3 in the paper). This means that an implementation using 64-bit
173integers could set *N=62* and compute the full transition matrix for 62 steps at once without any
174big integer arithmetic at all. This is the reason why this algorithm is efficient: it only needs
175to update the full-size *f*, *g*, *d*, and *e* numbers once every *N* steps.
176
177We still need functions to compute:
178
179```
180  [ out_f ] = (1/2^N * [ u,  v ]) * [ in_f ]
181  [ out_g ]   (        [ q,  r ])   [ in_g ]
182
183  [ out_d ] = (1/2^N * [ u,  v ]) * [ in_d ]  (mod M)
184  [ out_e ]   (        [ q,  r ])   [ in_e ]
185```
186
187Because the divsteps transformation only ever divides even numbers by two, the result of *t&thinsp;[f,g]* is always even. When *t* is a composition of *N* divsteps, it follows that the resulting *f*
188and *g* will be multiple of *2<sup>N</sup>*, and division by *2<sup>N</sup>* is simply shifting them down:
189
190```python
191def update_fg(f, g, t):
192    """Multiply matrix t/2^N with [f, g]."""
193    u, v, q, r = t
194    cf, cg = u*f + v*g, q*f + r*g
195    # (t / 2^N) should cleanly apply to [f,g] so the result of t*[f,g] should have N zero
196    # bottom bits.
197    assert cf % 2**N == 0
198    assert cg % 2**N == 0
199    return cf >> N, cg >> N
200```
201
202The same is not true for *d* and *e*, and we need an equivalent of the `div2` function for division by *2<sup>N</sup> mod M*.
203This is easy if we have precomputed *1/M mod 2<sup>N</sup>* (which always exists for odd *M*):
204
205```python
206def div2n(M, Mi, x):
207    """Compute x/2^N mod M, given Mi = 1/M mod 2^N."""
208    assert (M * Mi) % 2**N == 1
209    # Find a factor m such that m*M has the same bottom N bits as x. We want:
210    #     (m * M) mod 2^N = x mod 2^N
211    # <=> m mod 2^N = (x / M) mod 2^N
212    # <=> m mod 2^N = (x * Mi) mod 2^N
213    m = (Mi * x) % 2**N
214    # Subtract that multiple from x, cancelling its bottom N bits.
215    x -= m * M
216    # Now a clean division by 2^N is possible.
217    assert x % 2**N == 0
218    return (x >> N) % M
219
220def update_de(d, e, t, M, Mi):
221    """Multiply matrix t/2^N with [d, e], modulo M."""
222    u, v, q, r = t
223    cd, ce = u*d + v*e, q*d + r*e
224    return div2n(M, Mi, cd), div2n(M, Mi, ce)
225```
226
227With all of those, we can write a version of `modinv` that performs *N* divsteps at once:
228
229```python3
230def modinv(M, Mi, x):
231    """Compute the modular inverse of x mod M, given Mi=1/M mod 2^N."""
232    assert M & 1
233    delta, f, g, d, e = 1, M, x, 0, 1
234    while g != 0:
235        # Compute the delta and transition matrix t for the next N divsteps (this only needs
236        # (N+1)-bit signed integer arithmetic).
237        delta, t = divsteps_n_matrix(delta, f % 2**N, g % 2**N)
238        # Apply the transition matrix t to [f, g]:
239        f, g = update_fg(f, g, t)
240        # Apply the transition matrix t to [d, e]:
241        d, e = update_de(d, e, t, M, Mi)
242    return (d * f) % M
243```
244
245This means that in practice we'll always perform a multiple of *N* divsteps. This is not a problem
246because once *g=0*, further divsteps do not affect *f*, *g*, *d*, or *e* anymore (only *&delta;* keeps
247increasing). For variable time code such excess iterations will be mostly optimized away in later
248sections.
249
250