Discrete Module

The discrete module in SymPy implements methods to compute discrete transforms and convolutions of finite sequences.

This module contains functions which operate on discrete sequences.

Transforms - fft, ifft, ntt, intt, fwht, ifwht,

mobius_transform, inverse_mobius_transform

Convolutions - convolution, convolution_fft, convolution_ntt,

convolution_fwht, convolution_subset, covering_product, intersecting_product

Since the discrete transforms can be used to reduce the computational complexity of the discrete convolutions, the convolutions module makes use of the transforms module for efficient computation (notable for long input sequences).

Transforms

This section lists the methods which implement the basic transforms for discrete sequences.

Fast Fourier Transform

sympy.discrete.transforms.fft(seq, dps=None)

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.

This method should be used with default arguments only for short sequences as the complexity of expressions increases with the size of the sequence.

Parameters

seq : iterable

The sequence on which DFT is to be applied.

dps : Integer

Specifies the number of decimal digits for precision.

References

R86

https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm

R87

http://mathworld.wolfram.com/FastFourierTransform.html

Examples

>>> from sympy import fft, ifft

>>> fft([1, 2, 3, 4])
[10, -2 - 2*I, -2, -2 + 2*I]
>>> ifft(_)
[1, 2, 3, 4]

>>> ifft([1, 2, 3, 4])
[5/2, -1/2 + I/2, -1/2, -1/2 - I/2]
>>> fft(_)
[1, 2, 3, 4]

>>> 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]

sympy.discrete.transforms.ifft(seq, dps=None)

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.

This method should be used with default arguments only for short sequences as the complexity of expressions increases with the size of the sequence.

Parameters

seq : iterable

The sequence on which DFT is to be applied.

dps : Integer

Specifies the number of decimal digits for precision.

References

R88

https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm

R89

http://mathworld.wolfram.com/FastFourierTransform.html

Examples

>>> from sympy import fft, ifft

>>> fft([1, 2, 3, 4])
[10, -2 - 2*I, -2, -2 + 2*I]
>>> ifft(_)
[1, 2, 3, 4]

>>> ifft([1, 2, 3, 4])
[5/2, -1/2 + I/2, -1/2, -1/2 - I/2]
>>> fft(_)
[1, 2, 3, 4]

>>> 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]

Number Theoretic Transform

sympy.discrete.transforms.ntt(seq, prime)

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.

Parameters

seq : iterable

The sequence on which DFT is to be applied.

prime : Integer

Prime modulus of the form \((m 2^k + 1)\) to be used for performing NTT on the sequence.

References

R90

http://www.apfloat.org/ntt.html

R91

http://mathworld.wolfram.com/NumberTheoreticTransform.html

R92

https://en.wikipedia.org/wiki/Discrete_Fourier_transform_(general%29

Examples

>>> from sympy import ntt, intt
>>> ntt([1, 2, 3, 4], prime=3*2**8 + 1)
[10, 643, 767, 122]
>>> intt(_, 3*2**8 + 1)
[1, 2, 3, 4]
>>> intt([1, 2, 3, 4], prime=3*2**8 + 1)
[387, 415, 384, 353]
>>> ntt(_, prime=3*2**8 + 1)
[1, 2, 3, 4]

sympy.discrete.transforms.intt(seq, prime)

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.

Parameters

seq : iterable

The sequence on which DFT is to be applied.

prime : Integer

Prime modulus of the form \((m 2^k + 1)\) to be used for performing NTT on the sequence.

References

R93

http://www.apfloat.org/ntt.html

R94

http://mathworld.wolfram.com/NumberTheoreticTransform.html

R95

https://en.wikipedia.org/wiki/Discrete_Fourier_transform_(general%29

Examples

>>> from sympy import ntt, intt
>>> ntt([1, 2, 3, 4], prime=3*2**8 + 1)
[10, 643, 767, 122]
>>> intt(_, 3*2**8 + 1)
[1, 2, 3, 4]
>>> intt([1, 2, 3, 4], prime=3*2**8 + 1)
[387, 415, 384, 353]
>>> ntt(_, prime=3*2**8 + 1)
[1, 2, 3, 4]

Fast Walsh Hadamard Transform

sympy.discrete.transforms.fwht(seq)

Performs the Walsh Hadamard Transform (WHT), and uses Hadamard ordering for the sequence.

The sequence is automatically padded to the right with zeros, as the radix-2 FWHT requires the number of sample points to be a power of 2.

Parameters

seq : iterable

The sequence on which WHT is to be applied.

References

R96

https://en.wikipedia.org/wiki/Hadamard_transform

R97

