Re: log_int_rat

["Ruud H.G. van Tol" <rvtol@isolution.nl>]

On 2024-01-22 13:22, Bill Allombert wrote:
> 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)));

To me this is just another example of why to avoid floats and floor() 
and such,
when the goal is integers.

So I'm back at my original quest. :)
> ? a054414(n) = 1 + n + floor( n * log(2) / log(3/2) );
> Is there a cleaner way, similar to logint, to do that floor-expression?

The table() function (from Bill) made me think that a discrete a(n),
ideally without any slow while-loop, is doable,
so I will keep looking and trying in that direction.

-- Thanks, Ruud