update working with different float types

.gitignore

@@ -2,3 +2,4 @@
 /docs/Manifest.toml
 /docs/build/
 .DS_Store
+.vscode

Project.toml

@@ -6,6 +6,7 @@ version = "0.1.4"
 [deps]
 ComplexRegions = "c64915e2-6c82-11e9-38e9-1f159a780463"
 ComplexValues = "41a84b80-6cf2-11e9-379d-9df124847946"
+GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PyFormattedStrings = "5f89f4a4-a228-4886-b223-c468a82ed5b9"
 
@@ -28,10 +29,10 @@ julia = "1.9"
 [extras]
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ComplexRegions = "c64915e2-6c82-11e9-38e9-1f159a780463"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-docs = ["CairoMakie", "ComplexRegions", "DoubleFloats"]
-test = ["Test", "LinearAlgebra", "ComplexRegions"]
+docs = ["CairoMakie"]
+test = ["Test", "LinearAlgebra", "ComplexRegions", "DoubleFloats"]

src/RationalFunctionApproximation.jl

@@ -1,6 +1,6 @@
 module RationalFunctionApproximation
 
-using LinearAlgebra, ComplexRegions
+using LinearAlgebra, GenericLinearAlgebra, ComplexRegions
 using PyFormattedStrings
 
 export Barycentric, nodes, weights, degree, rewind,

src/aaa.jl

@@ -53,16 +53,17 @@ julia> r(1im * π / 2)
 See also [`approximate`](@ref) for approximating a function on a curve or region.
 """
 function aaa(z::AbstractVector{<:Number}, y::AbstractVector{<:Number};
-    max_degree = 150, float_type = Float64, tol = 1000*eps(float_type),
-    lookahead = 10, stats = false
+    max_degree = 150,
+    float_type = promote_type(typeof(float(1)), eltype(z), eltype(y)),
+    tol = 1000*eps(float_type),
+    lookahead = 10,
+    stats = false
     )
 
-    @assert float_type <: AbstractFloat
-    T = float_type
     fmax = norm(y, Inf)    # for scaling
     m = length(z)
     iteration = NamedTuple[]
-    err = T[]
+    err = float_type[]
     besterr, bestidx, best = Inf, NaN, nothing
 
     # Allocate space for Cauchy matrix, Loewner matrix, and residual
@@ -143,17 +144,20 @@ end
 
 function aaa(
     f::Function;
-    max_degree=150, float_type=Float64, tol=1000*eps(float_type),
-    refinement=3, lookahead=10, stats=false
+    max_degree=150,
+    float_type = promote_type(typeof(float(1)), real_type(f(11//23))),
+    tol=1000*eps(float_type),
+    refinement=3,
+    lookahead=10,
+    stats=false
     )
-    @assert float_type <: AbstractFloat
-    T = float_type
-    CT = Complex{T}
+
+    CT = Complex{float_type}
     # arrays for tracking convergence progress
-    err, nbad = T[], Int[]
-    nodes, vals, pol, weights = Vector{T}[], Vector{CT}[], Vector{CT}[], Vector{CT}[]
+    err, nbad = float_type[], Int[]
+    nodes, vals, pol, weights = Vector{float_type}[], Vector{CT}[], Vector{CT}[], Vector{CT}[]
 
-    S = [-one(T), one(T)]                       # initial nodes
+    S = [-one(float_type), one(float_type)]     # initial nodes
     fS = f.(S)
     besterr, bestm = Inf, NaN
     while true                                  # main loop

src/approximate.jl

@@ -60,17 +60,17 @@ end
 approximate(f::Function, d::ComplexCurve; kw...) = approximate(f, Path(d); kw...)
 approximate(f::Function, d::ComplexClosedCurve; kw...) = approximate(f, ClosedPath(d); kw...)
 
+# all methods end up here
 function approximate(f::Function, d::ComplexPath;
     max_degree = 150,
-    float_type = promote_type(typeof(float(1)), real_type(d)),
+    float_type = promote_type(real_type(d), typeof(float(1))),
     tol = 1000*eps(float_type),
     isbad = z->dist(z, d) < tol,
     refinement = 3,
     lookahead = 10,
     stats = false
     )
 
-    @assert float_type <: AbstractFloat
     err, nbad = float_type[], Int[]
     iteration = NamedTuple[]
 

src/methods.jl

@@ -70,10 +70,21 @@ roots(f::Approximation) = roots(f.fun)
 function roots(r::Barycentric)
     wf = r.w_times_f
     m = length(wf)
-    B = diagm( [0; ones(m)] )
+    ty = eltype(nodes(r))
+    B = diagm( [ty(0); ones(m)] )
     # Thanks to Daan Huybrechs:
     E = [0 transpose(wf); ones(m) diagm(nodes(r))]
-    return filter(isfinite, eigvals(E, B))
+    if ty in (Float64,)
+        return filter(isfinite, eigvals(E, B))
+    else
+        # super kludgy; thanks to Daan Huybrechs
+        μ = 11 - 17im  # shift since 0 may well be a root
+        E -= μ*B
+        EB = inv(E)*B
+        pp = eigvals(EB)
+        large_enough = abs.(pp) .> 1e-10
+        return μ .+ 1 ./ pp[large_enough]
+    end
 end
 
 """

