Coverage for local/lib/python2.7/site-packages/sage/rings/polynomial/convolution.py : 84%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
|
r""" Generic Convolution
Asymptotically fast convolution of lists over any commutative ring in which the multiply-by-two map is injective. (More precisely, if `x \in R`, and `x = 2^k*y` for some `k \geq 0`, we require that `R(x/2^k)` returns `y`.)
The main function to be exported is convolution().
EXAMPLES::
sage: convolution([1, 2, 3, 4, 5], [6, 7]) [6, 19, 32, 45, 58, 35]
The convolution function is reasonably fast, even though it is written in pure Python. For example, the following takes less than a second::
sage: v = convolution(list(range(1000)), list(range(1000)))
ALGORITHM: Converts the problem to multiplication in the ring `S[x]/(x^M - 1)`, where `S = R[y]/(y^K + 1)` (where `R` is the original base ring). Performs FFT with respect to the roots of unity `1, y, y^2, \ldots, y^{2K-1}` in `S`. The FFT/IFFT are accomplished with just additions and subtractions and rotating python lists. (I think this algorithm is essentially due to Schonhage, not completely sure.) The pointwise multiplications are handled recursively, switching to a classical algorithm at some point.
Complexity is O(n log(n) log(log(n))) additions/subtractions in R and O(n log(n)) multiplications in R.
AUTHORS:
- David Harvey (2007-07): first implementation
- William Stein: editing the docstrings for inclusion in Sage. """
################################################################################# # Copyright (C) 2007 William Stein <wstein@gmail.com> # David Harvey <dmharvey@math.harvard.edu> # # Distributed under the terms of the GNU General Public License (GPL) # # http://www.gnu.org/licenses/ #*****************************************************************************
from sage.structure.all import parent from math import log, ceil
#-------------------------------------------------------------------- # main exported routine
def convolution(L1, L2): """ Returns convolution of non-empty lists L1 and L2. L1 and L2 may have arbitrary lengths.
EXAMPLES::
sage: convolution([1, 2, 3], [4, 5, 6, 7]) [4, 13, 28, 34, 32, 21]
::
sage: R = Integers(47) sage: L1 = [R.random_element() for _ in range(1000)] sage: L2 = [R.random_element() for _ in range(3756)] sage: L3 = convolution(L1, L2) sage: L3[2000] == sum([L1[i] * L2[2000-i] for i in range(1000)]) True sage: len(L3) == 1000 + 3756 - 1 True """ raise ValueError("cannot compute convolution of empty lists")
else:
#-------------------------------------------------------------------- # naive convolution and negacyclic convolutions
def _convolution_naive(L1, L2): """ Returns convolution of non-empty lists L1 and L2, using naive algorithm. L1 and L2 may have arbitrary lengths.
EXAMPLES::
sage: from sage.rings.polynomial.convolution import _convolution_naive sage: _convolution_naive([2], [3]) [6] sage: _convolution_naive([2, 5], [3]) [6, 15] sage: _convolution_naive([2], [3, 6]) [6, 12] sage: _convolution_naive([1, 2, 3], [4, 5, 6, 7]) [4, 13, 28, 34, 32, 21] sage: _convolution_naive([4, 5, 6, 7], [1, 2, 3]) [4, 13, 28, 34, 32, 21] """
for i in range(max(0, k-m2+1), min(k+1, m1))]) \ for k in range(m1 + m2 - 1)]
def _negaconvolution_naive(L1, L2): """ Negacyclic convolution of L1 and L2, using naive algorithm. L1 and L2 must be the same length.
EXAMPLES::
sage: from sage.rings.polynomial.convolution import _negaconvolution_naive sage: from sage.rings.polynomial.convolution import _convolution_naive sage: _negaconvolution_naive([2], [3]) [6] sage: _convolution_naive([1, 2, 3], [3, 4, 5]) [3, 10, 22, 22, 15] sage: _negaconvolution_naive([1, 2, 3], [3, 4, 5]) [-19, -5, 22] """
sum([L1[i] * L2[N+j-i] for i in range(j+1, N)]) for j in range(N)]
#-------------------------------------------------------------------- # FFT/IFFT routines
def _forward_butterfly(L1, L2, r): r""" L1 and L2 are both lists of length K, and `0 \leq r \leq K`. They represent polynomials in `S = R[y]/(y^K + 1)`. This function returns `(L_1 + y^r L_2, L_1 - y^r L_2)`, as a list. """
[L1[i] + L2[i-r] for i in range(r, K)], \ [L1[i] + L2[i+K-r] for i in range(r)] + \ [L1[i] - L2[i-r] for i in range(r, K)]
def _inverse_butterfly(L1, L2, r): r""" L1 and L2 are both lists of length `K`, and `0 \leq r \leq K`. They represent polynomials in `S = R[y]/(y^K + 1)`. This function returns `(L_1 + L_2, y^{-r}*(L_1 - L_2))`, as a list. """
[L1[i] - L2[i] for i in range(r, K)] + \ [L2[i] - L1[i] for i in range(r)]
def _fft(L, K, start, depth, root): r""" L is a list of length `M = 2^m`, each entry a list of length `K = 2^k`.
This function only operates on the [start, start + D) portion of L, where `D = 2^\text{depth}`. This portion is interpreted as a polynomial in `S[x]/(x^D - y^(2*root))`, where `S = R[y]/(y^K + 1)`.
This function performs an inplace FFT, i.e. evaluates the polynomial at x = each D-th root of unity in S (namely the powers of `y^{2K/D}`), with results in bit-reversed order. """
# reduce mod (x^(D/2) - y^root) and mod (x^(D/2) + y^root) _forward_butterfly(L[start + i], L[start2 + i], root)
# recurse into each half
def _ifft(L, K, start, depth, root): r""" Inverse operation of ``_fft_trunc()`` (except that result is a factor of ``2^depth`` too big) """
# recurse into each half
# CRT together (x^(D/2) - y^root) and mod (x^(D/2) + y^root) _inverse_butterfly(L[start + i], L[start2 + i], root)
#-------------------------------------------------------------------- # splitting and recombining routines
def _split(L, m, k): """ Assumes L is a list of length `2^{m+k-1}`. Splits it into `2^m` lists of length `2^{k-1}`, returned as a list of lists. Each list is zero padded up to length `2^k`.
EXAMPLES::
sage: from sage.rings.polynomial.convolution import _split sage: _split([1, 2, 3, 4, 5, 6, 7, 8], 2, 2) [[1, 2, 0, 0], [3, 4, 0, 0], [5, 6, 0, 0], [7, 8, 0, 0]] sage: _split([1, 2, 3, 4, 5, 6, 7, 8], 1, 3) [[1, 2, 3, 4, 0, 0, 0, 0], [5, 6, 7, 8, 0, 0, 0, 0]] sage: _split([1, 2, 3, 4, 5, 6, 7, 8], 3, 1) [[1, 0], [2, 0], [3, 0], [4, 0], [5, 0], [6, 0], [7, 0], [8, 0]] """
def _combine(L, m, k): r""" Assumes L is a list of length `2^m`, each entry a list of length `2^k`. Combines together into a single list, effectively inverting ``_split()``, but overlaying coefficients, i.e. list #i gets added in starting at position `2^{k-1} i`. Note that the second half of the last list is ignored. """ [L[i+1][j] + L[i][j+half_K] \ for i in range(M-1) for j in range(half_K)]
def _nega_combine(L, m, k): r""" Same as ``_combine()``, but doesn't ignore the second half of the last list; instead it makes that piece wrap around negacyclically. """ [L[i+1][j] + L[i][j+half_K] \ for i in range(M-1) for j in range(half_K)]
#-------------------------------------------------------------------- # FFT-based convolution and negacyclic convolutions
def _negaconvolution(L1, L2, n): """ Negacyclic convolution of L1 and L2. L1 and L2 must both be length `2^n`. """ else:
def _negaconvolution_fft(L1, L2, n): """ Returns negacyclic convolution of lists L1 and L2, using FFT algorithm. L1 and L2 must both be length `2^n`, where `n \geq 3`. Assumes all entries of L1 and L2 belong to the same ring.
EXAMPLES::
sage: from sage.rings.polynomial.convolution import _negaconvolution_naive sage: from sage.rings.polynomial.convolution import _negaconvolution_fft sage: _negaconvolution_naive(list(range(8)), list(range(5, 13))) [-224, -234, -224, -192, -136, -54, 56, 196] sage: _negaconvolution_fft(list(range(8)), list(range(5, 13)), 3) [-224, -234, -224, -192, -136, -54, 56, 196]
::
sage: for n in range(3, 10): ....: L1 = [ZZ.random_element(100) for _ in range(1 << n)] ....: L2 = [ZZ.random_element(100) for _ in range(1 << n)] ....: assert _negaconvolution_naive(L1, L2) == _negaconvolution_fft(L1, L2, n) """
# split into 2^m pieces of 2^(k-1) coefficients each, with k as small # as possible, subject to m <= k (so that the ring of Fourier coefficients # has enough roots of unity)
# split inputs into polynomials
# fft each input
# pointwise multiply
# inverse fft
# combine back into a single list
# normalise
def _convolution_fft(L1, L2): """ Returns convolution of non-empty lists L1 and L2, using FFT algorithm. L1 and L2 may have arbitrary lengths `\geq 4`. Assumes all entries of L1 and L2 belong to the same ring.
EXAMPLES::
sage: from sage.rings.polynomial.convolution import _convolution_naive sage: from sage.rings.polynomial.convolution import _convolution_fft sage: _convolution_naive([1, 2, 3], [4, 5, 6]) [4, 13, 28, 27, 18] sage: _convolution_fft([1, 2, 3], [4, 5, 6]) [4, 13, 28, 27, 18]
::
sage: for len1 in range(4, 30): ....: for len2 in range(4, 30): ....: L1 = [ZZ.random_element(100) for _ in range(len1)] ....: L2 = [ZZ.random_element(100) for _ in range(len2)] ....: assert _convolution_naive(L1, L2) == _convolution_fft(L1, L2) """
# choose n so that output convolution length is 2^n
# split into 2^m pieces of 2^(k-1) coefficients each, with k as small # as possible, subject to m <= k + 1 (so that the ring of Fourier # coefficients has enough roots of unity)
# zero pad inputs up to length N
# split inputs into polynomials
# fft each input
# pointwise multiply
# inverse fft
# combine back into a single list
# normalise, and truncate to correct length
######## end of file |