1 2 3(* Copyright (c) 2011 Stefan Krah. All rights reserved. *) 4 5 6========================================================================== 7 Calculate (a * b) % p using special primes 8========================================================================== 9 10A description of the algorithm can be found in the apfloat manual by 11Tommila [1]. 12 13 14Definitions: 15------------ 16 17In the whole document, "==" stands for "is congruent with". 18 19Result of a * b in terms of high/low words: 20 21 (1) hi * 2**64 + lo = a * b 22 23Special primes: 24 25 (2) p = 2**64 - z + 1, where z = 2**n 26 27Single step modular reduction: 28 29 (3) R(hi, lo) = hi * z - hi + lo 30 31 32Strategy: 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 44The 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 59Maximum numbers of step b): 60--------------------------- 61 62# To avoid unnecessary formalism, define: 63 64def 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: 72hi=2**64-1; lo=2**64-1; z=2**32 73p1 = 2**64 - z + 1 74 75hi, lo = R(hi, lo, z) # First reduction 76hi, lo = R(hi, lo, z) # Second reduction 77hi * 2**64 + lo < 2 * p1 # True 78 79# For p2: 80hi=2**64-1; lo=2**64-1; z=2**34 81p2 = 2**64 - z + 1 82 83hi, lo = R(hi, lo, z) # First reduction 84hi, lo = R(hi, lo, z) # Second reduction 85hi, lo = R(hi, lo, z) # Third reduction 86hi * 2**64 + lo < 2 * p2 # True 87 88# For p3: 89hi=2**64-1; lo=2**64-1; z=2**40 90p3 = 2**64 - z + 1 91 92hi, lo = R(hi, lo, z) # First reduction 93hi, lo = R(hi, lo, z) # Second reduction 94hi, lo = R(hi, lo, z) # Third reduction 95hi * 2**64 + lo < 2 * p3 # True 96 97 98Step 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