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

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