• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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