From 6b046f0a2df05f54d02e04f8cabed5482e1aaeeb Mon Sep 17 00:00:00 2001 From: Joel Dahne Date: Fri, 19 Jan 2024 19:08:57 +0100 Subject: [PATCH 01/12] Fix typo in documentation --- src/integrate.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/integrate.jl b/src/integrate.jl index 5994ea1..88c8d70 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -86,7 +86,7 @@ end Compute an enclosure of the integral of `f on the interval `[a, b]`. -It uses a Guass-Legendre quadrature of order 2 through +It uses a Gauss-Legendre quadrature of order 2 through [`integrate_gauss_legendre`](@ref) together with bisection of the interval. On each subinterval the integral is computed using [`integrate_gauss_legendre`](@ref) and it checks if the result From a16368f83b3610165662db0fb0e33523807c0461 Mon Sep 17 00:00:00 2001 From: Joel Dahne Date: Fri, 19 Jan 2024 19:10:18 +0100 Subject: [PATCH 02/12] derivative_function: Define for polynomials and use in refine_root --- src/refine_root.jl | 2 +- src/utilities.jl | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/refine_root.jl b/src/refine_root.jl index ffd9622..a65f6d8 100644 --- a/src/refine_root.jl +++ b/src/refine_root.jl @@ -57,7 +57,7 @@ root. function refine_root( f, root::Arb; - df = f isa ArbPoly ? Arblib.derivative(f) : x -> f(ArbSeries((x, 1)))[1], + df = derivative_function(f), atol = 0, rtol = 4eps(one(root)), min_iterations = 1, diff --git a/src/utilities.jl b/src/utilities.jl index 39e20e7..102635b 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -481,6 +481,8 @@ The returned function accepts only `Arb`, `Acb`, `ArbSeries` and `AcbSeries` as input. For `Arb` and `ArbSeries` it calls `f` with `ArbSeries`. For `Acb` and `AcbSeries` it calls `f` with `AcbSeries`. +If `f` is a polynomial, then return the derivative of the polynomial. + **IMPROVE:** Avoid overflow in factorial function for large `n`. """ function derivative_function(f, n::Integer = 1) @@ -504,3 +506,5 @@ function derivative_function(f, n::Integer = 1) return res end end + +derivative_function(p::Union{ArbPoly,AcbPoly}, n::Integer = 1) = Arblib.derivative(p, n) From ff3d84288715690d85a169a6b6ed56f3e957a84f Mon Sep 17 00:00:00 2001 From: Joel Dahne Date: Fri, 19 Jan 2024 19:19:25 +0100 Subject: [PATCH 03/12] compose_zero!: Add documentation --- src/temporary.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/temporary.jl b/src/temporary.jl index e486122..1beb5e8 100644 --- a/src/temporary.jl +++ b/src/temporary.jl @@ -9,8 +9,12 @@ Return `true` if `p` is of the form `c + x`. iscpx(p::Union{ArbPoly,AcbPoly,ArbSeries,AcbSeries}) = length(Arblib.cstruct(p)) == 2 && isone(Arblib.ref(p, 1)) -function compose_zero!(res::T, p::T, q::T) where {T<:Union{ArbPoly,AcbPoly}} +""" + compose_zero!(res::T, p::T, q::T) where {T<:Union{ArbPoly,AcbPoly,ArbSeries,AcbSeries}} +In-place version of [`compose_zero`](@ref). +""" +function compose_zero!(res::T, p::T, q::T) where {T<:Union{ArbPoly,AcbPoly}} iscpx(q) && return Arblib.set!(res, p) if !iszero(q) && !iszero(Arblib.ref(q, 0)) From 161f747521d0d569bbc0be1acc752ed7832a2aee Mon Sep 17 00:00:00 2001 From: Joel Dahne Date: Fri, 19 Jan 2024 19:19:40 +0100 Subject: [PATCH 04/12] Update README --- README.md | 125 +++++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 106 insertions(+), 19 deletions(-) diff --git a/README.md b/README.md index 0599cdb..4d07617 100644 --- a/README.md +++ b/README.md @@ -1,21 +1,108 @@ # ArbExtras.jl - This package extends [Arblib](https://github.com/kalmarek/Arblib.jl) -with some extra functionality. It's mainly intended for personal use -in my research and should not be seen as a stable Julia package. It -will likely never be added to the Julia repository. It can be seen as -a follow up of [ArbTools](https://github.com/Joel-Dahne/ArbTools.jl) -but based on Arblib instead of -[Nemo](https://github.com/wbhart/Nemo.jl) for the Arb interface. - -## Installation -The package is not in the general Julia repository but can be -installed through the package manager with -``` julia -pkg> add https://github.com/Joel-Dahne/ArbExtras.jl -``` - -To see if it (some) things work correctly you can run the tests with -``` julia -pkg> test ArbExtras -``` +with some methods for enclosing roots, enclosing extrema and computing +integrals. + +The package started development during my PhD with the goal to give +reusable implementations of algorithms that were used in my research. +The type of methods that are implemented and which type of input they +are optimized for has been heavily influenced by the projects that I +worked on during that time. + +## Functionality +The main functionality provided by the package are methods for +isolating roots and enclosing extrema of univariate real functions. In +addition to this is provides a rudimentary integral routine and a few +minor things. + +### Enclosing roots +The package provides three different functions for isolating and +enclosing roots of univariate, real valued functions: + +- `isolate_roots(f, a::Arf, b::Arf)` is used for isolating all roots + of a function `f` on the interval `[a, b]`. +- `refine_root(f, root::Arb)` is used for refining an initial + enclosure `root` of a zero to `f`. It makes use of the interval + Newton method. +- `refine_root_bisection(f, a, b)` is used for refining an initial + enclosure `[a, b]` of a root of `f`. It requires that `f(a)` and + `f(b)` have different signs. It makes use if bisection, checking the + sign of the midpoint in each iteration. + +See the documentation of the respective methods for more details, in +particular the available keyword arguments. The first two methods +require access to the derivative of `f`, this is computed using +`ArbSeries` and for that reason `f` needs to support evaluation of +both `Arb` and `ArbSeries`. The `refine_root_bisection` doesn't +require access to the derivative and can therefore be used when this +is not available. + +### Enclosing extrema +The package provides three families of methods: + +- `extrema_polynomial`, `minimum_polynomial` and `maximum_polynomial` +- `extrema_series`, `minimum_series` and `maximum_series` +- `extrema_enclosure`, `minimum_enclosure` and `maximum_enclosure` + +As indicated by the names they compute either the minimum, the maximum +or both. + +- `extrema_polynomial(p::ArbPoly, a::Arf, b::Arb)` is used for + enclosing the extrema of a polynomial `p` on the interval `[a, b]`. + This is done by enclosing the roots of the derivative of `p`. +- `extrema_series(f, a::Arf, b::Arb; degree::Integer = 8)` is used for + enclosing the extrema of a function `f` on the interval `[a, b]`. + This is done by computing a Taylor series of the given degree. The + extrema of the series is computed using `extrema_polynomial` and + then an enclosure of the remainder term is added. +- `extrema_enclosure(f, a::Arf, b::Arb)` is also + used for enclosing the extrema of a function `f` on the interval + `[a, b]`. In this case it is done by iteratively bisecting the + interval and using `extrema_series` on each subinterval. The + bisection is continued until the result satisfies the required + tolerance. + +See the documentation of the respective methods for more details, in +particular the available keyword arguments. The `minimum_` and +`maximum_` versions have the same interface, the only difference being +that they return either the minimum or maximum instead of both. + +For `enclosure_series` it is also possible to call it like +`enclosure_series(f, x::Arb)`, in which case it computes an enclosure +of `f(x)` by computing the extrema and taking the union of them. + +The package also provides `bounded_by(f, a::Arf, b::Arf, C::Arf)` +which checks if the function `f` is bounded by `C` on the interval +`[a, b]`. It is similar to `maximum_enclosure` but only bisects enough +to be able to check that it is bounded by `C`. See also the +`ubound_tol` and `lbound_tol` arguments to `extrema_enclosure`. + +### Other functionality +In addition to the above the package also provides some smaller +things: + +- `integrate(f, a::Arb, b::Arb)` can be used to compute the integral + of `f` from `a` to `b` using a Gauss-Legendre quadrature of order 2, + combined with bisection. In general the integrator from Arb, + `Arblib.integrate` performs much better and should be preferred. + This method does however have the benefit that it does not require + evaluation on complex balls, instead it needs access to the fourth + derivative of `f`, which is computed using `ArbSeries`. So it can be + of use if there is no implementation of `f(::Acb)`, but there is one + for `f(::ArbSeries)`. +- `SpecialFunctions.besselj(ν::Arb, z::ArbSeries)` for computing + Taylor expansions of the Bessel function. +- Utility functions which are mainly for internal use, but could + be useful in other cases as well , see their documentation for + details. + - `bisect_interval`, `bisect_interval_recursive` and `bisect_intervals` + - `check_tolerance` + - `check_interval` + - `format_interval` + - `taylor_remainder` + - `enclosure_ubound`, `enclosure_lbound` and `enclosure_getinterval` + - `derivative_function` +- Some "temporary" functions that should possibly be added to + Arblib.jl or even the Arb library itself in the future + - `iscpx` + - `compose_zero` and `compose_zero!` From 199cb449f36d74d8f6dd277ca057e9d9963a0cce Mon Sep 17 00:00:00 2001 From: Joel Dahne Date: Fri, 19 Jan 2024 20:22:16 +0100 Subject: [PATCH 05/12] Add compat for PropCheck --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 953d6bf..5d057f1 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" [compat] Arblib = "1" +PropCheck = "0.10" Random = "1.6" SpecialFunctions = "1, 2" Test = "1.6" From 2cca2ff1c55dc2cf5d0655f5f76db8356a8768ec Mon Sep 17 00:00:00 2001 From: Joel Dahne Date: Sat, 20 Jan 2024 09:07:49 +0100 Subject: [PATCH 06/12] Add TagBot workflow --- .github/workflows/TagBot.yml | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 .github/workflows/TagBot.yml diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml new file mode 100644 index 0000000..f49313b --- /dev/null +++ b/.github/workflows/TagBot.yml @@ -0,0 +1,15 @@ +name: TagBot +on: + issue_comment: + types: + - created + workflow_dispatch: +jobs: + TagBot: + if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' + runs-on: ubuntu-latest + steps: + - uses: JuliaRegistries/TagBot@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} + ssh: ${{ secrets.DOCUMENTER_KEY }} From 90da9247400c231c985ade28d11e3a4c32872b1e Mon Sep 17 00:00:00 2001 From: Joel Dahne Date: Sat, 20 Jan 2024 09:34:28 +0100 Subject: [PATCH 07/12] check_tolerance: Add tests with ArbVector --- test/utilities.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/utilities.jl b/test/utilities.jl index 2084bda..4c697dc 100644 --- a/test/utilities.jl +++ b/test/utilities.jl @@ -102,6 +102,7 @@ @testset "x::AbstractVector" begin prop_check_tolerance_elementwise((x, atol, rtol)) = + ArbExtras.check_tolerance(ArbVector(x); atol, rtol) == ArbExtras.check_tolerance(x; atol, rtol) == all(xᵢ -> ArbExtras.check_tolerance(xᵢ; atol, rtol), x) From 9c258820b6940a33580c97da163e84ec54106ca0 Mon Sep 17 00:00:00 2001 From: Joel Dahne Date: Tue, 23 Jan 2024 19:38:14 +0100 Subject: [PATCH 08/12] TagBot: Generate a draft release --- .github/workflows/TagBot.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml index f49313b..0b0cab5 100644 --- a/.github/workflows/TagBot.yml +++ b/.github/workflows/TagBot.yml @@ -13,3 +13,4 @@ jobs: with: token: ${{ secrets.GITHUB_TOKEN }} ssh: ${{ secrets.DOCUMENTER_KEY }} + draft: true From e47b5fc1e8506b94536b7585192ba2fff1172ca8 Mon Sep 17 00:00:00 2001 From: Joel Dahne Date: Sun, 28 Jan 2024 13:38:06 +0100 Subject: [PATCH 09/12] check_tolerance: Use Arb for error also for vectors --- src/utilities.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/utilities.jl b/src/utilities.jl index 102635b..11af4a3 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -258,7 +258,7 @@ function check_tolerance( end isnothing(atol) && isnothing(rtol) && return true - error = Mag() + error = Arb() for i in eachindex(x) xᵢ = if x isa ArbVector @@ -269,7 +269,8 @@ function check_tolerance( Arblib.isexact(xᵢ) && continue - Arblib.mul_2exp!(error, Arblib.radref(xᵢ), 1) + Arblib.set!(error, Arblib.radref(xᵢ)) + Arblib.mul_2exp!(error, error, 1) !isnothing(atol) && !iszero(atol) && error <= atol && continue From 399f3e9d94366d9cf33645e5b24ece572849d4eb Mon Sep 17 00:00:00 2001 From: Joel Dahne Date: Sun, 28 Jan 2024 15:46:03 +0100 Subject: [PATCH 10/12] bisect_interval_recursive: Fix typo in throw method --- src/utilities.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utilities.jl b/src/utilities.jl index 11af4a3..95ec9f2 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -124,7 +124,7 @@ function bisect_interval_recursive( depth::Integer; log_midpoint::Bool = false, ) where {T<:Union{Arf,Arb}} - depth >= 0 || Throw(ArgumentError("depth needs to be non-negative, got $depth")) + depth >= 0 || throw(ArgumentError("depth needs to be non-negative, got $depth")) if T == Arf && log_midpoint buffer = Arb() From 1a58e369aee911e4ae9099a4e3d225534acff4ff Mon Sep 17 00:00:00 2001 From: Joel Dahne Date: Sun, 28 Jan 2024 15:50:21 +0100 Subject: [PATCH 11/12] Improve type stability issues from using Threads.@threads See https://discourse.julialang.org/t/type-instability-because-of-threads-boxing-variables/78395 In short, using Threads.@threads introduces an implicit closure and this causes issues with type stability if a binding is assigned to multiple times. Introducing an explicit let scope fixes these type stability issues. --- src/bounded_by.jl | 12 ++++++--- src/extrema_enclosure.jl | 54 ++++++++++++++++++++++++---------------- 2 files changed, 41 insertions(+), 25 deletions(-) diff --git a/src/bounded_by.jl b/src/bounded_by.jl index 752e09f..17ab593 100644 --- a/src/bounded_by.jl +++ b/src/bounded_by.jl @@ -90,8 +90,10 @@ function bounded_by( values = similar(intervals, Arb) if degree < 0 if threaded - Threads.@threads for i in eachindex(intervals) - values[i] = maybe_abs(f(Arb(intervals[i]))) + let intervals = intervals + Threads.@threads for i in eachindex(intervals) + values[i] = maybe_abs(f(Arb(intervals[i]))) + end end else for i in eachindex(intervals) @@ -100,8 +102,10 @@ function bounded_by( end else if threaded - Threads.@threads for i in eachindex(intervals) - values[i], _ = maximum_series(f, intervals[i]...; degree, abs_value) + let intervals = intervals + Threads.@threads for i in eachindex(intervals) + values[i], _ = maximum_series(f, intervals[i]...; degree, abs_value) + end end else for (i, (a, b)) in enumerate(intervals) diff --git a/src/extrema_enclosure.jl b/src/extrema_enclosure.jl index daaf12c..853a0b7 100644 --- a/src/extrema_enclosure.jl +++ b/src/extrema_enclosure.jl @@ -178,10 +178,12 @@ function extrema_enclosure( values_max = similar(intervals, Arb) if degree < 0 if threaded - Threads.@threads for i in eachindex(intervals) - v = f(Arb(intervals[i])) - abs_value && Arblib.abs!(v, v) - values_min[i] = values_max[i] = v + let intervals = intervals + Threads.@threads for i in eachindex(intervals) + v = f(Arb(intervals[i])) + abs_value && Arblib.abs!(v, v) + values_min[i] = values_max[i] = v + end end else for i in eachindex(intervals) @@ -193,9 +195,11 @@ function extrema_enclosure( else point_values = similar(values_min) if threaded - Threads.@threads for i in eachindex(intervals) - values_min[i], values_max[i], point_values[i] = - extrema_series(f, intervals[i]...; degree, abs_value) + let intervals = intervals + Threads.@threads for i in eachindex(intervals) + values_min[i], values_max[i], point_values[i] = + extrema_series(f, a, b; degree, abs_value) + end end else for (i, (a, b)) in enumerate(intervals) @@ -369,10 +373,12 @@ function minimum_enclosure( values = similar(intervals, Arb) if degree < 0 if threaded - Threads.@threads for i in eachindex(intervals) - v = f(Arb(intervals[i])) - abs_value && Arblib.abs!(v, v) - values[i] = v + let intervals = intervals + Threads.@threads for i in eachindex(intervals) + v = f(Arb(intervals[i])) + abs_value && Arblib.abs!(v, v) + values[i] = v + end end else for i in eachindex(intervals) @@ -384,9 +390,11 @@ function minimum_enclosure( else point_values = similar(values) if threaded - Threads.@threads for i in eachindex(intervals) - values[i], point_values[i] = - minimum_series(f, intervals[i]...; degree, abs_value) + let intervals = intervals + Threads.@threads for i in eachindex(intervals) + values[i], point_values[i] = + minimum_series(f, intervals[i]...; degree, abs_value) + end end else for (i, (a, b)) in enumerate(intervals) @@ -524,10 +532,12 @@ function maximum_enclosure( values = similar(intervals, Arb) if degree < 0 if threaded - Threads.@threads for i in eachindex(intervals) - v = f(Arb(intervals[i])) - abs_value && Arblib.abs!(v, v) - values[i] = v + let intervals = intervals + Threads.@threads for i in eachindex(intervals) + v = f(Arb(intervals[i])) + abs_value && Arblib.abs!(v, v) + values[i] = v + end end else for i in eachindex(intervals) @@ -539,9 +549,11 @@ function maximum_enclosure( else point_values = similar(values) if threaded - Threads.@threads for i in eachindex(intervals) - values[i], point_values[i] = - maximum_series(f, intervals[i]...; degree, abs_value) + let intervals = intervals + Threads.@threads for i in eachindex(intervals) + values[i], point_values[i] = + maximum_series(f, intervals[i]...; degree, abs_value) + end end else for (i, (a, b)) in enumerate(intervals) From 4cfbeee4b9f1acbf3207ec749366c31922045446 Mon Sep 17 00:00:00 2001 From: Joel Dahne Date: Sat, 20 Jan 2024 09:08:07 +0100 Subject: [PATCH 12/12] Set version to 1.0 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 5d057f1..ebf2700 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ArbExtras" uuid = "c87f5d39-a852-4149-8a88-e3a13a25afc6" authors = ["Joel Dahne "] -version = "0.1.10" +version = "1.0.0" [deps] Arblib = "fb37089c-8514-4489-9461-98f9c8763369"