Inside Bitcoin

251## 4. Avoiding modulus operations
252
253So far, there are two places where we compute a remainder of big numbers modulo *M*: at the end of
254`div2n` in every `update_de`, and at the very end of `modinv` after potentially negating *d* due to the
255sign of *f*. These are relatively expensive operations when done generically.
256
257To deal with the modulus operation in `div2n`, we simply stop requiring *d* and *e* to be in range
258*[0,M)* all the time. Let's start by inlining `div2n` into `update_de`, and dropping the modulus
259operation at the end:
260
261```python
262def update_de(d, e, t, M, Mi):
263    """Multiply matrix t/2^N with [d, e] mod M, given Mi=1/M mod 2^N."""
264    u, v, q, r = t
265    cd, ce = u*d + v*e, q*d + r*e
266    # Cancel out bottom N bits of cd and ce.
267    md = -((Mi * cd) % 2**N)
268    me = -((Mi * ce) % 2**N)
269    cd += md * M
270    ce += me * M
271    # And cleanly divide by 2**N.
272    return cd >> N, ce >> N
273```
274
275Let's look at bounds on the ranges of these numbers. It can be shown that *|u|+|v|* and *|q|+|r|*
276never exceed *2<sup>N</sup>* (see paragraph 8.3 in the paper), and thus a multiplication with *t* will have
277outputs whose absolute values are at most *2<sup>N</sup>* times the maximum absolute input value. In case the
278inputs *d* and *e* are in *(-M,M)*, which is certainly true for the initial values *d=0* and *e=1* assuming
279*M > 1*, the multiplication results in numbers in range *(-2<sup>N</sup>M,2<sup>N</sup>M)*. Subtracting less than *2<sup>N</sup>*
280times *M* to cancel out *N* bits brings that up to *(-2<sup>N+1</sup>M,2<sup>N</sup>M)*, and
281dividing by *2<sup>N</sup>* at the end takes it to *(-2M,M)*. Another application of `update_de` would take that
282to *(-3M,2M)*, and so forth. This progressive expansion of the variables' ranges can be
283counteracted by incrementing *d* and *e* by *M* whenever they're negative:
284
285```python
286    ...
287    if d < 0:
288        d += M
289    if e < 0:
290        e += M
291    cd, ce = u*d + v*e, q*d + r*e
292    # Cancel out bottom N bits of cd and ce.
293    ...
294```
295
296With inputs in *(-2M,M)*, they will first be shifted into range *(-M,M)*, which means that the
297output will again be in *(-2M,M)*, and this remains the case regardless of how many `update_de`
298invocations there are. In what follows, we will try to make this more efficient.
299
300Note that increasing *d* by *M* is equal to incrementing *cd* by *u&thinsp;M* and *ce* by *q&thinsp;M*. Similarly,
301increasing *e* by *M* is equal to incrementing *cd* by *v&thinsp;M* and *ce* by *r&thinsp;M*. So we could instead write:
302
303```python
304    ...
305    cd, ce = u*d + v*e, q*d + r*e
306    # Perform the equivalent of incrementing d, e by M when they're negative.
307    if d < 0:
308        cd += u*M
309        ce += q*M
310    if e < 0:
311        cd += v*M
312        ce += r*M
313    # Cancel out bottom N bits of cd and ce.
314    md = -((Mi * cd) % 2**N)
315    me = -((Mi * ce) % 2**N)
316    cd += md * M
317    ce += me * M
318    ...
319```
320
321Now note that we have two steps of corrections to *cd* and *ce* that add multiples of *M*: this
322increment, and the decrement that cancels out bottom bits. The second one depends on the first
323one, but they can still be efficiently combined by only computing the bottom bits of *cd* and *ce*
324at first, and using that to compute the final *md*, *me* values:
325
326```python
327def update_de(d, e, t, M, Mi):
328    """Multiply matrix t/2^N with [d, e], modulo M."""
329    u, v, q, r = t
330    md, me = 0, 0
331    # Compute what multiples of M to add to cd and ce.
332    if d < 0:
333        md += u
334        me += q
335    if e < 0:
336        md += v
337        me += r
338    # Compute bottom N bits of t*[d,e] + M*[md,me].
339    cd, ce = (u*d + v*e + md*M) % 2**N, (q*d + r*e + me*M) % 2**N
340    # Correct md and me such that the bottom N bits of t*[d,e] + M*[md,me] are zero.
341    md -= (Mi * cd) % 2**N
342    me -= (Mi * ce) % 2**N
343    # Do the full computation.
344    cd, ce = u*d + v*e + md*M, q*d + r*e + me*M
345    # And cleanly divide by 2**N.
346    return cd >> N, ce >> N
347```
348
349One last optimization: we can avoid the *md&thinsp;M* and *me&thinsp;M* multiplications in the bottom bits of *cd*
350and *ce* by moving them to the *md* and *me* correction:
351
352```python
353    ...
354    # Compute bottom N bits of t*[d,e].
355    cd, ce = (u*d + v*e) % 2**N, (q*d + r*e) % 2**N
356    # Correct md and me such that the bottom N bits of t*[d,e]+M*[md,me] are zero.
357    # Note that this is not the same as {md = (-Mi * cd) % 2**N} etc. That would also result in N
358    # zero bottom bits, but isn't guaranteed to be a reduction of [0,2^N) compared to the
359    # previous md and me values, and thus would violate our bounds analysis.
360    md -= (Mi*cd + md) % 2**N
361    me -= (Mi*ce + me) % 2**N
362    ...
363```
364
365The resulting function takes *d* and *e* in range *(-2M,M)* as inputs, and outputs values in the same
366range. That also means that the *d* value at the end of `modinv` will be in that range, while we want
367a result in *[0,M)*. To do that, we need a normalization function. It's easy to integrate the
368conditional negation of *d* (based on the sign of *f*) into it as well:
369
370```python
371def normalize(sign, v, M):
372    """Compute sign*v mod M, where v is in range (-2*M,M); output in [0,M)."""
373    assert sign == 1 or sign == -1
374    # v in (-2*M,M)
375    if v < 0:
376        v += M
377    # v in (-M,M). Now multiply v with sign (which can only be 1 or -1).
378    if sign == -1:
379        v = -v
380    # v in (-M,M)
381    if v < 0:
382        v += M
383    # v in [0,M)
384    return v
385```
386
387And calling it in `modinv` is simply:
388
389```python
390   ...
391   return normalize(f, d, M)
392```
393
394
395## 5. Constant-time operation
396
397The primary selling point of the algorithm is fast constant-time operation. What code flow still
398depends on the input data so far?
399
400- the number of iterations of the while *g &ne; 0* loop in `modinv`
401- the branches inside `divsteps_n_matrix`
402- the sign checks in `update_de`
403- the sign checks in `normalize`
404
405To make the while loop in `modinv` constant time it can be replaced with a constant number of
406iterations. The paper proves (Theorem 11.2) that *741* divsteps are sufficient for any *256*-bit
407inputs, and [safegcd-bounds](https://github.com/sipa/safegcd-bounds) shows that the slightly better bound *724* is
408sufficient even. Given that every loop iteration performs *N* divsteps, it will run a total of
409*&lceil;724/N&rceil;* times.
410
411To deal with the branches in `divsteps_n_matrix` we will replace them with constant-time bitwise
412operations (and hope the C compiler isn't smart enough to turn them back into branches; see
413`ctime_tests.c` for automated tests that this isn't the case). To do so, observe that a
414divstep can be written instead as (compare to the inner loop of `gcd` in section 1).
415
416```python
417    x = -f if delta > 0 else f         # set x equal to (input) -f or f
418    if g & 1:
419        g += x                         # set g to (input) g-f or g+f
420        if delta > 0:
421            delta = -delta
422            f += g                     # set f to (input) g (note that g was set to g-f before)
423    delta += 1
424    g >>= 1
425```
426
427To convert the above to bitwise operations, we rely on a trick to negate conditionally: per the
428definition of negative numbers in two's complement, (*-v == ~v + 1*) holds for every number *v*. As
429*-1* in two's complement is all *1* bits, bitflipping can be expressed as xor with *-1*. It follows
430that *-v == (v ^ -1) - (-1)*. Thus, if we have a variable *c* that takes on values *0* or *-1*, then
431*(v ^ c) - c* is *v* if *c=0* and *-v* if *c=-1*.
432
433Using this we can write:
434
435```python
436    x = -f if delta > 0 else f
437```
438
439in constant-time form as:
440
441```python
442    c1 = (-delta) >> 63
443    # Conditionally negate f based on c1:
444    x = (f ^ c1) - c1
445```
446
447To use that trick, we need a helper mask variable *c1* that resolves the condition *&delta;>0* to *-1*
448(if true) or *0* (if false). We compute *c1* using right shifting, which is equivalent to dividing by
449the specified power of *2* and rounding down (in Python, and also in C under the assumption of a typical two's complement system; see
450`assumptions.h` for tests that this is the case). Right shifting by *63* thus maps all
451numbers in range *[-2<sup>63</sup>,0)* to *-1*, and numbers in range *[0,2<sup>63</sup>)* to *0*.
452
453Using the facts that *x&0=0* and *x&(-1)=x* (on two's complement systems again), we can write:
454
455```python
456    if g & 1:
457        g += x
458```
459
460as:
461
462```python
463    # Compute c2=0 if g is even and c2=-1 if g is odd.
464    c2 = -(g & 1)
465    # This masks out x if g is even, and leaves x be if g is odd.
466    g += x & c2
467```
468
469Using the conditional negation trick again we can write:
470
471```python
472    if g & 1:
473        if delta > 0:
474            delta = -delta
475```
476
477as:
478
479```python
480    # Compute c3=-1 if g is odd and delta>0, and 0 otherwise.
481    c3 = c1 & c2
482    # Conditionally negate delta based on c3:
483    delta = (delta ^ c3) - c3
484```
485
486Finally:
487
488```python
489    if g & 1:
490        if delta > 0:
491            f += g
492```
493
494becomes:
495
496```python
497    f += g & c3
498```
499
500It turns out that this can be implemented more efficiently by applying the substitution