GitHub

Week 9 Lecture 1: Overview of FLINT and Arblib.jl

We have now covered many of the basic ideas in interval arithmetic and rigorous numerics. There are, however, still lots of things we haven't talked about. The goal of this week is to get a broad overview of all the capabilities that FLINT/Arblib.jl and IntervalArithmetic.jl have, starting with FLINT/Arblib.jl today.

Our primary references for today are the FLINT documentation for real and complex numbers as well as the Arblib.jl documentation.

FLINT vs Arblib.jl

Arblib.jl is a wrapper of the FLINT C-library. It doesn't implement any functionality by itself, but simply wraps the FLINT functions in a way that makes them convenient to use from Julia.

The Arblib.jl interface is split into two parts, a high level interface that works well with general Julia code, and a low level interface that gives more control. The low level interface wraps almost all of the functionality of FLINT, whereas the high level interface only wraps some parts of it.

Types

So far we have primarily made use of the Arb type, representing a ball with a midpoint and a radius, and the ArbSeries type, representing a truncated Taylor expansion with coefficients given by Arb balls. Let us take a look at the other relevant FLINT types.

There are 8 basic types in FLINT that are relevant for us. In the list below we give both the FLINT name for the type, as well as the name they have in Arblib.jl.

  • mag_t / Mag: Low-level type used to represent the radius of an Arb ball. It's a fixed (31-bit) precision floating point that only supports positive numbers. It has very limited functionality by itself.
  • arf_t / Arf: Low-level type used to represent the midpoint of an Arb ball. It's an arbitrary precision floating point. Similar to Mag it has very limited functionality by itself.
  • arb_t / Arb: This is, for us, the most fundamental type. It consists of an Arf midpoint and a Mag radius. It implements a lot of functionality and, unless otherwise specified, all operations are fully rigorous.
  • acf_t / Acf: This represents a complex number and consists of a pair of Arf values, representing the real and imaginary parts. Similar to Arf it has very limited functionality by itself.
  • acb_t / Acb: This represents a complex number and consists of a pair of Arb values, representing the real and imaginary parts. Note that this is not a proper complex ball, but instead represents a rectangular region. In many cases it is however referred to as a complex ball. Similar to Arb it implements a lot of functionality and, unless otherwise specified, all operations are fully rigorous.
  • arb_poly_t / ArbPoly: Represents a polynomial with coefficients given by Arb values. In particular this is also the type used to represent truncated Taylor series. Unless otherwise specified, operations are rigorous.
  • acb_poly_t / AcbPoly: Similar to ArbPoly, but the coefficients are Acb values.
  • arb_mat_t / ArbMatrix: Represents a matrix with coefficients given by Arb values. It implements a limited set of linear algebra routines, more about that later. Unless otherwise specified, operations are rigorous.
  • acb_mat_t / AcbMatrix: Similar to ArbMatrix, but the coefficients are Acb values.

The ArbSeries type (and also the AcbSeries type) in Arblib.jl does not have a direct correspondence in FLINT. It consists of an ArbPoly together with an integer specifying the degree of the expansion. In FLINT these values are kept separate, so functions will take the polynomial and the degree as two separate arguments.

FLINT also has a type nfloat_t for representing n-word floating points. This type is relatively new and has yet to be wrapped in Arblib.jl. In general the operations with nfloat_t are not rigorous, so it is also less relevant for computer-assisted proofs.

Here are some examples of constructing these types and performing basic arithmetic on them.

julia> using Arblib
julia> setprecision(Arb, 64)64
julia> Mag(1)1.0
julia> Mag(1) + Mag(3) # Note that this is rounded upwards!4.00000001
julia> Arf(1) / Arf(3)0.3333333333333333333
julia> sin(Arf(1)) # Most functions are not implemented for ArfERROR: MethodError: no method matching sin(::Arblib.Arf)
This error has been manually thrown, explicitly, so the method may exist but be intentionally marked as unimplemented.

Closest candidates are:
  sin(::Real)
   @ Base math.jl:1543
  sin(::Irrational{:π})
   @ Base mathconstants.jl:146
  sin(::ComplexF16)
   @ Base math.jl:1527
  ...
