From bfdb7677610dc8093e10175b75990274b7662373 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Richard?= Date: Wed, 29 Apr 2026 15:03:24 +0200 Subject: [PATCH 01/10] Complete out-of-domain tests --- test/roots.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/roots.jl b/test/roots.jl index a5a2ac8..6e9cc32 100644 --- a/test/roots.jl +++ b/test/roots.jl @@ -159,9 +159,19 @@ end end @testset "Out of domain" begin + s(x) = sqrt(x^2 - 4) - 2 # root: -6 and 6 + a(xy) = [asin(xy[1]), acos(xy[2]) - π/2] # root: (0, 0) + for contractor in newtonlike_methods @test length(roots(log, interval(-100, 2) ; contractor)) == 1 @test length(roots(log, interval(-100, -1) ; contractor)) == 0 + + @test length(roots(s, interval(-20, 20) ; contractor)) == 2 + @test length(roots(s, interval(0, 20) ; contractor)) == 1 + @test length(roots(s, interval(-1, 1) ; contractor)) == 0 + + @test length(roots(a, fill(interval(-20, 20), 2) ; contractor)) == 1 + @test length(roots(a, fill(interval(0.5, 20), 2) ; contractor)) == 0 end end From a3b8987903e7d5eef06938008bb998aa1b30b862 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Richard?= Date: Wed, 29 Apr 2026 15:37:49 +0200 Subject: [PATCH 02/10] Fix bisect_region to return com interval when bisecting infinite intervals --- Project.toml | 9 +-------- src/region.jl | 15 ++++++++++++++- test/Project.toml | 22 ++++++++++++++++++++++ test/roots.jl | 34 +++++++++++++++++++++++++--------- test/runtests.jl | 1 + 5 files changed, 63 insertions(+), 18 deletions(-) create mode 100644 test/Project.toml diff --git a/Project.toml b/Project.toml index 8ab743c..8f53ca0 100644 --- a/Project.toml +++ b/Project.toml @@ -18,11 +18,4 @@ InteractiveUtils = "1.10.0" IntervalArithmetic = "1.0.3" Reexport = "1" StaticArrays = "1" -julia = "1" - -[extras] -Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[targets] -test = ["Test", "Polynomials"] +julia = "1" \ No newline at end of file diff --git a/src/region.jl b/src/region.jl index 64ea4a1..630b708 100644 --- a/src/region.jl +++ b/src/region.jl @@ -29,7 +29,20 @@ mag_region(X::AbstractVector) = maximum(mag, X) mig_region(X::Interval) = mig(X) mig_region(X::AbstractVector) = minimum(mig, X) -bisect_region(X::Interval, α) = bisect(X, α) +function bisect_region(X::Interval, α) + X1, X2 = bisect(X, α) + + isbounded(X) && return X1, X2 + + if isbounded(X1) + X1 = IntervalArithmetic.setdecoration(X1, com) + end + + if isbounded(X2) + X2 = IntervalArithmetic.setdecoration(X2, com) + end + return X1, X2 +end function bisect_region(X::AbstractVector, α) X1 = copy(X) diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 0000000..5cae9d0 --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,22 @@ +[deps] +BranchAndPrune = "d3bc4f2e-91e6-11e9-365e-cd067da536ce" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253" +IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[compat] +BranchAndPrune = "0.2.1" +ForwardDiff = "0.10, 1" +InteractiveUtils = "1.10.0" +IntervalArithmetic = "1.0.3" +Reexport = "1" +StaticArrays = "1" +Suppressor = "0.2.8" +julia = "1" diff --git a/test/roots.jl b/test/roots.jl index 6e9cc32..4f3856f 100644 --- a/test/roots.jl +++ b/test/roots.jl @@ -163,22 +163,38 @@ end a(xy) = [asin(xy[1]), acos(xy[2]) - π/2] # root: (0, 0) for contractor in newtonlike_methods - @test length(roots(log, interval(-100, 2) ; contractor)) == 1 - @test length(roots(log, interval(-100, -1) ; contractor)) == 0 + @suppress @test length(roots(log, interval(-100, 2) ; contractor)) == 1 + @suppress @test length(roots(log, interval(-100, -1) ; contractor)) == 0 - @test length(roots(s, interval(-20, 20) ; contractor)) == 2 - @test length(roots(s, interval(0, 20) ; contractor)) == 1 - @test length(roots(s, interval(-1, 1) ; contractor)) == 0 + @suppress @test length(roots(s, interval(-20, 20) ; contractor)) == 2 + @suppress @test length(roots(s, interval(0, 20) ; contractor)) == 1 + @suppress @test length(roots(s, interval(-1, 1) ; contractor)) == 0 - @test length(roots(a, fill(interval(-20, 20), 2) ; contractor)) == 1 - @test length(roots(a, fill(interval(0.5, 20), 2) ; contractor)) == 0 + @suppress @test length(roots(a, fill(interval(-20, 20), 2) ; contractor)) == 1 + @suppress @test length(roots(a, fill(interval(0.5, 20), 2) ; contractor)) == 0 end end @testset "Infinite domain" begin + f(xy) = [sin(1/xy[1]), xy[1]*xy[2]] # largest roots: [1/π, 0] and [infinity, 0] + for contractor in newtonlike_methods - rts = roots(x -> x^2 - 2, interval(-Inf, Inf) ; contractor) - @test length(filter(isunique, rts)) == 2 + rts = @suppress roots(x -> x^2 - 2, interval(-Inf, Inf) ; contractor) + @test length(rts) == 2 + @test all(isbounded, [rt.region for rt in rts]) + @test all(rt.region.decoration == com for rt in rts) + @test all(isunique, rts) + + rts = @suppress roots(f, [interval(0.9/π, Inf), interval(-Inf, Inf)] ; contractor) + @test length(rts) == 2 + sort!(rts, by = mig) + @test isunique(rts[1]) + @test all(isbounded, rts[1].region) + @test minimum(X.decoration for X in rts[1].region) == com + + @test !isunique(rts[2]) + @test !all(isbounded, rts[2].region) + @test minimum(X.decoration for X in rts[2].region) == dac end end diff --git a/test/runtests.jl b/test/runtests.jl index 7367f5d..c1de6b8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,7 @@ using IntervalRootFinding using IntervalArithmetic.Symbols using BranchAndPrune +using Suppressor using Test newtonlike_methods = [Newton, Krawczyk] From 7cbde132709207f2835e65251af1546728b8ffc4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Richard?= Date: Thu, 30 Apr 2026 16:06:55 +0200 Subject: [PATCH 03/10] Revert decoration change for bisection --- src/region.jl | 15 +-------------- test/roots.jl | 6 +++--- test/runtests.jl | 3 ++- 3 files changed, 6 insertions(+), 18 deletions(-) diff --git a/src/region.jl b/src/region.jl index 630b708..64ea4a1 100644 --- a/src/region.jl +++ b/src/region.jl @@ -29,20 +29,7 @@ mag_region(X::AbstractVector) = maximum(mag, X) mig_region(X::Interval) = mig(X) mig_region(X::AbstractVector) = minimum(mig, X) -function bisect_region(X::Interval, α) - X1, X2 = bisect(X, α) - - isbounded(X) && return X1, X2 - - if isbounded(X1) - X1 = IntervalArithmetic.setdecoration(X1, com) - end - - if isbounded(X2) - X2 = IntervalArithmetic.setdecoration(X2, com) - end - return X1, X2 -end +bisect_region(X::Interval, α) = bisect(X, α) function bisect_region(X::AbstractVector, α) X1 = copy(X) diff --git a/test/roots.jl b/test/roots.jl index 4f3856f..0445fa6 100644 --- a/test/roots.jl +++ b/test/roots.jl @@ -104,7 +104,7 @@ end # Infinite interval X = [interval(-Inf, Inf), interval(-Inf, Inf)] - rts = roots(f, X ; contractor = Newton) + @suppress rts = roots(f, X ; contractor = Newton) @test length(rts) == 2 end @@ -182,7 +182,7 @@ end rts = @suppress roots(x -> x^2 - 2, interval(-Inf, Inf) ; contractor) @test length(rts) == 2 @test all(isbounded, [rt.region for rt in rts]) - @test all(rt.region.decoration == com for rt in rts) + @test all(rt.region.decoration == dac for rt in rts) @test all(isunique, rts) rts = @suppress roots(f, [interval(0.9/π, Inf), interval(-Inf, Inf)] ; contractor) @@ -190,7 +190,7 @@ end sort!(rts, by = mig) @test isunique(rts[1]) @test all(isbounded, rts[1].region) - @test minimum(X.decoration for X in rts[1].region) == com + @test minimum(X.decoration for X in rts[1].region) == dac @test !isunique(rts[2]) @test !all(isbounded, rts[2].region) diff --git a/test/runtests.jl b/test/runtests.jl index c1de6b8..48939b3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,7 @@ +using BranchAndPrune using IntervalRootFinding using IntervalArithmetic.Symbols -using BranchAndPrune +using ForwardDiff using Suppressor using Test From f7c956bca74e033bb0f6e11109589b95fca08251 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Richard?= Date: Thu, 30 Apr 2026 16:24:02 +0200 Subject: [PATCH 04/10] Add 1D example to NaN suite --- test/roots.jl | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/test/roots.jl b/test/roots.jl index 0445fa6..ea6e9b3 100644 --- a/test/roots.jl +++ b/test/roots.jl @@ -199,9 +199,17 @@ end end @testset "NaN return value" begin - f(xx) = ( (x, y) = xx; [log(y/x) + 3x, y - 2x] ) + f1(x) = x^2/x - 1 + f(xy) = [log(xy[2]/xy[1]) + 3xy[1], xy[2] - 2xy[1]] X = [interval(-100, 100), interval(-100, 100)] + for contractor in newtonlike_methods + rts = roots(f1, interval(-2, 2) ; contractor) + @test length(rts) == 2 + @test !isunique(rts[1]) + @test in_interval(0, root_region(rts[1])) + @test isunique(rts[2]) + rts = roots(f, X ; contractor) @test length(filter(isunique, rts)) == 1 @test length(filter(x -> all(in_interval.(0, x)), root_region.(rts))) == 1 From 97be102f3ac4a5ebdb6c61c696b750c7d03b4160 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Richard?= Date: Sun, 3 May 2026 19:34:04 +0200 Subject: [PATCH 05/10] Use tests from findroots.jl --- test/findroots.jl | 171 ---------------------------------------------- test/roots.jl | 71 ++++++++++++++++++- test/wilkinson.jl | 54 --------------- 3 files changed, 70 insertions(+), 226 deletions(-) delete mode 100644 test/findroots.jl delete mode 100644 test/wilkinson.jl diff --git a/test/findroots.jl b/test/findroots.jl deleted file mode 100644 index 602ff28..0000000 --- a/test/findroots.jl +++ /dev/null @@ -1,171 +0,0 @@ - -using IntervalArithmetic, IntervalRootFinding -using ForwardDiff -using Test - -const D = IntervalRootFinding.derivative - -include("wilkinson.jl") - - -function generate_wilkinson(n)#, T=BigFloat) # SLOW - - p = poly(collect(1:n)) - p_prime = polyder(p) - - coeffs = p.a #[convert(T, coeff) for coeff in p.a] - # coeffs_prime = [convert(T, coeff) for coeff in p_prime.a] - - function f(x) - total = coeffs[1] - for i in 2:length(coeffs) - total += coeffs[i]*x^(i-1) - end - total - end - - f -end - - -setprecision(Interval, Float64) -float_pi = @interval(pi) - -setprecision(Interval, 10000) -big_pi = @interval(pi) -# Using precision "only" 256 leads to overestimation of the true roots for `cos` -# i.e the Newton method gives more accurate results! - -half_pi = big_pi / 2 -three_halves_pi = 3*big_pi/2 - - -# Format: (function, derivative, lower_bound, upper_bound, [true_roots]) -function_list = [ - (sin, cos, -5, 5, [ -big_pi, @interval(0), big_pi ] ) , - (cos, x->D(cos, x), -7.5, 7.5, [ -three_halves_pi, -half_pi, half_pi, three_halves_pi ] ), - (W₃, x->D(W₃, x), -10, 10, [ @interval(1), @interval(2), @interval(3) ] ), - (W₇, x->D(W₇, x), -10, 10, [ @interval(i) for i in 1:7 ] ), - (x->exp(x)-2, y->D(x->exp(x),y), -20, 20, [log(@biginterval(2))] ), - (x->asin(sin(x)) - 0.1, y->D(x->asin(sin(x)),y), 0, 1, [@biginterval(0.1)]) - ] - - -@testset "Testing root finding" begin - - for rounding_type in (:wide, :narrow) - @testset "Interval rounding: $rounding_type" begin - # setrounding(Interval, rounding_type) - - for prec in ( (BigFloat,53), (BigFloat,256), (Float64,64) ) - @testset "Precision: $prec" begin - setprecision(Interval, prec) - - for method in (IntervalRootFinding.newton, IntervalRootFinding.krawczyk) - @testset "Method $method" begin - - for func in function_list - - f, f_prime, a_lower, a_upper, true_roots = func - a = @interval(a_lower, a_upper) - - @testset "Function $f; interval $a" begin - - for autodiff in (false, true) - - @testset "With autodiff=$autodiff" begin - - if autodiff - rts = find_roots(f, a, method) - else - rts = find_roots(f, f_prime, a, method) - end - - @test length(rts) == length(true_roots) - - for i in 1:length(rts) - root = rts[i] - - @test isa(root, Root) - @test is_unique(root) - @test true_roots[i] ⊆ root.interval - end - end - end - end - end - end - end - end - end - end - end -end - - - -setprecision(Interval, Float64) - -@testset "find_roots tests" begin - f(x) = x^2 - 2 - - rts = IntervalRootFinding.newton(f, @interval(10, 11)) - @test length(rts) == 0 - - - rts = IntervalRootFinding.find_roots_midpoint(f, -5, 5) - @test length(rts) == 3 - @test length(rts[1]) == 2 - - rts = find_roots(f, -5, 5) - @test length(rts) == 2 - - setprecision(Interval, 256) - - for method in (IntervalRootFinding.newton, IntervalRootFinding.krawczyk) - new_roots = method(f, rts) - @test length(new_roots) == length(rts) - - for (root1, root2) in zip(new_roots, rts) - @test root1 ⊆ root2 - end - end -end - -@testset "Multiple roots" begin - setprecision(Interval, Float64) - let - f(x) = (x-1) * (x^2 - 2)^3 * (x^3 - 2)^4 - - rts = IntervalRootFinding.newton(f, -5..5.1, maxlevel=1000) - - @test length(rts) == 4 - @test rts[1].status == :unknown - @test rts[1].interval == Interval(-1.4142135623730954, -1.414213562373095) - @test rts[3].interval == Interval(1.259921049894873, 1.2599210498948734) - - end -end - -@testset "Iterated logistic function fixed points" begin - ∘(f, g) = x -> f(g(x)) - iterate(f, n) = n == 1 ? f : f ∘ iterate(f, n-1) - - f(x) = 4x*(1-x) # logistic map - - for n in 1:10 - g = x -> iterate(f, n)(x) - x # look for fixed points of f^n - - rts = IntervalRootFinding.newton(g, 0..1) - @test length(rts) == 2^n - - rts = IntervalRootFinding.krawczyk(g, 0..1) - @test length(rts) == 2^n - - end -end - - - -# Example of a function with a double root at 0 from Burden & Faires, 9th ed, p.84 -# exp(x) - x - 1 diff --git a/test/roots.jl b/test/roots.jl index ea6e9b3..b4dd9b9 100644 --- a/test/roots.jl +++ b/test/roots.jl @@ -158,6 +158,31 @@ end end end +@testset "Multiple roots" begin + f(x) = (x-1) * (x^2 - 2)^3 * (x^3 - 2)^4 + g(x) = exp(x) - x - 1 # Double root at 0 from Burden & Faires, 9th ed, p.84 + + for contractor in newtonlike_methods + rts = roots(f, interval(-5, 5) ; contractor) + + @test length(rts) == 4 + @test rts[1].status == :unknown + @test in_interval(-sqrt(2), rts[1].region) + + @test rts[2].status == :unique + @test in_interval(1, rts[2].region) + + @test rts[3].status == :unknown + @test in_interval(cbrt(2), rts[3].region) + + @test rts[4].status == :unknown + @test in_interval(sqrt(2), rts[4].region) + + rts = roots(g, interval(-5, 5) ; contractor) + @test all(rt.status == :unknown for rt in rts) + @test in_interval(0, hull([rt.region for rt in rts]...)) + end +end @testset "Out of domain" begin s(x) = sqrt(x^2 - 4) - 2 # root: -6 and 6 a(xy) = [asin(xy[1]), acos(xy[2]) - π/2] # root: (0, 0) @@ -446,4 +471,48 @@ end g(x) = (x < 1 ? x : error()) rts = roots(g, interval(-10, 10) ; ignored_errors = [ErrorException, IntervalArithmetic.InconclusiveBooleanOperation]) @test any(isunique, rts) -end \ No newline at end of file +end + +# Wilkinson-type polynomial defined by its roots: +# Wn(x) = (x-1)⋅(x-2)⋅ ⋯ ⋅ (x-n) +W3(x) = prod(x .- (1:3)) +W7(x) = prod(x .- (1:7)) + +# Format: (function, derivative, lower_bound, upper_bound, [true_roots]) +function_list = [ + (sin, cos, -5, 5, [-big(π), 0, big(π)]) , + (cos, x -> -sin(x), -7.5, 7.5, [-3big(π)/2, -big(π)/2, big(π)/2, 3big(π)/2]), + (W3, nothing, -10, 10, [1, 2, 3]), + (W7, nothing, -10, 10, collect(1:7)), + (x -> exp(x) - 2, exp, -20, 20, [log(big(2))]), + (x -> asin(sin(x)) - big"0.1", Returns(1.0), 0, 1, [big"0.1"]) +] + +@testset "Various precisions" begin + @testset "$T" for T in [Float64, BigFloat] + @testset "Function $(func[1])" for func in function_list + @testset "$contractor" for contractor in newtonlike_methods + f, f_prime, a_lower, a_upper, true_roots = func + a = interval(T, a_lower, a_upper) + + @testset "autodiff = $autodiff" for autodiff in (false, true) + isnothing(f_prime) && continue + + if autodiff + rts = roots(f, a ; contractor) + else + rts = roots(f, a ; contractor, derivative = f_prime) + end + + @test length(rts) == length(true_roots) + + for (root, true_root) in zip(rts, true_roots) + @test isa(root, Root) + @test isunique(root) + @test in_interval(true_root, root.region) + end + end + end + end + end +end diff --git a/test/wilkinson.jl b/test/wilkinson.jl deleted file mode 100644 index 0ef6b8d..0000000 --- a/test/wilkinson.jl +++ /dev/null @@ -1,54 +0,0 @@ -# Wilkinson-type polynomial defined by its roots: -# wₙ(x) = (x-1)⋅(x-2)⋅ ⋯ ⋅ (x-n) - -using Polynomials - - -# function generate_wilkinson(n)#, T=BigFloat) # SLOW - -# p = poly(collect(1:n)) -# p_prime = polyder(p) - -# coeffs = p.a #[convert(T, coeff) for coeff in p.a] -# # coeffs_prime = [convert(T, coeff) for coeff in p_prime.a] - -# function f(x) -# total = coeffs[1] -# for i in 2:length(coeffs) -# total += coeffs[i]*x^(i-1) -# end -# total -# end - -# f -# end - -function generate_wilkinson_horner(n) - p = poly(collect(1:n)) # make polynomial with roots 1 to n - - coeffs = p.a - - expr = :($(coeffs[end])) - - for i in length(coeffs)-1:-1:1 - expr = :($expr*x + $(coeffs[i])) - end - - expr -end - -function subscriptify(n::Int) - subscript_digits = [c for c in "₀₁₂₃₄₅₆₇₈₉"] - dig = reverse(digits(n)) - join([subscript_digits[i+1] for i in dig]) -end - -# Generate Wilkinson functions W₃ etc.: -for n in 1:15 - - fn_name = Symbol("W", subscriptify(n)) - - expr = generate_wilkinson_horner(n) - - @eval $(fn_name)(x) = $(expr) -end From 19305af2853e53ed53b16b8e2aa38207cf1806ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Richard?= Date: Sun, 3 May 2026 19:45:48 +0200 Subject: [PATCH 06/10] Mini clean up of import in tests --- test/roots.jl | 11 ++++------- test/runtests.jl | 1 + 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/test/roots.jl b/test/roots.jl index b4dd9b9..787dc35 100644 --- a/test/roots.jl +++ b/test/roots.jl @@ -1,6 +1,3 @@ -using IntervalArithmetic, IntervalRootFinding, StaticArrays, ForwardDiff -using Test - function all_unique(rts) all(root_status.(rts) .== :unique) end @@ -28,6 +25,10 @@ function test_newtonlike(f, derivative, X, contractor, nsol, tol=1e-10) @test sum(roots_dist.(rts, roots(f, X ; contractor, derivative))) < tol end +function in_solution_set(point, solution_intervals) + return any(map(Y -> in_region(point, Y), solution_intervals)) +end + @testset "1D roots" begin # Default rts = roots(sin, interval(-5, 5)) @@ -79,10 +80,6 @@ end end -function in_solution_set(point, solution_intervals) - return any(map(Y -> in_region(point, Y), solution_intervals)) -end - @testset "2D roots" begin f(x, y) = [x^2 + y^2 - 1, y - 2x] f(X) = f(X...) diff --git a/test/runtests.jl b/test/runtests.jl index 48939b3..3cb8f68 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,6 +2,7 @@ using BranchAndPrune using IntervalRootFinding using IntervalArithmetic.Symbols using ForwardDiff +using StaticArrays using Suppressor using Test From a5c55491c863ce52892ecf57058b546b2c699ba2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Richard?= Date: Sun, 3 May 2026 19:51:04 +0200 Subject: [PATCH 07/10] Add test with itnerval parameters --- test/newton1d.jl | 86 ------------------------------------------------ test/roots.jl | 11 +++++++ 2 files changed, 11 insertions(+), 86 deletions(-) delete mode 100644 test/newton1d.jl diff --git a/test/newton1d.jl b/test/newton1d.jl deleted file mode 100644 index 2c11c0b..0000000 --- a/test/newton1d.jl +++ /dev/null @@ -1,86 +0,0 @@ - -using IntervalArithmetic, IntervalRootFinding -using ForwardDiff -using Test - -const D = IntervalRootFinding.derivative - -setprecision(Interval, Float64) -float_pi = @interval(pi) - -setprecision(Interval, 10000) -big_pi = @interval(pi) -# Using precision "only" 256 leads to overestimation of the true roots for `cos` -# i.e the Newton method gives more accurate results! - -half_pi = big_pi / 2 -three_halves_pi = 3*big_pi/2 - -@testset "Testing newton1d" begin - - f(x) = exp(x^2) - cos(x) #double root - f′(x) = 2*x*exp(x^2) + sin(x) - f1(x) = x^4 - 10x^3 + 35x^2 - 50x + 24 #four unique roots - f1′(x) = 4x^3 - 30x^2 + 70x - 50 - f2(x) = 4567x^2 - 9134x + 4567 #double root - f2′(x) = 9134x - 9134 - f3(x) = (x^2 - 2)^2 #two double roots - f3′(x) = 4x * (x^2 - 2) - f4(x) = sin(x) - x #triple root at 0 - f4′(x) = cos(x) - 1 - f5(x) = (x^2 - 1)^4 * (x^2 - 2)^4 #two quadruple roots - f5′(x) = 8x * (-3 + 2x^2) * (2 - 3x^2 + x^4)^3 - for autodiff in (false, true) - if autodiff - rts1 = newton1d(sin, -5..5) - rts2 = newton1d(f, -∞..∞) - rts3 = newton1d(f1, -10..10) - rts4 = newton1d(f2, -10..11) - rts5 = newton1d(f3, -10..10) - rts6 = newton1d(f4, -10..10) - rts7 = newton1d(f5, -10..10) - - else - rts1 = newton1d(sin, cos, -5..5) - rts2 = newton1d(f, f′, -∞..∞) - rts3 = newton1d(f1, f1′, -10..10) - rts4 = newton1d(f2, f2′, -10..11) - rts5 = newton1d(f3, f3′, -10..10) - rts6 = newton1d(f4, f4′, -10..10) - rts7 = newton1d(f5, f5′, -10..10, reltol=0) - end - - @test length(rts1) == 3 - L = [-pi, 0, pi] - for i = 1:length(rts1) - @test L[i] in rts1[i].interval && :unique == rts1[i].status - end - - @test length(rts2) == 1 - @test (0..0) == rts2[1].interval && :unknown == rts2[1].status - - @test length(rts3) == 4 - L = [1, 2, 3, 4] - for i = 1:length(rts3) - @test L[i] in rts3[i].interval && :unique == rts3[i].status - end - - @test length(rts4) == 1 - @test 1 in rts4[1].interval && :unknown == rts4[1].status - - L1 = [-sqrt(2), sqrt(2)] - for i = 1:length(rts5) - @test L1[i] in rts5[i].interval && :unknown == rts5[i].status - end - - @test length(rts6) == 1 - @test 0 in rts6[1].interval && :unknown == rts6[1].status - - @test length(rts7) == 4 - L = [-sqrt(2), -1, 1, sqrt(2)] - for i = 1:length(rts7) - @test L[i] in rts7[i].interval && :unknown == rts7[i].status - end - - end -end diff --git a/test/roots.jl b/test/roots.jl index 787dc35..3cde8e1 100644 --- a/test/roots.jl +++ b/test/roots.jl @@ -513,3 +513,14 @@ function_list = [ end end end + + +@testset "Interval parameters" begin + a = interval(1.9, 2.1) + for contractor in newtonlike_methods + rts = roots(x -> x^2 - a^2, interval(0, 10) ; contractor) + + @test length(rts) == 1 + @test issubset_interval(a, rts[1].region) + end +end \ No newline at end of file From b41f8ef007d7e0da551bc8e6144196c68c385572 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Richard?= Date: Sun, 3 May 2026 19:56:59 +0200 Subject: [PATCH 08/10] Add double reciprocal from #131 --- test/roots.jl | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/test/roots.jl b/test/roots.jl index 3cde8e1..c56bfbc 100644 --- a/test/roots.jl +++ b/test/roots.jl @@ -76,6 +76,14 @@ end @test length(rts) == 1 @test root_status(only(rts)) == :unknown @test in_interval(0, root_region(only(rts))) + + # Double reciprocal from #131 + rr(x) = 1/(1/(1 + x)*(1 - x)) + rts = roots(rr, interval(-10, 10) ; contractor) + @test length(rts) == 2 + @test all(root_status(rt) == :unknown for rt in rts) + @test in_interval(-1, root_region(rts[1])) + @test in_interval(1, root_region(rts[2])) end end @@ -482,7 +490,7 @@ function_list = [ (W3, nothing, -10, 10, [1, 2, 3]), (W7, nothing, -10, 10, collect(1:7)), (x -> exp(x) - 2, exp, -20, 20, [log(big(2))]), - (x -> asin(sin(x)) - big"0.1", Returns(1.0), 0, 1, [big"0.1"]) + (x -> asin(sin(x)) - big"0.1", Returns(1.0), 0, 1, [big"0.1"]), ] @testset "Various precisions" begin From 0a454c77a8fa64126dc0f51c8e1112e211346e9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Richard?= Date: Sun, 3 May 2026 20:20:56 +0200 Subject: [PATCH 09/10] Greater clean up --- test/multidimensional.jl | 74 ++++++++++++++++++++++++ test/quadratic.jl | 21 ------- test/roots.jl | 116 +------------------------------------- test/runtests.jl | 39 +++++++++++-- test/stationary_points.jl | 4 -- test/test_smiley.jl | 30 +++------- 6 files changed, 121 insertions(+), 163 deletions(-) create mode 100644 test/multidimensional.jl delete mode 100644 test/quadratic.jl diff --git a/test/multidimensional.jl b/test/multidimensional.jl new file mode 100644 index 0000000..be0bf16 --- /dev/null +++ b/test/multidimensional.jl @@ -0,0 +1,74 @@ +@testset "2D roots" begin + f(x, y) = [x^2 + y^2 - 1, y - 2x] + f(X) = f(X...) + X = [interval(-6, 6), interval(-6, 6)] + + # Bisection + rts = roots(f, X ; contractor = Bisection, abstol = 1e-3) + exact_sol = [sqrt(1/5), 2sqrt(1/5)] + @test in_solution_set(exact_sol, root_region.(rts)) + @test in_solution_set(-exact_sol, root_region.(rts)) + + for contractor in newtonlike_methods + deriv = xx -> ForwardDiff.jacobian(f, xx) + test_newtonlike(f, deriv, X, contractor, 2) + rts = roots(f, X ; contractor) + @test in_solution_set(exact_sol, root_region.(rts)) + @test in_solution_set(-exact_sol, root_region.(rts)) + end + + # Infinite interval + X = [interval(-Inf, Inf), interval(-Inf, Inf)] + @suppress rts = roots(f, X ; contractor = Newton) + @test length(rts) == 2 +end + + +# From R docs: +# https://www.rdocumentation.org/packages/pracma/versions/1.9.9/topics/broyden + +@testset "3D roots" begin + function g(x) + (x1, x2, x3) = x + SVector( x1^2 + x2^2 + x3^2 - 1, + x1^2 + x3^2 - 0.25, + x1^2 + x2^2 - 4x3 + ) + end + dg = xx -> ForwardDiff.jacobian(g, xx) + + X = interval(-5, 5) + XX = [X, X, X] + + for contractor in newtonlike_methods + rts = roots(g, XX ; contractor) + @test length(rts) == 4 + @test all(isunique, rts) + @test all(rts .== roots(g, XX ; contractor, derivative = dg)) + end +end + +@testset "10D roots" begin + include("../examples/10_dimensional.jl") + + for contractor in newtonlike_methods + rts = roots(f10d, X10d_large ; contractor, max_iteration = 50_000) + @test any(X -> all(in_interval.(sol10d, X.region)), rts) + + rts = roots(f10d, X10d_close ; contractor, max_iteration = 1_000_000) + @test length(rts) == 1 + @test isunique(rts[1]) + end +end + +@testset "Dimension mismatch" begin + f21(xy) = [xy[1]^2 - 2] + f23(xy) = [xy[1]^2 - 2, xy[2]^2 - 3, xy[1] + xy[2]] + + X = [interval(0, 5), interval(0, 5)] + + for contractor in newtonlike_methods + @test_throws DimensionMismatch roots(f21, X ; contractor) + @test_throws DimensionMismatch roots(f23, X ; contractor) + end +end diff --git a/test/quadratic.jl b/test/quadratic.jl deleted file mode 100644 index 5e4dda0..0000000 --- a/test/quadratic.jl +++ /dev/null @@ -1,21 +0,0 @@ -using IntervalArithmetic, IntervalRootFinding - -struct Quad{T} - a::Interval{T} - b::Interval{T} - c::Interval{T} - x::Vector{Interval{T}} -end - - -@testset "Quadratic Interval Equations" begin - rts = Quad{Float64}[] - push!(rts, Quad(interval(1, 5), interval(2, 4), interval(0, 2), [interval(-4, -2), interval(-0, -0)])) - push!(rts, Quad(interval(-1, 1), interval(-2, 2), interval(-1, 1), [interval(-∞, -1), interval(-1, -1), interval(-1, ∞)])) - push!(rts, Quad(interval(1, 2), interval(1, 3), interval(2, 6), [interval(-2, -1)])) - push!(rts, Quad(interval(1, 1), interval(2, 2), interval(1, 1), [interval(-1, -1), interval(-1, -1)])) - - for i in 1:length(rts) - @test quadratic_roots(rts[i].a, rts[i].b, rts[i].c) == rts[i].x - end -end diff --git a/test/roots.jl b/test/roots.jl index c56bfbc..c0b9c93 100644 --- a/test/roots.jl +++ b/test/roots.jl @@ -1,39 +1,8 @@ -function all_unique(rts) - all(root_status.(rts) .== :unique) -end - -function roots_dist(rt1::Root{<:Interval}, rt2::Root{<:Interval}) - d = dist(root_region(rt1), root_region(rt2)) - return sum(d) -end - -function roots_dist(rt1::Root{<:AbstractVector}, rt2::Root{<:AbstractVector}) - return sum(dist.(root_region(rt1), root_region(rt2))) -end - -function roots_dist(rt1::Root{Complex{<:Interval}}, rt2::Root{Complex{<:Interval}}) - dreal = dist(real(root_region(rt1)), real(root_region(rt2))) - dimag = dist(imag(root_region(rt1)), imag(root_region(rt2))) - - return sum(dreal) + sum(dimag) -end - -function test_newtonlike(f, derivative, X, contractor, nsol, tol=1e-10) - rts = roots(f, X ; contractor) - @test length(rts) == nsol - @test all_unique(rts) - @test sum(roots_dist.(rts, roots(f, X ; contractor, derivative))) < tol -end - -function in_solution_set(point, solution_intervals) - return any(map(Y -> in_region(point, Y), solution_intervals)) -end - -@testset "1D roots" begin +@testset "Simple tests" begin # Default rts = roots(sin, interval(-5, 5)) @test length(rts) == 3 - @test all_unique(rts) + @test all(isunique, rts) # Bisection rts = roots(sin, interval(-5, 6) ; contractor = Bisection, abstol = 1e-3) @@ -50,10 +19,6 @@ end test_newtonlike(sin, cos, interval(-5, 5), contractor, 3) - # Infinite interval - rts = roots(x -> x^2 - 2, interval(-Inf, Inf) ; contractor) - @test length(rts) == 2 - # abs g(p) = abs(1 / ( (1+p)^30 ) * 10_000 - 100) g2(p) = (1 / ( (1+p)^30 ) * 10_000 - 100)^2 @@ -87,82 +52,6 @@ end end end - -@testset "2D roots" begin - f(x, y) = [x^2 + y^2 - 1, y - 2x] - f(X) = f(X...) - X = [interval(-6, 6), interval(-6, 6)] - - # Bisection - rts = roots(f, X ; contractor = Bisection, abstol = 1e-3) - exact_sol = [sqrt(1/5), 2sqrt(1/5)] - @test in_solution_set(exact_sol, root_region.(rts)) - @test in_solution_set(-exact_sol, root_region.(rts)) - - for contractor in newtonlike_methods - deriv = xx -> ForwardDiff.jacobian(f, xx) - test_newtonlike(f, deriv, X, contractor, 2) - rts = roots(f, X ; contractor) - @test in_solution_set(exact_sol, root_region.(rts)) - @test in_solution_set(-exact_sol, root_region.(rts)) - end - - # Infinite interval - X = [interval(-Inf, Inf), interval(-Inf, Inf)] - @suppress rts = roots(f, X ; contractor = Newton) - @test length(rts) == 2 -end - - -# From R docs: -# https://www.rdocumentation.org/packages/pracma/versions/1.9.9/topics/broyden - -@testset "3D roots" begin - function g(x) - (x1, x2, x3) = x - SVector( x1^2 + x2^2 + x3^2 - 1, - x1^2 + x3^2 - 0.25, - x1^2 + x2^2 - 4x3 - ) - end - dg = xx -> ForwardDiff.jacobian(g, xx) - - X = interval(-5, 5) - XX = [X, X, X] - - for contractor in newtonlike_methods - rts = roots(g, XX ; contractor) - @test length(rts) == 4 - @test all_unique(rts) - @test all(rts .== roots(g, XX ; contractor, derivative = dg)) - end -end - -@testset "10D roots" begin - include("../examples/10_dimensional.jl") - - for contractor in newtonlike_methods - rts = roots(f10d, X10d_large ; contractor, max_iteration = 50_000) - @test any(X -> all(in_interval.(sol10d, X.region)), rts) - - rts = roots(f10d, X10d_close ; contractor, max_iteration = 1_000_000) - @test length(rts) == 1 - @test isunique(rts[1]) - end -end - -@testset "Dimension mismatch" begin - f21(xy) = [xy[1]^2 - 2] - f23(xy) = [xy[1]^2 - 2, xy[2]^2 - 3, xy[1] + xy[2]] - - X = [interval(0, 5), interval(0, 5)] - - for contractor in newtonlike_methods - @test_throws DimensionMismatch roots(f21, X ; contractor) - @test_throws DimensionMismatch roots(f23, X ; contractor) - end -end - @testset "Multiple roots" begin f(x) = (x-1) * (x^2 - 2)^3 * (x^3 - 2)^4 g(x) = exp(x) - x - 1 # Double root at 0 from Burden & Faires, 9th ed, p.84 @@ -206,6 +95,7 @@ end end @testset "Infinite domain" begin + # Infinite interval f(xy) = [sin(1/xy[1]), xy[1]*xy[2]] # largest roots: [1/π, 0] and [infinity, 0] for contractor in newtonlike_methods diff --git a/test/runtests.jl b/test/runtests.jl index 3cb8f68..19d0e91 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,12 +8,43 @@ using Test newtonlike_methods = [Newton, Krawczyk] +# Compute the distance between set of roots +function roots_dist(rt1::Root{<:Interval}, rt2::Root{<:Interval}) + d = dist(root_region(rt1), root_region(rt2)) + return sum(d) +end + +function roots_dist(rt1::Root{<:AbstractVector}, rt2::Root{<:AbstractVector}) + return sum(dist.(root_region(rt1), root_region(rt2))) +end + +function roots_dist(rt1::Root{Complex{<:Interval}}, rt2::Root{Complex{<:Interval}}) + dreal = dist(real(root_region(rt1)), real(root_region(rt2))) + dimag = dist(imag(root_region(rt1)), imag(root_region(rt2))) + + return sum(dreal) + sum(dimag) +end + +# Test that a contractor that use derivatives +# Check that all return roots are unique, and that giving the derivative explictly +# makes no difference +function test_newtonlike(f, derivative, X, contractor, nsol, tol=1e-10) + rts = roots(f, X ; contractor) + @test length(rts) == nsol + @test all(isunique, rts) + @test sum(roots_dist.(rts, roots(f, X ; contractor, derivative))) < tol +end + +# Check that a given point is in one of the given intervals +function in_solution_set(point, solution_intervals) + return any(map(Y -> in_region(point, Y), solution_intervals)) +end + include("roots.jl") include("stationary_points.jl") + +include("multidimensional.jl") include("svectors.jl") + include("test_smiley.jl") include("linear_eq.jl") - -# include("newton1d.jl") -# include("quadratic.jl") -# include("slopes.jl") diff --git a/test/stationary_points.jl b/test/stationary_points.jl index ce74265..e46301c 100644 --- a/test/stationary_points.jl +++ b/test/stationary_points.jl @@ -1,7 +1,3 @@ -using Test - -import ForwardDiff - @testset "Stationary points" begin f(xx) = ( (x, y) = xx; sin(x) * sin(y) ) gradf = xx -> ForwardDiff.gradient(f, xx) diff --git a/test/test_smiley.jl b/test/test_smiley.jl index 749d08c..bd1cc45 100644 --- a/test/test_smiley.jl +++ b/test/test_smiley.jl @@ -1,33 +1,23 @@ include("../examples/smiley_examples.jl") -using Test -using IntervalArithmetic, IntervalRootFinding using .SmileyExample22, .SmileyExample52, .SmileyExample54, .SmileyExample55 -function test_all_unique(xs) - for x in xs - @test x.status == :unique - end - return nothing -end - const abstol = 1e-6 -for contractor in (Newton, Krawczyk) # NOTE: Bisection method performs badly in all examples - -@info("Testing method $contractor") @testset "$(SmileyExample22.title)" begin - roots_found = roots(SmileyExample22.f, SmileyExample22.region ; contractor, abstol) - @test length(roots_found) == 8 - test_all_unique(roots_found) - # no reference data for roots given + for contractor in (Newton, Krawczyk) + roots_found = roots(SmileyExample22.f, SmileyExample22.region ; contractor, abstol) + @test length(roots_found) == 8 + @test all(isunique, roots_found) + # no reference data for roots given + end end -for example in (SmileyExample52, SmileyExample54)#, SmileyExample55) - @testset "$(example.title)" begin +@testset "$(example.title)" for example in (SmileyExample52, SmileyExample54) + for contractor in (Newton, Krawczyk) roots_found = roots(example.f, example.region ; contractor, abstol) @test length(roots_found) == length(example.known_roots) - test_all_unique(roots_found) + @test all(isunique, roots_found) for rf in roots_found # check there is exactly one known root for each found root @test 1 == sum(example.known_roots) do rk @@ -36,5 +26,3 @@ for example in (SmileyExample52, SmileyExample54)#, SmileyExample55) end end end - -end # method loop From 33725782900c928c861a1d183c530b0700df32d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Richard?= Date: Mon, 11 May 2026 22:01:35 +0200 Subject: [PATCH 10/10] Test that examples are running --- test/runtests.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 19d0e91..f5e6282 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -48,3 +48,10 @@ include("svectors.jl") include("test_smiley.jl") include("linear_eq.jl") + +@testset "Example files" begin + for file in readdir("../examples/", join = true) + include(file) + @test true + end +end \ No newline at end of file