Coverage for /usr/lib/python3/dist-packages/sympy/discrete/convolutions.py: 12%
103 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2Convolution (using **FFT**, **NTT**, **FWHT**), Subset Convolution,
3Covering Product, Intersecting Product
4"""
6from sympy.core import S, sympify
7from sympy.core.function import expand_mul
8from sympy.discrete.transforms import (
9 fft, ifft, ntt, intt, fwht, ifwht,
10 mobius_transform, inverse_mobius_transform)
11from sympy.utilities.iterables import iterable
12from sympy.utilities.misc import as_int
15def convolution(a, b, cycle=0, dps=None, prime=None, dyadic=None, subset=None):
16 """
17 Performs convolution by determining the type of desired
18 convolution using hints.
20 Exactly one of ``dps``, ``prime``, ``dyadic``, ``subset`` arguments
21 should be specified explicitly for identifying the type of convolution,
22 and the argument ``cycle`` can be specified optionally.
24 For the default arguments, linear convolution is performed using **FFT**.
26 Parameters
27 ==========
29 a, b : iterables
30 The sequences for which convolution is performed.
31 cycle : Integer
32 Specifies the length for doing cyclic convolution.
33 dps : Integer
34 Specifies the number of decimal digits for precision for
35 performing **FFT** on the sequence.
36 prime : Integer
37 Prime modulus of the form `(m 2^k + 1)` to be used for
38 performing **NTT** on the sequence.
39 dyadic : bool
40 Identifies the convolution type as dyadic (*bitwise-XOR*)
41 convolution, which is performed using **FWHT**.
42 subset : bool
43 Identifies the convolution type as subset convolution.
45 Examples
46 ========
48 >>> from sympy import convolution, symbols, S, I
49 >>> u, v, w, x, y, z = symbols('u v w x y z')
51 >>> convolution([1 + 2*I, 4 + 3*I], [S(5)/4, 6], dps=3)
52 [1.25 + 2.5*I, 11.0 + 15.8*I, 24.0 + 18.0*I]
53 >>> convolution([1, 2, 3], [4, 5, 6], cycle=3)
54 [31, 31, 28]
56 >>> convolution([111, 777], [888, 444], prime=19*2**10 + 1)
57 [1283, 19351, 14219]
58 >>> convolution([111, 777], [888, 444], prime=19*2**10 + 1, cycle=2)
59 [15502, 19351]
61 >>> convolution([u, v], [x, y, z], dyadic=True)
62 [u*x + v*y, u*y + v*x, u*z, v*z]
63 >>> convolution([u, v], [x, y, z], dyadic=True, cycle=2)
64 [u*x + u*z + v*y, u*y + v*x + v*z]
66 >>> convolution([u, v, w], [x, y, z], subset=True)
67 [u*x, u*y + v*x, u*z + w*x, v*z + w*y]
68 >>> convolution([u, v, w], [x, y, z], subset=True, cycle=3)
69 [u*x + v*z + w*y, u*y + v*x, u*z + w*x]
71 """
73 c = as_int(cycle)
74 if c < 0:
75 raise ValueError("The length for cyclic convolution "
76 "must be non-negative")
78 dyadic = True if dyadic else None
79 subset = True if subset else None
80 if sum(x is not None for x in (prime, dps, dyadic, subset)) > 1:
81 raise TypeError("Ambiguity in determining the type of convolution")
83 if prime is not None:
84 ls = convolution_ntt(a, b, prime=prime)
85 return ls if not c else [sum(ls[i::c]) % prime for i in range(c)]
87 if dyadic:
88 ls = convolution_fwht(a, b)
89 elif subset:
90 ls = convolution_subset(a, b)
91 else:
92 ls = convolution_fft(a, b, dps=dps)
94 return ls if not c else [sum(ls[i::c]) for i in range(c)]
97#----------------------------------------------------------------------------#
98# #
99# Convolution for Complex domain #
100# #
101#----------------------------------------------------------------------------#
103def convolution_fft(a, b, dps=None):
104 """
105 Performs linear convolution using Fast Fourier Transform.
107 Parameters
108 ==========
110 a, b : iterables
111 The sequences for which convolution is performed.
112 dps : Integer
113 Specifies the number of decimal digits for precision.
115 Examples
116 ========
118 >>> from sympy import S, I
119 >>> from sympy.discrete.convolutions import convolution_fft
121 >>> convolution_fft([2, 3], [4, 5])
122 [8, 22, 15]
123 >>> convolution_fft([2, 5], [6, 7, 3])
124 [12, 44, 41, 15]
125 >>> convolution_fft([1 + 2*I, 4 + 3*I], [S(5)/4, 6])
126 [5/4 + 5*I/2, 11 + 63*I/4, 24 + 18*I]
128 References
129 ==========
131 .. [1] https://en.wikipedia.org/wiki/Convolution_theorem
132 .. [2] https://en.wikipedia.org/wiki/Discrete_Fourier_transform_(general%29
134 """
136 a, b = a[:], b[:]
137 n = m = len(a) + len(b) - 1 # convolution size
139 if n > 0 and n&(n - 1): # not a power of 2
140 n = 2**n.bit_length()
142 # padding with zeros
143 a += [S.Zero]*(n - len(a))
144 b += [S.Zero]*(n - len(b))
146 a, b = fft(a, dps), fft(b, dps)
147 a = [expand_mul(x*y) for x, y in zip(a, b)]
148 a = ifft(a, dps)[:m]
150 return a
153#----------------------------------------------------------------------------#
154# #
155# Convolution for GF(p) #
156# #
157#----------------------------------------------------------------------------#
159def convolution_ntt(a, b, prime):
160 """
161 Performs linear convolution using Number Theoretic Transform.
163 Parameters
164 ==========
166 a, b : iterables
167 The sequences for which convolution is performed.
168 prime : Integer
169 Prime modulus of the form `(m 2^k + 1)` to be used for performing
170 **NTT** on the sequence.
172 Examples
173 ========
175 >>> from sympy.discrete.convolutions import convolution_ntt
176 >>> convolution_ntt([2, 3], [4, 5], prime=19*2**10 + 1)
177 [8, 22, 15]
178 >>> convolution_ntt([2, 5], [6, 7, 3], prime=19*2**10 + 1)
179 [12, 44, 41, 15]
180 >>> convolution_ntt([333, 555], [222, 666], prime=19*2**10 + 1)
181 [15555, 14219, 19404]
183 References
184 ==========
186 .. [1] https://en.wikipedia.org/wiki/Convolution_theorem
187 .. [2] https://en.wikipedia.org/wiki/Discrete_Fourier_transform_(general%29
189 """
191 a, b, p = a[:], b[:], as_int(prime)
192 n = m = len(a) + len(b) - 1 # convolution size
194 if n > 0 and n&(n - 1): # not a power of 2
195 n = 2**n.bit_length()
197 # padding with zeros
198 a += [0]*(n - len(a))
199 b += [0]*(n - len(b))
201 a, b = ntt(a, p), ntt(b, p)
202 a = [x*y % p for x, y in zip(a, b)]
203 a = intt(a, p)[:m]
205 return a
208#----------------------------------------------------------------------------#
209# #
210# Convolution for 2**n-group #
211# #
212#----------------------------------------------------------------------------#
214def convolution_fwht(a, b):
215 """
216 Performs dyadic (*bitwise-XOR*) convolution using Fast Walsh Hadamard
217 Transform.
219 The convolution is automatically padded to the right with zeros, as the
220 *radix-2 FWHT* requires the number of sample points to be a power of 2.
222 Parameters
223 ==========
225 a, b : iterables
226 The sequences for which convolution is performed.
228 Examples
229 ========
231 >>> from sympy import symbols, S, I
232 >>> from sympy.discrete.convolutions import convolution_fwht
234 >>> u, v, x, y = symbols('u v x y')
235 >>> convolution_fwht([u, v], [x, y])
236 [u*x + v*y, u*y + v*x]
238 >>> convolution_fwht([2, 3], [4, 5])
239 [23, 22]
240 >>> convolution_fwht([2, 5 + 4*I, 7], [6*I, 7, 3 + 4*I])
241 [56 + 68*I, -10 + 30*I, 6 + 50*I, 48 + 32*I]
243 >>> convolution_fwht([S(33)/7, S(55)/6, S(7)/4], [S(2)/3, 5])
244 [2057/42, 1870/63, 7/6, 35/4]
246 References
247 ==========
249 .. [1] https://www.radioeng.cz/fulltexts/2002/02_03_40_42.pdf
250 .. [2] https://en.wikipedia.org/wiki/Hadamard_transform
252 """
254 if not a or not b:
255 return []
257 a, b = a[:], b[:]
258 n = max(len(a), len(b))
260 if n&(n - 1): # not a power of 2
261 n = 2**n.bit_length()
263 # padding with zeros
264 a += [S.Zero]*(n - len(a))
265 b += [S.Zero]*(n - len(b))
267 a, b = fwht(a), fwht(b)
268 a = [expand_mul(x*y) for x, y in zip(a, b)]
269 a = ifwht(a)
271 return a
274#----------------------------------------------------------------------------#
275# #
276# Subset Convolution #
277# #
278#----------------------------------------------------------------------------#
280def convolution_subset(a, b):
281 """
282 Performs Subset Convolution of given sequences.
284 The indices of each argument, considered as bit strings, correspond to
285 subsets of a finite set.
287 The sequence is automatically padded to the right with zeros, as the
288 definition of subset based on bitmasks (indices) requires the size of
289 sequence to be a power of 2.
291 Parameters
292 ==========
294 a, b : iterables
295 The sequences for which convolution is performed.
297 Examples
298 ========
300 >>> from sympy import symbols, S
301 >>> from sympy.discrete.convolutions import convolution_subset
302 >>> u, v, x, y, z = symbols('u v x y z')
304 >>> convolution_subset([u, v], [x, y])
305 [u*x, u*y + v*x]
306 >>> convolution_subset([u, v, x], [y, z])
307 [u*y, u*z + v*y, x*y, x*z]
309 >>> convolution_subset([1, S(2)/3], [3, 4])
310 [3, 6]
311 >>> convolution_subset([1, 3, S(5)/7], [7])
312 [7, 21, 5, 0]
314 References
315 ==========
317 .. [1] https://people.csail.mit.edu/rrw/presentations/subset-conv.pdf
319 """
321 if not a or not b:
322 return []
324 if not iterable(a) or not iterable(b):
325 raise TypeError("Expected a sequence of coefficients for convolution")
327 a = [sympify(arg) for arg in a]
328 b = [sympify(arg) for arg in b]
329 n = max(len(a), len(b))
331 if n&(n - 1): # not a power of 2
332 n = 2**n.bit_length()
334 # padding with zeros
335 a += [S.Zero]*(n - len(a))
336 b += [S.Zero]*(n - len(b))
338 c = [S.Zero]*n
340 for mask in range(n):
341 smask = mask
342 while smask > 0:
343 c[mask] += expand_mul(a[smask] * b[mask^smask])
344 smask = (smask - 1)&mask
346 c[mask] += expand_mul(a[smask] * b[mask^smask])
348 return c
351#----------------------------------------------------------------------------#
352# #
353# Covering Product #
354# #
355#----------------------------------------------------------------------------#
357def covering_product(a, b):
358 """
359 Returns the covering product of given sequences.
361 The indices of each argument, considered as bit strings, correspond to
362 subsets of a finite set.
364 The covering product of given sequences is a sequence which contains
365 the sum of products of the elements of the given sequences grouped by
366 the *bitwise-OR* of the corresponding indices.
368 The sequence is automatically padded to the right with zeros, as the
369 definition of subset based on bitmasks (indices) requires the size of
370 sequence to be a power of 2.
372 Parameters
373 ==========
375 a, b : iterables
376 The sequences for which covering product is to be obtained.
378 Examples
379 ========
381 >>> from sympy import symbols, S, I, covering_product
382 >>> u, v, x, y, z = symbols('u v x y z')
384 >>> covering_product([u, v], [x, y])
385 [u*x, u*y + v*x + v*y]
386 >>> covering_product([u, v, x], [y, z])
387 [u*y, u*z + v*y + v*z, x*y, x*z]
389 >>> covering_product([1, S(2)/3], [3, 4 + 5*I])
390 [3, 26/3 + 25*I/3]
391 >>> covering_product([1, 3, S(5)/7], [7, 8])
392 [7, 53, 5, 40/7]
394 References
395 ==========
397 .. [1] https://people.csail.mit.edu/rrw/presentations/subset-conv.pdf
399 """
401 if not a or not b:
402 return []
404 a, b = a[:], b[:]
405 n = max(len(a), len(b))
407 if n&(n - 1): # not a power of 2
408 n = 2**n.bit_length()
410 # padding with zeros
411 a += [S.Zero]*(n - len(a))
412 b += [S.Zero]*(n - len(b))
414 a, b = mobius_transform(a), mobius_transform(b)
415 a = [expand_mul(x*y) for x, y in zip(a, b)]
416 a = inverse_mobius_transform(a)
418 return a
421#----------------------------------------------------------------------------#
422# #
423# Intersecting Product #
424# #
425#----------------------------------------------------------------------------#
427def intersecting_product(a, b):
428 """
429 Returns the intersecting product of given sequences.
431 The indices of each argument, considered as bit strings, correspond to
432 subsets of a finite set.
434 The intersecting product of given sequences is the sequence which
435 contains the sum of products of the elements of the given sequences
436 grouped by the *bitwise-AND* of the corresponding indices.
438 The sequence is automatically padded to the right with zeros, as the
439 definition of subset based on bitmasks (indices) requires the size of
440 sequence to be a power of 2.
442 Parameters
443 ==========
445 a, b : iterables
446 The sequences for which intersecting product is to be obtained.
448 Examples
449 ========
451 >>> from sympy import symbols, S, I, intersecting_product
452 >>> u, v, x, y, z = symbols('u v x y z')
454 >>> intersecting_product([u, v], [x, y])
455 [u*x + u*y + v*x, v*y]
456 >>> intersecting_product([u, v, x], [y, z])
457 [u*y + u*z + v*y + x*y + x*z, v*z, 0, 0]
459 >>> intersecting_product([1, S(2)/3], [3, 4 + 5*I])
460 [9 + 5*I, 8/3 + 10*I/3]
461 >>> intersecting_product([1, 3, S(5)/7], [7, 8])
462 [327/7, 24, 0, 0]
464 References
465 ==========
467 .. [1] https://people.csail.mit.edu/rrw/presentations/subset-conv.pdf
469 """
471 if not a or not b:
472 return []
474 a, b = a[:], b[:]
475 n = max(len(a), len(b))
477 if n&(n - 1): # not a power of 2
478 n = 2**n.bit_length()
480 # padding with zeros
481 a += [S.Zero]*(n - len(a))
482 b += [S.Zero]*(n - len(b))
484 a, b = mobius_transform(a, subset=False), mobius_transform(b, subset=False)
485 a = [expand_mul(x*y) for x, y in zip(a, b)]
486 a = inverse_mobius_transform(a, subset=False)
488 return a