c 20082020 stefan krah redistribution and use in source and binary forms with or without modification are permitted provided that the following conditions are met 1 redistributions of source code must retain the above notice this list of conditions and the following disclaimer 2 redistributions in binary form must reproduce the above notice this list of conditions and the following disclaimer in the documentation andor other materials provided with the distribution this software is provided by the and contributors as is and any express or implied warranties including but not limited to the implied warranties of merchantability and fitness for a particular purpose are disclaimed in no event shall the or contributors be liable for any direct indirect incidental special exemplary or consequential damages including but not limited to procurement of substitute goods or services loss of use data or profits or business interruption however caused and on any theory of liability whether in contract strict liability or tort including negligence or otherwise arising in any way out of the use of this software even if advised of the possibility of such damage this file lists and checks some of the constants and limits used in libmpdec s number theoretic transform at the end of the file there is an example function for the plain dft transform number theoretic transforms are done in subfields of fp pi are the primes di pi 1 are highly composite and wi are the respective primitive roots of fp the strategy is to convolute two coefficients modulo all three primes then use the chinese remainder theorem on the three result arrays to recover the result in the usual base radix form primitive roots verify primitive roots for a prime field r is a primitive root if and only if for all prime factors f of p1 rp1f 1 mod p check that the factorization of p1 is correct f is the list of factors of p1 e lists the number of occurrences of each factor x 1 for y z in zipf e x yz return x def isprimitiverootr p factors exponents constants and limits for the 64bit version primes p1 p2 and p3 p1 highly composite the transform length d is variable and must divide d p1 since all d are divisible by 3 232 transform lengths can be 2n or 3 2n where n 32 prime factors of p1 and their exponents maximum transform length for 2n above that only 3 231 or 3 232 are possible limits in the terminology of pollard s paper primitive roots of fp1 fp2 and fp3 the primitive roots are correct constants and limits for the 32bit version primes p1 p2 and p3 p1 highly composite all d p1 are divisible by 3 225 allowing for transform lengths up to 3 225 words prime factors of p1 and their exponents maximum transform length for 2n above that only 3 224 or 3 225 are possible limits in the terminology of pollard s paper primitive roots of fp1 fp2 and fp3 the primitive roots are correct example transform using a single prime perform a transform on the elements of lst lenlst must be 2n or 3 2n where n 25 this is the slow dft p 2113929217 prime d lenlst transform length dprime powd p2 p inverse of d xi p1d w 5 primitive root of fp r poww xi p primitive root of the subfield rprime poww p1xi p inverse of r if dir 1 forward transform a lst input array a 0 d transformed values for i in ranged s 0 for j in ranged s aj powr ij p ai s p return a elif dir 1 backward transform a lst input array a 0 d transformed values for j in ranged s 0 for i in ranged s ai powrprime ij p aj dprime s p return a def nttconvolutea b example two arrays representing 21 and 81 in littleendian c 2008 2020 stefan krah redistribution and use in source and binary forms with or without modification are permitted provided that the following conditions are met 1 redistributions of source code must retain the above notice this list of conditions and the following disclaimer 2 redistributions in binary form must reproduce the above notice this list of conditions and the following disclaimer in the documentation and or other materials provided with the distribution this software is provided by the and contributors as is and any express or implied warranties including but not limited to the implied warranties of merchantability and fitness for a particular purpose are disclaimed in no event shall the or contributors be liable for any direct indirect incidental special exemplary or consequential damages including but not limited to procurement of substitute goods or services loss of use data or profits or business interruption however caused and on any theory of liability whether in contract strict liability or tort including negligence or otherwise arising in any way out of the use of this software even if advised of the possibility of such damage this file lists and checks some of the constants and limits used in libmpdec s number theoretic transform at the end of the file there is an example function for the plain dft transform number theoretic transforms are done in subfields of f p p i are the primes d i p i 1 are highly composite and w i are the respective primitive roots of f p the strategy is to convolute two coefficients modulo all three primes then use the chinese remainder theorem on the three result arrays to recover the result in the usual base radix form primitive roots verify primitive roots for a prime field r is a primitive root if and only if for all prime factors f of p 1 r p 1 f 1 mod p check that the factorization of p 1 is correct f is the list of factors of p 1 e lists the number of occurrences of each factor check if r is a primitive root of f p constants and limits for the 64 bit version primes p1 p2 and p3 p 1 highly composite the transform length d is variable and must divide d p 1 since all d are divisible by 3 2 32 transform lengths can be 2 n or 3 2 n where n 32 prime factors of p 1 and their exponents maximum transform length for 2 n above that only 3 2 31 or 3 2 32 are possible limits in the terminology of pollard s paper maximum length of the smaller array maximum value per single word primitive roots of f p1 f p2 and f p3 the primitive roots are correct constants and limits for the 32 bit version primes p1 p2 and p3 p 1 highly composite all d p 1 are divisible by 3 2 25 allowing for transform lengths up to 3 2 25 words prime factors of p 1 and their exponents maximum transform length for 2 n above that only 3 2 24 or 3 2 25 are possible limits in the terminology of pollard s paper maximum length of the smaller array maximum value per single word primitive roots of f p1 f p2 and f p3 the primitive roots are correct example transform using a single prime perform a transform on the elements of lst len lst must be 2 n or 3 2 n where n 25 this is the slow dft prime transform length inverse of d primitive root of f p primitive root of the subfield inverse of r forward transform input array transformed values backward transform input array transformed values convolute arrays a and b example two arrays representing 21 and 81 in little endian
def prod(F, E):
    
    x = 1
    for y, z in zip(F, E):
        x *= y**z
    return x

