(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 .
Mike
Sent from my iPad
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/
`
|