https://en.wikipedia.org/wiki/Fast_Walsh%E2%80%93Hadamard_transform

Examples

>>> from sympy import fwht, ifwht
>>> fwht([4, 2, 2, 0, 0, 2, -2, 0])
[8, 0, 8, 0, 8, 8, 0, 0]
>>> ifwht(_)
[4, 2, 2, 0, 0, 2, -2, 0]

>>> ifwht([19, -1, 11, -9, -7, 13, -15, 5])
[2, 0, 4, 0, 3, 10, 0, 0]
>>> fwht(_)
[19, -1, 11, -9, -7, 13, -15, 5]

sympy.discrete.transforms.ifwht(seq)

Performs the Walsh Hadamard Transform (WHT), and uses Hadamard ordering for the sequence.

The sequence is automatically padded to the right with zeros, as the radix-2 FWHT requires the number of sample points to be a power of 2.

Parameters

seq : iterable

The sequence on which WHT is to be applied.

References

R98

https://en.wikipedia.org/wiki/Hadamard_transform

R99

https://en.wikipedia.org/wiki/Fast_Walsh%E2%80%93Hadamard_transform

Examples

>>> from sympy import fwht, ifwht
>>> fwht([4, 2, 2, 0, 0, 2, -2, 0])
[8, 0, 8, 0, 8, 8, 0, 0]
>>> ifwht(_)
[4, 2, 2, 0, 0, 2, -2, 0]

>>> ifwht([19, -1, 11, -9, -7, 13, -15, 5])
[2, 0, 4, 0, 3, 10, 0, 0]
>>> fwht(_)
[19, -1, 11, -9, -7, 13, -15, 5]

Möbius Transform

sympy.discrete.transforms.mobius_transform(seq, subset=True)

Performs the Möbius Transform for subset lattice with indices of sequence as bitmasks.

The indices of each argument, considered as bit strings, correspond to subsets of a finite set.

The sequence is automatically padded to the right with zeros, as the definition of subset/superset based on bitmasks (indices) requires the size of sequence to be a power of 2.

Parameters

seq : iterable

The sequence on which Möbius Transform is to be applied.

subset : bool

Specifies if Möbius Transform is applied by enumerating subsets or supersets of the given set.

References

R100

https://en.wikipedia.org/wiki/M%C3%B6bius_inversion_formula

R101

https://people.csail.mit.edu/rrw/presentations/subset-conv.pdf

R102

https://arxiv.org/pdf/1211.0189.pdf

Examples

>>> from sympy import symbols
>>> from sympy import mobius_transform, inverse_mobius_transform
>>> x, y, z = symbols('x y z')

>>> mobius_transform([x, y, z])
[x, x + y, x + z, x + y + z]
>>> inverse_mobius_transform(_)
[x, y, z, 0]

>>> mobius_transform([x, y, z], subset=False)
[x + y + z, y, z, 0]
>>> inverse_mobius_transform(_, subset=False)
[x, y, z, 0]

>>> mobius_transform([1, 2, 3, 4])
[1, 3, 4, 10]
>>> inverse_mobius_transform(_)
[1, 2, 3, 4]
>>> mobius_transform([1, 2, 3, 4], subset=False)
[10, 6, 7, 4]
>>> inverse_mobius_transform(_, subset=False)
[1, 2, 3, 4]

sympy.discrete.transforms.inverse_mobius_transform(seq, subset=True)

Performs the Möbius Transform for subset lattice with indices of sequence as bitmasks.

The indices of each argument, considered as bit strings, correspond to subsets of a finite set.

The sequence is automatically padded to the right with zeros, as the definition of subset/superset based on bitmasks (indices) requires the size of sequence to be a power of 2.

Parameters

seq : iterable

The sequence on which Möbius Transform is to be applied.

subset : bool

Specifies if Möbius Transform is applied by enumerating subsets or supersets of the given set.

References

R103

https://en.wikipedia.org/wiki/M%C3%B6bius_inversion_formula

R104

https://people.csail.mit.edu/rrw/presentations/subset-conv.pdf

R105

https://arxiv.org/pdf/1211.0189.pdf

Examples

>>> from sympy import symbols
>>> from sympy import mobius_transform, inverse_mobius_transform
>>> x, y, z = symbols('x y z')

>>> mobius_transform([x, y, z])
[x, x + y, x + z, x + y + z]
>>> inverse_mobius_transform(_)
[x, y, z, 0]

>>> mobius_transform([x, y, z], subset=False)
[x + y + z, y, z, 0]
>>> inverse_mobius_transform(_, subset=False)
[x, y, z, 0]

