Re: guarantees on digits of floating point values

[Karim Belabas <Karim.Belabas@u-bordeaux.fr>]

* Ruud H.G. van Tol [2025-11-01 22:07]:
> 
> As a test, I checked 1000021 digits of exp(1),
> 
> my(N=10^6+21); localprec(N); exp(1) * 10^(N-1)\1
> 
> against https://apod.nasa.gov/htmltest/gifcity/e.1mil
> and all digits matched. :)
> 
> Is there a guarantee in the documentation about such?
> Similar to how all primes up to 2^64 are guaranteed.

We don't guarantee correct rounding (as per IEEE 754 floating point
standard) so the notion of "correct" digits is moot for us. We don't
have interval arithmetic either.

The general guarantee is that relative bit accuracy should be respected
(barring bugs) for an atomic function:
1) the input is assumed exact (with infinitely many trailing zeroes)
2) if x != 0 is the mathematically correct result, y is the returned result
and b is the bit accuracy
then we must have |x - y| / |x| < 2^(-b).
3) if x = 0 we should have |y| < 2^(-b).

For chains of computation (e.g., x + y + z or sqrt(cos(x))), compound
atomic errors.

A few functions (e.g., L functions or inverse Mellin transforms)
guarantee an *absolute* accuracy instead, so that |x - y| < 2^(-b); this
is stated in their documentation.

Numerical sums and integral *strive* for this but are heuristic in
nature since the summand/integrand is handled as a black box: no
theoretical information is available to the integration/summation
routine except evaluations at a finite number of points.

Cheers,

    K.B.
-- 
Pr. Karim Belabas, U. Bordeaux, Vice-président en charge du Numérique
Institut de Mathématiques de Bordeaux UMR 5251 - (+33) 05 40 00 29 77
http://www.math.u-bordeaux.fr/~kbelabas/