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 = 3177A071521 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.