1# 2# Copyright (c) 2008-2016 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# 57def 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 65def 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 82RADIX = 10**19 83 84# Primes P1, P2 and P3: 85P = [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). 90D = [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: 95F = [(2,3,5,17,257,65537), (2,3,7,11,31,151,331), (2,3,5,7,13,17,241)] 96E = [(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. 100MPD_MAXTRANSFORM_2N = 2**32 101 102 103# Limits in the terminology of Pollard's paper: 104m2 = (MPD_MAXTRANSFORM_2N * 3) // 2 # Maximum length of the smaller array. 105M1 = M2 = RADIX-1 # Maximum value per single word. 106L = m2 * M1 * M2 107P[0] * P[1] * P[2] > 2 * L 108 109 110# Primitive roots of F(P1), F(P2) and F(P3): 111w = [7, 10, 19] 112 113# The primitive roots are correct: 114for 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 123RADIX = 10**9 124 125# Primes P1, P2 and P3: 126P = [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. 130D = [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: 135F = [(2,3,7), (2,3,5), (2,3)] 136E = [(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. 140MPD_MAXTRANSFORM_2N = 2**25 141 142 143# Limits in the terminology of Pollard's paper: 144m2 = (MPD_MAXTRANSFORM_2N * 3) // 2 # Maximum length of the smaller array. 145M1 = M2 = RADIX-1 # Maximum value per single word. 146L = m2 * M1 * M2 147P[0] * P[1] * P[2] > 2 * L 148 149 150# Primitive roots of F(P1), F(P2) and F(P3): 151w = [5, 31, 13] 152 153# The primitive roots are correct: 154for 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 163def 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 192def 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: 204a = [1, 2, 0, 0] 205b = [1, 8, 0, 0] 206 207assert(ntt_convolute(a, b) == [1, 10, 16, 0]) 208assert(21 * 81 == (1*10**0 + 10*10**1 + 16*10**2 + 0*10**3)) 209