Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

allow endpoints evaluation for tiny initial intervals #99

Merged
merged 1 commit into from
Jan 13, 2024
Merged

allow endpoints evaluation for tiny initial intervals #99

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 4 additions & 9 deletions src/adapt.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@ function do_quadgk(f::F, s::NTuple{N,T}, n, atol, rtol, maxevals, nrm, segbuf) w
else
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
Expand Down Expand Up @@ -176,15 +175,11 @@ function handle_infinities(workfunc, f::InplaceIntegrand, s)
return workfunc(f, s, identity)
end

function check_endpoint_roundoff(a, b, x; throw_error::Bool=false)
throw_error && a == b && return false # don't throw on empty intervals
function check_endpoint_roundoff(a, b, x)
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)"))
eval_at_a = a == a + (1+x[1])*c
eval_at_b = b == a + (1-x[1])*c
return eval_at_a || eval_at_b
end

# Gauss-Kronrod quadrature of f from a to b to c...
Expand Down
1 change: 0 additions & 1 deletion src/batch.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,6 @@ 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
Expand Down
21 changes: 2 additions & 19 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,8 +27,9 @@ using QuadGK, LinearAlgebra, Test
# order=1 (issue #66)
@test quadgk_count(x -> 1, 0, 1, order=1) == (1.0, 0.0, 3)

# empty intervals (issue #97)
# empty or nearly empty intervals (issue #97)
@test quadgk(x -> 1, 0,0) == quadgk(x -> 1, 0,0,0) == (0.0,0.0)
@test quadgk(x -> 1, 1,nextfloat(1.0)) == (eps(),0.0)
end

module Test19626
Expand Down Expand Up @@ -348,22 +349,4 @@ end
@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
Loading