Bill Allombert on Thu, 06 Nov 2025 23:38:37 +0100


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: guarantees on digits of floating point values


On Sat, Nov 01, 2025 at 08:11:06PM +0100, Ruud H.G. van Tol wrote:
> 
> 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. :)

You can avoid floating points entirely and do the computation with
rational numbers. The formula PARI/GP uses is

/* binary spliiting */
exprat_rec(a,b)=
{
  if(a==b,[0,1],
     a==b-1,[1,b],
     my(c=(a+b)\2,[u,v]=exprat_rec(a,c),[u1,v1]=exprat_rec(c,b));[u*v1+u1,v*v1])
}
/* compute sum(i=0,n,1/i!) as a rational number */
exprat(n)=my([p,q]=exprat_rec(0,n));1+p/q

We have
sum(i=0,n,1/i!) < exp(1) < sum(i=0,n,1/i!) + 1/(n*n!)

so for n = 205022, B = 10^(10^6+1)

floor(B*(exprat(n))) < B * exp(1) < ceil(B*(exprat(n) + 1/(n*n!)))

? r=exprat(n);a=floor(B*r);b=ceil(B*(r+1/(n*n!)));
? b-a
%59 = 1
? a%10^20
%60 = 37981764476942281888
? b%10^20
%61 = 37981764476942281889

so to 10^6 digits, it should end with 3798176447694228188

Cheers,
Bill