Home | History | Annotate | Download | only in literature 
      1 
      2 
      3 (* Copyright (c) 2011 Stefan Krah. All rights reserved. *)
      4 
      5 
      6 ==========================================================================
      7                 Calculate (a * b) % p using special primes
      8 ==========================================================================
      9 
     10 A description of the algorithm can be found in the apfloat manual by
     11 Tommila [1].
     12 
     13 
     14 Definitions:
     15 ------------
     16 
     17 In the whole document, "==" stands for "is congruent with".
     18 
     19 Result of a * b in terms of high/low words:
     20 
     21    (1) hi * 2**64 + lo = a * b
     22 
     23 Special primes:
     24 
     25    (2) p = 2**64 - z + 1, where z = 2**n
     26 
     27 Single step modular reduction:
     28 
     29    (3) R(hi, lo) = hi * z - hi + lo
     30 
     31 
     32 Strategy:
     33 ---------
     34 
     35    a) Set (hi, lo) to the result of a * b.
     36 
     37    b) Set (hi', lo') to the result of R(hi, lo).
     38 
     39    c) Repeat step b) until 0 <= hi' * 2**64 + lo' < 2*p.
     40 
     41    d) If the result is less than p, return lo'. Otherwise return lo' - p.
     42 
     43 
     44 The reduction step b) preserves congruence:
     45 -------------------------------------------
     46 
     47     hi * 2**64 + lo == hi * z - hi + lo   (mod p)
     48 
     49     Proof:
     50     ~~~~~~
     51 
     52        hi * 2**64 + lo = (2**64 - z + 1) * hi + z * hi - hi + lo
     53 
     54                        = p * hi               + z * hi - hi + lo
     55 
     56                        == z * hi - hi + lo   (mod p)
     57 
     58 
     59 Maximum numbers of step b):
     60 ---------------------------
     61 
     62 # To avoid unnecessary formalism, define:
     63 
     64 def R(hi, lo, z):
     65      return divmod(hi * z - hi + lo, 2**64)
     66 
     67 # For simplicity, assume hi=2**64-1, lo=2**64-1 after the
     68 # initial multiplication a * b. This is of course impossible
     69 # but certainly covers all cases.
     70 
     71 # Then, for p1:
     72 hi=2**64-1; lo=2**64-1; z=2**32
     73 p1 = 2**64 - z + 1
     74 
     75 hi, lo = R(hi, lo, z)    # First reduction
     76 hi, lo = R(hi, lo, z)    # Second reduction
     77 hi * 2**64 + lo < 2 * p1 # True
     78 
     79 # For p2:
     80 hi=2**64-1; lo=2**64-1; z=2**34
     81 p2 = 2**64 - z + 1
     82 
     83 hi, lo = R(hi, lo, z)    # First reduction
     84 hi, lo = R(hi, lo, z)    # Second reduction
     85 hi, lo = R(hi, lo, z)    # Third reduction
     86 hi * 2**64 + lo < 2 * p2 # True
     87 
     88 # For p3:
     89 hi=2**64-1; lo=2**64-1; z=2**40
     90 p3 = 2**64 - z + 1
     91 
     92 hi, lo = R(hi, lo, z)    # First reduction
     93 hi, lo = R(hi, lo, z)    # Second reduction
     94 hi, lo = R(hi, lo, z)    # Third reduction
     95 hi * 2**64 + lo < 2 * p3 # True
     96 
     97 
     98 Step d) preserves congruence and yields a result < p:
     99 -----------------------------------------------------
    100 
    101    Case hi = 0:
    102 
    103        Case lo < p: trivial.
    104 
    105        Case lo >= p:
    106 
    107           lo == lo - p   (mod p)             # result is congruent
    108 
    109           p <= lo < 2*p  ->  0 <= lo - p < p # result is in the correct range
    110 
    111    Case hi = 1:
    112 
    113        p < 2**64 /\ 2**64 + lo < 2*p  ->  lo < p  # lo is always less than p
    114 
    115        2**64 + lo == 2**64 + (lo - p)   (mod p)   # result is congruent
    116 
    117                   = lo - p   # exactly the same value as the previous RHS
    118                              # in uint64_t arithmetic.
    119 
    120        p < 2**64 + lo < 2*p  ->  0 < 2**64 + (lo - p) < p  # correct range
    121 
    122 
    123 
    124 [1]  http://www.apfloat.org/apfloat/2.40/apfloat.pdf
    125 
    126 
    127 
    128