From a2b0940dcf00d45e882e32a611ac10dbf2ced114 Mon Sep 17 00:00:00 2001 From: KirilBangachev <51961981+KirilBangachev@users.noreply.github.com> Date: Mon, 2 Sep 2019 15:10:32 +0300 Subject: [PATCH 1/5] Add radix2 FFT Created a dynamic implementation of the radix - 2 Fast Fourier Transform for fast polynomial multiplication. Reference: https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm#The_radix-2_DIT_case --- maths/radix2_FFT.py | 202 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 202 insertions(+) create mode 100644 maths/radix2_FFT.py diff --git a/maths/radix2_FFT.py b/maths/radix2_FFT.py new file mode 100644 index 000000000000..5717026de6db --- /dev/null +++ b/maths/radix2_FFT.py @@ -0,0 +1,202 @@ +""" +Fast Polynomial Multiplication using radix-2 fast Fourier Transform. + +Reference: +https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm#The_radix-2_DIT_case + +For polynomials of degree m and n the algorithms has complexity +O(n*logn + m*logm) + +The main part of the algorithm is split in two parts: + 1) __DFT: We compute the discrete fourier transform (DFT) of A and B using a + bottom-up dynamic approach - + 2) __multiply: Once we obtain the DFT of A*B, we can similarly + invert it to obtain A*B + +The class FFT takes two polynomials A and B with complex coefficients as arguments; +The two polynomials should be represented as a sequence of coefficients starting +from the free term. Thus, for instance x + 2*x^3 could be represented as +[0,1,0,2] or (0,1,0,2). The constructor adds some zeros at the end so that the +polynomials have the same length which is a power of 2 at least the length of +their product. + +The unit tests demonstrate how the class can be used. +""" + +import mpmath # for roots of unity +import numpy as np + + +class FFT: + def __init__(self, polyA=[0], polyB=[0]): + # Input as list + self.polyA = list(polyA)[:] + self.polyB = list(polyB)[:] + + # Remove leading zero coefficients + while self.polyA[-1] == 0: + self.polyA.pop() + self.len_A = len(self.polyA) + + while self.polyB[-1] == 0: + self.polyB.pop() + self.len_B = len(self.polyB) + + # Add 0 to make lengths equal a power of 2 + self.C_max_length = int( + 2 + ** np.ceil( + np.log2( + len(self.polyA) + len(self.polyB) - 1 + ) + ) + ) + + while len(self.polyA) < self.C_max_length: + self.polyA.append(0) + while len(self.polyB) < self.C_max_length: + self.polyB.append(0) + # A complex root used for the fourier transform + self.root = complex( + mpmath.root(x=1, n=self.C_max_length, k=1) + ) + + # The product + self.product = self.__multiply() + + # Discrete fourier transform of A and B + def __DFT(self, which): + if which == "A": + dft = [[x] for x in self.polyA] + else: + dft = [[x] for x in self.polyB] + # Corner case + if len(dft) <= 1: + return dft[0] + # + next_ncol = self.C_max_length // 2 + while next_ncol > 0: + new_dft = [[] for i in range(next_ncol)] + root = self.root ** next_ncol + + # First half of next step + current_root = 1 + for j in range( + self.C_max_length // (next_ncol * 2) + ): + for i in range(next_ncol): + new_dft[i].append( + dft[i][j] + + current_root + * dft[i + next_ncol][j] + ) + current_root *= root + # Second half of next step + current_root = 1 + for j in range( + self.C_max_length // (next_ncol * 2) + ): + for i in range(next_ncol): + new_dft[i].append( + dft[i][j] + - current_root + * dft[i + next_ncol][j] + ) + current_root *= root + # Update + dft = new_dft + next_ncol = next_ncol // 2 + return dft[0] + + # multiply the DFTs of A and B and find A*B + def __multiply(self): + dftA = self.__DFT("A") + dftB = self.__DFT("B") + inverseC = [ + [ + dftA[i] * dftB[i] + for i in range(self.C_max_length) + ] + ] + del dftA + del dftB + + # Corner Case + if len(inverseC[0]) <= 1: + return inverseC[0] + # Inverse DFT + next_ncol = 2 + while next_ncol <= self.C_max_length: + new_inverseC = [[] for i in range(next_ncol)] + root = self.root ** (next_ncol // 2) + current_root = 1 + # First half of next step + for j in range(self.C_max_length // next_ncol): + for i in range(next_ncol // 2): + # Even positions + new_inverseC[i].append( + ( + inverseC[i][j] + + inverseC[i][ + j + + self.C_max_length + // next_ncol + ] + ) + / 2 + ) + # Odd positions + new_inverseC[i + next_ncol // 2].append( + ( + inverseC[i][j] + - inverseC[i][ + j + + self.C_max_length + // next_ncol + ] + ) + / (2 * current_root) + ) + current_root *= root + # Update + inverseC = new_inverseC + next_ncol *= 2 + # Unpack + inverseC = [ + round(x[0].real, 8) + round(x[0].imag, 8) * 1j + for x in inverseC + ] + + # Remove leading 0's + while inverseC[-1] == 0: + inverseC.pop() + return inverseC + + # Overwrite __str__ for print(); Shows A, B and A*B + def __str__(self): + A = "A = " + B = "B = " + C = "A*B = " + for i in range(self.len_A): + A += str(self.polyA[i]) + "*x^" + str(i) + " + " + for i in range(self.len_B): + B += str(self.polyB[i]) + "*x^" + str(i) + " + " + for i in range(self.len_B + self.len_A - 1): + C += ( + str(self.product[i]) + + "*x^" + + str(i) + + " + " + ) + return A + "\n \n" + B + "\n \n" + C + + +# Unit tests +if __name__ == "__main__": + + A = [0, 1, 0, 2] # x+2x^3 + B = (2, 3, 4, 0) # 2+3x+4x^2 + x = FFT(A, B) + print(x.product) # 2x + 3x^2 + 8x^3 + 4x^4 + 6x^5, + # as [(-0+0j), (2+0j), (3+0j), (8+0j), (6+0j), (8+0j)] + print(x) From 9737b436c6be56c55d3a10dc8f4b9bc0584e3049 Mon Sep 17 00:00:00 2001 From: KirilBangachev <51961981+KirilBangachev@users.noreply.github.com> Date: Mon, 2 Sep 2019 15:18:57 +0300 Subject: [PATCH 2/5] Rename radix2_FFT.py to radix2_fft.py --- maths/{radix2_FFT.py => radix2_fft.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename maths/{radix2_FFT.py => radix2_fft.py} (100%) diff --git a/maths/radix2_FFT.py b/maths/radix2_fft.py similarity index 100% rename from maths/radix2_FFT.py rename to maths/radix2_fft.py From 0e8774252fad58e2f02f141c2b4cc24b9bdfc23d Mon Sep 17 00:00:00 2001 From: KirilBangachev <51961981+KirilBangachev@users.noreply.github.com> Date: Tue, 3 Sep 2019 11:51:47 +0300 Subject: [PATCH 3/5] Update radix2_fft printing Improved the printing method with f.prefix and String.join() --- maths/radix2_fft.py | 39 ++++++++++++++++++++++----------------- 1 file changed, 22 insertions(+), 17 deletions(-) diff --git a/maths/radix2_fft.py b/maths/radix2_fft.py index 5717026de6db..f10e1fbf6a62 100644 --- a/maths/radix2_fft.py +++ b/maths/radix2_fft.py @@ -174,20 +174,25 @@ def __multiply(self): # Overwrite __str__ for print(); Shows A, B and A*B def __str__(self): - A = "A = " - B = "B = " - C = "A*B = " - for i in range(self.len_A): - A += str(self.polyA[i]) + "*x^" + str(i) + " + " - for i in range(self.len_B): - B += str(self.polyB[i]) + "*x^" + str(i) + " + " - for i in range(self.len_B + self.len_A - 1): - C += ( - str(self.product[i]) - + "*x^" - + str(i) - + " + " - ) + A = "A = " + " + ".join( + [ + f"{self.polyA[i]}*x^{i}" + for i in range(self.len_A) + ] + ) + B = "B = " + " + ".join( + [ + f"{self.polyB[i]}*x^{i}" + for i in range(self.len_B) + ] + ) + C = "A*B = " + " + ".join( + [ + f"{self.product[i]}*x^{i}" + for i in range(len(self.product)) + ] + ) + return A + "\n \n" + B + "\n \n" + C @@ -197,6 +202,6 @@ def __str__(self): A = [0, 1, 0, 2] # x+2x^3 B = (2, 3, 4, 0) # 2+3x+4x^2 x = FFT(A, B) - print(x.product) # 2x + 3x^2 + 8x^3 + 4x^4 + 6x^5, - # as [(-0+0j), (2+0j), (3+0j), (8+0j), (6+0j), (8+0j)] - print(x) + print(x.product) # 2x + 3x^2 + 8x^3 + 4x^4 + 6x^5, + # as [(-0+0j), (2+0j), (3+0j), (8+0j), (6+0j), (8+0j)] + print(x) From 611f6aad9e6466f9ee73d641414b56c734934d43 Mon Sep 17 00:00:00 2001 From: KirilBangachev <51961981+KirilBangachev@users.noreply.github.com> Date: Tue, 3 Sep 2019 17:23:36 +0300 Subject: [PATCH 4/5] __str__ method update --- maths/radix2_fft.py | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/maths/radix2_fft.py b/maths/radix2_fft.py index f10e1fbf6a62..40009aeec2b0 100644 --- a/maths/radix2_fft.py +++ b/maths/radix2_fft.py @@ -175,25 +175,23 @@ def __multiply(self): # Overwrite __str__ for print(); Shows A, B and A*B def __str__(self): A = "A = " + " + ".join( - [ - f"{self.polyA[i]}*x^{i}" - for i in range(self.len_A) - ] + f"{coef}*x^{i}" + for coef, i in enumerate( + self.polyA[: self.len_A] + ) ) B = "B = " + " + ".join( - [ - f"{self.polyB[i]}*x^{i}" - for i in range(self.len_B) - ] + f"{coef}*x^{i}" + for coef, i in enumerate( + self.polyB[: self.len_B] + ) ) C = "A*B = " + " + ".join( - [ - f"{self.product[i]}*x^{i}" - for i in range(len(self.product)) - ] + f"{coef}*x^{i}" + for coef, i in enumerate(self.product) ) - return A + "\n \n" + B + "\n \n" + C + return "\n\n".join((A, B, C)) # Unit tests From d932f15a001d93fb0898e8ad6ad64b653c6cd397 Mon Sep 17 00:00:00 2001 From: KirilBangachev <51961981+KirilBangachev@users.noreply.github.com> Date: Thu, 5 Sep 2019 13:26:48 +0300 Subject: [PATCH 5/5] Turned the tests into doctests --- maths/radix2_fft.py | 73 ++++++++++++++++++++++++++++----------------- 1 file changed, 45 insertions(+), 28 deletions(-) diff --git a/maths/radix2_fft.py b/maths/radix2_fft.py index 40009aeec2b0..c7ffe96528b4 100644 --- a/maths/radix2_fft.py +++ b/maths/radix2_fft.py @@ -1,26 +1,5 @@ """ Fast Polynomial Multiplication using radix-2 fast Fourier Transform. - -Reference: -https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm#The_radix-2_DIT_case - -For polynomials of degree m and n the algorithms has complexity -O(n*logn + m*logm) - -The main part of the algorithm is split in two parts: - 1) __DFT: We compute the discrete fourier transform (DFT) of A and B using a - bottom-up dynamic approach - - 2) __multiply: Once we obtain the DFT of A*B, we can similarly - invert it to obtain A*B - -The class FFT takes two polynomials A and B with complex coefficients as arguments; -The two polynomials should be represented as a sequence of coefficients starting -from the free term. Thus, for instance x + 2*x^3 could be represented as -[0,1,0,2] or (0,1,0,2). The constructor adds some zeros at the end so that the -polynomials have the same length which is a power of 2 at least the length of -their product. - -The unit tests demonstrate how the class can be used. """ import mpmath # for roots of unity @@ -28,6 +7,48 @@ class FFT: + """ + Fast Polynomial Multiplication using radix-2 fast Fourier Transform. + + Reference: + https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm#The_radix-2_DIT_case + + For polynomials of degree m and n the algorithms has complexity + O(n*logn + m*logm) + + The main part of the algorithm is split in two parts: + 1) __DFT: We compute the discrete fourier transform (DFT) of A and B using a + bottom-up dynamic approach - + 2) __multiply: Once we obtain the DFT of A*B, we can similarly + invert it to obtain A*B + + The class FFT takes two polynomials A and B with complex coefficients as arguments; + The two polynomials should be represented as a sequence of coefficients starting + from the free term. Thus, for instance x + 2*x^3 could be represented as + [0,1,0,2] or (0,1,0,2). The constructor adds some zeros at the end so that the + polynomials have the same length which is a power of 2 at least the length of + their product. + + Example: + + Create two polynomials as sequences + >>> A = [0, 1, 0, 2] # x+2x^3 + >>> B = (2, 3, 4, 0) # 2+3x+4x^2 + + Create an FFT object with them + >>> x = FFT(A, B) + + Print product + >>> print(x.product) # 2x + 3x^2 + 8x^3 + 4x^4 + 6x^5 + [(-0+0j), (2+0j), (3+0j), (8+0j), (6+0j), (8+0j)] + + __str__ test + >>> print(x) + A = 0*x^0 + 1*x^1 + 2*x^0 + 3*x^2 + B = 0*x^2 + 1*x^3 + 2*x^4 + 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) + """ + def __init__(self, polyA=[0], polyB=[0]): # Input as list self.polyA = list(polyA)[:] @@ -191,15 +212,11 @@ def __str__(self): for coef, i in enumerate(self.product) ) - return "\n\n".join((A, B, C)) + return "\n".join((A, B, C)) # Unit tests if __name__ == "__main__": + import doctest - A = [0, 1, 0, 2] # x+2x^3 - B = (2, 3, 4, 0) # 2+3x+4x^2 - x = FFT(A, B) - print(x.product) # 2x + 3x^2 + 8x^3 + 4x^4 + 6x^5, - # as [(-0+0j), (2+0j), (3+0j), (8+0j), (6+0j), (8+0j)] - print(x) + doctest.testmod()