|
1 | 1 | """
|
2 | 2 | Fast Polynomial Multiplication using radix-2 fast Fourier Transform.
|
3 |
| -
|
4 |
| -Reference: |
5 |
| -https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm#The_radix-2_DIT_case |
6 |
| -
|
7 |
| -For polynomials of degree m and n the algorithms has complexity |
8 |
| -O(n*logn + m*logm) |
9 |
| -
|
10 |
| -The main part of the algorithm is split in two parts: |
11 |
| - 1) __DFT: We compute the discrete fourier transform (DFT) of A and B using a |
12 |
| - bottom-up dynamic approach - |
13 |
| - 2) __multiply: Once we obtain the DFT of A*B, we can similarly |
14 |
| - invert it to obtain A*B |
15 |
| -
|
16 |
| -The class FFT takes two polynomials A and B with complex coefficients as arguments; |
17 |
| -The two polynomials should be represented as a sequence of coefficients starting |
18 |
| -from the free term. Thus, for instance x + 2*x^3 could be represented as |
19 |
| -[0,1,0,2] or (0,1,0,2). The constructor adds some zeros at the end so that the |
20 |
| -polynomials have the same length which is a power of 2 at least the length of |
21 |
| -their product. |
22 |
| -
|
23 |
| -The unit tests demonstrate how the class can be used. |
24 | 3 | """
|
25 | 4 |
|
26 | 5 | import mpmath # for roots of unity
|
27 | 6 | import numpy as np
|
28 | 7 |
|
29 | 8 |
|
30 | 9 | class FFT:
|
| 10 | + """ |
| 11 | + Fast Polynomial Multiplication using radix-2 fast Fourier Transform. |
| 12 | +
|
| 13 | + Reference: |
| 14 | + https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm#The_radix-2_DIT_case |
| 15 | + |
| 16 | + For polynomials of degree m and n the algorithms has complexity |
| 17 | + O(n*logn + m*logm) |
| 18 | + |
| 19 | + The main part of the algorithm is split in two parts: |
| 20 | + 1) __DFT: We compute the discrete fourier transform (DFT) of A and B using a |
| 21 | + bottom-up dynamic approach - |
| 22 | + 2) __multiply: Once we obtain the DFT of A*B, we can similarly |
| 23 | + invert it to obtain A*B |
| 24 | +
|
| 25 | + The class FFT takes two polynomials A and B with complex coefficients as arguments; |
| 26 | + The two polynomials should be represented as a sequence of coefficients starting |
| 27 | + from the free term. Thus, for instance x + 2*x^3 could be represented as |
| 28 | + [0,1,0,2] or (0,1,0,2). The constructor adds some zeros at the end so that the |
| 29 | + polynomials have the same length which is a power of 2 at least the length of |
| 30 | + their product. |
| 31 | + |
| 32 | + Example: |
| 33 | + |
| 34 | + Create two polynomials as sequences |
| 35 | + >>> A = [0, 1, 0, 2] # x+2x^3 |
| 36 | + >>> B = (2, 3, 4, 0) # 2+3x+4x^2 |
| 37 | + |
| 38 | + Create an FFT object with them |
| 39 | + >>> x = FFT(A, B) |
| 40 | + |
| 41 | + Print product |
| 42 | + >>> print(x.product) # 2x + 3x^2 + 8x^3 + 4x^4 + 6x^5 |
| 43 | + [(-0+0j), (2+0j), (3+0j), (8+0j), (6+0j), (8+0j)] |
| 44 | + |
| 45 | + __str__ test |
| 46 | + >>> print(x) |
| 47 | + A = 0*x^0 + 1*x^1 + 2*x^0 + 3*x^2 |
| 48 | + B = 0*x^2 + 1*x^3 + 2*x^4 |
| 49 | + A*B = 0*x^(-0+0j) + 1*x^(2+0j) + 2*x^(3+0j) + 3*x^(8+0j) + 4*x^(6+0j) + 5*x^(8+0j) |
| 50 | + """ |
| 51 | + |
31 | 52 | def __init__(self, polyA=[0], polyB=[0]):
|
32 | 53 | # Input as list
|
33 | 54 | self.polyA = list(polyA)[:]
|
@@ -191,15 +212,11 @@ def __str__(self):
|
191 | 212 | for coef, i in enumerate(self.product)
|
192 | 213 | )
|
193 | 214 |
|
194 |
| - return "\n\n".join((A, B, C)) |
| 215 | + return "\n".join((A, B, C)) |
195 | 216 |
|
196 | 217 |
|
197 | 218 | # Unit tests
|
198 | 219 | if __name__ == "__main__":
|
| 220 | + import doctest |
199 | 221 |
|
200 |
| - A = [0, 1, 0, 2] # x+2x^3 |
201 |
| - B = (2, 3, 4, 0) # 2+3x+4x^2 |
202 |
| - x = FFT(A, B) |
203 |
| - print(x.product) # 2x + 3x^2 + 8x^3 + 4x^4 + 6x^5, |
204 |
| - # as [(-0+0j), (2+0j), (3+0j), (8+0j), (6+0j), (8+0j)] |
205 |
| - print(x) |
| 222 | + doctest.testmod() |
0 commit comments