Michael Day on Sun, 27 Nov 2022 10:50:48 +0100

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

Re: A071521(n)

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



On 26/11/2022 10:29, Ruud H.G. van Tol wrote:

"A071521 Number of 3-smooth numbers <= n"

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.