fredrikj.net / blog /

# Arb example programs, rigorous calculus

*September 26, 2013*

I have starting adding example programs to Arb. Right now the following programs are available:

- pi.c – computes $\pi$ to any precision (this is mandatory to have…)
- hilbert_matrix.c – computes determinants of huge Hilbert matrices with adaptive precision
- keiper_li.c – computes the Keiper-Li coefficients (see the previous post on testing Li’s criterion and my new paper)
- real_roots.c – rigorously isolates all single roots of a real analytic function on an interval

See the documentation page for more detailed descriptions and some sample output.

I will talk a bit more about the last example program. It is is based on the module fmprb_calc which I’ve just added to Arb. Quoting from the documentation page:

This module provides functions for operations of calculus over the real numbers (intended to include root-finding, optimization, integration, and so on). It is planned that the module will include two types of algorithms:

- Interval algorithms that give provably correct results. An example would be numerical integration on an interval by dividing the interval into small balls and evaluating the function on each ball, giving rigorous upper and lower bounds.
- Conventional numerical algorithms that use heuristics to estimate the accuracy of a result, without guaranteeing that it is correct. An example would be numerical integration based on pointwise evaluation, where the error is estimated by comparing the results with two different sets of evaluation points. Ball arithmetic then still tracks the accuracy of the function evaluations.

Right now, the *fmprb_calc* module contains a few algorithms for rigorously finding the roots of real analytic functions. I hope to find the time to add numerical integration soon. Expect a complex analysis module to appear in the future too (the complex counterpart should eventually be much more powerful, as it will be able to use the omnipotent Cauchy integral formula and other magic tricks).

For now, using the root-finding routines takes a bit of work. It’s quite difficult to write a black-box algorithm that does what the user expects: you’ll always have some tradeoffs (for example, if convergence seems slow, should it fail early and give the user the chance to try something else, or should it try hard to solve the problem but possibly waste a lot of time; should it be optimized for one type of input or the other?). So my philosophy here is to implement several algorithms that give precise control over what goes on, and let the user connect them in the most appropriate way for a given situation. The real_roots program shows one way to do this. I’m lazily copying its description from the documentation page:

This program isolates the roots of a function on the interval $(a,b)$ (where $a$ and $b$ are input as double-precision literals) using the routines in the *fmprb_calc* module. The program takes the following arguments:

real_roots function a b [-refine d] [-verbose] [-maxdepth n] [-maxeval n] [-maxfound n] [-prec n]

The following functions (specified by an integer code) are implemented:

* 0 – $Z(x)$ (Riemann-Siegel Z-function)

* 1 – $\sin(x)$

* 2 – $\sin(x^2)$

* 3 – $\sin(1/x)$

The following options are available:

`-refine d`: If provided, after isolating the roots, attempt to refine the roots to*d*digits of accuracy using a few bisection steps followed by Newton’s method with adaptive precision, and then print them.`-verbose`: Print more information.`-maxdepth n`: Stop searching after*n*recursive subdivisions.`-maxeval n`: Stop searching after approximately*n*function evaluations (the actual number evaluations will be a small multiple of this).`-maxfound n`: Stop searching after having found*n*isolated roots.`-prec n`: Working precision to use for the root isolation.

With *function* 0, the program isolates roots of the Riemann zeta function on the critical line, and guarantees that no roots are missed (there are more efficient ways to do this, but it is a nice example):

> build/examples/real_roots 0 0.0 50.0 -verbose interval: 25 +/- 25 maxdepth = 30, maxeval = 100000, maxfound = 100000, low_prec = 30 found isolated root in: 14.12353515625 +/- 0.012207 found isolated root in: 21.0205078125 +/- 0.024414 found isolated root in: 25.0244140625 +/- 0.024414 found isolated root in: 30.43212890625 +/- 0.012207 found isolated root in: 32.9345703125 +/- 0.024414 found isolated root in: 37.5732421875 +/- 0.024414 found isolated root in: 40.9423828125 +/- 0.024414 found isolated root in: 43.32275390625 +/- 0.012207 found isolated root in: 48.01025390625 +/- 0.012207 found isolated root in: 49.76806640625 +/- 0.012207 --------------------------------------------------------------- Found roots: 10 Subintervals possibly containing undetected roots: 0 Function evaluations: 3425 cpu/wall(s): 1.22 1.229 virt/peak/res/peak(MB): 20.63 20.66 2.23 2.23

Find just one root and refine it to approximately 75 digits:

