Rietveld Code Review Tool
Help | Bug tracker | Discussion group | Source code | Sign in
(13617)

Unified Diff: Modules/_decimal/libmpdec/literature/fnt.py

Issue 7652: Merge C version of decimal into py3k.
Patch Set: Created 7 years, 7 months ago
Use n/p to move between diff chunks; N/P to move between comments. Please Sign in to add in-line comments.
Jump to:
View side-by-side diff with in-line comments
Download patch
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Modules/_decimal/libmpdec/literature/fnt.py Sat Mar 10 18:12:20 2012 +0100
@@ -0,0 +1,212 @@
+#
+# Copyright (c) 2008-2010 Stefan Krah. All rights reserved.
+#
+# 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 copyright
+# notice, this list of conditions and the following disclaimer.
+#
+# 2. Redistributions in binary form must reproduce the above copyright
+# 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 AUTHOR 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 AUTHOR 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).
+#
+def prod(F, E):
+ """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."""
+ x = 1
+ for y, z in zip(F, E):
+ x *= y**z
+ return x
+
+def is_primitive_root(r, p, factors, exponents):
+ """Check if r is a primitive root of F(p)."""
+ 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
+
+
+# =================================================================
+# Constants and limits for the 64-bit version
+# =================================================================
+
+RADIX = 10**19
+
+# Primes P1, P2 and P3:
+P = [2**64-2**32+1, 2**64-2**34+1, 2**64-2**40+1]
+
+# 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).
+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)]
+
+# Prime factors of P-1 and their exponents:
+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)]
+
+# Maximum transform length for 2**n. Above that only 3 * 2**31
+# or 3 * 2**32 are possible.
+MPD_MAXTRANSFORM_2N = 2**32
+
+
+# Limits in the terminology of Pollard's paper:
+m2 = (MPD_MAXTRANSFORM_2N * 3) // 2 # Maximum length of the smaller array.
+M1 = M2 = RADIX-1 # Maximum value per single word.
+L = m2 * M1 * M2
+P[0] * P[1] * P[2] > 2 * L
+
+
+# Primitive roots of F(P1), F(P2) and F(P3):
+w = [7, 10, 19]
+
+# The primitive roots are correct:
+for i in range(3):
+ if not is_primitive_root(w[i], P[i], F[i], E[i]):
+ print("FAIL")
+
+
+# =================================================================
+# Constants and limits for the 32-bit version
+# =================================================================
+
+RADIX = 10**9
+
+# Primes P1, P2 and P3:
+P = [2113929217, 2013265921, 1811939329]
+
+# P-1, highly composite. All D = P-1 are divisible by 3 * 2**25,
+# allowing for transform lengths up to 3 * 2**25 words.
+D = [2**25 * 3**2 * 7,
+ 2**27 * 3 * 5,
+ 2**26 * 3**3]
+
+# Prime factors of P-1 and their exponents:
+F = [(2,3,7), (2,3,5), (2,3)]
+E = [(25,2,1), (27,1,1), (26,3)]
+
+# Maximum transform length for 2**n. Above that only 3 * 2**24 or
+# 3 * 2**25 are possible.
+MPD_MAXTRANSFORM_2N = 2**25
+
+
+# Limits in the terminology of Pollard's paper:
+m2 = (MPD_MAXTRANSFORM_2N * 3) // 2 # Maximum length of the smaller array.
+M1 = M2 = RADIX-1 # Maximum value per single word.
+L = m2 * M1 * M2
+P[0] * P[1] * P[2] > 2 * L
+
+
+# Primitive roots of F(P1), F(P2) and F(P3):
+w = [5, 31, 13]
+
+# The primitive roots are correct:
+for i in range(3):
+ if not is_primitive_root(w[i], P[i], F[i], E[i]):
+ print("FAIL")
+
+
+# ======================================================================
+# Example transform using a single prime
+# ======================================================================
+
+def ntt(lst, dir):
+ """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."""
+ p = 2113929217 # prime
+ d = len(lst) # transform length
+ d_prime = pow(d, (p-2), p) # inverse of d
+ xi = (p-1)//d
+ w = 5 # primitive root of F(p)
+ r = pow(w, xi, p) # primitive root of the subfield
+ r_prime = pow(w, (p-1-xi), p) # inverse of r
+ if dir == 1: # forward transform
+ a = lst # input array
+ A = [0] * d # transformed values
+ 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: # backward transform
+ A = lst # input array
+ a = [0] * d # transformed values
+ 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):
+ """convolute arrays a and 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
+
+
+# Example: Two arrays representing 21 and 81 in little-endian:
+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))
+
+
+
+
« no previous file with comments | « Modules/_decimal/libmpdec/literature/bignum.txt ('k') | Modules/_decimal/libmpdec/literature/matrix-transform.txt » ('j') | no next file with comments »

RSS Feeds Recent Issues | This issue
This is Rietveld 894c83f36cb7+