| 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