> build/examples/real_roots 0 0.0 50.0 -maxfound 1 -refine 75 interval: 25 +/- 25 maxdepth = 30, maxeval = 100000, maxfound = 1, low_prec = 30 refined root: 14.134725141734693790457251983562470270784257115699243175685567460149963429809 +/- 8.4532e-81 --------------------------------------------------------------- Found roots: 1 Subintervals possibly containing undetected roots: 8 Function evaluations: 992 cpu/wall(s): 0.41 0.415 virt/peak/res/peak(MB): 20.76 20.76 2.23 2.23

Find roots of $\sin(x^2)$ on $(0,100)$. The algorithm cannot isolate the root at $x = 0$ (it is at the endpoint of the interval, and in any case a root of multiplicity higher than one). The failure is reported:

> build/examples/real_roots 2 0 100 interval: 50 +/- 50 maxdepth = 30, maxeval = 100000, maxfound = 100000, low_prec = 30 --------------------------------------------------------------- Found roots: 3183 Subintervals possibly containing undetected roots: 1 Function evaluations: 34058 cpu/wall(s): 0.26 0.263 virt/peak/res/peak(MB): 20.73 20.76 1.72 1.72

This does not miss any roots:

> build/examples/real_roots 2 1 100 interval: 50.5 +/- 49.5 maxdepth = 30, maxeval = 100000, maxfound = 100000, low_prec = 30 --------------------------------------------------------------- Found roots: 3183 Subintervals possibly containing undetected roots: 0 Function evaluations: 34039 cpu/wall(s): 0.26 0.266 virt/peak/res/peak(MB): 20.73 20.76 1.70 1.70

Looking for roots of $\sin(1/x)$ on $(0,1)$, the algorithm finds many roots, but will never find all of them since there are infinitely many:

> build/examples/real_roots 3 0.0 1.0 interval: 0.5 +/- 0.5 maxdepth = 30, maxeval = 100000, maxfound = 100000, low_prec = 30 --------------------------------------------------------------- Found roots: 10198 Subintervals possibly containing undetected roots: 24695 Function evaluations: 202587 cpu/wall(s): 1.73 1.731 virt/peak/res/peak(MB): 21.84 22.89 2.76 2.76

Arbitrary-precision ball arithmetic is not ideal for root-isolation: most of the time double precision interval arithmetic would be much faster. But it does work, even if you really try to push it. For example, if you increase the default termination bounds, it will find all the 318309 roots of $\sin(x^2)$ on the interval $1 < x < 1000$ (proving that no others exist) in 25 seconds:

> build/examples/real_roots 2 1.0 1000.0 -maxeval 10000000 -maxfound 1000000 -maxdepth 40 -prec 53 interval: 500.5 +/- 499.5 maxdepth = 40, maxeval = 10000000, maxfound = 1000000, low_prec = 53 ————————————————————— Found roots: 318309 Subintervals possibly containing undetected roots: 0 Function evaluations: 3230804 cpu/wall(s): 25.23 25.256 virt/peak/res/peak(MB): 31.56 38.64 12.54 12.54

Finally, very high precision is no problem (computing a million digits of $\sqrt{\pi}$ as a root of $\sin(x^2)$ finishes in less than 20 seconds):

> build/examples/real_roots 2 1.0 2.0 -refine 1000000 -verbose interval: 1.5 +/- 0.5 maxdepth = 30, maxeval = 100000, maxfound = 100000, low_prec = 30 found isolated root in: 1.875 +/- 0.125 after bisection 1: 1.76953125 +/- 0.0039062 after bisection 2: 1.7723388671875 +/- 0.00012207 newton initial accuracy: 12 newton step: wp = 18 + 10 = 28 newton step: wp = 31 + 10 = 41 newton step: wp = 56 + 10 = 66 newton step: wp = 107 + 10 = 117 newton step: wp = 208 + 10 = 218 newton step: wp = 411 + 10 = 421 newton step: wp = 817 + 10 = 827 newton step: wp = 1628 + 10 = 1638 newton step: wp = 3250 + 10 = 3260 newton step: wp = 6494 + 10 = 6504 newton step: wp = 12982 + 10 = 12992 newton step: wp = 25958 + 10 = 25968 newton step: wp = 51911 + 10 = 51921 newton step: wp = 103816 + 10 = 103826 newton step: wp = 207626 + 10 = 207636 newton step: wp = 415247 + 10 = 415257 newton step: wp = 830489 + 10 = 830499 newton step: wp = 1660973 + 10 = 1660983 newton step: wp = 3321941 + 10 = 3321951 1.77245385[omitting huge string of digits]3729443 +/- 3.9829e-1000007 --------------------------------------------------------------- Found roots: 1 Subintervals possibly containing undetected roots: 0 Function evaluations: 46 cpu/wall(s): 18.3 18.322 virt/peak/res/peak(MB): 35.12 48.81 16.30 29.85