From 99975e18b0733c08e109b02b65f497f7f9de6866 Mon Sep 17 00:00:00 2001 From: sidhantnagpal Date: Fri, 18 May 2018 12:26:57 +0530 Subject: [PATCH 01/11] Add discrete module, fft --- setup.py | 2 + sympy/discrete/__init__.py | 9 ++++ sympy/discrete/tests/__init__.py | 0 sympy/discrete/tests/test_transforms.py | 25 +++++++++ sympy/discrete/transforms.py | 68 +++++++++++++++++++++++++ 5 files changed, 104 insertions(+) create mode 100644 sympy/discrete/__init__.py create mode 100644 sympy/discrete/tests/__init__.py create mode 100644 sympy/discrete/tests/test_transforms.py create mode 100644 sympy/discrete/transforms.py diff --git a/setup.py b/setup.py index 196a6ba655ce..60e8c06f64c1 100755 --- a/setup.py +++ b/setup.py @@ -89,6 +89,7 @@ 'sympy.core', 'sympy.core.benchmarks', 'sympy.crypto', + 'sympy.discrete', 'sympy.deprecated', 'sympy.diffgeom', 'sympy.external', @@ -310,6 +311,7 @@ def run(self): 'sympy.concrete.tests', 'sympy.core.tests', 'sympy.crypto.tests', + 'sympy.discrete.tests', 'sympy.deprecated.tests', 'sympy.diffgeom.tests', 'sympy.external.tests', diff --git a/sympy/discrete/__init__.py b/sympy/discrete/__init__.py new file mode 100644 index 000000000000..71891af2167a --- /dev/null +++ b/sympy/discrete/__init__.py @@ -0,0 +1,9 @@ +"""A module containing discrete functions. + +Transforms - fft, ifft, ntt, intt, fwht, ifwht, fmt, fzt +Convolution - conv, conv_xor, conv_and, conv_or, conv_sub, conv_sup +Recurrence Evaluation - reval_hcc +""" + + +from .transforms import (fft, ifft) diff --git a/sympy/discrete/tests/__init__.py b/sympy/discrete/tests/__init__.py new file mode 100644 index 000000000000..e69de29bb2d1 diff --git a/sympy/discrete/tests/test_transforms.py b/sympy/discrete/tests/test_transforms.py new file mode 100644 index 000000000000..14c445c741f5 --- /dev/null +++ b/sympy/discrete/tests/test_transforms.py @@ -0,0 +1,25 @@ +from __future__ import print_function, division + +from sympy import sqrt +from sympy.core import S, Symbol, I +from sympy.core.compatibility import range +from sympy.discrete import fft, ifft +from sympy.utilities.pytest import raises + +def test_fft_ifft(): + ls = list(range(6)) + fls = [15, -7*sqrt(2)/2 - 4 - sqrt(2)*I/2 + 2*I, 2 + 3*I, \ + -4 + 7*sqrt(2)/2 - 2*I - sqrt(2)*I/2, -3, \ + -4 + 7*sqrt(2)/2 + sqrt(2)*I/2 + 2*I, \ + 2 - 3*I, -7*sqrt(2)/2 - 4 - 2*I + sqrt(2)*I/2] + assert fft(ls) == fls + assert ifft(fls) == ls + [S.Zero]*2 + + ls = [1 + 2*I, 3 + 4*I, 5 + 6*I] + ifls = [S(9)/4 + 3*I, -7*I/4, S(3)/4 + I, -2 - I/4] + assert ifft(ls) == ifls + assert fft(ifls) == ls + [S.Zero] + + x = Symbol('x', real=True) + raises(TypeError, lambda: fft(x)) + raises(ValueError, lambda: ifft([x, 2*x, 3*x**2, 4*x**3])) diff --git a/sympy/discrete/transforms.py b/sympy/discrete/transforms.py new file mode 100644 index 000000000000..f76fd11f93e5 --- /dev/null +++ b/sympy/discrete/transforms.py @@ -0,0 +1,68 @@ +""" +Discrete Fourier Transform, Walsh Hadamard Transform, +Number Theoretic Transform, Zeta Transform, Mobius Transform +""" +from __future__ import print_function, division + +from sympy.core import S, Symbol, sympify +from sympy.core.compatibility import as_int, range +from sympy.core.function import expand, expand_mul +from sympy.core.numbers import pi, I +from sympy.functions.elementary.exponential import exp +from sympy.ntheory import isprime, primitive_root +from sympy.utilities.iterables import ibin + +def _fourier_transform(seq, inverse=False): + """ + Performs the DFT in complex field, by padding zeros as + radix 2 FFT, requires the number of sample points to be + a power of 2. + """ + + if not hasattr(seq, '__iter__'): + raise TypeError("Expected a sequence of numeric coefficients " + + "for Fourier Transform") + + a = sympify(seq) + if any(x.has(Symbol) for x in a): + raise ValueError("Expected non-symbolic coefficients") + + n = len(a) + if n < 2: + return + + while n&(n - 1): + n += n&-n + + a += [S.Zero]*(n - len(a)) + b = n.bit_length() - 1 + for i in range(1, n): + j = int(ibin(i, b, str=True)[::-1], 2) + if i < j: + a[i], a[j] = a[j], a[i] + + ang = -2*pi/n if inverse else 2*pi/n + w = [exp(I*ang*i).expand(complex=True) for i in range(n // 2)] + + h = 2 + while h <= n: + hf, ut = h // 2, n // h + for i in range(0, n, h): + for j in range(hf): + u, v = a[i + j], expand_mul(a[i + j + hf]*w[ut * j]) + a[i + j], a[i + j + hf] = u + v, u - v + h *= 2 + + if inverse: + for i in range(n): + a[i] /= n + + return a + + +def fft(seq): + return _fourier_transform(seq) + + +def ifft(seq): + return _fourier_transform(seq, inverse=True) From 5d2bef1f7bf5db725113226fa67006dc3deaaf8e Mon Sep 17 00:00:00 2001 From: sidhantnagpal Date: Fri, 18 May 2018 13:16:14 +0530 Subject: [PATCH 02/11] Add discrete to __init__ --- setup.py | 4 ++-- sympy/__init__.py | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index 60e8c06f64c1..a4327e32eed1 100755 --- a/setup.py +++ b/setup.py @@ -89,9 +89,9 @@ 'sympy.core', 'sympy.core.benchmarks', 'sympy.crypto', - 'sympy.discrete', 'sympy.deprecated', 'sympy.diffgeom', + 'sympy.discrete', 'sympy.external', 'sympy.functions', 'sympy.functions.combinatorial', @@ -311,9 +311,9 @@ def run(self): 'sympy.concrete.tests', 'sympy.core.tests', 'sympy.crypto.tests', - 'sympy.discrete.tests', 'sympy.deprecated.tests', 'sympy.diffgeom.tests', + 'sympy.discrete.tests', 'sympy.external.tests', 'sympy.functions.combinatorial.tests', 'sympy.functions.elementary.tests', diff --git a/sympy/__init__.py b/sympy/__init__.py index 21936f61b96a..982b6fa66c88 100644 --- a/sympy/__init__.py +++ b/sympy/__init__.py @@ -62,6 +62,7 @@ def __sympy_debug(): from .functions import * from .ntheory import * from .concrete import * +from .discrete import * from .simplify import * from .sets import * from .solvers import * From 8773627278ebf18e77b34b2cde5adac307316c18 Mon Sep 17 00:00:00 2001 From: sidhantnagpal Date: Sat, 19 May 2018 09:02:24 +0530 Subject: [PATCH 03/11] add docstring, tests --- sympy/discrete/tests/test_transforms.py | 10 ++- sympy/discrete/transforms.py | 85 +++++++++++++++++++++---- 2 files changed, 78 insertions(+), 17 deletions(-) diff --git a/sympy/discrete/tests/test_transforms.py b/sympy/discrete/tests/test_transforms.py index 14c445c741f5..121fe6312b5c 100644 --- a/sympy/discrete/tests/test_transforms.py +++ b/sympy/discrete/tests/test_transforms.py @@ -7,16 +7,20 @@ from sympy.utilities.pytest import raises def test_fft_ifft(): + assert all(tf(ls) == ls for tf in (fft, ifft) for ls in ([], [S(5)/3])) + ls = list(range(6)) - fls = [15, -7*sqrt(2)/2 - 4 - sqrt(2)*I/2 + 2*I, 2 + 3*I, \ - -4 + 7*sqrt(2)/2 - 2*I - sqrt(2)*I/2, -3, \ - -4 + 7*sqrt(2)/2 + sqrt(2)*I/2 + 2*I, \ + fls = [15, -7*sqrt(2)/2 - 4 - sqrt(2)*I/2 + 2*I, 2 + 3*I, + -4 + 7*sqrt(2)/2 - 2*I - sqrt(2)*I/2, -3, + -4 + 7*sqrt(2)/2 + sqrt(2)*I/2 + 2*I, 2 - 3*I, -7*sqrt(2)/2 - 4 - 2*I + sqrt(2)*I/2] + assert fft(ls) == fls assert ifft(fls) == ls + [S.Zero]*2 ls = [1 + 2*I, 3 + 4*I, 5 + 6*I] ifls = [S(9)/4 + 3*I, -7*I/4, S(3)/4 + I, -2 - I/4] + assert ifft(ls) == ifls assert fft(ifls) == ls + [S.Zero] diff --git a/sympy/discrete/transforms.py b/sympy/discrete/transforms.py index f76fd11f93e5..95c5d01dd23f 100644 --- a/sympy/discrete/transforms.py +++ b/sympy/discrete/transforms.py @@ -5,44 +5,41 @@ from __future__ import print_function, division from sympy.core import S, Symbol, sympify -from sympy.core.compatibility import as_int, range +from sympy.core.compatibility import as_int, range, iterable from sympy.core.function import expand, expand_mul from sympy.core.numbers import pi, I -from sympy.functions.elementary.exponential import exp +from sympy.functions.elementary.trigonometric import sin, cos from sympy.ntheory import isprime, primitive_root from sympy.utilities.iterables import ibin def _fourier_transform(seq, inverse=False): - """ - Performs the DFT in complex field, by padding zeros as - radix 2 FFT, requires the number of sample points to be - a power of 2. - """ + """Utility function for the discrete Fourier transform (DFT)""" - if not hasattr(seq, '__iter__'): + if not iterable(seq): raise TypeError("Expected a sequence of numeric coefficients " + "for Fourier Transform") - a = sympify(seq) + a = [sympify(arg) for arg in seq] if any(x.has(Symbol) for x in a): raise ValueError("Expected non-symbolic coefficients") n = len(a) if n < 2: - return + return a - while n&(n - 1): - n += n&-n + b = n.bit_length() - 1 + if n&(n - 1): # not a power of 2 + b += 1 + n = 2**b a += [S.Zero]*(n - len(a)) - b = n.bit_length() - 1 for i in range(1, n): j = int(ibin(i, b, str=True)[::-1], 2) if i < j: a[i], a[j] = a[j], a[i] ang = -2*pi/n if inverse else 2*pi/n - w = [exp(I*ang*i).expand(complex=True) for i in range(n // 2)] + w = [cos(ang*i) + I*sin(ang*i) for i in range(n // 2)] h = 2 while h <= n: @@ -61,8 +58,68 @@ def _fourier_transform(seq, inverse=False): def fft(seq): + r""" + Performs the discrete Fourier transform (DFT) in the complex domain. + + The sequence is automatically padded to the right with zeros, as the + radix 2 FFT requires the number of sample points to be a power of 2. + + Examples + ======== + + >>> from sympy import fft, ifft + >>> fft([1, 2, 3, 4]) + [10, -2 - 2*I, -2, -2 + 2*I] + >>> ifft(_) + [1, 2, 3, 4] + >>> fft([5]) + [5] + + References + ========== + + .. [1] https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm + .. [2] http://mathworld.wolfram.com/FastFourierTransform.html + + See Also + ======== + + ifft + + """ + return _fourier_transform(seq) def ifft(seq): + r""" + Performs the inverse discrete Fourier transform (IDFT) in the complex domain. + + The sequence is automatically padded to the right with zeros, as the + radix 2 FFT requires the number of sample points to be a power of 2. + + Examples + ======== + + >>> from sympy import fft, ifft + >>> ifft([1, 2, 3, 4]) + [5/2, -1/2 + I/2, -1/2, -1/2 - I/2] + >>> fft(_) + [1, 2, 3, 4] + >>> ifft([7]) + [7] + + References + ========== + + .. [1] https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm + .. [2] http://mathworld.wolfram.com/FastFourierTransform.html + + See Also + ======== + + fft + + """ + return _fourier_transform(seq, inverse=True) From a2698c3aa5081d0730d261c73dfe0666703c8173 Mon Sep 17 00:00:00 2001 From: sidhantnagpal Date: Sun, 20 May 2018 15:55:47 +0530 Subject: [PATCH 04/11] Support symbolic, numeric fft --- sympy/discrete/__init__.py | 2 +- sympy/discrete/tests/test_transforms.py | 4 +- sympy/discrete/transforms.py | 75 ++++++++++++------------- 3 files changed, 40 insertions(+), 41 deletions(-) diff --git a/sympy/discrete/__init__.py b/sympy/discrete/__init__.py index 71891af2167a..7410739b8094 100644 --- a/sympy/discrete/__init__.py +++ b/sympy/discrete/__init__.py @@ -1,6 +1,6 @@ """A module containing discrete functions. -Transforms - fft, ifft, ntt, intt, fwht, ifwht, fmt, fzt +Transforms - fft, ifft, ntt, intt, fwht, ifwht, fzt, ifzt, fmt, ifmt Convolution - conv, conv_xor, conv_and, conv_or, conv_sub, conv_sup Recurrence Evaluation - reval_hcc """ diff --git a/sympy/discrete/tests/test_transforms.py b/sympy/discrete/tests/test_transforms.py index 121fe6312b5c..d6654373330e 100644 --- a/sympy/discrete/tests/test_transforms.py +++ b/sympy/discrete/tests/test_transforms.py @@ -6,8 +6,10 @@ from sympy.discrete import fft, ifft from sympy.utilities.pytest import raises + def test_fft_ifft(): - assert all(tf(ls) == ls for tf in (fft, ifft) for ls in ([], [S(5)/3])) + assert all(tf(ls) == ls for tf in (fft, ifft) \ + for ls in ([], [S(5)/3])) ls = list(range(6)) fls = [15, -7*sqrt(2)/2 - 4 - sqrt(2)*I/2 + 2*I, 2 + 3*I, diff --git a/sympy/discrete/transforms.py b/sympy/discrete/transforms.py index 95c5d01dd23f..ff45aa65d5bf 100644 --- a/sympy/discrete/transforms.py +++ b/sympy/discrete/transforms.py @@ -8,12 +8,20 @@ from sympy.core.compatibility import as_int, range, iterable from sympy.core.function import expand, expand_mul from sympy.core.numbers import pi, I +from sympy.functions.elementary.exponential import exp from sympy.functions.elementary.trigonometric import sin, cos from sympy.ntheory import isprime, primitive_root from sympy.utilities.iterables import ibin -def _fourier_transform(seq, inverse=False): - """Utility function for the discrete Fourier transform (DFT)""" + +#----------------------------------------------------------------------------# +# # +# Discrete Fourier Transform # +# # +#----------------------------------------------------------------------------# + +def _fourier_transform(seq, symbolic, inverse=False): + """Utility function for the Discrete Fourier Transform (DFT)""" if not iterable(seq): raise TypeError("Expected a sequence of numeric coefficients " + @@ -39,7 +47,11 @@ def _fourier_transform(seq, inverse=False): a[i], a[j] = a[j], a[i] ang = -2*pi/n if inverse else 2*pi/n - w = [cos(ang*i) + I*sin(ang*i) for i in range(n // 2)] + w = None + if symbolic: + w = [cos(ang*i) + I*sin(ang*i) for i in range(n // 2)] + else: + w = [exp(I*ang*i).evalf() for i in range(n // 2)] h = 2 while h <= n: @@ -57,55 +69,39 @@ def _fourier_transform(seq, inverse=False): return a -def fft(seq): +def fft(seq, symbolic=True): r""" - Performs the discrete Fourier transform (DFT) in the complex domain. + Performs the Discrete Fourier Transform (DFT) in the complex domain. The sequence is automatically padded to the right with zeros, as the radix 2 FFT requires the number of sample points to be a power of 2. + Parameters + ========== + + seq : iterable + The sequence on which DFT is to be applied. + symbolic : bool + Determines if DFT is to be performed using symbolic values or + numerical values. + Examples ======== >>> from sympy import fft, ifft + >>> fft([1, 2, 3, 4]) [10, -2 - 2*I, -2, -2 + 2*I] >>> ifft(_) [1, 2, 3, 4] - >>> fft([5]) - [5] - - References - ========== - - .. [1] https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm - .. [2] http://mathworld.wolfram.com/FastFourierTransform.html - - See Also - ======== - - ifft - - """ - return _fourier_transform(seq) - - -def ifft(seq): - r""" - Performs the inverse discrete Fourier transform (IDFT) in the complex domain. - - The sequence is automatically padded to the right with zeros, as the - radix 2 FFT requires the number of sample points to be a power of 2. - - Examples - ======== - - >>> from sympy import fft, ifft >>> ifft([1, 2, 3, 4]) [5/2, -1/2 + I/2, -1/2, -1/2 - I/2] >>> fft(_) [1, 2, 3, 4] + + >>> fft([5]) + [5] >>> ifft([7]) [7] @@ -115,11 +111,12 @@ def ifft(seq): .. [1] https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm .. [2] http://mathworld.wolfram.com/FastFourierTransform.html - See Also - ======== + """ - fft + return _fourier_transform(seq, symbolic=symbolic) - """ - return _fourier_transform(seq, inverse=True) +def ifft(seq, symbolic=True): + return _fourier_transform(seq, symbolic=symbolic, inverse=True) + +ifft.__doc__ = fft.__doc__ From 5fa296b5f942d2ce3fed46a59464ffa534991545 Mon Sep 17 00:00:00 2001 From: sidhantnagpal Date: Wed, 23 May 2018 06:54:25 +0530 Subject: [PATCH 05/11] Add NTT, use evalf for ang in FFT --- sympy/discrete/__init__.py | 2 +- sympy/discrete/tests/test_transforms.py | 28 +++++- sympy/discrete/transforms.py | 119 ++++++++++++++++++++++-- 3 files changed, 140 insertions(+), 9 deletions(-) diff --git a/sympy/discrete/__init__.py b/sympy/discrete/__init__.py index 7410739b8094..39b406ab7269 100644 --- a/sympy/discrete/__init__.py +++ b/sympy/discrete/__init__.py @@ -6,4 +6,4 @@ """ -from .transforms import (fft, ifft) +from .transforms import (fft, ifft, ntt, intt) diff --git a/sympy/discrete/tests/test_transforms.py b/sympy/discrete/tests/test_transforms.py index d6654373330e..bde733858fd9 100644 --- a/sympy/discrete/tests/test_transforms.py +++ b/sympy/discrete/tests/test_transforms.py @@ -3,7 +3,7 @@ from sympy import sqrt from sympy.core import S, Symbol, I from sympy.core.compatibility import range -from sympy.discrete import fft, ifft +from sympy.discrete import fft, ifft, ntt, intt from sympy.utilities.pytest import raises @@ -29,3 +29,29 @@ def test_fft_ifft(): x = Symbol('x', real=True) raises(TypeError, lambda: fft(x)) raises(ValueError, lambda: ifft([x, 2*x, 3*x**2, 4*x**3])) + + +def test_ntt_intt(): + # prime modulo of the form (m*2**k + 1), sequence length should be a divisor of 2**k + p = 7*17*2**23 + 1 + q = 2*500000003 + 1 # prime modulo not of the required form (for length > 1) + r = 2*3*5*7 # composite modulo + + assert all(tf(ls, p) == ls for tf in (ntt, intt) \ + for ls in ([], [5])) + + ls = list(range(6)) + nls = [15, 801133602, 738493201, 334102277, 998244350, 849020224, 259751156, 12232587] + assert ntt(ls, p) == nls + assert intt(nls, p) == ls + [0]*2 + + ls = [1 + 2*I, 3 + 4*I, 5 + 6*I] + x = Symbol('x', integer=True) + raises(TypeError, lambda: ntt(x, p)) + raises(ValueError, lambda: intt([x, 2*x, 3*x**2, 4*x**3], p)) + raises(ValueError, lambda: intt(ls, p)) + raises(ValueError, lambda: ntt([1.2, 2.1, 3.5], p)) + raises(ValueError, lambda: ntt([3, 5, 6], q)) + raises(ValueError, lambda: ntt([4, 5, 7], r)) + + assert ntt([1.0, 2.0, 3.0], p) == ntt([1, 2, 3], p) diff --git a/sympy/discrete/transforms.py b/sympy/discrete/transforms.py index ff45aa65d5bf..e7c0b86b1221 100644 --- a/sympy/discrete/transforms.py +++ b/sympy/discrete/transforms.py @@ -1,6 +1,6 @@ """ -Discrete Fourier Transform, Walsh Hadamard Transform, -Number Theoretic Transform, Zeta Transform, Mobius Transform +Discrete Fourier Transform, Number Theoretic Transform, +Walsh Hadamard Transform, Zeta Transform, Mobius Transform """ from __future__ import print_function, division @@ -47,11 +47,11 @@ def _fourier_transform(seq, symbolic, inverse=False): a[i], a[j] = a[j], a[i] ang = -2*pi/n if inverse else 2*pi/n - w = None - if symbolic: - w = [cos(ang*i) + I*sin(ang*i) for i in range(n // 2)] - else: - w = [exp(I*ang*i).evalf() for i in range(n // 2)] + + if not symbolic: + ang = ang.evalf() + + w = [cos(ang*i) + I*sin(ang*i) for i in range(n // 2)] h = 2 while h <= n: @@ -120,3 +120,108 @@ def ifft(seq, symbolic=True): return _fourier_transform(seq, symbolic=symbolic, inverse=True) ifft.__doc__ = fft.__doc__ + + +#----------------------------------------------------------------------------# +# # +# Number Theoretic Transform # +# # +#----------------------------------------------------------------------------# + +def _number_theoretic_transform(seq, q, inverse=False): + """Utility function for the Number Theoretic transform (NTT)""" + + if not iterable(seq): + raise TypeError("Expected a sequence of integer coefficients " + + "for Number Theoretic Transform") + + q = as_int(q) + if isprime(q) == False: + raise ValueError("Expected prime modulo for " + + "Number Theoretic Transform") + + a = [as_int(x) for x in seq] + + n = len(a) + if n < 1: + return a + + b = n.bit_length() - 1 + if n&(n - 1): + b += 1 + n = 2**b + + if (q - 1) % n: + raise ValueError("Expected prime modulo of the form (m*2**k + 1)") + + a += [0]*(n - len(a)) + for i in range(1, n): + j = int(ibin(i, b, str=True)[::-1], 2) + if i < j: + a[i], a[j] = a[j] % q, a[i] % q + + pr = primitive_root(q) + + rt = pow(pr, (q - 1) // n, q) + if inverse: + rt = pow(rt, q - 2, q) + + w = [1]*(n // 2) + for i in range(1, n // 2): + w[i] = w[i - 1] * rt % q + + h = 2 + while h <= n: + hf, ut = h // 2, n // h + for i in range(0, n, h): + for j in range(hf): + u, v = a[i + j], a[i + j + hf]*w[ut * j] + a[i + j], a[i + j + hf] = (u + v) % q, (u - v) % q + h *= 2 + + if inverse: + rv = pow(n, q - 2, q) + for i in range(n): + a[i] = a[i]*rv % q + + return a + + +def ntt(seq, q): + r""" + Performs the Number Theoretic Transform (NTT), which specializes the + Discrete Fourier Transform (DFT) over quotient ring Z/pZ for prime p + instead of complex numbers C. + + + The sequence is automatically padded to the right with zeros, as the + radix 2 NTT requires the number of sample points to be a power of 2. + + Examples + ======== + + >>> from sympy import ntt, intt + >>> ntt([1, 2, 3, 4], 3*2**8 + 1) + [10, 643, 767, 122] + >>> intt(_, 3*2**8 + 1) + [1, 2, 3, 4] + >>> intt([1, 2, 3, 4], 3*2**8 + 1) + [387, 415, 384, 353] + >>> ntt(_, 3*2**8 + 1) + [1, 2, 3, 4] + + References + ========== + + .. [1] http://www.apfloat.org/ntt.html + .. [2] http://mathworld.wolfram.com/NumberTheoreticTransform.html + + """ + + return _number_theoretic_transform(seq, q) + + +def intt(seq, q): + return _number_theoretic_transform(seq, q, inverse=True) + +intt.__doc__ = ntt.__doc__ From b63a021c7bd2c4aa2d3a17c0941c73b6f6c13ed3 Mon Sep 17 00:00:00 2001 From: sidhantnagpal Date: Wed, 23 May 2018 07:53:44 +0530 Subject: [PATCH 06/11] update comments for tests --- sympy/discrete/tests/test_transforms.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/sympy/discrete/tests/test_transforms.py b/sympy/discrete/tests/test_transforms.py index bde733858fd9..0fc1aa497162 100644 --- a/sympy/discrete/tests/test_transforms.py +++ b/sympy/discrete/tests/test_transforms.py @@ -32,21 +32,25 @@ def test_fft_ifft(): def test_ntt_intt(): - # prime modulo of the form (m*2**k + 1), sequence length should be a divisor of 2**k + # prime modulo of the form (m*2**k + 1), sequence length should + # be a divisor of 2**k p = 7*17*2**23 + 1 - q = 2*500000003 + 1 # prime modulo not of the required form (for length > 1) + q = 2*500000003 + 1 # prime modulo only for (length = 1) r = 2*3*5*7 # composite modulo assert all(tf(ls, p) == ls for tf in (ntt, intt) \ for ls in ([], [5])) ls = list(range(6)) - nls = [15, 801133602, 738493201, 334102277, 998244350, 849020224, 259751156, 12232587] + nls = [15, 801133602, 738493201, 334102277, 998244350, 849020224, \ + 259751156, 12232587] + assert ntt(ls, p) == nls assert intt(nls, p) == ls + [0]*2 ls = [1 + 2*I, 3 + 4*I, 5 + 6*I] x = Symbol('x', integer=True) + raises(TypeError, lambda: ntt(x, p)) raises(ValueError, lambda: intt([x, 2*x, 3*x**2, 4*x**3], p)) raises(ValueError, lambda: intt(ls, p)) From 7635b868bfd04dc5ba38bd9834bb884cdef3fa86 Mon Sep 17 00:00:00 2001 From: sidhantnagpal Date: Wed, 23 May 2018 13:06:52 +0530 Subject: [PATCH 07/11] improve documentation --- sympy/discrete/transforms.py | 41 ++++++++++++++++++++++-------------- 1 file changed, 25 insertions(+), 16 deletions(-) diff --git a/sympy/discrete/transforms.py b/sympy/discrete/transforms.py index e7c0b86b1221..de628f32edfb 100644 --- a/sympy/discrete/transforms.py +++ b/sympy/discrete/transforms.py @@ -128,15 +128,15 @@ def ifft(seq, symbolic=True): # # #----------------------------------------------------------------------------# -def _number_theoretic_transform(seq, q, inverse=False): +def _number_theoretic_transform(seq, p, inverse=False): """Utility function for the Number Theoretic transform (NTT)""" if not iterable(seq): raise TypeError("Expected a sequence of integer coefficients " + "for Number Theoretic Transform") - q = as_int(q) - if isprime(q) == False: + p = as_int(p) + if isprime(p) == False: raise ValueError("Expected prime modulo for " + "Number Theoretic Transform") @@ -151,24 +151,24 @@ def _number_theoretic_transform(seq, q, inverse=False): b += 1 n = 2**b - if (q - 1) % n: + if (p - 1) % n: raise ValueError("Expected prime modulo of the form (m*2**k + 1)") a += [0]*(n - len(a)) for i in range(1, n): j = int(ibin(i, b, str=True)[::-1], 2) if i < j: - a[i], a[j] = a[j] % q, a[i] % q + a[i], a[j] = a[j] % p, a[i] % p - pr = primitive_root(q) + pr = primitive_root(p) - rt = pow(pr, (q - 1) // n, q) + rt = pow(pr, (p - 1) // n, p) if inverse: - rt = pow(rt, q - 2, q) + rt = pow(rt, p - 2, p) w = [1]*(n // 2) for i in range(1, n // 2): - w[i] = w[i - 1] * rt % q + w[i] = w[i - 1] * rt % p h = 2 while h <= n: @@ -176,18 +176,18 @@ def _number_theoretic_transform(seq, q, inverse=False): for i in range(0, n, h): for j in range(hf): u, v = a[i + j], a[i + j + hf]*w[ut * j] - a[i + j], a[i + j + hf] = (u + v) % q, (u - v) % q + a[i + j], a[i + j + hf] = (u + v) % p, (u - v) % p h *= 2 if inverse: - rv = pow(n, q - 2, q) + rv = pow(n, p - 2, p) for i in range(n): - a[i] = a[i]*rv % q + a[i] = a[i]*rv % p return a -def ntt(seq, q): +def ntt(seq, p): r""" Performs the Number Theoretic Transform (NTT), which specializes the Discrete Fourier Transform (DFT) over quotient ring Z/pZ for prime p @@ -197,6 +197,15 @@ def ntt(seq, q): The sequence is automatically padded to the right with zeros, as the radix 2 NTT requires the number of sample points to be a power of 2. + Parameters + ========== + + seq : iterable + The sequence on which DFT is to be applied. + p : Integer + Prime modulus of the form (m*2**k + 1) to be used for performing + NTT on the sequence. + Examples ======== @@ -218,10 +227,10 @@ def ntt(seq, q): """ - return _number_theoretic_transform(seq, q) + return _number_theoretic_transform(seq, p) -def intt(seq, q): - return _number_theoretic_transform(seq, q, inverse=True) +def intt(seq, p): + return _number_theoretic_transform(seq, p, inverse=True) intt.__doc__ = ntt.__doc__ From 14d03bf1399084617430619acdc7718a1cb92420 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9Csidhantnagpal=E2=80=9D?= Date: Wed, 23 May 2018 16:16:12 +0530 Subject: [PATCH 08/11] remove backward slash, fix comments --- sympy/discrete/tests/test_transforms.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/sympy/discrete/tests/test_transforms.py b/sympy/discrete/tests/test_transforms.py index 0fc1aa497162..d35fa8f5c363 100644 --- a/sympy/discrete/tests/test_transforms.py +++ b/sympy/discrete/tests/test_transforms.py @@ -8,8 +8,8 @@ def test_fft_ifft(): - assert all(tf(ls) == ls for tf in (fft, ifft) \ - for ls in ([], [S(5)/3])) + assert all(tf(ls) == ls for tf in (fft, ifft) + for ls in ([], [S(5)/3])) ls = list(range(6)) fls = [15, -7*sqrt(2)/2 - 4 - sqrt(2)*I/2 + 2*I, 2 + 3*I, @@ -32,17 +32,17 @@ def test_fft_ifft(): def test_ntt_intt(): - # prime modulo of the form (m*2**k + 1), sequence length should + # prime modulus of the form (m*2**k + 1), sequence length should # be a divisor of 2**k p = 7*17*2**23 + 1 - q = 2*500000003 + 1 # prime modulo only for (length = 1) - r = 2*3*5*7 # composite modulo + q = 2*500000003 + 1 # prime modulus (valid only for unity length) + r = 2*3*5*7 # composite modulus - assert all(tf(ls, p) == ls for tf in (ntt, intt) \ + assert all(tf(ls, p) == ls for tf in (ntt, intt) for ls in ([], [5])) ls = list(range(6)) - nls = [15, 801133602, 738493201, 334102277, 998244350, 849020224, \ + nls = [15, 801133602, 738493201, 334102277, 998244350, 849020224, 259751156, 12232587] assert ntt(ls, p) == nls From c77648fb612b1bcc89554aba3a2da54a7bd83462 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9Csidhantnagpal=E2=80=9D?= Date: Wed, 23 May 2018 20:17:56 +0530 Subject: [PATCH 09/11] Add decimal places arg for dft --- sympy/discrete/transforms.py | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/sympy/discrete/transforms.py b/sympy/discrete/transforms.py index de628f32edfb..b6d0670c1433 100644 --- a/sympy/discrete/transforms.py +++ b/sympy/discrete/transforms.py @@ -20,7 +20,7 @@ # # #----------------------------------------------------------------------------# -def _fourier_transform(seq, symbolic, inverse=False): +def _fourier_transform(seq, symbolic, dps, inverse=False): """Utility function for the Discrete Fourier Transform (DFT)""" if not iterable(seq): @@ -49,7 +49,9 @@ def _fourier_transform(seq, symbolic, inverse=False): ang = -2*pi/n if inverse else 2*pi/n if not symbolic: - ang = ang.evalf() + if dps is None: + dps = 15 + ang = ang.evalf(dps + 10) w = [cos(ang*i) + I*sin(ang*i) for i in range(n // 2)] @@ -63,13 +65,13 @@ def _fourier_transform(seq, symbolic, inverse=False): h *= 2 if inverse: - for i in range(n): - a[i] /= n + a = [x/n for x in a] if symbolic else \ + [(x/n).evalf(dps) for x in a] return a -def fft(seq, symbolic=True): +def fft(seq, symbolic=True, dps=None): r""" Performs the Discrete Fourier Transform (DFT) in the complex domain. @@ -84,6 +86,8 @@ def fft(seq, symbolic=True): symbolic : bool Determines if DFT is to be performed using symbolic values or numerical values. + dps : Integer + Specifies the number of decimal digits for precision. Examples ======== @@ -100,6 +104,11 @@ def fft(seq, symbolic=True): >>> fft(_) [1, 2, 3, 4] + >>> ifft([1, 7, 3, 4], symbolic=False) + [3.75, -0.5 - 0.75*I, -1.75, -0.5 + 0.75*I] + >>> fft(_) + [1.0, 7.0, 3.0, 4.0] + >>> fft([5]) [5] >>> ifft([7]) @@ -113,11 +122,11 @@ def fft(seq, symbolic=True): """ - return _fourier_transform(seq, symbolic=symbolic) + return _fourier_transform(seq, symbolic=symbolic, dps=dps) -def ifft(seq, symbolic=True): - return _fourier_transform(seq, symbolic=symbolic, inverse=True) +def ifft(seq, symbolic=True, dps=None): + return _fourier_transform(seq, symbolic=symbolic, dps=dps, inverse=True) ifft.__doc__ = fft.__doc__ From 568bfc33900af3197da6da1b217cf46ee667aabd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E2=80=9Csidhantnagpal=E2=80=9D?= Date: Thu, 24 May 2018 07:14:48 +0530 Subject: [PATCH 10/11] use only dps arg for dft --- sympy/discrete/transforms.py | 27 +++++++++++---------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/sympy/discrete/transforms.py b/sympy/discrete/transforms.py index b6d0670c1433..7c9baef8b10d 100644 --- a/sympy/discrete/transforms.py +++ b/sympy/discrete/transforms.py @@ -20,7 +20,7 @@ # # #----------------------------------------------------------------------------# -def _fourier_transform(seq, symbolic, dps, inverse=False): +def _fourier_transform(seq, dps, inverse=False): """Utility function for the Discrete Fourier Transform (DFT)""" if not iterable(seq): @@ -48,9 +48,7 @@ def _fourier_transform(seq, symbolic, dps, inverse=False): ang = -2*pi/n if inverse else 2*pi/n - if not symbolic: - if dps is None: - dps = 15 + if dps is not None: ang = ang.evalf(dps + 10) w = [cos(ang*i) + I*sin(ang*i) for i in range(n // 2)] @@ -65,13 +63,13 @@ def _fourier_transform(seq, symbolic, dps, inverse=False): h *= 2 if inverse: - a = [x/n for x in a] if symbolic else \ - [(x/n).evalf(dps) for x in a] + a = [(x/n).evalf(dps) for x in a] if dps is not None \ + else [x/n for x in a] return a -def fft(seq, symbolic=True, dps=None): +def fft(seq, dps=None): r""" Performs the Discrete Fourier Transform (DFT) in the complex domain. @@ -83,9 +81,6 @@ def fft(seq, symbolic=True, dps=None): seq : iterable The sequence on which DFT is to be applied. - symbolic : bool - Determines if DFT is to be performed using symbolic values or - numerical values. dps : Integer Specifies the number of decimal digits for precision. @@ -104,7 +99,7 @@ def fft(seq, symbolic=True, dps=None): >>> fft(_) [1, 2, 3, 4] - >>> ifft([1, 7, 3, 4], symbolic=False) + >>> ifft([1, 7, 3, 4], dps=15) [3.75, -0.5 - 0.75*I, -1.75, -0.5 + 0.75*I] >>> fft(_) [1.0, 7.0, 3.0, 4.0] @@ -122,11 +117,11 @@ def fft(seq, symbolic=True, dps=None): """ - return _fourier_transform(seq, symbolic=symbolic, dps=dps) + return _fourier_transform(seq, dps=dps) -def ifft(seq, symbolic=True, dps=None): - return _fourier_transform(seq, symbolic=symbolic, dps=dps, inverse=True) +def ifft(seq, dps=None): + return _fourier_transform(seq, dps=dps, inverse=True) ifft.__doc__ = fft.__doc__ @@ -146,7 +141,7 @@ def _number_theoretic_transform(seq, p, inverse=False): p = as_int(p) if isprime(p) == False: - raise ValueError("Expected prime modulo for " + + raise ValueError("Expected prime modulus for " + "Number Theoretic Transform") a = [as_int(x) for x in seq] @@ -161,7 +156,7 @@ def _number_theoretic_transform(seq, p, inverse=False): n = 2**b if (p - 1) % n: - raise ValueError("Expected prime modulo of the form (m*2**k + 1)") + raise ValueError("Expected prime modulus of the form (m*2**k + 1)") a += [0]*(n - len(a)) for i in range(1, n): From a962db89b612a4a03c59eab60a38faa54a9060d7 Mon Sep 17 00:00:00 2001 From: Sidhant Nagpal Date: Thu, 24 May 2018 11:54:52 +0530 Subject: [PATCH 11/11] improv dps, documentation --- sympy/discrete/tests/test_transforms.py | 6 +++--- sympy/discrete/transforms.py | 3 ++- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/sympy/discrete/tests/test_transforms.py b/sympy/discrete/tests/test_transforms.py index d35fa8f5c363..d053903a4909 100644 --- a/sympy/discrete/tests/test_transforms.py +++ b/sympy/discrete/tests/test_transforms.py @@ -32,10 +32,10 @@ def test_fft_ifft(): def test_ntt_intt(): - # prime modulus of the form (m*2**k + 1), sequence length should - # be a divisor of 2**k + # prime moduli of the form (m*2**k + 1), sequence length + # should be a divisor of 2**k p = 7*17*2**23 + 1 - q = 2*500000003 + 1 # prime modulus (valid only for unity length) + q = 2*500000003 + 1 # only for sequences of length 1 or 2 r = 2*3*5*7 # composite modulus assert all(tf(ls, p) == ls for tf in (ntt, intt) diff --git a/sympy/discrete/transforms.py b/sympy/discrete/transforms.py index 7c9baef8b10d..1c72cc9a9b82 100644 --- a/sympy/discrete/transforms.py +++ b/sympy/discrete/transforms.py @@ -49,7 +49,7 @@ def _fourier_transform(seq, dps, inverse=False): ang = -2*pi/n if inverse else 2*pi/n if dps is not None: - ang = ang.evalf(dps + 10) + ang = ang.evalf(dps + 2) w = [cos(ang*i) + I*sin(ang*i) for i in range(n // 2)] @@ -228,6 +228,7 @@ def ntt(seq, p): .. [1] http://www.apfloat.org/ntt.html .. [2] http://mathworld.wolfram.com/NumberTheoreticTransform.html + .. [3] https://en.wikipedia.org/wiki/Discrete_Fourier_transform_(general) """