diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml new file mode 100644 index 00000000..4bf42c6b --- /dev/null +++ b/.github/workflows/documentation.yml @@ -0,0 +1,23 @@ +name: Documentation +on: + push: + branches: + - master + tags: '*' + pull_request: + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - uses: julia-actions/setup-julia@latest + with: + version: '1.6' + - name: Install dependencies + run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' + - name: Build and deploy + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # If authenticating with SSH deploy key + run: julia --project=docs/ docs/make.jl diff --git a/.gitignore b/.gitignore index ba2144b3..c4fbd212 100644 --- a/.gitignore +++ b/.gitignore @@ -2,4 +2,5 @@ deps/usr deps/build_* deps/build.log deps/deps.jl +docs/build Manifest.toml diff --git a/Project.toml b/Project.toml index e52407e0..dd71e835 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,9 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" + +[extras] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] @@ -17,3 +20,6 @@ Arb_jll = "~200.2300" FLINT_jll = "~200.900.000" SpecialFunctions = "1.0, 2" julia = "1.6" + +[targets] +test = ["Documenter", "Test"] diff --git a/README.md b/README.md index 87c7bf56..a6009bfa 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,8 @@ # Arblib.jl +[![][docs-stable-img]][docs-stable-url] +[![][docs-dev-img]][docs-dev-url] + This package is a thin, efficient wrapper around [Arb](http://arblib.org) - a C library for arbitrary-precision ball arithmetic. The package is currently in early development. More features and @@ -254,3 +257,8 @@ Enabling a threaded version of flint can be done by setting the environment variable `NEMO_THREADED=1`. Note that this should be set before `Arblib.jl` is loaded. To set the actual number of threads, use `Arblib.flint_set_num_threads($numberofthreads)`. + +[docs-dev-img]: https://img.shields.io/badge/docs-dev-blue.svg +[docs-dev-url]: https://Kalmarek.github.io/Arblib.jl/dev/ +[docs-stable-img]: https://img.shields.io/badge/docs-stable-blue.svg +[docs-stable-url]: https://Kalmarek.github.io/Arblib.jl/stable diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 00000000..208ecd62 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,7 @@ +[deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" + +[compat] +BenchmarkTools = "1" +Documenter = "1" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 00000000..2ffdd466 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,27 @@ +using Documenter, Arblib + +DocMeta.setdocmeta!(Arblib, :DocTestSetup, :(using Arblib); recursive = true) + +makedocs( + sitename = "Arblib.jl", + modules = [Arblib], + pages = [ + "index.md", + "Low level wrapper" => [ + "Types" => "wrapper-types.md", + "Methods" => "wrapper-methods.md", + "Floating point wrapper" => "wrapper-fpwrap.md", + ], + "High level interface" => [ + "Types" => "interface-types.md", + "Ball methods" => "interface-ball.md", + "Integration" => "interface-integration.md", + "Series" => "interface-series.md", + "Mutable arithmetic" => "interface-mutable.md", + ], + "Rigorous numerics" => "rigorous.md", + ], + warnonly = [:missing_docs], +) + +deploydocs(repo = "github.com/kalmarek/Arblib.jl") diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 00000000..a496b4fd --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,20 @@ +# Arblib.jl Documentation + +This package is a wrapper around [Arb](http://arblib.org) - a C +library for arbitrary-precision ball arithmetic. Other wrappers of Arb +for Julia include [Nemo](https://github.com/Nemocas/Nemo.jl) and +[ArbNumerics.jl](https://github.com/JeffreySarnoff/ArbNumerics.jl). + +The **goal** of Arblib.jl is to supply a **low lever wrapper** of the +methods in Arb as well as a **high level interface**. The low level +wrapper should allow for writing methods using mutability and with +performance very close to that of those written in C. The high level +interface should make it easy to use in generic Julia code, similarly +to how `BigFloat` is a wrapper around the MPFR library. In addition it +should be possible to seamlessly switch between the high level +interface and the low level wrapper when needed. + +The above goals can be put into contrast with Nemo, whose high level +interface is made for use in the +[AbstractAlgebra.jl](https://github.com/Nemocas/AbstractAlgebra.jl) +universe and not general Julia code. diff --git a/docs/src/interface-ball.md b/docs/src/interface-ball.md new file mode 100644 index 00000000..51b0d3d4 --- /dev/null +++ b/docs/src/interface-ball.md @@ -0,0 +1,37 @@ +# Ball methods +The following methods are useful for explicitly dealing with the ball +representation of `Arb` and related values. + +## Construction +For constructing balls the methods below are useful. Note that there +is no `setinterval` method, this is instead accomplished with `Arb((a, +b))` for constructing a ball containing the interval ``[a, b]``. + +``` @docs +setball +add_error +``` + +## Destruction +For extracting information about the ball representation the following +methods are useful. + +``` @docs +radius +midpoint +lbound +ubound +abs_lbound +abs_ubound +getinterval +getball +``` + +## Union and intersection +The `Base.union` and `Base.intersect` methods are overloaded to +compute the union and intersection of balls. + +``` @docs +union(::Arb, ::Arb) +intersect(::Arb, ::Arb) +``` diff --git a/docs/src/interface-integration.md b/docs/src/interface-integration.md new file mode 100644 index 00000000..37bf57e2 --- /dev/null +++ b/docs/src/interface-integration.md @@ -0,0 +1,4 @@ +``` @docs +Arblib.integrate +Arblib.integrate! +``` diff --git a/docs/src/interface-mutable.md b/docs/src/interface-mutable.md new file mode 100644 index 00000000..e6b756eb --- /dev/null +++ b/docs/src/interface-mutable.md @@ -0,0 +1,85 @@ +# Mutable arithmetic +The high level interface can be combined with the low level wrapper to +allow for efficient computations using mutable arithmetic. + +In the future it would be nice to have an interface to +[MutableArithmetics.jl](https://github.com/jump-dev/MutableArithmetics.jl), +see [#118](https://github.com/kalmarek/Arblib.jl/issues/118). + +The following methods are useful for mutating part of a value + +``` @docs +Arblib.radref +Arblib.midref +Arblib.realref +Arblib.imagref +Arblib.ref +``` + +## Examples + +Compare computing ``\sqrt{x^2 + y^2}`` using mutable arithmetic with +the default. + +``` @repl +using Arblib, BenchmarkTools + +x = Arb(1 // 3) +y = Arb(1 // 5) +res = zero(x) + +f(x, y) = sqrt(x^2 + y^2) +f!(res, x, y) = begin + Arblib.sqr!(res, x) + Arblib.fma!(res, res, y, y) + return Arblib.sqrt!(res, res) +end + +@benchmark f($x, $y) samples=10000 evals=500 +@benchmark f!($res, $x, $y) samples=10000 evals=500 +``` + +Set the radius of the real part of an `Acb`. + +``` @repl +using Arblib + +z = Acb(1, 2) +Arblib.set!(Arblib.radref(Arblib.realref(z)), 1e-10) +z +``` + +Compare a naive implementation of polynomial evaluation using +mutable arithmetic with one not using using it. + +``` @repl +using Arblib, BenchmarkTools + +p = ArbPoly(1:10) +x = Arb(1 // 3) +res = zero(x) + +function eval(p, x) + res = zero(x) + xi = one(x) + for i in 0:Arblib.degree(p) + res += p[i] * xi + xi *= x + end + return res +end + +function eval!(res, p, x) + Arblib.zero!(res) + xi = one(x) + for i in 0:Arblib.degree(p) + Arblib.addmul!(res, Arblib.ref(p, i), xi) + Arblib.mul!(xi, xi, x) + end + return res +end + +@benchmark eval($p, $x) samples = 10000 evals = 30 +@benchmark eval!($res, $p, $x) samples = 10000 evals = 30 +@benchmark $p($x) samples = 10000 evals = 30 # Arb implementation for reference +``` diff --git a/docs/src/interface-series.md b/docs/src/interface-series.md new file mode 100644 index 00000000..f43d96e2 --- /dev/null +++ b/docs/src/interface-series.md @@ -0,0 +1,25 @@ +# Series +Taylor series arithmetic allows for the computation of truncated +Taylor series of functions and is a form of higher order automatic +differentiation. See e.g. +[`TaylorSeries.jl`](https://github.com/JuliaDiff/TaylorSeries.jl) and +[`TaylorDiff.jl`](https://github.com/JuliaDiff/TaylorDiff.jl) for +implementations of Taylor series in Julia. + +The Arb library has good support for computing with polynomials as +Taylor expansions. The types [`ArbSeries`](@ref) and +[`AcbSeries`](@ref) are intended to make this easy to use from Julia. +They are given by an [`ArbPoly`](@ref)/[`AcbPoly`](@ref) together with +the length of the expansion. + +## Example +```@repl 1 +using Arblib + +x0 = Arb(1 // 3, prec = 64) +x = ArbSeries((x0, 1), degree = 5) + +sin(x) + +sin(x)^2 + cos(x)^2 +``` diff --git a/docs/src/interface-types.md b/docs/src/interface-types.md new file mode 100644 index 00000000..1f311c4c --- /dev/null +++ b/docs/src/interface-types.md @@ -0,0 +1,61 @@ +# Types + +The package defines a number of types for the high level interface. + +## Basic +These types directly map to corresponding Arb types. + +``` @docs +Mag +Arf +Arb +Acb +ArbVector +AcbVector +ArbPoly +AcbPoly +ArbMatrix +AcbMatrix +``` + +## Series +The package defines two series types, which are wrapper for the +polynomial types with a specified degree. + +``` @docs +ArbSeries +AcbSeries +``` + +## Ref +In addition to these there are a number of `Ref` types, which allow +for non-allocating access in a number of cases. + +``` @docs +MagRef +ArfRef +ArbRef +AcbRef +ArbRefVector +AcbRefVector +ArbRefMatrix +AcbRefMatrix +``` + +## Correspondence between types +We have the following table for the correspondence with between the +[Low level wrapper types](wrapper-types.md) and the high level +interface types. + +| Arb | Wrapper | High level | Ref | +|----------|------------------|-------------|----------------| +| `mag_t` | `mag_struct` | `Mag` | `MagRef` | +| `arf_t` | `arf_struct` | `Arf` | `ArfRef` | +| `arb_t` | `arb_struct` | `Arb` | `ArbRef` | +| `acb_t` | `acb_struct` | `Acb` | `AcbRef` | +| `arb_t*` | `arb_vec_struct` | `ArbVector` | `ArbRefVector` | +| `acb_t*` | `acb_vec_struct` | `AcbVector` | `AcbRefVector` | +| `arb_poly_t` | `arb_poly_struct` | `ArbPoly` or `ArbSeries` | | +| `acb_poly_t` | `acb_poly_struct` | `AcbPoly` or `AcbSeries` | | +| `arb_mat_t` | `arb_mat_struct` | `ArbMatrix` | `ArbRefMatrix` | +| `acb_mat_t` | `acb_mat_struct` | `AcbMatrix` | `AcbRefMatrix` | diff --git a/docs/src/rigorous.md b/docs/src/rigorous.md new file mode 100644 index 00000000..2ce6c042 --- /dev/null +++ b/docs/src/rigorous.md @@ -0,0 +1,48 @@ +# Rigorous numerics +Arb is made for rigorous numerics and any functions which do not +produce rigorous results are clearly marked as such. This is not the +case with Julia in general and you therefore have to be careful when +interacting with the ecosystem if you want your results to be +completely rigorous. Below we discuss some things to be extra careful +with. + +## Implicit promotion +Julia automatically promotes types in many cases and in particular you +have to watch out for temporary non-rigorous values. For example +`2(π * Arb(1 // 3))` is okay, but not `2π * Arb(1 // 3)` + +``` repl +x = 2(π * Arb(1 // 3)) +y = 2π * Arb(1 // 3) +Arblib.overlaps(x, y) +``` + +## Non-rigorous algorithms +Standard numerical algorithms typically return (hopefully good) +approximations. These algorithms can then not directly be used in +rigorous numerical computations unless the error can be bounded. + +For example Julias built in methods for solving linear systems doesn't +produce rigorous results. Instead you would have to use the solves +provided by Arb, such as `Arblib.solve!`. + +Other examples would include integration and solving of differential +equations. + +## Implementation details +In some cases the implementation in Julia implicitly makes certain +assumptions to improve performance and this can lead to issues. + +For example, prior to Julia version 1.8 the `minimum` and `maximum` +methods in Julia checked for `NaN` results (on which is short fuses) +using `x == x`, which works for most numerical types but not for `Arb` +(`x == x` is only true if the radius is zero). See + and in particular + for more details. +Since Julia version 1.8 the `minimum` and `maximum` methods work +correctly for `Arb`, for earlier versions of Julia it only works +correctly in some cases. + +These types of problems are the hardest to find since they are not +clear from the documentation but you have to read the implementation, +`@which` and `@less` are your friends in these cases. diff --git a/docs/src/wrapper-fpwrap.md b/docs/src/wrapper-fpwrap.md new file mode 100644 index 00000000..76b3ae08 --- /dev/null +++ b/docs/src/wrapper-fpwrap.md @@ -0,0 +1,33 @@ +# Floating point wrapper + +The functions in the +[`arb_fpwrap`](https://www.arblib.org/arb_fpwrap.html) module of Arb +are handled different from other functions. The methods for them are +generated using + +``` @docs +Arblib.ArbCall.@arbfpwrapcall_str +Arblib.ArbCall.ArbFPWrapFunction +Arblib.ArbCall.jlcode(af::Arblib.ArbCall.ArbFPWrapFunction) +``` + +The name for the generated method is given by removing the `arb` +prefix and the `double` or `cdouble` in the middle of the name. + +The `flag` argument that the C functions take are split into several +keyword arguments in Julia. For the `double` type this is +`correct_rounding::Bool = false` and `work_limit::Integer = 8`. For +the `cdouble` type it also includes `accurate_parts::Bool = false`. +The default values correspond to setting the flag to 0 in C. + +The C functions return an `int` flag indicating the status. If the +return flag is `FPWRAP_SUCCESS` the computed value is returned. If the +return flag is `FPWRAP_UNABLE` it throws an error if the keyword +argument `error_on_failure` is true and returns `NaN` otherwise. The +default value for `error_on_failure` is handled by the following two +methods + +``` @docs +Arblib.ArbCall.fpwrap_error_on_failure_default +Arblib.ArbCall.fpwrap_error_on_failure_default_set +``` diff --git a/docs/src/wrapper-methods.md b/docs/src/wrapper-methods.md new file mode 100644 index 00000000..a51b68af --- /dev/null +++ b/docs/src/wrapper-methods.md @@ -0,0 +1,118 @@ +# Methods + +The methods for the low level wrapper are automatically generated by +parsing the Arb documentation. This is handled by the `Arblib.ArbCall` +module. + +## Parsing +The parsing is handled by + +``` @docs +Arblib.ArbCall.parse_and_generate_arbdoc +``` + +Note that the parsing is done ahead of time and the generated files in +`src/arbcalls/` are added to git. As a user of the package you +therefore typically don't need to care about this step. + +## Generated methods +The automatic generation of the methods is handled by + +``` @docs +Arblib.ArbCall.@arbcall_str +Arblib.ArbCall.ArbFunction +Arblib.ArbCall.jlcode(af::Arblib.ArbCall.ArbFunction) +``` + +The main things to understand is how the name of the generated +function is determined, how the arguments are handled and the return +value. + +### Naming +The name of the Arb function is "Juliafied" using the following +guidelines: + +1. Prefixes and suffixes a function name which only refer to the type + of input are removed since Julia has multiple dispatch to deal with + this problem. +2. Functions which modify the first argument get an `!` appended to + the name. + +The implementation is based on heuristics for determining when part of +the function name is referring to the type or when the function +modifies the argument. This works well for the majority of functions +but gives a few odd cases. + +### Arguments +The arguments of the function are represented by + +``` @docs +Arblib.ArbCall.Carg +``` + +For the generated method the allowed types for the argument is +determined by + +``` @docs +Arblib.ArbCall.jltype +Arblib.ArbCall.ctype +``` + +Some arguments are automatically converted to keyword arguments. +1. For functions which take a precision argument this argument becomes + a `prec::Integer` keyword argument which is by default set to the + precision of the first argument (if applicable). +2. For functions which take a rounding mode argument this argument + becomes a `rnd::Union{Arblib.arb_rnd,RoundingMode}` keyword + argument which is by default set to `RoundNearest`. +3. For functions which takes a flag argument this argument becomes a + `flag::Integer` keyword argument which is by default set to `0`. +4. For functions which takes an argument giving the length of the + vector preceding the argument this argument becomes a keyword + argument which is by default set to the length of the preceding + vector. In this case the name of the keyword argument is the same + as the argument name in the function declaration. + +As with the naming the implementation is based on heuristics for +determining when an argument is supposed to be a certain kind of +keyword argument. + +### Return value +The returned value is determined in the following way + +1. For functions which have the C function has return type `void` and + modify the first argument the generated method returns the first + argument. This is follows the normal convention in Julia. +2. For predicates, for which the C function returns `int`, the return + value is converted to a `Bool`. +3. Otherwise the return type is the same as for the C function. + +### Examples + +For example Arb declares the following functions + +1. `void arb_zero(arb_t x)` +2. `slong arb_rel_error_bits(const arb_t x)` +3. `int arb_is_zero(const arb_t x)` +4. `void arb_add(arb_t z, const arb_t x, const arb_t y, slong prec)` +5. `void arb_add_arf(arb_t z, const arb_t x, const arf_t y, slong prec)` +6. `void arb_add_ui(arb_t z, const arb_t x, ulong y, slong prec)` +7. `void arb_add_si(arb_t z, const arb_t x, slong y, slong prec)` +8. `void arb_sin(arb_t s, const arb_t x, slong prec)` +9. `void arb_cos(arb_t c, const arb_t x, slong prec)` +10. `void arb_sin_cos(arb_t s, arb_t c, const arb_t x, slong prec)` +11. `int arf_add(arf_t res, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd)` + +For which the following methods are generated + +1. `zero!(x::ArbLike)::ArbLike` +2. `rel_error_bits(x::ArbLike)::Int` +3. `is_zero(x::ArbLike)::Bool` +4. `add!(z::ArbLike, x::ArbLike, y::ArbLike; prec::Integer = _precision(z))::ArbLike` +5. `add!(z::ArbLike, x::ArbLike, y::ArfLike; prec::Integer = _precision(z))::ArbLike` +6. `add!(z::ArbLike, x::ArbLike, y::Unsigned; prec::Integer = _precision(z))::ArbLike` +7. `add!(z::ArbLike, x::ArbLike, y::Integer; prec::Integer = _precision(z))::ArbLike` +8. `sin!(s::ArbLike, x::ArbLike; prec::Integer = _precision(s))::ArbLike` +9. `cos!(c::ArbLike, x::ArbLike; prec::Integer = _precision(c))::ArbLike` +10. `sin_cos!(s::ArbLike, c::ArbLike, x::ArbLike, prec::Integer = _precision(s))::ArbLike` +11. `add!(res::ArfLike, x::ArfLike, y::ArfLike; prec::Integer = _precision(res), rnd::Union{Arblib.arb_rnd, RoundingMode} = RoundNearest)::Int32` diff --git a/docs/src/wrapper-types.md b/docs/src/wrapper-types.md new file mode 100644 index 00000000..363f6ad1 --- /dev/null +++ b/docs/src/wrapper-types.md @@ -0,0 +1,50 @@ +# Types + +The package defines the following types which map directly to Arb +types with the corresponding name. In most cases you should not use +these types directly but use the types from the [High level +interface](interface-types.md). + +``` @docs +Arblib.mag_struct +Arblib.arf_struct +Arblib.arb_struct +Arblib.acb_struct +Arblib.arb_vec_struct +Arblib.acb_vec_struct +Arblib.arb_poly_struct +Arblib.acb_poly_struct +Arblib.arb_mat_struct +Arblib.acb_mat_struct +Arblib.calc_integrate_opt_struct +``` + +For each low-level type there is a union of types that can be +interpreted as the low-level type. These are the types that can be +given directly as arguments to the low-level methods. Below you find +these union-types. + +``` @docs +Arblib.MagLike +Arblib.ArfLike +Arblib.ArbLike +Arblib.AcbLike +Arblib.ArbVectorLike +Arblib.AcbVectorLike +Arblib.ArbMatrixLike +Arblib.AcbMatrixLike +Arblib.ArbPolyLike +Arblib.AcbPolyLike +``` + +Note that the `show` method is overloaded for these union types, this +is to make method declarations easier to read in the REPL. As an +example we can print the methods for `Arblib.sin!` and we see that it +prints `ArbLike` instead of the much longer +`Union{Ptr{Arblib.arb_struct}, Arb, ArbRef, Arblib.arb_struct}`. + +``` @repl +using Arblib + +methods(Arblib.sin!) +``` diff --git a/src/ArbCall/ArbCall.jl b/src/ArbCall/ArbCall.jl index 7d74ce9d..69380fc9 100644 --- a/src/ArbCall/ArbCall.jl +++ b/src/ArbCall/ArbCall.jl @@ -1,3 +1,9 @@ +""" + ArbCall + +This module handles the automatic generation of methods from parsing +the Arb documentation. +""" module ArbCall using ..Arblib diff --git a/src/ArbCall/ArbFPWrapFunction.jl b/src/ArbCall/ArbFPWrapFunction.jl index 60cc0a7f..7311afa7 100644 --- a/src/ArbCall/ArbFPWrapFunction.jl +++ b/src/ArbCall/ArbFPWrapFunction.jl @@ -129,6 +129,12 @@ function jlargs(af::ArbFPWrapFunction) return args, kwargs end +""" + jlcode(af::ArbFPWrapFunction, jl_fname = jlfname(af)) + +Generate the Julia code for calling a function from the `fpwrap` +module of Arb from Julia. +""" function jlcode(af::ArbFPWrapFunction, jl_fname = jlfname(af)) T = basetype(af) cargs = arguments(af) @@ -171,6 +177,26 @@ function jlcode(af::ArbFPWrapFunction, jl_fname = jlfname(af)) return func end +""" + @arbfpwrapcall_str str + +Parse a string as an [`Arblib.ArbCall.ArbFPWrapFunction`](@ref), +generate the code for a corresponding method with +[`Arblib.ArbCall.jlcode`](@ref) and evaluate the code. + +For example +``` +arbfpwrapcall"int arb_fpwrap_double_exp(double * res, double x, int flags)" +``` +defines the method +``` +fpwrap_exp( + x::Union{Float16, Float32, Float64}; + error_on_failure::Bool = Arblib.ArbCall.fpwrap_error_on_failure_default(), + correct_rounding::Bool = false, work_limit::Integer = 8, +) +``` +""" macro arbfpwrapcall_str(str) af = ArbFPWrapFunction(str) return esc(jlcode(af)) diff --git a/src/ArbCall/ArbFunction.jl b/src/ArbCall/ArbFunction.jl index 3e8cae1c..4b6b9de6 100644 --- a/src/ArbCall/ArbFunction.jl +++ b/src/ArbCall/ArbFunction.jl @@ -1,8 +1,25 @@ +""" + ArbFunction{T}(fname::String, args::Vector{Carg}) + +Struct representing a C-function in the Arb library. +""" struct ArbFunction{T} fname::String args::Vector{Carg} end +""" + ArbFunction(str::AbstractString) + +Parse a string as an ArbFunction. The string should be as a function +declaration in C. Note that note all types of function declarations +are supported, only those which are relevant for the Arb library. + +```jldoctest +julia> Arblib.ArbCall.ArbFunction("void arb_zero(arb_t x)") +Arblib.ArbCall.ArbFunction{Nothing}("arb_zero", Arblib.ArbCall.Carg[Arblib.ArbCall.Carg{Arb}(:x, false)]) +``` +""" function ArbFunction(str::AbstractString) m = match(r"(?\w+(\s\*)?)\s+(?[\w_]+)\((?.*)\)", str) isnothing(m) && @@ -119,6 +136,11 @@ function jlargs(af::ArbFunction; argument_detection::Bool = true) return args, kwargs end +""" + jlcode(af::ArbFunction, jl_fname = jlfname(af)) + +Generate the Julia code for calling the Arb function from Julia. +""" function jlcode(af::ArbFunction, jl_fname = jlfname(af)) jl_args, jl_kwargs = jlargs(af, argument_detection = true) jl_full_args, _ = jlargs(af, argument_detection = false) @@ -156,6 +178,19 @@ function jlcode(af::ArbFunction, jl_fname = jlfname(af)) end end +""" + @arbcall_str str + +Parse a string as an [`Arblib.ArbCall.ArbFunction`](@ref), generate +the code for a corresponding method with +[`Arblib.ArbCall.jlcode`](@ref) and evaluate the code. + +For example +``` +arbcall"void arb_zero(arb_t x)" +``` +defines the method `zero!(x::ArbLike)`. +""" macro arbcall_str(str) af = ArbFunction(str) return esc(jlcode(af)) diff --git a/src/ArbCall/Carg.jl b/src/ArbCall/Carg.jl index 716402ec..cef969fc 100644 --- a/src/ArbCall/Carg.jl +++ b/src/ArbCall/Carg.jl @@ -48,12 +48,12 @@ end jltype(ca::Carg{T}) The most general Julia type for which we allow automatic conversion to -the [`ctype`](@ref) of `ca`. +the [`Arblib.ArbCall.ctype`](@ref) of `ca`. These conversations should be done without any loss of information, for example for floating point numbers we only allow conversion from types with lower precision. In general the conversion is done using -[`Base.cconvert`](@ref). +`Base.cconvert`. """ jltype(ca::Carg) = rawtype(ca) jltype(::Carg{Cint}) = Integer diff --git a/src/ArbCall/parse.jl b/src/ArbCall/parse.jl index 635cdf3b..dfe3c3e7 100644 --- a/src/ArbCall/parse.jl +++ b/src/ArbCall/parse.jl @@ -73,9 +73,9 @@ end generate_file(filename, title, sections; verbose, manual_overrides) Given a title and a list of sections, as returned by -[`parse_arbdoc`](@ref), return a string with the Julia code to load -all of those functions. If a file name is given then also write the -string to that file. +[`Arblib.ArbCall.parse_arbdoc`](@ref), return a string with the Julia +code to load all of those functions. If a file name is given then also +write the string to that file. For each function it checks that it is correctly parsed. Those that do not parse due to types which are not supported and likely never will diff --git a/src/arb_types.jl b/src/arb_types.jl index 4121abe8..1fb275c8 100644 --- a/src/arb_types.jl +++ b/src/arb_types.jl @@ -1,3 +1,6 @@ +""" + arf_struct +""" mutable struct arf_struct exponent::UInt # fmpz size::UInt # mp_size_t @@ -19,6 +22,9 @@ mutable struct arf_struct end end +""" + mag_struct +""" mutable struct mag_struct exponent::UInt # fmpz mantissa::UInt # mp_limb_t @@ -38,6 +44,9 @@ mutable struct mag_struct end end +""" + arb_struct +""" mutable struct arb_struct # ┌ arf_struct (midpoint) exponent::UInt # │ fmpz @@ -58,6 +67,9 @@ mutable struct arb_struct end end +""" + acb_struct +""" mutable struct acb_struct # ┌ arb_struct (real) exponent_r::UInt # │ fmpz @@ -83,6 +95,9 @@ mutable struct acb_struct end end +""" + arb_vec_struct +""" mutable struct arb_vec_struct entries::Ptr{arb_struct} n::Int @@ -94,6 +109,9 @@ mutable struct arb_vec_struct end end +""" + acb_vec_struct +""" mutable struct acb_vec_struct entries::Ptr{acb_struct} n::Int @@ -105,6 +123,9 @@ mutable struct acb_vec_struct end end +""" + arb_poly_struct +""" mutable struct arb_poly_struct coeffs::Ptr{arb_struct} length::Int @@ -118,6 +139,9 @@ mutable struct arb_poly_struct end end +""" + acb_poly_struct +""" mutable struct acb_poly_struct coeffs::Ptr{acb_struct} length::Int @@ -131,6 +155,9 @@ mutable struct acb_poly_struct end end +""" + arb_mat_struct +""" mutable struct arb_mat_struct entries::Ptr{arb_struct} r::Int @@ -145,6 +172,9 @@ mutable struct arb_mat_struct end end +""" + acb_mat_struct +""" mutable struct acb_mat_struct entries::Ptr{acb_struct} r::Int @@ -199,6 +229,9 @@ function Base.deepcopy_internal( return y end +""" + calc_integrate_opt_struct +""" mutable struct calc_integrate_opt_struct deg_limit::Int eval_limit::Int diff --git a/src/calc_integrate.jl b/src/calc_integrate.jl index d89ca49b..270d90cb 100644 --- a/src/calc_integrate.jl +++ b/src/calc_integrate.jl @@ -165,7 +165,9 @@ julia> Arblib.integrate!(Arblib.inv!, Acb(0), Acb(1, -5), Acb(1, 5)) # Integrate [+/- 2.02e-75] + [2.74680153389003172172254385288992229730199919179940161793956671182574846633 +/- 2.83e-75]im julia> # Integrate √z from 1 - 5im to 1 + 5im, taking into account the branch cut at (-∞, 0] + julia> f! = (res, z; analytic = false) -> Arblib.sqrt_analytic!(res, z, analytic); + julia> Arblib.integrate!(f!, Acb(0), Acb(1, -5), Acb(1, 10), check_analytic = true, prec = 64) [-9.0064084416559764 +/- 6.53e-17] + [23.8636067095598007 +/- 6.98e-17]im ``` @@ -298,6 +300,8 @@ paper > https://doi.org/10.1007/978-3-319-96418-8 > https://arxiv.org/abs/1802.07942 +See also: [`integrate!`](@ref). + # Examples ```jldoctest julia> Arblib.integrate(sin, 0, 10) # Integrate sin from 0 to 10 @@ -307,6 +311,7 @@ julia> Arblib.integrate(z -> 1/z, Acb(1, -5), Acb(1, 5)) # Integrate 1/z from 1 [+/- 2.02e-75] + [2.74680153389003172172254385288992229730199919179940161793956671182574846633 +/- 2.83e-75]im julia> # Integrate √z from 1 - 5im to 1 + 5im, taking into account the branch cut at (-∞, 0] + julia> f = (z; analytic = false) -> begin if analytic && Arblib.contains_nonpositive(real(z)) return Acb(NaN, prec = precision(z)) @@ -314,6 +319,7 @@ julia> f = (z; analytic = false) -> begin return sqrt(z) end end; + julia> Arblib.integrate(f, Acb(1, -5), Acb(1, 10), check_analytic = true, prec = 64) [-9.0064084416559764 +/- 7.40e-17] + [23.8636067095598007 +/- 9.03e-17]im ``` diff --git a/src/interval.jl b/src/interval.jl index 27556499..de428124 100644 --- a/src/interval.jl +++ b/src/interval.jl @@ -178,7 +178,7 @@ Returns a copy of `x` with the absolute value of `err` added to the radius. For complex `x` it adds the error to both the real and imaginary parts. For matrices it adds it elementwise. -See also [`set_ball`](@ref). +See also [`setball`](@ref). """ add_error(x::Union{ArbOrRef,AcbOrRef}, err::Union{MagOrRef,ArfOrRef,ArbOrRef}) = add_error!(copy(x), err) diff --git a/src/poly.jl b/src/poly.jl index db510f88..46dd8392 100644 --- a/src/poly.jl +++ b/src/poly.jl @@ -93,12 +93,12 @@ as the degree of the underlying polynomial (in case higher order coefficients are zero) but the way they are constructed ensures that the coefficients will always be initialised. -!!! Note: If you use this to change the coefficient in a way so that the degree - of the polynomial might change you need to normalise the polynomial - afterwards to make sure that Arb recognises the possibly new degree of - the polynomial. If the new degree is the same or lower this can be - done using [`Arblib.normalise!`](@ref). If the new degree is higher - you need to manually set the correct value of +!!! Note: If you use this to change the coefficient in a way so that + the degree of the polynomial might change you need to normalise + the polynomial afterwards to make sure that Arb recognises the + possibly new degree of the polynomial. If the new degree is the + same or lower this can be done using `Arblib.normalise!`. If the + new degree is higher you need to manually set the correct value of `Arblib.cstruct(p).length` to be one higher than the new degree. """ Base.@propagate_inbounds function ref(p::Union{ArbPoly,ArbSeries}, i::Integer) diff --git a/src/ref.jl b/src/ref.jl index 1e758370..a26fe8c6 100644 --- a/src/ref.jl +++ b/src/ref.jl @@ -46,19 +46,41 @@ Base.getindex(x::ArfRef) = Arf(x) Base.getindex(x::ArbRef) = Arb(x) Base.getindex(x::AcbRef) = Acb(x) +""" + midref(x::ArbLike, prec = precision(x)) + +Return an `ArfRef` referencing the midpoint of `x`. +""" function midref(x::ArbLike, prec = precision(x)) mid_ptr = ccall(@libarb(arb_mid_ptr), Ptr{arf_struct}, (Ref{arb_struct},), x) ArfRef(mid_ptr, prec, parentstruct(x)) end + +""" + radref(x::ArbLike, prec = precision(x)) + +Return a `MagRef` referencing the radius of `x`. +""" function radref(x::ArbLike) rad_ptr = ccall(@libarb(arb_rad_ptr), Ptr{mag_struct}, (Ref{arb_struct},), x) MagRef(rad_ptr, parentstruct(x)) end +""" + realref(z::AcbLike, prec = precision(z)) + +Return an `ArbRef` referencing the real part of `x`. +""" function realref(z::AcbLike; prec = precision(z)) real_ptr = ccall(@libarb(acb_real_ptr), Ptr{arb_struct}, (Ref{acb_struct},), z) ArbRef(real_ptr, prec, parentstruct(z)) end + +""" + imagref(z::AcbLike, prec = precision(z)) + +Return an `ArbRef` referencing the imaginary part of `x`. +""" function imagref(z::AcbLike; prec = precision(z)) real_ptr = ccall(@libarb(acb_imag_ptr), Ptr{arb_struct}, (Ref{acb_struct},), z) ArbRef(real_ptr, prec, parentstruct(z)) diff --git a/src/setters.jl b/src/setters.jl index d3a4d2d4..9f7aac71 100644 --- a/src/setters.jl +++ b/src/setters.jl @@ -86,7 +86,7 @@ end function set!(res::ArbLike, (a, b)::Tuple{<:Real,<:Real}; prec::Integer = precision(res)) # This is not strictly required to check since the union will give - # an enclosure anyway. But since thus method is designed for a <= + # an enclosure anyway. But since this method is designed for a <= # b adding this check could catch some bugs. a > b && throw(ArgumentError("must have a <= b, got a = $a and b = $b")) if !(a isa ArbLike) diff --git a/src/show.jl b/src/show.jl index 97bbc979..7d283fd6 100644 --- a/src/show.jl +++ b/src/show.jl @@ -128,6 +128,8 @@ for T in [ :AcbVectorLike, :ArbMatrixLike, :AcbMatrixLike, + :ArbPolyLike, + :AcbPolyLike, ] @eval Base.show(io::IO, ::Type{$T}) = print(io, $(QuoteNode(T))) end diff --git a/src/types.jl b/src/types.jl index 344ea830..8bdb4f44 100644 --- a/src/types.jl +++ b/src/types.jl @@ -273,16 +273,48 @@ const MagOrRef = Union{Mag,MagRef} const ArfOrRef = Union{Arf,ArfRef} const ArbOrRef = Union{Arb,ArbRef} const AcbOrRef = Union{Acb,AcbRef} + +""" + MagLike = Union{Mag,MagRef,mag_struct,Ptr{mag_struct}} +""" const MagLike = Union{Mag,MagRef,cstructtype(Mag),Ptr{cstructtype(Mag)}} +""" + ArfLike = Union{Arf,ArfRef,arf_struct,Ptr{arf_struct}}} +""" const ArfLike = Union{Arf,ArfRef,cstructtype(Arf),Ptr{cstructtype(Arf)}} +""" + ArbLike = Union{Arb,ArbRef,arb_struct,Ptr{arb_struct}}} +""" const ArbLike = Union{Arb,ArbRef,cstructtype(Arb),Ptr{cstructtype(Arb)}} +""" + AcbLike = Union{Acb,AcbRef,acb_struct,Ptr{acb_struct}}} +""" const AcbLike = Union{Acb,AcbRef,cstructtype(Acb),Ptr{cstructtype(Acb)}} +""" + ArbVectorLike = Union{ArbVector,ArbRefVector,arb_vec_struct} +""" const ArbVectorLike = Union{ArbVector,ArbRefVector,cstructtype(ArbVector)} +""" + AcbVectorLike = Union{AcbVector,AcbRefVector,acb_vec_struct} +""" const AcbVectorLike = Union{AcbVector,AcbRefVector,cstructtype(AcbVector)} -const ArbMatrixLike = Union{ArbMatrix,ArbRefMatrix,cstructtype(ArbMatrix)} -const AcbMatrixLike = Union{AcbMatrix,AcbRefMatrix,cstructtype(AcbMatrix)} +""" + ArbPolyLike = Union{ArbPoly,ArbSeries,arb_poly_struct} +""" const ArbPolyLike = Union{ArbPoly,ArbSeries,cstructtype(ArbPoly)} +""" + AcbPolyLike = Union{AcbPoly,AcbSeries,acb_poly_struct} +""" const AcbPolyLike = Union{AcbPoly,AcbSeries,cstructtype(AcbPoly)} +""" + ArbMatrixLike = Union{ArbMatrix,ArbRefMatrix,arb_mat_struct)} +""" +const ArbMatrixLike = Union{ArbMatrix,ArbRefMatrix,cstructtype(ArbMatrix)} +""" + AcbMatrixLike = Union{AcbMatrix,AcbRefMatrix,acb_mat_struct} +""" +const AcbMatrixLike = Union{AcbMatrix,AcbRefMatrix,cstructtype(AcbMatrix)} + const ArbTypes = Union{ Mag, MagRef, diff --git a/test/runtests.jl b/test/runtests.jl index 151d9273..9e83e6a4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,8 +1,13 @@ ENV["NEMO_THREADED"] = 1 using Arblib, Test, LinearAlgebra, Random, Serialization, SpecialFunctions +using Documenter + +DocMeta.setdocmeta!(Arblib, :DocTestSetup, :(using Arblib); recursive = true) @testset "Arblib" begin + doctest(Arblib) + include("ArbCall/runtests.jl") include("arb_types.jl")