Re: log_int_rat

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Mon, Jan 22, 2024 at 11:56:29AM +0100, Ruud H.G. van Tol wrote:
> 
> On 2024-01-22 11:25, Ruud H.G. van Tol wrote:
> > [...]
> > On the non-discrete code branch:
> > 
> >   a(n) = localbitprec(logint(n+2,2)); 1 + n + floor(n*log(2)/log(3/2));
> > 
> > with test:
> > 
> >   {default(parisizemax,2^30); my(v=alist(2^20)); for(n=1,#v,
> > (v[n]!=a(n-1)) && print(n);break);}
> > 
> > I received this report:
> > 
> > > that code wrong result a(83690) = 226760 (in a 32-bit pari)
> > > with no hint that it didn't know right from wrong.
> > 
> > 
> > ? my(v=alist(83700)); v[83691]
> > cpu time = 485 ms, real time = 490 ms.
> > % 226759
> 
> My current guess is that it needs something more like
>   localbitprec(2*n+1); ...

No, the issue run deeper than that.

? 83690*log(2)/log(3/2)
%19 = 143068.9999732032502851373630

To get the correct floor() you need to 'skip' the four 9 digits.

Now consider this other example
? \p200
   realprecision = 202 significant digits (200 digits displayed)
? 282896829528374*log(2)/log(3/2)
%22 = 483615324366282.99999999999999972480016707721405476228331340155888020025561304044368677831069845019580206096595446347432023644929708447688723054552955187462035282662255837003805169825570885964555199

Here you have to skip fifteen 9 digits.
So unless you have a formula for the maximal number of 9 or zeros that can
appear, you cannot get a sufficient precision.

Instead I offer this (written while waiting for the water to boil, so be careful)

safedigits(f)=
{
  my(z=f(),e=exponent(z),k=3,r);
  while(1,
    my(K=10^k,y=round(K*z,&r));
    if (r<0 && (y+1)%K>1, return(y\K));
    k*=2;
    localbitprec(e+4*k);z=f());
}
a(n) = localbitprec(logint(n+2,2)); 1+n+safedigits(()->(n*log(2)/log(3/2)));

Cheers,
Bill