test/interval.jl

@@ -1,46 +1,58 @@
-pts = 10 .^ range(-15,0,500); pts = [-reverse(pts); 0; pts]
-approx(f; kw...) = approximate(f, unit_interval; kw...)
 
-@testset "Basic functions" begin
-    f = x -> abs(x + 0.5 + 0.01im); @test pass(f, approx(f), pts, atol=5e-13)
-    f = x -> sin(1/(1.05-x)); @test pass(f, approx(f), pts, atol=2e-13)
-    f = x -> exp(-1/(x^2)); @test pass(f, approx(f), pts, rtol=4e-13)
-    f = x -> exp(-100x^2); @test pass(f, approx(f), pts, rtol=2e-13)
-    f = x -> exp(-10/(1.2-x)); @test pass(f, approx(f), pts, rtol=1e-12)
-    f = x -> 1/(1+exp(100*(x+.5))); @test pass(f, approx(f), pts, atol=2e-13)
-    f = x -> sin(100x) * exp(-10x^2); @test pass(f, approx(f), pts, atol=1e-11)
-    f = x -> abs(x);  @test pass(f, approx(f), pts, atol=1e-8)
-    f = x -> abs(x - 0.95);  @test pass(f, approx(f), pts, atol=1e-6)
+@testset "Basic functions in $T" for T in (Float64, Double64)
+    pts = T(10) .^ range(T(-15), T(0), 500); pts = [-reverse(pts); 0; pts]
+    tol = 2000*eps(T)
+    approx(f; kw...) = approximate(f, Segment{T}(-1, 1); kw...)
+    f = x -> abs(x + 1//2 + 1im//100); @test pass(f, approx(f), pts; rtol=tol)
+    f = x -> sin(1 / (21//20 - x)); @test pass(f, approx(f), pts; rtol=tol)
+    f = x -> exp(-1 / x^2); @test pass(f, approx(f), pts; rtol=tol)
+    f = x -> exp(-100x^2); @test pass(f, approx(f), pts; rtol=tol)
+    f = x -> exp(-10 / (6//5 - x)); @test pass(f, approx(f), pts; rtol=tol)
+    f = x -> 1/(1 + exp(100*(x + 1//2))); @test pass(f, approx(f), pts; rtol=tol)
+    f = x -> sin(100x) * exp(-10x^2); @test pass(f, approx(f), pts; rtol=tol)
+    f = x -> tanh(100*(x - 1//5)); @test pass(f, approx(f), pts, rtol=tol)
+    f = x -> tanh(100x); @test pass(f, approx(f), pts, rtol=tol)
+    f = x -> exp(x); @test pass(f, approx(f), pts, rtol=tol)
+    f = x -> cis(x); @test pass(f, approx(f), pts, rtol=tol)
 end
 
-@testset "Low-accuracy" begin
+@testset "Low accuracy" begin
+    pts = 10.0 .^ range(-15, 0, 500); pts = [-reverse(pts); 0; pts]
+    approx(f; kw...) = approximate(f, unit_interval; kw...)
     f = x -> exp(3x);
-    @test !pass(f, approximate(f, unit_interval, tol=1e-4), pts, atol=1e-8)
+    r = approximate(f, unit_interval, tol=1e-5)
+    @test !pass(f, r, pts, atol=1e-7)
+    @test pass(f, r, pts, atol=5e-5)
+    f = x -> abs(x);  @test pass(f, approx(f), pts, atol=1e-8)
+    f = x -> abs(x - 0.95);  @test pass(f, approx(f), pts, atol=1e-6)
 end
 
-@testset "Poles, zeros, residues" begin
+@testset "Poles, zeros, residues" for T in (Float64,)
+    setprecision(80)
+    approx(f; kw...) = approximate(f, Segment{T}(-1, 1); kw...)
     f = z -> (z+1) * (z+2) / ((z+3) * (z+4))
     r = approx(f)
     pol = poles(r)
     zer = roots(r)
-    @test isapprox(sum(pol+zer), -10, atol=1e-12)
+    @test isapprox(sum(pol+zer), -10, rtol=5000*eps(T))
 
     f = x -> 2/(3 + x) + 5/(x - 2im);  r = approx(f)
-    @test isapprox(prod(values(residues(r))), 10, atol=1e-8)
+    @test isapprox(prod(values(residues(r))), 10, rtol=sqrt(eps(T)))
 
     f = x -> sinpi(10x);  r = approx(f);
-    @test isapprox(sort(abs.(roots(r)))[19], 0.9, atol=1e-12)
+    zer = roots(r)
+    miss = minimum(abs, zer .- 9//10)
+    @test miss < 1e3*eps(T)
 
     f = z -> (z - (3 + 3im))/(z + 2);  r = approx(f)
     pol,zer = poles(r), roots(r)
-    @test isapprox(pol[1]*zer[1], -6-6im, atol=1e-12)
+    @test isapprox(pol[1]*zer[1], -6-6im, rtol=5000*eps(T))
 end
 
-@testset "Vertical scaling" begin
-    f = x -> 1e100*sin(x); @test pass(f, approx(f), pts, rtol=2e-13)
-    f = x -> 1e-100*cos(x); @test pass(f, approx(f), pts, rtol=2e-13)
-    # f = x -> 1e100*sin(x); @test pass(f, approx(f, =5, lawson=20), pts, rtol=1e-6)
-    # f = x -> 1e-100*cos(x); @test pass(f, approx(f, =5, lawson=20), pts, rtol=1e-6)
+@testset "Vertical scaling in $T" for T in (Float64, Double64)
+    pts = T(10) .^ range(T(-15), T(0), 500); pts = [-reverse(pts); 0; pts]
+    f = x -> T(10)^50*sin(x); @test pass(f, approx(f), pts, rtol=2000*eps(T))
+    f = x -> T(10)^(-50)*cos(x); @test pass(f, approx(f), pts, rtol=2000*eps(T))
 end
 
 # @testset "Lawson" begin
@@ -72,17 +84,14 @@ end
     f = x -> 1/(1+1im*x); @test pass(f, approx(f, max_degree=3), pts, atol=2e-13)
     f = x -> 1/(3+x+x^2); @test pass(f, approx(f, max_degree=2), pts, atol=2e-13)
     f = x -> 1/(1.01+x^3); @test pass(f, approx(f, max_degree=3), pts, atol=2e-13)
-    f = x -> tanh(100x); @test pass(f, approx(f), pts, atol=2e-13)
-    f = x -> tanh(100*(x-.2)); @test pass(f, approx(f), pts, atol=2e-13)
-    f = x -> exp(x); @test pass(f, approx(f, tol=1e-13), pts, atol=2e-13)
-    f = x -> cis(x); @test pass(f, approx(f), pts, atol=2e-13)
     f = x -> cis(x); r = approx(f);
     @test sum(@. abs(r(pts))-1)/length(pts) < 2e-13
     f = x -> exp(exp(x))/(x - 0.2im); r = approx(f);
     @test minimum(abs.(poles(r) .- .2im)) < 1e-12
 end
 
 @testset "Other intervals" begin
+    pts = 10 .^ range(-15, 0, 500); pts = [-reverse(pts); 0; pts]
     zz(a,b) = (pts .+ 1)*(b-a)/2 .+ a
     for (a,b) in ((-2,3), (0,1), (-0.01,0) , (-1e4, 2e6))
         for f in ( x -> abs(x + 0.5 + 0.01im), x -> sin(1/(1.05im-x)), x -> exp(-10/(1.2*(b+1)-x)) )

test/runtests.jl

@@ -1,4 +1,4 @@
-using RationalFunctionApproximation, Test, LinearAlgebra, ComplexRegions
+using RationalFunctionApproximation, Test, LinearAlgebra, ComplexRegions, DoubleFloats
 
 # ii = 1im * range(-300,400,500)
 # zz = cispi.(2 * (0:500) / 500)