Issue 519

Project.toml

@@ -2,28 +2,29 @@ name = "Polynomials"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 license = "MIT"
 author = "JuliaMath"
-version = "4.0.1"
+version = "4.0.2"
 
 [deps]
-ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
-MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 
 [weakdeps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 
 [extensions]
 PolynomialsChainRulesCoreExt = "ChainRulesCore"
+PolynomialsFFTWExt = "FFTW"
 PolynomialsMakieCoreExt = "MakieCore"
 PolynomialsMutableArithmeticsExt = "MutableArithmetics"
 
 [compat]
 ChainRulesCore = "1"
+FFTW = "1"
 MakieCore = "0.6"
 MutableArithmetics = "1"
 RecipesBase = "0.7, 0.8, 1"
@@ -35,6 +36,7 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
 DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
@@ -44,4 +46,4 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "ChainRulesCore", "DualNumbers", "LinearAlgebra", "SparseArrays", "OffsetArrays", "SpecialFunctions", "Test"]
+test = ["Aqua", "ChainRulesCore", "DualNumbers", "FFTW", "LinearAlgebra", "MakieCore", "MutableArithmetics", "SparseArrays", "OffsetArrays", "SpecialFunctions", "Test"]

ext/PolynomialsFFTWExt.jl

@@ -0,0 +1,24 @@
+module PolynomialsFFTWExt
+
+using Polynomials
+import Polynomials: MutableDensePolynomial, StandardBasis, Pad
+import FFTW
+import FFTW: fft, ifft
+function Polynomials.poly_multiplication_fft(p::P, q::Q) where {B <: StandardBasis,X,
+                                                                T <: AbstractFloat, P<:MutableDensePolynomial{B,T,X},
+                                                                S <: AbstractFloat, Q<:MutableDensePolynomial{B,S,X}}
+
+    N = 1 + degree(p) + degree(q)
+    as = Pad(p.coeffs, N)
+    bs = Pad(q.coeffs, N)
+    us = fft(as)
+    vs = fft(bs)
+    cs = ifft(us .* vs)
+    map!(real, cs, cs)
+    MutableDensePolynomial{B, eltype(cs), X}(cs)
+
+end
+
+
+
+end

src/Polynomials.jl

@@ -53,9 +53,7 @@ include("legacy/misc.jl")
 include("legacy/Poly.jl")
 
 if !isdefined(Base, :get_extension)
-    include("../ext/PolynomialsChainRulesCoreExt.jl")
     include("../ext/PolynomialsMakieCoreExt.jl")
-    include("../ext/PolynomialsMutableArithmeticsExt.jl")
 end
 
 include("precompiles.jl")

src/polynomials/standard-basis/standard-basis.jl

@@ -865,10 +865,31 @@ function practical_polynomial_composition(f::StandardBasisPolynomial, g::Standar
 end
 
 ## issue #519 polynomial multiplication via FFT
+##
+## This implements [Cooley-Tukey](https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm) radix-2
 ## cf. http://www.cs.toronto.edu/~denisp/csc373/docs/tutorial3-adv-writeup.pdf
-## Compute recursive_fft
-## assumes length(as) = 2^k for some k
-## ωₙ is exp(-2pi*im/n) or Cyclotomics.E(n), the latter slower but non-lossy
+##
+## This is **much slower** than that of FFTW.jl (and more restrictive). However it does allow for exact computation
+## using `Cyclotomics.jl`, or with `Mods.jl`.
+## This assumes length(as) = 2^k for some k
+## ωₙ is an nth root of unity, for example `exp(-2pi*im/n)` (also available with `sincos(2pi/n)`) for floating point
+## or Cyclotomics.E(n), the latter much slower but non-lossy.
+##
+## Should implement NTT https://www.nayuki.io/page/number-theoretic-transform-integer-dft to close #519
+
+struct Pad{T} <: AbstractVector{T}
+    a::Vector{T}
+    n::Int
+end
+Base.length(a::Pad) = a.n
+Base.size(a::Pad) = (a.n,)
+function Base.getindex(a::Pad, i)
+    u = length(a.a)
+    i ≤ u && return a.a[i]
+    return zero(first(a.a))
+end
+
+
 function recursive_fft(as, ωₙ = nothing)
     n = length(as)
     N = 2^ceil(Int, log2(n))
@@ -886,7 +907,7 @@ function inverse_fft(as, ωₙ=nothing)
     recursive_fft(as, conj(ω)) / n
 end
 
-# note: can write version for big coefficients, but still allocates a bit
+# note: could write version for big coefficients, but still allocates a bit
 function recursive_fft!(ys, as, ωₙ)
 
     n = length(as)
@@ -915,22 +936,19 @@ end
 
 # This *should* be faster -- (O(nlog(n)), but this version is definitely not so.
 # when `ωₙ = Cyclotomics.E` and T,S are integer, this can be exact
+# using `FFTW.jl` over `Float64` types is much better and is
+# implemented in an extension
 function poly_multiplication_fft(p::P, q::Q, ωₙ=nothing) where {T,P<:StandardBasisPolynomial{T},
                                                                 S,Q<:StandardBasisPolynomial{S}}
     as, bs = coeffs0(p), coeffs0(q)
     n = maximum(length, (as, bs))
     N = 2^ceil(Int, log2(n))
-
-    as′ = zeros(promote_type(T,S), 2N)
-    copy!(view(as′, 1:length(as)), as)
+    as′ = Pad(as, 2N)
+    bs′ = Pad(bs, 2N)
 
     ω = something(ωₙ, n -> exp(-2im*pi/n))(2N)
     âs = recursive_fft(as′, ω)
-
-    as′ .= 0
-    copy!(view(as′, 1:length(bs)), bs)
-    b̂s = recursive_fft(as′, ω)
-
+    b̂s = recursive_fft(bs′, ω)
     âb̂s = âs .* b̂s
 
     PP = promote_type(P,Q)

test/Project.toml

@@ -1,7 +1,10 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"

test/runtests.jl

@@ -13,6 +13,7 @@ using OffsetArrays
 @testset "Rational functions" begin include("rational-functions.jl") end
 @testset "Poly, Pade (compatibility)" begin include("Poly.jl") end
 if VERSION >= v"1.9.0-"
-    @testset "MutableArithmetics" begin include("mutable-arithmetics.jl") end
     @testset "Aqua" begin include("aqua.jl") end
+    @testset "MutableArithmetics" begin include("mutable-arithmetics.jl") end
+    @testset "Extensions" begin include("test-extensions.jl") end
 end

test/test-extensions.jl

@@ -0,0 +1,3 @@
+using FFTW
+using MakieCore
+using ChainRulesCore