OLD | NEW |
(Empty) | |
| 1 # |
| 2 # Copyright (c) 2008-2010 Stefan Krah. All rights reserved. |
| 3 # |
| 4 # Redistribution and use in source and binary forms, with or without |
| 5 # modification, are permitted provided that the following conditions |
| 6 # are met: |
| 7 # |
| 8 # 1. Redistributions of source code must retain the above copyright |
| 9 # notice, this list of conditions and the following disclaimer. |
| 10 # |
| 11 # 2. Redistributions in binary form must reproduce the above copyright |
| 12 # notice, this list of conditions and the following disclaimer in the |
| 13 # documentation and/or other materials provided with the distribution. |
| 14 # |
| 15 # THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND |
| 16 # ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
| 17 # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
| 18 # ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE |
| 19 # FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL |
| 20 # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS |
| 21 # OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |
| 22 # HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT |
| 23 # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY |
| 24 # OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF |
| 25 # SUCH DAMAGE. |
| 26 # |
| 27 |
| 28 |
| 29 ###################################################################### |
| 30 # This file lists and checks some of the constants and limits used # |
| 31 # in libmpdec's Number Theoretic Transform. At the end of the file # |
| 32 # there is an example function for the plain DFT transform. # |
| 33 ###################################################################### |
| 34 |
| 35 |
| 36 # |
| 37 # Number theoretic transforms are done in subfields of F(p). P[i] |
| 38 # are the primes, D[i] = P[i] - 1 are highly composite and w[i] |
| 39 # are the respective primitive roots of F(p). |
| 40 # |
| 41 # The strategy is to convolute two coefficients modulo all three |
| 42 # primes, then use the Chinese Remainder Theorem on the three |
| 43 # result arrays to recover the result in the usual base RADIX |
| 44 # form. |
| 45 # |
| 46 |
| 47 # ====================================================================== |
| 48 # Primitive roots |
| 49 # ====================================================================== |
| 50 |
| 51 # |
| 52 # Verify primitive roots: |
| 53 # |
| 54 # For a prime field, r is a primitive root if and only if for all prime |
| 55 # factors f of p-1, r**((p-1)/f) =/= 1 (mod p). |
| 56 # |
| 57 def prod(F, E): |
| 58 """Check that the factorization of P-1 is correct. F is the list of |
| 59 factors of P-1, E lists the number of occurrences of each factor.""" |
| 60 x = 1 |
| 61 for y, z in zip(F, E): |
| 62 x *= y**z |
| 63 return x |
| 64 |
| 65 def is_primitive_root(r, p, factors, exponents): |
| 66 """Check if r is a primitive root of F(p).""" |
| 67 if p != prod(factors, exponents) + 1: |
| 68 return False |
| 69 for f in factors: |
| 70 q, control = divmod(p-1, f) |
| 71 if control != 0: |
| 72 return False |
| 73 if pow(r, q, p) == 1: |
| 74 return False |
| 75 return True |
| 76 |
| 77 |
| 78 # ================================================================= |
| 79 # Constants and limits for the 64-bit version |
| 80 # ================================================================= |
| 81 |
| 82 RADIX = 10**19 |
| 83 |
| 84 # Primes P1, P2 and P3: |
| 85 P = [2**64-2**32+1, 2**64-2**34+1, 2**64-2**40+1] |
| 86 |
| 87 # P-1, highly composite. The transform length d is variable and |
| 88 # must divide D = P-1. Since all D are divisible by 3 * 2**32, |
| 89 # transform lengths can be 2**n or 3 * 2**n (where n <= 32). |
| 90 D = [2**32 * 3 * (5 * 17 * 257 * 65537), |
| 91 2**34 * 3**2 * (7 * 11 * 31 * 151 * 331), |
| 92 2**40 * 3**2 * (5 * 7 * 13 * 17 * 241)] |
| 93 |
| 94 # Prime factors of P-1 and their exponents: |
| 95 F = [(2,3,5,17,257,65537), (2,3,7,11,31,151,331), (2,3,5,7,13,17,241)] |
| 96 E = [(32,1,1,1,1,1), (34,2,1,1,1,1,1), (40,2,1,1,1,1,1)] |
| 97 |
| 98 # Maximum transform length for 2**n. Above that only 3 * 2**31 |
| 99 # or 3 * 2**32 are possible. |
| 100 MPD_MAXTRANSFORM_2N = 2**32 |
| 101 |
| 102 |
| 103 # Limits in the terminology of Pollard's paper: |
| 104 m2 = (MPD_MAXTRANSFORM_2N * 3) // 2 # Maximum length of the smaller array. |
| 105 M1 = M2 = RADIX-1 # Maximum value per single word. |
| 106 L = m2 * M1 * M2 |
| 107 P[0] * P[1] * P[2] > 2 * L |
| 108 |
| 109 |
| 110 # Primitive roots of F(P1), F(P2) and F(P3): |
| 111 w = [7, 10, 19] |
| 112 |
| 113 # The primitive roots are correct: |
| 114 for i in range(3): |
| 115 if not is_primitive_root(w[i], P[i], F[i], E[i]): |
| 116 print("FAIL") |
| 117 |
| 118 |
| 119 # ================================================================= |
| 120 # Constants and limits for the 32-bit version |
| 121 # ================================================================= |
| 122 |
| 123 RADIX = 10**9 |
| 124 |
| 125 # Primes P1, P2 and P3: |
| 126 P = [2113929217, 2013265921, 1811939329] |
| 127 |
| 128 # P-1, highly composite. All D = P-1 are divisible by 3 * 2**25, |
| 129 # allowing for transform lengths up to 3 * 2**25 words. |
| 130 D = [2**25 * 3**2 * 7, |
| 131 2**27 * 3 * 5, |
| 132 2**26 * 3**3] |
| 133 |
| 134 # Prime factors of P-1 and their exponents: |
| 135 F = [(2,3,7), (2,3,5), (2,3)] |
| 136 E = [(25,2,1), (27,1,1), (26,3)] |
| 137 |
| 138 # Maximum transform length for 2**n. Above that only 3 * 2**24 or |
| 139 # 3 * 2**25 are possible. |
| 140 MPD_MAXTRANSFORM_2N = 2**25 |
| 141 |
| 142 |
| 143 # Limits in the terminology of Pollard's paper: |
| 144 m2 = (MPD_MAXTRANSFORM_2N * 3) // 2 # Maximum length of the smaller array. |
| 145 M1 = M2 = RADIX-1 # Maximum value per single word. |
| 146 L = m2 * M1 * M2 |
| 147 P[0] * P[1] * P[2] > 2 * L |
| 148 |
| 149 |
| 150 # Primitive roots of F(P1), F(P2) and F(P3): |
| 151 w = [5, 31, 13] |
| 152 |
| 153 # The primitive roots are correct: |
| 154 for i in range(3): |
| 155 if not is_primitive_root(w[i], P[i], F[i], E[i]): |
| 156 print("FAIL") |
| 157 |
| 158 |
| 159 # ====================================================================== |
| 160 # Example transform using a single prime |
| 161 # ====================================================================== |
| 162 |
| 163 def ntt(lst, dir): |
| 164 """Perform a transform on the elements of lst. len(lst) must |
| 165 be 2**n or 3 * 2**n, where n <= 25. This is the slow DFT.""" |
| 166 p = 2113929217 # prime |
| 167 d = len(lst) # transform length |
| 168 d_prime = pow(d, (p-2), p) # inverse of d |
| 169 xi = (p-1)//d |
| 170 w = 5 # primitive root of F(p) |
| 171 r = pow(w, xi, p) # primitive root of the subfield |
| 172 r_prime = pow(w, (p-1-xi), p) # inverse of r |
| 173 if dir == 1: # forward transform |
| 174 a = lst # input array |
| 175 A = [0] * d # transformed values |
| 176 for i in range(d): |
| 177 s = 0 |
| 178 for j in range(d): |
| 179 s += a[j] * pow(r, i*j, p) |
| 180 A[i] = s % p |
| 181 return A |
| 182 elif dir == -1: # backward transform |
| 183 A = lst # input array |
| 184 a = [0] * d # transformed values |
| 185 for j in range(d): |
| 186 s = 0 |
| 187 for i in range(d): |
| 188 s += A[i] * pow(r_prime, i*j, p) |
| 189 a[j] = (d_prime * s) % p |
| 190 return a |
| 191 |
| 192 def ntt_convolute(a, b): |
| 193 """convolute arrays a and b.""" |
| 194 assert(len(a) == len(b)) |
| 195 x = ntt(a, 1) |
| 196 y = ntt(b, 1) |
| 197 for i in range(len(a)): |
| 198 y[i] = y[i] * x[i] |
| 199 r = ntt(y, -1) |
| 200 return r |
| 201 |
| 202 |
| 203 # Example: Two arrays representing 21 and 81 in little-endian: |
| 204 a = [1, 2, 0, 0] |
| 205 b = [1, 8, 0, 0] |
| 206 |
| 207 assert(ntt_convolute(a, b) == [1, 10, 16, 0]) |
| 208 assert(21 * 81 == (1*10**0 + 10*10**1 + 16*10**2 + 0*10**3)) |
| 209 |
| 210 |
| 211 |
| 212 |
OLD | NEW |