julia> Arb(1 // 3) # Arb we have already worked with a lot![0.3333333333333333333 +/- 4.24e-20]
julia> Acb(1, π)1.0 + [3.141592653589793239 +/- 5.96e-19]im
julia> exp(Acb(1, π))[-2.718281828459045235 +/- 5.91e-19] + [+/- 4.32e-19]im
julia> ArbPoly([1, 2, 3])1.0 + 2.0⋅x + 3.0⋅x^2
julia> ArbPoly([1, 2, 3])^51.0 + 10.0⋅x + 55.0⋅x^2 + 200.0⋅x^3 + 530.0⋅x^4 + 1052.0⋅x^5 + 1590.0⋅x^6 + 1800.0⋅x^7 + 1485.0⋅x^8 + 810.0⋅x^9 + 243.0⋅x^10
julia> AcbPoly([Acb(1, 2), Acb(3, 4), Acb(5, 6)])(1.0 + 2.0im) + (3.0 + 4.0im)⋅x + (5.0 + 6.0im)⋅x^2
julia> AcbPoly([Acb(1, 1 // 3), Acb(3, 4), Acb(5, 6)])^2([0.888888888888888889 +/- 1.96e-19] + [0.6666666666666666666 +/- 8.48e-20]im) + ([3.333333333333333333 +/- 6.95e-19] + [10.0000000000000000 +/- 1.90e-18]im)⋅x + ([-1.000000000000000000 +/- 7.59e-19] + [39.33333333333333333 +/- 6.79e-18]im)⋅x^2 + (-18.0 + 76.0im)⋅x^3 + (-11.0 + 60.0im)⋅x^4
julia> ArbMatrix(2, 3) # Zero matrix of size 3 x 32×3 Arblib.ArbMatrix:
 0  0  0
 0  0  0
julia> ArbMatrix([1 2 3; 4 5 6; 7 8 9])3×3 Arblib.ArbMatrix:
 1.0  2.0  3.0
 4.0  5.0  6.0
 7.0  8.0  9.0
julia> ArbMatrix([1 2 3; 4 5 6; 7 8 9])^23×3 Arblib.ArbMatrix:
  30.0   36.0   42.0
  66.0   81.0   96.0
 102.0  126.0  150.0
julia> AcbMatrix(2, 3) # Zero matrix of size 3 x 32×3 Arblib.AcbMatrix:
 0  0  0
 0  0  0
julia> AcbMatrix([1 2 (3 + 2im); 4 5 6; 7 8 9])3×3 Arblib.AcbMatrix:
 1.0  2.0  3.0 + 2.0im
 4.0  5.0          6.0
 7.0  8.0          9.0
julia> AcbMatrix([1 2 (3 + 2im); 4 5 6; 7 8 9])^23×3 Arblib.AcbMatrix:
 30.0 + 14.0im  36.0 + 16.0im   42.0 + 20.0im
          66.0           81.0    96.0 + 8.0im
         102.0          126.0  150.0 + 14.0im

FLINT documentation

For each of the types discussed above, there is an associated page in the FLINT documentation. For example, take a look at the following pages:

Functionality

Let us take a look at some of the functionality that FLINT implements. Some of the most important groups are:

  1. Special functions, including support for Taylor arithmetic
  2. Linear algebra routines
  3. Polynomial routines
  4. Rigorous integration
  5. Basic root finding
  6. FFT

Special functions

This is the largest class of functionality and the documentation is spread out over many different pages. For example, the $\zeta$ function is described in the documentation for arb.h and acb.h, whereas hypergeometric functions have their own documentation in acb_hypgeom.h. There are also separate modules for elliptic integrals (acb_elliptic.h), modular forms (acb_modular.h) theta functions (acb_theta.h) and Dirichlet functions (acb_dirichlet.h).

In many cases, there are implementations for direct evaluation of the value as well as computing series expansions. For some functions, evaluation of the series expansion is however not implemented.

The high level Arblib.jl interface wraps all functions that are also in the Julia package SpecialFunctions.jl. FLINT does however implement many more special functions, these can only be accessed through the low level interface.

julia> using SpecialFunctions
julia> gamma(Arb(0.5)) # This is accessible in the high level interface[1.772453850905516027 +/- 3.11e-19]
julia> gamma(AcbSeries((1 + 2im, 1)))([0.1519040026700361374 +/- 5.54e-20] + [0.01980488016185498197 +/- 5.26e-21]im) + ([0.082390881509719160 +/- 3.00e-19] + [0.2147883123157484927 +/- 7.91e-20]im)⋅x + 𝒪(x^2)
julia> Arblib.dirichlet_lerch_phi!(Acb(), Acb(1), Acb(2), Acb(3)) # This needs the low level interface[0.3949340668482264365 +/- 2.82e-20]

Linear algebra routines

These are documented in arb_mat.h and acb_mat.h. It implements solvers for linear systems (which also allows for computation of inverses and determinants) and computation of eigenvalues.

Most of these routines are accessible from the high level interface, though the low level interface can give more control. For eigenvalues, only AcbMatrix is supported.

julia> using LinearAlgebra
julia> M = ArbMatrix([1 2; 3 4])2×2 Arblib.ArbMatrix:
 1.0  2.0
 3.0  4.0
julia> v = ArbMatrix([1, 2])2×1 Arblib.ArbMatrix:
 1.0
 2.0
julia> M \ v2×1 Arblib.ArbMatrix:
 [+/- 2.17e-19]
     [0.500000000000000000 +/- 1.63e-19]
julia> det(M)-2.0
julia> A = AcbMatrix(M)2×2 Arblib.AcbMatrix:
 1.0  2.0
 3.0  4.0
julia> eigvals(A)2-element Arblib.AcbVector:
 [-0.37228132326901433 +/- 1.71e-18] + [+/- 1.36e-18]im
  [5.37228132326901433 +/- 4.00e-18] + [+/- 1.36e-18]im

Polynomial routines

These are documented in arb_poly.h and acb_poly.h. It implements evaluation (including multipoint evaluation), composition, interpolation and root finding.

Some of these are available from the high level interface, others require the low level interface.

julia> p = ArbPoly([1, 2, 3, 4])1.0 + 2.0⋅x + 3.0⋅x^2 + 4.0⋅x^3
julia> p(Arb(1 // 3))[2.148148148148148148 +/- 4.63e-19]
julia> Arblib.compose(p, p)10.0 + 40.0⋅x + 120.0⋅x^2 + 292.0⋅x^3 + 519.0⋅x^4 + 768.0⋅x^5 + 924.0⋅x^6 + 816.0⋅x^7 + 576.0⋅x^8 + 256.0⋅x^9
julia> Arblib.evaluate_vec_fast!(ArbVector(4), p, ArbVector([7, 8, 9, 10]))4-element Arblib.ArbVector:
 1534.0
 2257.0
 3178.0
 4321.0

Rigorous integration

This is documented in acb_calc.h. It implements rigorous integration of holomorphic functions. It can deal with non-holomorphic functions as well, but with much slower convergence in the regions where they are not holomorphic.

This is available through the high level interface.

julia> Arblib.integrate(x -> sin(exp(x)), Arb(0), Arb(8))[0.62502040893251 +/- 3.75e-15]
julia> Arblib.integrate(x -> besselj(Acb(0), x)^2 + sin(2x), Acb(1 + im), Acb(7im))[-300653.10315846613 +/- 3.67e-12] + [15481.714432735356 +/- 6.41e-13]im

Handling functions which are not meromorphic requires a bit of extra work. The documentation of Arblib.integrate has more details.

Basic root finding

This is documented in arb_calc.h. It implements root finding using bisection and Newton iterations.

This is not directly accessible through either the high level or low level interface. It requires passing functions as arguments, which needs extra work in Julia. For my use case I primarily use similar methods implemented in ArbExtras.jl.

FFT

This is documented in acb_dft.h. Most of the documentation is very abstract, talking about discrete Fourier transforms over general Abelian groups. But it does implement the usual FFT with acb_dft.

This is only available from the low level interface. Currently the wrapping doesn't support precomputations.

julia> n = 1616
julia> xs = range(Arb(0), 2Arb(π), n + 1)[1:end-1]16-element LinRange{Arblib.Arb, Int64}:
 0, [0.39{…14 digits…}49 ± 4.09e-19], …, [5.89{…14 digits…}22 ± 7.24e-19]
julia> ys = exp.(im * xs) + 3exp.(2im * xs) + 2im * exp.(3im * xs)16-element Vector{Complex{Arblib.Arb}}:
                                                          4.0 + 2.0im
  [1.19744081104835586 +/- 6.10e-18] + [3.26937064065491189 +/- 2.22e-18]im
 [-0.70710678118654750 +/- 4.40e-18] + [2.29289321881345251 +/- 6.31e-18]im
 [-0.97327004646437323 +/- 7.13e-18] + [1.19744081104835571 +/- 7.03e-18]im
  [-0.99999999999999998 +/- 3.03e-18] + [1.0000000000000000 +/- 5.90e-18]im
 [-1.73863691119455256 +/- 7.82e-18] + [0.65031825397421795 +/- 4.53e-18]im
  [-2.1213203435596424 +/- 4.16e-17] - [0.87867965644035807 +/- 6.23e-18]im
 [-0.65031825397421743 +/- 6.11e-18] - [2.50400377592473244 +/- 8.30e-18]im
   [1.9999999999999997 +/- 4.14e-17] - [2.00000000000000022 +/- 5.70e-18]im
   [3.04519987607092953 +/- 9.17e-18] + [0.9732700464643721 +/- 1.51e-17]im
   [0.7071067811865496 +/- 1.49e-17] + [3.70710678118654675 +/- 7.76e-18]im
   [-3.2693706406549092 +/- 3.25e-17] + [3.0451998760709314 +/- 5.12e-17]im
  [-4.99999999999999974 +/- 4.13e-18] - [1.0000000000000031 +/- 4.68e-17]im
   [-2.5040037759247305 +/- 3.04e-17] - [4.8929589410935039 +/- 2.01e-17]im
    [2.1213203435596435 +/- 3.54e-17] - [5.1213203435596421 +/- 5.43e-17]im
   [4.89295894109350284 +/- 8.53e-18] - [1.7386369111945528 +/- 1.66e-17]im
julia> Arblib.dft!(AcbVector(n), AcbVector(ys))16-element Arblib.AcbVector:
             [9e-15 +/- 5.37e-16] + [-4e-15 +/- 4.69e-16]im
 [16.000000000000001 +/- 2.57e-16] + [8e-15 +/- 5.77e-16]im
       [47.999999999999997 +/- 2.35e-16] + [+/- 4.12e-16]im
       [+/- 1.76e-16] + [32.000000000000003 +/- 4.62e-16]im
            [-7e-15 +/- 1.61e-16] + [-2e-15 +/- 3.49e-16]im
             [5e-15 +/- 1.76e-16] + [-8e-15 +/- 5.23e-16]im
              [4e-15 +/- 3.88e-16] + [5e-15 +/- 3.82e-16]im
             [-3e-15 +/- 3.19e-16] + [3e-15 +/- 5.60e-16]im
            [-2e-15 +/- 2.31e-16] + [-4e-15 +/- 3.39e-16]im
                    [4e-15 +/- 5.98e-16] + [+/- 2.12e-16]im
                    [+/- 8.75e-16] + [4e-15 +/- 1.98e-16]im
            [-4e-15 +/- 4.85e-16] + [-1e-15 +/- 4.83e-16]im
                   [+/- 8.13e-16] + [-3e-15 +/- 1.52e-16]im
              [4e-15 +/- 2.05e-16] + [3e-15 +/- 5.40e-16]im
             [-4e-15 +/- 4.52e-16] + [3e-15 +/- 5.05e-16]im
            [-4e-15 +/- 1.96e-16] + [-5e-15 +/- 2.99e-16]im