Fourier{BigFloat}

[ioannisPApapadopoulos]

Is there a way to work with Fourier and BigFloat? It goes via FFTW which is causing some issues.

using ClassicalOrthogonalPolynomials
F = Fourier{BigFloat}()
F[:, 1:10] \ sin.(axes(F,1))

ERROR: MethodError: no method matching mul!(::Vector{BigFloat}, ::FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}, ::Vector{BigFloat}, ::Bool, ::Bool)

Closest candidates are:
  mul!(::StridedVecOrMat, ::SparseArrays.AbstractSparseMatrixCSC, ::Union{LinearAlgebra.Adjoint{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.LowerTriangular, LinearAlgebra.Transpose{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedVector, StridedMatrix, BitMatrix, BitVector}, ::Number, ::Number)
   @ SparseArrays C:\Users\john.papad\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\SparseArrays\src\linalg.jl:30
  mul!(::StridedVecOrMat, ::LinearAlgebra.Adjoint{<:Any, <:SparseArrays.AbstractSparseMatrixCSC}, ::Union{LinearAlgebra.Adjoint{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.LowerTriangular, LinearAlgebra.Transpose{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedVector, StridedMatrix, BitMatrix, BitVector}, ::Number, ::Number)
   @ SparseArrays C:\Users\john.papad\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\SparseArrays\src\linalg.jl:56
  mul!(::StridedVecOrMat, ::LinearAlgebra.Transpose{<:Any, <:SparseArrays.AbstractSparseMatrixCSC}, ::Union{LinearAlgebra.Adjoint{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.LowerTriangular, LinearAlgebra.Transpose{<:Any, <:Union{LinearAlgebra.LowerTriangular, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedMatrix, BitMatrix}}, LinearAlgebra.UnitLowerTriangular, LinearAlgebra.UnitUpperTriangular, LinearAlgebra.UpperTriangular, StridedVector, StridedMatrix, BitMatrix, BitVector}, ::Number, ::Number)
   @ SparseArrays C:\Users\john.papad\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\SparseArrays\src\linalg.jl:56
  ...

Stacktrace:
  [1] mul!
    @ C:\Users\john.papad\AppData\Local\Programs\Julia-1.9.2\share\julia\stdlib\v1.9\LinearAlgebra\src\matmul.jl:276 [inlined]
  [2] mul!(ret::Vector{BigFloat}, F::ClassicalOrthogonalPolynomials.ShuffledR2HC{BigFloat, FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}}, bin::Vector{Float64})
    @ ClassicalOrthogonalPolynomials C:\Users\john.papad\.julia\packages\ClassicalOrthogonalPolynomials\KZ9ds\src\classical\fourier.jl:128
  [3] *(F::ClassicalOrthogonalPolynomials.ShuffledR2HC{BigFloat, FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}}, b::Vector{Float64})
    @ ClassicalOrthogonalPolynomials C:\Users\john.papad\.julia\packages\ClassicalOrthogonalPolynomials\KZ9ds\src\classical\fourier.jl:133
  [4] \(a::ContinuumArrays.TransformFactorization{BigFloat, Vector{Float64}, ClassicalOrthogonalPolynomials.ShuffledR2HC{BigFloat, FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}}}, b::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
    @ ContinuumArrays C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\plans.jl:26
  [5] \(a::ContinuumArrays.ProjectionFactorization{BigFloat, ContinuumArrays.TransformFactorization{BigFloat, Vector{Float64}, ClassicalOrthogonalPolynomials.ShuffledR2HC{BigFloat, FFTW.r2rFFTWPlan{Float64, Vector{Int32}, false, 1, Int64}}}, UnitRange{Int64}}, b::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
    @ ContinuumArrays C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:213
  [6] transform_ldiv_size(#unused#::Tuple{Infinities.InfiniteCardinal{1}, Int64}, A::QuasiArrays.SubQuasiArray{BigFloat, 2, Fourier{BigFloat}, Tuple{Inclusion{Float64, DomainSets.RealNumbers}, UnitRange{Int64}}, false}, B::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
    @ ContinuumArrays C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:260
  [7] transform_ldiv(A::QuasiArrays.SubQuasiArray{BigFloat, 2, Fourier{BigFloat}, Tuple{Inclusion{Float64, DomainSets.RealNumbers}, UnitRange{Int64}}, false}, B::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
    @ ContinuumArrays C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:261
  [8] basis_ldiv_size
    @ C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:310 [inlined]
  [9] copy
    @ C:\Users\john.papad\.julia\packages\ContinuumArrays\6Nbuj\src\bases\bases.jl:302 [inlined]
 [10] #materialize#13
    @ C:\Users\john.papad\.julia\packages\ArrayLayouts\tD7jV\src\ldiv.jl:22 [inlined]
 [11] materialize
    @ C:\Users\john.papad\.julia\packages\ArrayLayouts\tD7jV\src\ldiv.jl:22 [inlined]
 [12] #ldiv#20
    @ C:\Users\john.papad\.julia\packages\ArrayLayouts\tD7jV\src\ldiv.jl:98 [inlined]
 [13] ldiv
    @ C:\Users\john.papad\.julia\packages\ArrayLayouts\tD7jV\src\ldiv.jl:98 [inlined]
 [14] \(A::QuasiArrays.SubQuasiArray{BigFloat, 2, Fourier{BigFloat}, Tuple{Inclusion{Float64, DomainSets.RealNumbers}, UnitRange{Int64}}, false}, B::QuasiArrays.BroadcastQuasiVector{Float64, typeof(f), Tuple{Inclusion{Float64, DomainSets.RealNumbers}}})
    @ QuasiArrays C:\Users\john.papad\.julia\packages\QuasiArrays\GYEdf\src\matmul.jl:34
 [15] top-level scope
    @ Untitled-1:6

[dlfivefifty]

GenericFFT.jl

try loading that it might just work

[TSGut]

I thought GenericFFT was a mirror of the Fourier functionality in FastTransforms?

[MikaelSlevinsky]

Did you try loading GenericFFT? Here's an alternative that works in ApproxFun

julia> using ApproxFun, GenericFFT

julia> f = (r,θ)-> sin(100*r*cos(θ))
#19 (generic function with 1 method)

julia> S = Chebyshev(big"0.1"..big"1.0")⊗Fourier(big"0.0"..big"2.0"*π)
Chebyshev(0.1000000000000000000000000000000000000000000000000000000000000000000000000000002 .. 1.0) ⊗ Fourier(【0.0,6.283185307179586476925286766559005768394338798750211641949889184615632812572396❫)

julia> F = ProductFun(LowRankFun(f, S; gridx = 256, gridy = 256))
ProductFun on Chebyshev(0.1000000000000000000000000000000000000000000000000000000000000000000000000000002 .. 1.0) ⊗ Fourier(【0.0,6.283185307179586476925286766559005768394338798750211641949889184615632812572396❫)

julia> coefficients(F); # get the array of coefficients

julia> f(big"0.123", big"0.456")
-0.9988661226402205869689991698845071591411740075338145723676043325127213580466635

julia> F(big"0.123", big"0.456")
-0.998866122640220586968999169884507159141174007533814572367604332512721358046603

[ioannisPApapadopoulos]

GenericFFT.jl gives the same error when using ClassicalOrthogonalPolynomials.jl but the ApproxFun code does work. Thanks @MikaelSlevinsky

[dlfivefifty]

Odd as ApproxFunFourier.jl uses the exact same transform (FFTW.plan_r2r!(x, FFTW.HC2R))

So I'm confused how the ApproxFun code is working

[dlfivefifty]

Ah I see we special cased the transform for BigFloat to call fft:

https://github.com/JuliaApproximation/ApproxFun.jl/blob/6a8fdd95e6147547d964d76900f65dcb47f803b0/src/Extras/fftGeneric.jl#L47