Question: is multipliciation of polynomials implemented with FFT

[mathprofessor]

How is multiplication of Polynomials implemented? With the basic O(n^2) algorithm or the Fast Fourier Transform O(n log n) algorithm?

[jverzani]

Here is a nice reference https://faculty.sites.iastate.edu/jia/files/inline-files/polymultiply.pdf

[jverzani]

Sorry, I misread the question. Multiplication is currently done with the basic algorithm. It would be nice to have it done more efficiently.

[mathprofessor]

It’s currently beyond my technical skills, but it would be great if somebody could bind the C++ library number theory library NLT so Julia could use it. It does multiplication very fast. I did implement the fast Fourier transform in Julia, but my experiment said it was growing quadratically even though I know it is really n log n. Maybe the recursive calls are too much and not efficient.

On Tue, Aug 15, 2023 at 7:54 AM john verzani ***@***.***> wrote: Sorry, I misread the question. Multiplication is currently done with the basic algorithm. It would be nice to have it done more efficiently. — Reply to this email directly, view it on GitHub <#519 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ASZZERESLXBC6CGMGVDL5LTXVNPP5ANCNFSM6AAAAAA3MJO2HM> . You are receiving this because you authored the thread.Message ID: ***@***.***>

-- Joshua Friedman PhD ***@***.*** http://www.math.sunysb.edu/~joshua

[jverzani]

I had hopes this would be reasonably quick, but it doesn't seem to be. Maybe you can play around?

using Polynomials, Test

recursive_fft(as, ys=zeros(Complex{Float64}, length(as)), s=1) = recursive_fft!(ys,as,s)

function recursive_fft!(ys, as, s=1)

    n = length(as)
    if isone(n)
        ys[1] = as[1]
        return ys
    end

    o = 1
    eidx = o .+ range(0, n-2, step=2)
    oidx = o .+ range(1, n-1, step=2)

    n2 = n ÷ 2

    ye = recursive_fft!(view(ys, 1:n2), view(as, eidx), s)
    yo = recursive_fft!(view(ys, (n2+1):n), view(as, oidx), s)

    ωₙ = exp(sign(s) * im * 2pi/n)
    ω = one(ωₙ)

    @inbounds for k ∈ o .+ range(0, n2 - 1)
        yₑ, yₒ = ye[k], yo[k]
        ys[k     ] = yₑ + ω * yₒ
        ys[k + n2] = yₑ - ω * yₒ
        ω = ω * ωₙ
    end
    return ys
end
inverse_fft(as, ys=zeros(Complex{Float64}, length(as))) = recursive_fft!(ys,as,-1) / length(as)


function poly_mul(as::AbstractVector{T}, bs::AbstractVector{S}) where {T,S}
    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)
    âs = recursive_fft(as′)

    as′ .= 0
    copy!(view(as′, 1:length(bs)), bs)
    b̂s = recursive_fft(as′)

    âb̂s = âs .* b̂s

    inverse_fft(âb̂s)
end

P = Polynomial
@testset for n ∈ (4,8,16)
    as = rand(n)
    @test norm(P(poly_mul(as, as)) - P(as)*P(as)) <= 1e-10
end

[mathprofessor]

Thanks John I will take a look at your code. I also gave it a try. Below is the FFT for exact rationals. Unfortunately, as the length of the Vector grows, the performance seems to be quadratic and not O(n log n). using Polynomials using Cyclotomics function FFT(a::Vector{Cyclotomic},n::Int64,wn::Cyclotomic) if n == 1 return a end w = 1 a0 = Vector{Cyclotomic}() a1 = Vector{Cyclotomic}() for i in 1:2:n-1 push!(a0,a[i]) push!(a1,a[i+1]) end nov2 = n÷2 a0h = FFT(a0,nov2,wn*wn) a1h = FFT(a1,nov2,wn*wn) ah = Vector{Cyclotomic}(undef,n) for k in 1:nov2 ah[k] = a0h[k] + w*a1h[k] ah[k+nov2] = a0h[k] - w*a1h[k] w = w*wn end return ah end function main() N = 2^11; a = Vector{Cyclotomic}(undef,N); for i in 1:N a[i] = BigInt(i)*E(N)^0; end ah = FFT(a,N,E(N)); b = FFT(ah,N,conj(E(N))); #println(b // N) end @time main()

