Fix array evaluation at nodes

src/barycentric.jl

@@ -167,12 +167,19 @@ function evaluate!(u::AbstractArray, r::Barycentric, z::AbstractArray{<:Number})
     num = @. r.w_times_f[1] / (z - r.nodes[1])
     den = @. r.weights[1] / (z - r.nodes[1])
     safe = trues(size(z))
-    idx = isinf.(num)
+    # Infinite arguments take the same limit as the scalar method.
+    idx = isinf.(z)
+    safe[idx] .= false
+    u[idx] .= sum(r.w_times_f) / sum(r.weights)
+    # A point at a node interpolates exactly. Find it by comparison rather than by
+    # inspecting the quotient: 1/(0+0im) is NaN+NaN*im rather than Inf, and a node whose
+    # value is zero makes the numerator 0/0.
+    idx = z .== r.nodes[1]
     safe[idx] .= false
     u[idx] .= r.values[1]
     @inbounds for k in 2:length(r.nodes)
         c = @. 1 / (z - r.nodes[k])
-        idx = isinf.(c)
+        idx = z .== r.nodes[k]
         safe[idx] .= false
         u[idx] .= r.values[k]
         axpy!(r.w_times_f[k], c, num)

src/thiele.jl

@@ -157,6 +157,8 @@ function evaluate!(t::AbstractArray, r::Thiele, z::AbstractArray{<:Number},
             a .= t
         end
         t .= r.weights[1] .+ b ./ a
+        # retroactively correct at infinite inputs
+        t[isinf.(z)] .= r(Inf)
     end
     return t
 end

test/RFATests.jl

@@ -3,7 +3,17 @@ module RFATests
 using RationalFunctionApproximation, ReTest, ComplexRegions, DoubleFloats, IntervalSets, Logging
 const RFA = RationalFunctionApproximation
 
-pass(f, r, z; kw...) = isapprox(f.(z), r.(z), norm=u->maximum(abs, u); kw...)
+# Equivalent to `isapprox(f.(z), r.(z), norm=u->maximum(abs, u); kw...)`, written out
+# because Julia 1.13.0-rc3's `isapprox` for arrays ignores the `norm` keyword
+# (JuliaLang/LinearAlgebra.jl#1675).
+function pass(f, r, z; atol::Real=0, rtol::Real=-1)
+    fz, rz = f.(z), r.(z)
+    if rtol < 0    # same default as `isapprox`
+        rtol = atol > 0 ? 0 : sqrt(eps(float(real(promote_type(eltype(fz), eltype(rz))))))
+    end
+    nrm(u) = maximum(abs, u)
+    return nrm(fz - rz) <= max(atol, rtol * max(nrm(fz), nrm(rz)))
+end
 logger = Logging.SimpleLogger(stderr, Logging.Error)
 global_logger(logger)
 

test/circle.jl

@@ -13,9 +13,24 @@
     end
 
     @testset "Unit circle for $method" for method in (Barycentric, Thiele)
-        f = z -> abs(z-1im);  @test pass(f, approximate(f, UC, Barycentric()), pts, rtol=2e-10)
-        f = z -> tan(π*z);  @test pass(f, approximate(f, UC, Barycentric()), pts, rtol=2e-13)
-        f = z -> tanh(100z); @test pass(f, approximate(f, UC, Barycentric()), pts, rtol=2e-13)
+        f = z -> abs(z-1im);  @test pass(f, approximate(f, UC, method(); stagnation=20), pts, rtol=2e-10)
+        f = z -> tan(π*z);  @test pass(f, approximate(f, UC, method()), pts, rtol=2e-13)
+        f = z -> tanh(100z); @test pass(f, approximate(f, UC, method()), pts, rtol=2e-12)
+    end
+
+    @testset "Array evaluation for Barycentric" begin
+        f = z -> sin(10z) * exp(-z^2)
+        r = approximate(f, UC, Barycentric())
+        @test isapprox(f.(pts), r(pts), norm=u->maximum(abs, u), rtol=2e-11)
+        # array evaluation at a node must interpolate, on or off the real line
+        z = nodes(r)
+        @test r(z) ≈ values(r)
+        @test r([z[1], 0.5]) ≈ [values(r)[1], r(0.5)]
+        # a node whose value is zero makes the numerator 0/0
+        s = Barycentric([0.0, 1.0, 2.0], [0.0, 2.0, 4.0], [1.0, -2.0, 1.0])
+        @test s([0.0, 1.0]) ≈ [0.0, 2.0]
+        # infinite arguments take the same limit as the scalar method
+        @test r([Inf]) ≈ [r(Inf)]
     end
 
     @testset "Array evaluation for Thiele" begin