Merge pull request #34 from Joel-Dahne/v1-release

Release v1.0

.github/workflows/TagBot.yml

@@ -0,0 +1,16 @@
+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 }}
+          draft: true

Project.toml

@@ -1,14 +1,15 @@
 name = "ArbExtras"
 uuid = "c87f5d39-a852-4149-8a88-e3a13a25afc6"
 authors = ["Joel Dahne <[email protected]>"]
-version = "0.1.10"
+version = "1.0.0"
 
 [deps]
 Arblib = "fb37089c-8514-4489-9461-98f9c8763369"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 Arblib = "1"
+PropCheck = "0.10"
 Random = "1.6"
 SpecialFunctions = "1, 2"
 Test = "1.6"

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!`

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)

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)

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

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,

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))

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()
@@ -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
 
@@ -481,6 +482,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 +507,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)

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)