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