Re: A071521(n)

[Karim Belabas <Karim.Belabas@math.u-bordeaux.fr>]

Hi,

  I didn't check the definition of A071521 in OEIS but I guess the
problem is the unexpected (if one is used to C) semantic of the 'a op= b'
operator in PARI/GP, whose value is A = a op b (and not 'a' !). I.e. in 

  A071521(n)= my(t=1);sum(k=1,logint(n,3),logint(n\t*=3,2)+1)

when k = 1, we set t = 3 and compute logint(n \ 3, 2), where it was probably 
expected to compute logint(n \ 1, 2). If so, a corrected version would be

  A071521bis(n) = my(t=1/3); sum(k=1,logint(n,3),logint(n\t*=3,2)+1) + 1;

(the final + 1 is for k = 0). We then have

(11:15) gp > A071521bis(2000000000000000000000000000000)
%1 = 3278

While I'm at it, here's a simpler and slightly faster version of a2():

  a3(n) = sum(k=0, logint(n, 2), logint(n >> k, 3) + 1)

(11:15) gp > a3(2000000000000000000000000000000)
%2 = 3278

A071521bis remains a little faster of course, so this is of marginal interest.
(Good for debugging, maybe.)

Cheers,

      K.B.

* Michael Day [2022-11-27 10:49]:
> I mostly use Pari GP when I need better precision than I can get with J,  so
> I'm no expert!
> 
> I had a look at this to see if I could get a nice answer in J,  but noticed
> that
> Charles' function has a pleasing variant which replaces the multiply by 3
> with a doubling by addition,  so I've called it a2:
> 
> a2(n)=my(t=1/2); sum(k=0, logint(n, 2), t+=t; logint(n\t, 3)+1)
> 
> So instead of a loop on powers of 3, we've got rather more loops on
> powers of 2.
> 
> Unfortunately, on the way,  I found what seems to be a problem with your
> suggestion (Ruud van Tol);  perhaps there's a typo?
> 
> (18:04) gp > a(2000000000000000000000000000000)
> %28 = 3278
> (18:04) gp > a2(2000000000000000000000000000000)
> %29 = 3278
> (18:04) gp > A071521(2000000000000000000000000000000)
> %30 = 3177
> 
> A071521 does seem a tad faster than a,  while a2 is a little slower than a, 
> as to be expected.
> 
> (I'm relieved to say that my J function,  n3sm,  works ok for this argument:
>    n3sm 2000000000000000000000000000000
> 3278
> )
> 
> Cheers,
> 
> Mike
> 
> On 26/11/2022 10:29, Ruud H.G. van Tol wrote:
> > 
> > "A071521 Number of 3-smooth numbers <= n"
> > https://oeis.org/A071521
> > 
> > The (to me very educational) PARI-code
> > from Charles, on that page, looks to me
> > like a variant of a formula shown in
> > the "Raphael Schumacher" link.
> > 
> > I would now probably write it as:
> > A071521(n)= my(t=1);sum(k=1,logint(n,3),logint(n\t*=3,2)+1)
> > 
> > (just trying to shave off another millisecond,
> > though I was blown away by its speed already:
> > a 3^k variant takes 40x as long)
> > 
> > -- Ruud
> > 
> > P.S. As n could be any ridiculous integer,
> > a "powers()" approach isn't applicable.
> > A "Power(b,e)=exp(e*log(b))" approach might be feasible,
> > but that probably just takes longer,
> > and will add precision treats.
> > 
> 
> 
> 

    K.B.
--
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/
`