>>> mobius_transform([1, 2, 3, 4])
[1, 3, 4, 10]
>>> inverse_mobius_transform(_)
[1, 2, 3, 4]
>>> mobius_transform([1, 2, 3, 4], subset=False)
[10, 6, 7, 4]
>>> inverse_mobius_transform(_, subset=False)
[1, 2, 3, 4]

Convolutions

This section lists the methods which implement the basic convolutions for discrete sequences.

Convolution

This is a general method for calculating the convolution of discrete sequences, which internally calls one of the methods convolution_fft, convolution_ntt, convolution_fwht, or convolution_subset.

sympy.discrete.convolutions.convolution(a, b, cycle=0, dps=None, prime=None, dyadic=None, subset=None)

Performs convolution by determining the type of desired convolution using hints.

Exactly one of dps, prime, dyadic, subset arguments should be specified explicitly for identifying the type of convolution, and the argument cycle can be specified optionally.

For the default arguments, linear convolution is performed using FFT.

Parameters

a, b : iterables

The sequences for which convolution is performed.

cycle : Integer

Specifies the length for doing cyclic convolution.

dps : Integer

Specifies the number of decimal digits for precision for performing FFT on the sequence.

prime : Integer

Prime modulus of the form \((m 2^k + 1)\) to be used for performing NTT on the sequence.

dyadic : bool

Identifies the convolution type as dyadic (bitwise-XOR) convolution, which is performed using FWHT.

subset : bool

Identifies the convolution type as subset convolution.

Examples

>>> from sympy import convolution, symbols, S, I
>>> u, v, w, x, y, z = symbols('u v w x y z')

>>> convolution([1 + 2*I, 4 + 3*I], [S(5)/4, 6], dps=3)
[1.25 + 2.5*I, 11.0 + 15.8*I, 24.0 + 18.0*I]
>>> convolution([1, 2, 3], [4, 5, 6], cycle=3)
[31, 31, 28]

>>> convolution([111, 777], [888, 444], prime=19*2**10 + 1)
[1283, 19351, 14219]
>>> convolution([111, 777], [888, 444], prime=19*2**10 + 1, cycle=2)
[15502, 19351]

>>> convolution([u, v], [x, y, z], dyadic=True)
[u*x + v*y, u*y + v*x, u*z, v*z]
>>> convolution([u, v], [x, y, z], dyadic=True, cycle=2)
[u*x + u*z + v*y, u*y + v*x + v*z]

>>> convolution([u, v, w], [x, y, z], subset=True)
[u*x, u*y + v*x, u*z + w*x, v*z + w*y]
>>> convolution([u, v, w], [x, y, z], subset=True, cycle=3)
[u*x + v*z + w*y, u*y + v*x, u*z + w*x]

Convolution using Fast Fourier Transform

sympy.discrete.convolutions.convolution_fft(a, b, dps=None)

Performs linear convolution using Fast Fourier Transform.

Parameters

a, b : iterables

The sequences for which convolution is performed.

dps : Integer

Specifies the number of decimal digits for precision.

References

R106

https://en.wikipedia.org/wiki/Convolution_theorem

R107

https://en.wikipedia.org/wiki/Discrete_Fourier_transform_(general%29

Examples

>>> from sympy import S, I
>>> from sympy.discrete.convolutions import convolution_fft

>>> convolution_fft([2, 3], [4, 5])
[8, 22, 15]
>>> convolution_fft([2, 5], [6, 7, 3])
[12, 44, 41, 15]
>>> convolution_fft([1 + 2*I, 4 + 3*I], [S(5)/4, 6])
[5/4 + 5*I/2, 11 + 63*I/4, 24 + 18*I]

Convolution using Number Theoretic Transform

sympy.discrete.convolutions.convolution_ntt(a, b, prime)

Performs linear convolution using Number Theoretic Transform.

Parameters

a, b : iterables

The sequences for which convolution is performed.

prime : Integer

Prime modulus of the form \((m 2^k + 1)\) to be used for performing NTT on the sequence.

References

R108

https://en.wikipedia.org/wiki/Convolution_theorem

R109

https://en.wikipedia.org/wiki/Discrete_Fourier_transform_(general%29

Examples

>>> from sympy.discrete.convolutions import convolution_ntt
>>> convolution_ntt([2, 3], [4, 5], prime=19*2**10 + 1)
[8, 22, 15]
>>> convolution_ntt([2, 5], [6, 7, 3], prime=19*2**10 + 1)
[12, 44, 41, 15]
>>> convolution_ntt([333, 555], [222, 666], prime=19*2**10 + 1)
[15555, 14219, 19404]

Convolution using Fast Walsh Hadamard Transform

sympy.discrete.convolutions.convolution_fwht(a, b)

