ยปCore Development>Code coverage>Modules/_decimal/libmpdec/literature/fnt.py

Python code coverage for Modules/_decimal/libmpdec/literature/fnt.py

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