From be37955b5d57476a17f90d9bb5c38a8d3be1fc02 Mon Sep 17 00:00:00 2001 From: john verzani Date: Fri, 25 Aug 2023 17:18:05 -0400 Subject: [PATCH] Issue 519 (#532) * better job with 519, but still not what is wanted... * more comments on #519 * did extensions get copied incorrectly? * version bump * extension testing * add extensions to test/Project.toml --- Project.toml | 10 +++-- ext/PolynomialsFFTWExt.jl | 24 +++++++++++ src/Polynomials.jl | 2 - .../standard-basis/standard-basis.jl | 42 +++++++++++++------ test/Project.toml | 3 ++ test/runtests.jl | 3 +- test/test-extensions.jl | 3 ++ 7 files changed, 68 insertions(+), 19 deletions(-) create mode 100644 ext/PolynomialsFFTWExt.jl create mode 100644 test/test-extensions.jl diff --git a/Project.toml b/Project.toml index 5f42ab01..7e92cdef 100644 --- a/Project.toml +++ b/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"] diff --git a/ext/PolynomialsFFTWExt.jl b/ext/PolynomialsFFTWExt.jl new file mode 100644 index 00000000..0797f86a --- /dev/null +++ b/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 diff --git a/src/Polynomials.jl b/src/Polynomials.jl index d7008de9..3fbf66ec 100644 --- a/src/Polynomials.jl +++ b/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") diff --git a/src/polynomials/standard-basis/standard-basis.jl b/src/polynomials/standard-basis/standard-basis.jl index 6b5a5dfc..7d55c242 100644 --- a/src/polynomials/standard-basis/standard-basis.jl +++ b/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) diff --git a/test/Project.toml b/test/Project.toml index 01312580..5fd62c26 100644 --- a/test/Project.toml +++ b/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" diff --git a/test/runtests.jl b/test/runtests.jl index 0502bc75..6fdde551 100644 --- a/test/runtests.jl +++ b/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 diff --git a/test/test-extensions.jl b/test/test-extensions.jl new file mode 100644 index 00000000..be23523f --- /dev/null +++ b/test/test-extensions.jl @@ -0,0 +1,3 @@ +using FFTW +using MakieCore +using ChainRulesCore