1 #
2 # Copyright (c) 2008-2020 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))