def is_primitive_root(r, p, factors, exponents):
    
    if p != prod(factors, exponents) + 1:
        return False
    for f in factors:
        q, control = divmod(p-1, f)
        if control != 0:
            return False
        if pow(r, q, p) == 1:
            return False
    return True






RADIX = 10**19


P = [2**64-2**32+1, 2**64-2**34+1, 2**64-2**40+1]




D = [2**32 * 3    * (5 * 17 * 257 * 65537),
     2**34 * 3**2 * (7 * 11 * 31 * 151 * 331),
     2**40 * 3**2 * (5 * 7 * 13 * 17 * 241)]


F = [(2,3,5,17,257,65537), (2,3,7,11,31,151,331), (2,3,5,7,13,17,241)]
E = [(32,1,1,1,1,1), (34,2,1,1,1,1,1), (40,2,1,1,1,1,1)]



MPD_MAXTRANSFORM_2N = 2**32



m2 = (MPD_MAXTRANSFORM_2N * 3) // 2 
M1 = M2 = RADIX-1                   
L = m2 * M1 * M2
P[0] * P[1] * P[2] > 2 * L



w = [7, 10, 19]


for i in range(3):
    if not is_primitive_root(w[i], P[i], F[i], E[i]):
        print("FAIL")






RADIX = 10**9


P = [2113929217, 2013265921, 1811939329]



D = [2**25 * 3**2 * 7,
     2**27 * 3    * 5,
     2**26 * 3**3]


F = [(2,3,7), (2,3,5), (2,3)]
E = [(25,2,1), (27,1,1), (26,3)]



MPD_MAXTRANSFORM_2N = 2**25



m2 = (MPD_MAXTRANSFORM_2N * 3) // 2 
M1 = M2 = RADIX-1                   
L = m2 * M1 * M2
P[0] * P[1] * P[2] > 2 * L



w = [5, 31, 13]


for i in range(3):
    if not is_primitive_root(w[i], P[i], F[i], E[i]):
        print("FAIL")






def ntt(lst, dir):
    
    p = 2113929217             
    d = len(lst)               
    d_prime = pow(d, (p-2), p) 
    xi = (p-1)//d
    w = 5                         
    r = pow(w, xi, p)             
    r_prime = pow(w, (p-1-xi), p) 
    if dir == 1:      
        a = lst       
        A = [0] * d   
        for i in range(d):
            s = 0
            for j in range(d):
                s += a[j] * pow(r, i*j, p)
            A[i] = s % p
        return A
    elif dir == -1: 
        A = lst     
        a = [0] * d 
        for j in range(d):
            s = 0
            for i in range(d):
                s += A[i] * pow(r_prime, i*j, p)
            a[j] = (d_prime * s) % p
        return a

def ntt_convolute(a, b):
    
    assert(len(a) == len(b))
    x = ntt(a, 1)
    y = ntt(b, 1)
    for i in range(len(a)):
        y[i] = y[i] * x[i]
    r = ntt(y, -1)
    return r



a = [1, 2, 0, 0]
b = [1, 8, 0, 0]

assert(ntt_convolute(a, b) == [1,        10,        16,        0])
assert(21 * 81             == (1*10**0 + 10*10**1 + 16*10**2 + 0*10**3))