Michael Day on Mon, 28 Nov 2022 00:58:01 +0100

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

Re: A071521(n)

(Resending, with slight edit, to include users' group - I suggest you (KB) offer a3 to the OEIS site!)

Beautiful! a3 is so elegant.  My effort at a J-equivalent,  doing the sum of a vector of
log-base-3 of n right-shifted by the vector 0,1,2...,logint(n,2) works fine unless n is so
large as to become a float!

It’s not so satisfyingly concise.  Here’s a candy code version;  the one-liner, not shown here, 
is rather ugly:

a3 =: monad define

shifts  =. i.1+2<.@^.y NB. vector of positions to shift

rshift  =. 33 b.~-     NB. verb to do the shift

logint3 =. 3 & (<.@^.) NB. J-verb equivalent of Pari GP fn

+/ 1 + logint3 y rshift shifts NB. do it


I recall that in Dyalog APL ,
   q <- a op <- b     
results in q = a, not q = a op b .


Sent from my iPad

On 27 Nov 2022, at 10:21, Karim Belabas <Karim.Belabas@math.u-bordeaux.fr> wrote:


 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.)



* 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
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.

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