…

On Tue, Aug 15, 2023 at 6:40 PM john verzani ***@***.***> wrote: I had hopes this would be reasonably quick, but it doesn't seem to be. Maybe you can play around? using Polynomials, Test recursive_fft(as, ys=zeros(Complex{Float64}, length(as)), s=1) = recursive_fft!(ys,as,s) function recursive_fft!(ys, as, s=1) n = length(as) if isone(n) ys[1] = as[1] return ys end o = 1 eidx = o .+ range(0, n-2, step=2) oidx = o .+ range(1, n-1, step=2) n2 = n ÷ 2 ye = recursive_fft!(view(ys, 1:n2), view(as, eidx), s) yo = recursive_fft!(view(ys, (n2+1):n), view(as, oidx), s) ωₙ = exp(sign(s) * im * 2pi/n) ω = one(ωₙ) @inbounds for k ∈ o .+ range(0, n2 - 1) yₑ, yₒ = ye[k], yo[k] ys[k ] = yₑ + ω * yₒ ys[k + n2] = yₑ - ω * yₒ ω = ω * ωₙ end return ys end inverse_fft(as, ys=zeros(Complex{Float64}, length(as))) = recursive_fft!(ys,as,-1) / length(as) function poly_mul(as::AbstractVector{T}, bs::AbstractVector{S}) where {T,S} 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) âs = recursive_fft(as′) as′ .= 0 copy!(view(as′, 1:length(bs)), bs) b̂s = recursive_fft(as′) âb̂s = âs .* b̂s inverse_fft(âb̂s) end P = Polynomial @testset for n ∈ (4,8,16) as = rand(n) @test norm(P(poly_mul(as, as)) - P(as)*P(as)) <= 1e-10 end — Reply to this email directly, view it on GitHub <#519 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ASZZERA6I45S7ODQDKGLGYTXVP3ENANCNFSM6AAAAAA3MJO2HM> . You are receiving this because you authored the thread.Message ID: ***@***.***>

-- Joshua Friedman PhD ***@***.*** http://www.math.sunysb.edu/~joshua

[jverzani]

I like how Cyclotomics is exact, but I can't seem to get the tricks from MutableArithmetics to work to eliminate all the allocations of using big numbers. Here is basically your function, with a slightly different flavor:

function recursive_fft(as::Vector{T}, inverse=false) where {T}
    n = length(as)
    N = 2^ceil(Int, log2(n))
    @assert n == N
    R = typeof(as[1]*E(N))
    ys = Vector{R}(undef, N)
    recursive_fft!(ys, as, inverse)
    ys
end

inverse_fft(as) = recursive_fft(as, true) // length(as)

function recursive_fft!(ys, as, inverse=false)

    n = length(as)
    isone(n) && (ys[1] = as[1]; return ys)

    o = 1
    eidx = o .+ range(0, n-2, step=2)
    oidx = o .+ range(1, n-1, step=2)

    n2 = n ÷ 2

    ye = recursive_fft!(view(ys, 1:n2), view(as, eidx), inverse)
    yo = recursive_fft!(view(ys, (n2+1):n), view(as, oidx), inverse)

    ωₙ = E(n)
    inverse && (ωₙ = ωₙ^(n-1)) # inv(ωₙ) changes type)
    ω = one(ωₙ)

    @inbounds for k ∈ o .+ range(0, n2 - 1)
        yₑ, yₒ = ye[k], yo[k]
        ys[k     ] = yₑ + ω * yₒ
        ys[k + n2] = yₑ - ω * yₒ
        ω *= ωₙ
    end
    return ys
end

We can see it is working as expected and the time:

julia> as = collect(big.(1:2^11));

julia> @time b = inverse_fft(recursive_fft(a));
 29.592997 seconds (188.26 M allocations: 4.266 GiB, 45.06% gc time)

julia> maximum(abs, b-as)
0.0

I'd bet there is a good way to do this, but it seems the Cyclotomic exactness means each omega is really a sparse vector with N terms and that makes it a bit of work to get working with the BigInt tricks in MutableArithmetics.

I thought maybe trying with ArbNumerics might help, but no luck.