Performs dyadic (bitwise-XOR) convolution using Fast Walsh Hadamard Transform.

The convolution is automatically padded to the right with zeros, as the radix-2 FWHT requires the number of sample points to be a power of 2.

Parameters

a, b : iterables

The sequences for which convolution is performed.

References

R110

https://www.radioeng.cz/fulltexts/2002/02_03_40_42.pdf

R111

https://en.wikipedia.org/wiki/Hadamard_transform

Examples

>>> from sympy import symbols, S, I
>>> from sympy.discrete.convolutions import convolution_fwht

>>> u, v, x, y = symbols('u v x y')
>>> convolution_fwht([u, v], [x, y])
[u*x + v*y, u*y + v*x]

>>> convolution_fwht([2, 3], [4, 5])
[23, 22]
>>> convolution_fwht([2, 5 + 4*I, 7], [6*I, 7, 3 + 4*I])
[56 + 68*I, -10 + 30*I, 6 + 50*I, 48 + 32*I]

>>> convolution_fwht([S(33)/7, S(55)/6, S(7)/4], [S(2)/3, 5])
[2057/42, 1870/63, 7/6, 35/4]

Subset Convolution

sympy.discrete.convolutions.convolution_subset(a, b)

Performs Subset Convolution of given sequences.

The indices of each argument, considered as bit strings, correspond to subsets of a finite set.

The sequence is automatically padded to the right with zeros, as the definition of subset based on bitmasks (indices) requires the size of sequence to be a power of 2.

Parameters

a, b : iterables

The sequences for which convolution is performed.

References

R112

https://people.csail.mit.edu/rrw/presentations/subset-conv.pdf

Examples

>>> from sympy import symbols, S, I
>>> from sympy.discrete.convolutions import convolution_subset
>>> u, v, x, y, z = symbols('u v x y z')

>>> convolution_subset([u, v], [x, y])
[u*x, u*y + v*x]
>>> convolution_subset([u, v, x], [y, z])
[u*y, u*z + v*y, x*y, x*z]

>>> convolution_subset([1, S(2)/3], [3, 4])
[3, 6]
>>> convolution_subset([1, 3, S(5)/7], [7])
[7, 21, 5, 0]

Covering Product

sympy.discrete.convolutions.covering_product(a, b)

Returns the covering product of given sequences.

The indices of each argument, considered as bit strings, correspond to subsets of a finite set.

The covering product of given sequences is a sequence which contains the sum of products of the elements of the given sequences grouped by the bitwise-OR of the corresponding indices.

The sequence is automatically padded to the right with zeros, as the definition of subset based on bitmasks (indices) requires the size of sequence to be a power of 2.

Parameters

a, b : iterables

The sequences for which covering product is to be obtained.

References

R113

https://people.csail.mit.edu/rrw/presentations/subset-conv.pdf

Examples

>>> from sympy import symbols, S, I, covering_product
>>> u, v, x, y, z = symbols('u v x y z')

>>> covering_product([u, v], [x, y])
[u*x, u*y + v*x + v*y]
>>> covering_product([u, v, x], [y, z])
[u*y, u*z + v*y + v*z, x*y, x*z]

>>> covering_product([1, S(2)/3], [3, 4 + 5*I])
[3, 26/3 + 25*I/3]
>>> covering_product([1, 3, S(5)/7], [7, 8])
[7, 53, 5, 40/7]

Intersecting Product

sympy.discrete.convolutions.intersecting_product(a, b)

Returns the intersecting product of given sequences.

The indices of each argument, considered as bit strings, correspond to subsets of a finite set.

The intersecting product of given sequences is the sequence which contains the sum of products of the elements of the given sequences grouped by the bitwise-AND of the corresponding indices.

The sequence is automatically padded to the right with zeros, as the definition of subset based on bitmasks (indices) requires the size of sequence to be a power of 2.

Parameters

a, b : iterables

The sequences for which intersecting product is to be obtained.

References

R114

https://people.csail.mit.edu/rrw/presentations/subset-conv.pdf

Examples

>>> from sympy import symbols, S, I, intersecting_product
>>> u, v, x, y, z = symbols('u v x y z')

>>> intersecting_product([u, v], [x, y])
[u*x + u*y + v*x, v*y]
>>> intersecting_product([u, v, x], [y, z])
[u*y + u*z + v*y + x*y + x*z, v*z, 0, 0]

>>> intersecting_product([1, S(2)/3], [3, 4 + 5*I])
[9 + 5*I, 8/3 + 10*I/3]
>>> intersecting_product([1, 3, S(5)/7], [7, 8])
[327/7, 24, 0, 0]