Change clenshaw!/horner! to mutate first argument (#255)

* Change clenshaw!/horner! to mutate first argument

* Update libfasttransformstests.jl

* Update Project.toml

* Update make.jl

* Update sphere.jl

* Update make.jl

* Move semiclassical.jl to ApproxFunExamples

Author: dlfivefifty

Project.toml

@@ -1,6 +1,7 @@
 name = "FastTransforms"
 uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
-version = "0.16.8"
+version = "0.17"
+
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
@@ -15,7 +16,6 @@ LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 RecurrenceRelationships = "807425ed-42ea-44d6-a357-6771516d7b2c"
-Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
@@ -29,8 +29,7 @@ FastTransforms_jll = "0.6.2"
 FillArrays = "0.9, 0.10, 0.11, 0.12, 0.13, 1"
 GenericFFT = "0.1"
 LazyArrays = "2.2"
-RecurrenceRelationships = "0.1"
-Reexport = "0.2, 1.0"
+RecurrenceRelationships = "0.2"
 SpecialFunctions = "0.10, 1, 2"
 ToeplitzMatrices = "0.7.1, 0.8"
 julia = "1.7"

docs/Project.toml

@@ -1,5 +1,4 @@
 [deps]
-ApproxFun = "28f2ccd6-bb30-5033-b560-165f7b14dc2f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"

docs/make.jl

@@ -13,7 +13,6 @@ examples = [
     "halfrange.jl",
     "nonlocaldiffusion.jl",
     "padua.jl",
-    "semiclassical.jl",
     "sphere.jl",
     "spinweighted.jl",
     "subspaceangles.jl",
@@ -48,7 +47,6 @@ makedocs(
                         "generated/halfrange.md",
                         "generated/nonlocaldiffusion.md",
                         "generated/padua.md",
-                        "generated/semiclassical.md",
                         "generated/sphere.md",
                         "generated/spinweighted.md",
                         "generated/subspaceangles.md",

examples/semiclassical.jl

@@ -61,7 +61,7 @@ C = [k/(k+1) for k in 0:N]
 c = zeros(N); c[N] = 1
 pts = vec([z(θ, φ)⋅y for θ in θ, φ in φ])
 phi0 = ones(N*M)
-F = reshape(FastTransforms.clenshaw!(c, A, B, C, pts, phi0, zeros(N*M)), N, M)
+F = reshape(FastTransforms.clenshaw!(zeros(N*M), c, A, B, C, pts, phi0), N, M)
 
 # We superpose a surface plot of $f$ on top of the grid:
 X = [sinpi(θ)*cospi(φ) for θ in θ, φ in φ]
@@ -91,7 +91,7 @@ U = threshold!(P\V, 400*eps())
 nrm1 = norm(U)
 
 # Similarly, on the tensor product grid, our function samples are:
-Pnxy = FastTransforms.clenshaw!(c, A, B, C, [x⋅y], [1.0], [0.0])[1]
+Pnxy = FastTransforms.clenshaw!([0.0], c, A, B, C, [x⋅y], [1.0])[1]
 F = [(F[n, m] - Pnxy)/(z(θ[n], φ[m])⋅y - x⋅y) for n in 1:N, m in 1:M]
 
 # We superpose a surface plot of $f$ on top of the grid:
@@ -108,7 +108,7 @@ U = threshold!(P\V, 400*eps())
 
 # Finally, the Legendre polynomial $P_n(z\cdot x)$ is aligned with the grid:
 pts = vec([z(θ, φ)⋅x for θ in θ, φ in φ])
-F = reshape(FastTransforms.clenshaw!(c, A, B, C, pts, phi0, zeros(N*M)), N, M)
+F = reshape(FastTransforms.clenshaw!(zeros(N*M), c, A, B, C, pts, phi0), N, M)
 
 # We superpose a surface plot of $f$ on top of the grid:
 scatter3d(vec(X), vec(Y), vec(Z); markersize=1.25, markercolor=:violetred)

examples/sphere.jl

@@ -1,11 +1,11 @@
 module FastTransforms
 
 using ArrayLayouts, BandedMatrices, FastGaussQuadrature, FillArrays, LazyArrays, LinearAlgebra,
-      Reexport, SpecialFunctions, ToeplitzMatrices, RecurrenceRelationships
+      SpecialFunctions, ToeplitzMatrices, RecurrenceRelationships
 
-@reexport using AbstractFFTs
-@reexport using FFTW
-@reexport using GenericFFT
+using AbstractFFTs
+using FFTW
+using GenericFFT
 
 import Base: convert, unsafe_convert, eltype, ndims, adjoint, transpose, show,
              *, \, inv, length, size, view, getindex, tail, OneTo
@@ -34,11 +34,8 @@ import LinearAlgebra: cholesky, issymmetric, isposdef, mul!, lmul!, ldiv!
 
 import GenericFFT: interlace # imported in downstream packages
 
-import RecurrenceRelationships: clenshaw!, check_clenshaw_recurrences
+import RecurrenceRelationships: check_clenshaw_recurrences
 
-const _forwardrecurrence! = RecurrenceRelationships.forwardrecurrence!
-const _clenshaw_next = RecurrenceRelationships.clenshaw_next
-const _forwardrecurrence_next = RecurrenceRelationships.forwardrecurrence_next
 
 export leg2cheb, cheb2leg, ultra2ultra, jac2jac,
        lag2lag, jac2ultra, ultra2jac, jac2cheb,

src/FastTransforms.jl

@@ -49,13 +49,13 @@ function renew!(x::AbstractArray{BigFloat})
     return x
 end
 
-function horner!(c::StridedVector{Float64}, x::Vector{Float64}, f::Vector{Float64})
+function horner!(f::Vector{Float64}, c::StridedVector{Float64}, x::Vector{Float64})
     @assert length(x) == length(f)
     ccall((:ft_horner, libfasttransforms), Cvoid, (Cint, Ptr{Float64}, Cint, Cint, Ptr{Float64}, Ptr{Float64}), length(c), c, stride(c, 1), length(x), x, f)
     f
 end
 
-function horner!(c::StridedVector{Float32}, x::Vector{Float32}, f::Vector{Float32})
+function horner!(f::Vector{Float32}, c::StridedVector{Float32}, x::Vector{Float32})
     @assert length(x) == length(f)
     ccall((:ft_hornerf, libfasttransforms), Cvoid, (Cint, Ptr{Float32}, Cint, Cint, Ptr{Float32}, Ptr{Float32}), length(c), c, stride(c, 1), length(x), x, f)
     f
@@ -69,27 +69,27 @@ function check_clenshaw_points(x, f)
     length(x) == length(f) || throw(ArgumentError("Dimensions must match"))
 end
 
-function clenshaw!(c::StridedVector{Float64}, x::Vector{Float64}, f::Vector{Float64})
+function clenshaw!(f::Vector{Float64}, c::StridedVector{Float64}, x::Vector{Float64})
     @boundscheck check_clenshaw_points(x, f)
     ccall((:ft_clenshaw, libfasttransforms), Cvoid, (Cint, Ptr{Float64}, Cint, Cint, Ptr{Float64}, Ptr{Float64}), length(c), c, stride(c, 1), length(x), x, f)
     f
 end
 
-function clenshaw!(c::StridedVector{Float32}, x::Vector{Float32}, f::Vector{Float32})
+function clenshaw!(f::Vector{Float32}, c::StridedVector{Float32}, x::Vector{Float32})
     @boundscheck check_clenshaw_points(x, f)
     ccall((:ft_clenshawf, libfasttransforms), Cvoid, (Cint, Ptr{Float32}, Cint, Cint, Ptr{Float32}, Ptr{Float32}), length(c), c, stride(c, 1), length(x), x, f)
     f
 end
 
-function clenshaw!(c::StridedVector{Float64}, A::Vector{Float64}, B::Vector{Float64}, C::Vector{Float64}, x::Vector{Float64}, ϕ₀::Vector{Float64}, f::Vector{Float64})
+function clenshaw!(f::Vector{Float64}, c::StridedVector{Float64}, A::Vector{Float64}, B::Vector{Float64}, C::Vector{Float64}, x::Vector{Float64}, ϕ₀::Vector{Float64})
     N = length(c)
     @boundscheck check_clenshaw_recurrences(N, A, B, C)
     @boundscheck check_clenshaw_points(x, ϕ₀, f)
     ccall((:ft_orthogonal_polynomial_clenshaw, libfasttransforms), Cvoid, (Cint, Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}), N, c, stride(c, 1), A, B, C, length(x), x, ϕ₀, f)
     f
 end
 
-function clenshaw!(c::StridedVector{Float32}, A::Vector{Float32}, B::Vector{Float32}, C::Vector{Float32}, x::Vector{Float32}, ϕ₀::Vector{Float32}, f::Vector{Float32})
+function clenshaw!(f::Vector{Float32}, c::StridedVector{Float32}, A::Vector{Float32}, B::Vector{Float32}, C::Vector{Float32}, x::Vector{Float32}, ϕ₀::Vector{Float32})
     N = length(c)
     @boundscheck check_clenshaw_recurrences(N, A, B, C)
     @boundscheck check_clenshaw_points(x, ϕ₀, f)

src/libfasttransforms.jl

@@ -8,18 +8,18 @@ FastTransforms.ft_set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))
         c = one(T) ./ (1:n)
         x = collect(-1 .+ 2*(0:n-1)/T(n))
         f = similar(x)
-        @test FastTransforms.horner!(c, x, f) == f
+        @test FastTransforms.horner!(f, c, x) == f
         fd = T[sum(c[k]*x^(k-1) for k in 1:length(c)) for x in x]
         @test f ≈ fd
-        @test FastTransforms.clenshaw!(c, x, f) == f
+        @test FastTransforms.clenshaw!(f, c, x) == f
         fd = T[sum(c[k]*cos((k-1)*acos(x)) for k in 1:length(c)) for x in x]
         @test f ≈ fd
         A = T[(2k+one(T))/(k+one(T)) for k in 0:length(c)-1]
         B = T[zero(T) for k in 0:length(c)-1]
         C = T[k/(k+one(T)) for k in 0:length(c)]
         phi0 = ones(T, length(x))
         c = FastTransforms.lib_cheb2leg(c)
-        @test FastTransforms.clenshaw!(c, A, B, C, x, phi0, f) == f
+        @test FastTransforms.clenshaw!(f, c, A, B, C, x, phi0) == f
         @test f ≈ fd
     end