Fix cheb2leg accuracy

Author: dlfivefifty

src/toeplitzhankel.jl

@@ -102,6 +102,36 @@ function hankel_partialchol(v::Vector{T}) where T
     C[:,1:r]
 end
 
+# cholesky for D .* H .* D'
+function hankel_partialchol(v::Vector, D::AbstractVector)
+    T = promote_type(eltype(v), eltype(D))
+    # Assumes positive definite
+    σ = T[]
+    n = isempty(v) ? 0 : (length(v)+2) ÷ 2
+    C = Matrix{T}(undef, n, n)
+    d = v[1:2:end] .* D.^2 # diag of D .* H .* D'
+    @assert length(v) ≥ 2n-1
+    reltol = maximum(abs,d)*eps(T)*log(n)
+    r = 0
+    for k = 1:n
+        mx,idx = findmax(d)
+        if mx ≤ reltol break end
+        push!(σ, inv(mx))
+        C[:,k] .= view(v,idx:n+idx-1) .*D.*D[idx]
+        for j = 1:k-1
+            nCjidxσj = -C[idx,j]*σ[j]
+            LinearAlgebra.axpy!(nCjidxσj, view(C,:,j), view(C,:,k))
+        end
+        @inbounds for p=1:n
+            d[p] -= C[p,k]^2/mx
+        end
+        r += 1
+    end
+    for k=1:length(σ) rmul!(view(C,:,k), sqrt(σ[k])) end
+    C[:,1:r]
+end
+
+
 
 # Diagonally-scaled Toeplitz∘Hankel polynomial transforms
 
@@ -209,9 +239,10 @@ function _cheb2legTH_TLC(::Type{S}, mn, d) where S
         t[1:2:end] = Λ.(0:one(S̃):div(n-2,2), -half(S̃), one(S̃))
     end
     h = Λ.(1:half(S̃):n-1, zero(S̃), 3half(S̃))
-    DL = (3half(S̃):n-half(S̃))
-    DR = -(one(S̃):n-one(S̃))./4
-    C = hankel_partialchol(h)
+    D = 1:one(S):n-1
+    DL = (3half(S̃):n-half(S̃)) ./ D
+    DR = -(one(S̃):n-one(S̃)) ./ (4 .* D)
+    C = hankel_partialchol(h, D)
     T = plan_uppertoeplitz!(t, (_sub_dim_by_one(d, mn...)..., size(C,2)), d)
     T, DL .* C, DR .* C
 end

test/toeplitzhankeltests.jl

@@ -4,7 +4,7 @@ import FastTransforms: th_leg2cheb, th_cheb2leg, th_ultra2ultra,th_jac2jac, th_l
                         plan_th_cheb2leg!, plan_th_leg2cheb!
 
 @testset "ToeplitzHankel" begin
-    for x in ([1.0], [1.0,2,3,4,5], [1.0+im,2-3im,3+4im,4-5im,5+10im])
+    for x in ([1.0], [1.0,2,3,4,5], [1.0+im,2-3im,3+4im,4-5im,5+10im], collect(1.0:1000))
         @test th_leg2cheb(x) ≈ lib_leg2cheb(x)
         @test th_cheb2leg(x) ≈ lib_cheb2leg(x)
         @test th_leg2chebu(x) ≈ lib_ultra2ultra(x, 0.5, 1.0)
@@ -14,6 +14,9 @@ import FastTransforms: th_leg2cheb, th_cheb2leg, th_ultra2ultra,th_jac2jac, th_l
 
         @test all(th_leg2cheb(x) .=== leg2cheb(x))
         @test all(th_cheb2leg(x) .=== cheb2leg(x))
+
+        @test th_cheb2leg(th_leg2cheb(x)) ≈ x atol=1E-9
+        @test th_leg2leg(th_cheb2cheb(x)) ≈ x atol=1E-11
     end
 
     for X in (randn(5,4), randn(5,4) + im*randn(5,4))
@@ -37,4 +40,13 @@ import FastTransforms: th_leg2cheb, th_cheb2leg, th_ultra2ultra,th_jac2jac, th_l
         @test leg2cheb(x) ≈ lib_leg2cheb(x)
         @test cheb2leg(x) ≈ lib_cheb2leg(x)
     end
+
+    @testset "jishnub example" begin
+        x = chebyshevpoints(4096);
+        f = x -> cospi(1000x);  
+        y = f.(x);
+        v = cheb2leg(chebyshevtransform(y))
+        @test norm(v - cheb2leg(leg2cheb(v)), Inf) ≤ 1E-13
+        @test norm(v - cheb2leg(leg2cheb(v)))/norm(v) ≤ 1E-14
+    end
 end