Karim Belabas on Wed, 29 Nov 2023 21:30:33 +0100


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

Re: log_int_rat 

* Bill Allombert [2023-11-29 21:18]:
> On Wed, Nov 29, 2023 at 07:44:11PM +0100, Ruud H.G. van Tol wrote:
> > 
> > ? a054414(n) = 1 + n + floor( n * log(2) / log(3/2) );
> > 
> > ? [ a054414(n) |n<-[0..20]]
> > % [1, 3, 6, 9, 11, 14, 17, 19, 22, 25, 28, 30, 33, 36, 38, 41, 44, 47, 49,
> > 52, 55]
> > 
> > Is there a cleaner way, similar to logint, to do that floor-expression?
> 
> Using logint(2^n,3/2) would be much slower even if it worked, because that would
> force you to compute 2^n.
> 
> Really, the only defect of your function is that you do not take the accuracy
> into account so
> 
> ? a054414(2^1000)
>   *** floor: precision too low in truncr (precision loss in truncation).
> 
> but this is easily fixed:
> 
> a054414(n) = localbitprec(logint(n+2,2));1 + n + floor( n * log(2) / log(3/2) );
> 
> ? a054414(2^1000)
> %8 = 29032646699514618648207963388262277036465203935797085941659905401800597040332207967488444244921882872990640968909754332584729580377856680178651318617855116668095773049095433042501905586769687829752260479426676858529835781762375806913040422676074996080564646620120844538289072080537948895898995900461439
>   ***   last result computed in 0 ms.

Another defect is that you're computing log(2)/log(3/2) over and over again.

? for (n=1,10^6,a054414(n))
time = 3,265 ms.

? N = 10^6; /* max index */
? localbitprec(logint(N,2)); L = log(2)/log(3/2);
? a054414(n) = 1 + n + floor(n * L);
? for (n=1, N, a054414(n))
time = 203 ms.

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/