Arb: Efficient Arbitrary-Precision Midpoint-Radius Interval Arithmetic (2016)
arxiv.org
arxiv.org
Some improvements not covered in the 2016 preprint include much faster code for arithmetic and matrix multiplication (https://arxiv.org/abs/1901.04289, https://fredrikj.net/blog/2018/07/high-precision-linear-alge...), and code for rigorous arbitrary-precision numerical integration (https://arxiv.org/abs/1802.07942, http://fredrikj.net/blog/2017/11/new-rigorous-numerical-inte...).
<link rel="alternate" type="application/rss+xml"
title="RSS" href="http://www.example.net/feed.rss" />
<atom:link href="https://fredrikj.net/rss.xml" rel="self" type="application/rss+xml"/>A good technical introduction to ball arithmetic is this paper by Joris van der Hoeven: https://www.texmacs.org/joris/ball/ball.html
So that sent me chasing interval arithmetic's history -- https://en.wikipedia.org/wiki/Interval_arithmetic -- and that page is also silent as to where did it come from. But academic.ru has the (surprising to me) answer that none other than Archimedes (not a surprise) used it:
"Interval arithmetic is not a completely new phenomenon in mathematics; it has appeared several times under different names in the course of history.
For example Archimedes calculated lower and upper bounds 223/71 < π < 22/7 in the 3rd century BC. Actual calculations with intervals has neither been as popular as other numerical techniques, nor been completely forgotten."
d(x, y) = |sum_i((x_i - y_i)^n)|^(1/n)
Of course in one dimension, all of these Lp metrics are identical.
I’m amused by your subdivision of space note, given Figure 1 :). There are clearly some problems and transforms that are better on N-spheres (N-balls?) and others on rectangles / N-cubes. Do you have a deeper intuition for which? (The x^2 example was simple and cute).
Years ago [1], we applied interval arithmetic to tracing groups of rays (and compared to geometric bounding via planes/frustums). I’d be curious to think through an equivalent with your midpoint ball arithmetic, but I feel like it would “need” to be parameterized as a ball of origins (easy) and then something else for the cone of directions —- maybe theta/phi clusters or “cluster of points on the unit sphere” (but converting back and forth is more expensive than the gains).
"No general way exists to predict how many extra digits will have to be carried to compute a transcendental expression and round it correctly to some preassigned number of digits. Even the fact (if true) that a finite number of extra digits will ultimately suffice may be a deep theorem."
My guess is they aim for convergence and it works well in practice but it's still unknown if that's guaranteed to be correct... is that right?
k * x * (1-x)
It is possible to do grid scans of parameters and control the interval such thst you know you used enough digits that you can trust the result.
The table-makers dilemma only applies if you want a specific number of decimal places, correctly rounded. The midpoint-radius representation gives you a value and an upper bound for the distance to the true value.
You may find that the radius is too big to get the number of correctly rounded decimal places you want, but that doesn't make the result incorrect, just vague.
That's true, but if you can do this, then you can also increase the precision until the radius shrinks to something sharp, right? Which means you should be able to round correctly?
Of course if you're not a table maker trying to make a fixed-width table with correctly rounded values, you can just increase your precision to 10⁻¹⁰⁰ instead of worrying about the rounding error of a 10-digit representation.
Somehow that's not how I've been interpreting the statement until now. Rather, the way I read it, we don't know how to predict the amount of work it takes to guarantee correct rounding to any number of digits (even if we allow for extra digits beyond those we really want)—not merely to a fixed number of digits. Or to put it another way, it sounded as if for pathological inputs even guaranteeing >= N digits of accuracy isn't possible with a predictable amount of effort.
What led me to believe this is a couple of different statements from Kahan [1]:
A desideratum generally acknowledged but not always achieved is: "In error by rather less than an ulp if by more than half."
Reputable math libraries compute elementary transcendental functions mostly within slightly more than half an ulp and almost always well within one ulp.
Both of these suggested to me that errors can occasionally exceed half an ulp by a non-negligible amount. So, if we were rounding to the nearest integer, this would imply occasionally not only rounding arbitrary overshoots like 0.5 + 1E-100 to zero (whose difficulty you're pointing out), but in fact even occasionally even producing 0 when the output is 0.7. That's far worse. Is this interpretation wrong? If it is, then why should errors sometimes not be "well within one ulp"?
Yes, but in many cases that's not desirable. E.g. if you want to compute things in hardware where adding an extra digit means adding more transistors to your chip. (And then why not always use that extra digit instead of rounding to one digit less?)
> mostly within slightly more than half an ulp
The reason this doesn't say "at most half an ulp" but allows for "slightly more" is the table-makers dilemma. (If you keep reading after the part you quoted, you'll see that he's talking about precisely the case where a value is close to x.5)
> almost always well within one ulp
This is a different case, concerning functions that are expensive to compute to a high degree of precision (e.g. due to high sensitivity to rounding of intermediate results) but "reputable" libraries will still try really hard to give good results.
In terms of the midpoint-radius representation, the first is about not knowing how much smaller the radius needs to be to achieve correct rounding to a predetermined number of digits, whereas the second is about needing to wait for a long (but usually predictable) time before the radius drops below a predetermined target.
I guess you refer only to the computation of transcendentals sin, exp, etc. (because even in native floating-point formats, most other individual operations are guaranteed to be rounded correctly).
I haven't checked their docs carefully enough, but if it is like MPFR, the way they do it is by trading time guarantee for precision guarantee. So they can't predict how long each operation will take (let alone guarantee that it will go fast), so that they can guarantee that they meet your required precision.
This is why Arb has a lot more functions than MPFR: it's extremely difficult to implement correctly rounded functions as required by MPFR; Arb's contract (returning an error bound) is far easier to satisfy and works just as well 99% of the time.
From a software development point of view, the problem with correct rounding is that it doesn't compose: given correctly rounded functions f and g and correctly rounded input x, f(g(x)) is not (in general) correctly rounded. In interval arithmetic, the corresponding composition law (the inclusion property) does hold; this means that you can compose arbitrarily complicated functions without thinking.
It's unfortunate that this isn't something that was standardized and more widely available a long time ago. By contrast, IEEE binary128 and binary256 are practically useless due to lack of hardware support.
I agree this package is not new conceptually, but it has a convenient metaprogramming to generate code for arithmetic for any [hardware ieee number type] x N.
It also lists the complexity as a function of N, most operations scale as N^3, so yes, for larger N it might fail to be fast.
LGPL
This is from TAOCP 4A part 1:
> Dudeney's Digital Century puzzle.) There are many curious ways to obtain the number 100 by inserting arithmetical operators and possibly also parentheses into the sequence 123456789. For example, 100 = 1 + 2 x 3 + 4 x 5 - 6 + 7 + 8 x 9 = (1 + 2 - 3 - 4 ) x (5-6-7-8-9) = ((1/((2 + 3)/4 - 5 + 6)) x 7 + 8) x 9.
Of course, instead of 100 we can put something else, add exponentiation, decimal point and an operator that reinterprets the number as if written in base X (12345_6789 = 2124962081444303 in base 10).
I also like how Knuth adds the number of possible trees in a grammar that generates these expressions without redundant parentheses :D