Re: A071521(n)

[Michael Day <mike_liz.day@tiscali.co.uk>]

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.