avoid integrand evaluations at endpoints (#91)

* avoid integrand evaluations at endpoints

* improve error message

* remove type instability from evalrule and refactor

* increase order in issue 86 test

* fix test on CI

Author: lxvm

src/adapt.jl

@@ -10,7 +10,11 @@ function do_quadgk(f::F, s::NTuple{N,T}, n, atol, rtol, maxevals, nrm, segbuf) w
     if f isa BatchIntegrand
         segs = evalrules(f, s, x,w,gw, nrm)
     else
-        segs = ntuple(i -> evalrule(f, s[i],s[i+1], x,w,gw, nrm), Val{N-1}())
+        segs = ntuple(Val{N-1}()) do i
+            a, b = s[i], s[i+1]
+            check_endpoint_roundoff(a, b, x, throw_error=true)
+            evalrule(f, a,b, x,w,gw, nrm)
+        end
     end
     if f isa InplaceIntegrand
         I = f.I .= segs[1].I
@@ -57,8 +61,16 @@ end
 function refine(f::F, segs::Vector{T}, I, E, numevals, x,w,gw,n, atol, rtol, maxevals, nrm) where {F, T}
     s = heappop!(segs, Reverse)
     mid = (s.a + s.b) / 2
+
+    # early return if integrand evaluated at endpoints
+    if check_endpoint_roundoff(s.a, mid, x) || check_endpoint_roundoff(mid, s.b, x)
+        heappush!(segs, s, Reverse)
+        return segs
+    end
+
     s1 = evalrule(f, s.a, mid, x,w,gw, nrm)
     s2 = evalrule(f, mid, s.b, x,w,gw, nrm)
+
     if f isa InplaceIntegrand
         I .= (I .- s.I) .+ s1.I .+ s2.I
     else
@@ -164,6 +176,16 @@ function handle_infinities(workfunc, f::InplaceIntegrand, s)
     return workfunc(f, s, identity)
 end
 
+function check_endpoint_roundoff(a, b, x; throw_error::Bool=false)
+    c = convert(eltype(x), 0.5) * (b-a)
+    (eval_at_a = a == a + (1+x[1])*c) && throw_error && throw_endpoint_error(a, a, b)
+    (eval_at_b = b == a + (1-x[1])*c) && throw_error && throw_endpoint_error(b, a, b)
+    eval_at_a || eval_at_b
+end
+function throw_endpoint_error(x, a, b)
+    throw(DomainError(x, "roundoff error detected near endpoint of the initial interval ($a, $b)"))
+end
+
 # Gauss-Kronrod quadrature of f from a to b to c...
 
 """

src/batch.jl

@@ -69,6 +69,7 @@ function evalrules(f::BatchIntegrand, s::NTuple{N}, x,w,gw, nrm) where {N}
     resize!(f.y, n)
     for i in 1:(N-1)    # fill buffer with evaluation points
         a = s[i]; b = s[i+1]
+        check_endpoint_roundoff(a, b, x, throw_error=true)
         c = convert(eltype(x), 0.5) * (b-a)
         o = (i-1)*m
         f.x[l+o] = a + c
@@ -110,10 +111,12 @@ function refine(f::BatchIntegrand, segs::Vector{T}, I, E, numevals, x,w,gw,n, at
         s = segs[len-i+1]
         mid = (s.a+s.b)/2
         for (j,a,b) in ((2,s.a,mid), (1,mid,s.b))
+            check_endpoint_roundoff(a, b, x) && return segs
             c = convert(eltype(x), 0.5) * (b-a)
             o = (2i-j)*m
             f.x[l+o] = a + c
             for k in 1:l-1
+                # early return if integrand evaluated at endpoints
                 f.x[k+o] = a + (1 + x[k]) * c
                 f.x[m+1-k+o] = a + (1 - x[k]) * c
             end

test/runtests.jl

@@ -338,3 +338,29 @@ end
     A = Test89()
     @test quadgk(A, 0, 10) == (50, 0)
 end
+
+# issue 86: nodes roundoff to endpoints
+@testset "issue 86" begin
+    I, = quadgk(y->10*y^9/(y^10)^1.1,1,Inf)[1]
+    @test quadgk(x->1/x^1.1,1,Inf)[1] ≈ I rtol=0.05
+    @test quadgk!((y,x)-> y .= 1/x^1.1,[0.0],1,Inf)[1][1] ≈ I rtol=0.05
+    @test quadgk(BatchIntegrand{Float64}((y,x)-> @.(y = 1/x^1.1)),1,Inf)[1] ≈ I rtol=0.05
+
+    # if order < 85, there is also a DomainError, but due to overflow of the change of variables
+    errmsg = "roundoff error detected near endpoint of the initial interval"
+    a = Float16(1)
+    b = Float16(Inf)
+    for (routine, args) in (
+        (quadgk,  (x -> x,)),
+        (quadgk!, ((y,x) -> x, Float16[1])),
+        (quadgk,  (BatchIntegrand{Float16}((y,x) -> y .= x),)),
+    )
+        # need Julia 1.8 for @test_throws with the error message
+        try
+            routine(args..., a, b; order=85)
+            error()
+        catch e
+            @test startswith(e.msg, errmsg)
+        end
+    end
+end