login
A306458
a(n) = A001694(n)/A007947(A001694(n)), the powerful numbers divided by their squarefree kernel.
4
1, 2, 4, 3, 8, 5, 9, 16, 6, 7, 32, 12, 27, 10, 18, 11, 25, 64, 24, 13, 14, 20, 36, 15, 81, 128, 48, 17, 54, 49, 19, 28, 40, 72, 21, 22, 50, 256, 23, 96, 125, 108, 45, 26, 243, 56, 80, 29, 144, 30, 31, 44, 162, 100, 512, 33, 75, 192, 34, 35, 216, 63, 121, 52
OFFSET
1,2
COMMENTS
A permutation of the positive integers.
LINKS
FORMULA
A064549(a(n)) = A001694(n).
MAPLE
N:= 10^4: # to get terms corresponding to powerful numbers <= N
rad:= n -> convert(numtheory:-factorset(n), `*`):
S:= {1}:
p:= 1:
do
p:= nextprime(p);
if p^2 > N then break fi;
S:= S union map(t -> seq(t*p^i, i=2..floor(log[p](N/t))), select(`<=`, S, N/p^2));
od:
map(t -> t/rad(t), sort(convert(S, list))); # Robert Israel, Mar 20 2019
MATHEMATICA
p=Join[{1}, Select[ Range@ 12500, Min@ FactorInteger[#][[All, 2]] > 1 &]]; rad[n_] := Times @@ (First@# & /@ FactorInteger@ n); p/(rad/@p) (* after Harvey P. Dale at A001694 and Robert G. Wilson v at A007947 *)
PROG
(PARI) apply(x->(x/factorback(factorint(x)[, 1])), select(x->ispowerful(x), vector(1600, k, k))) \\ Michel Marcus, Feb 17 2019
(Python)
from math import isqrt, prod
from sympy import mobius, integer_nthroot, primefactors
def A306458(n):
def squarefreepi(n): return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+1)))
def bisection(f, kmin=0, kmax=1):
while f(kmax) > kmax: kmax <<= 1
while kmax-kmin > 1:
kmid = kmax+kmin>>1
if f(kmid) <= kmid:
kmax = kmid
else:
kmin = kmid
return kmax
def f(x):
c, l = n+x-squarefreepi(integer_nthroot(x, 3)[0]), 0
j = isqrt(x)
while j>1:
k2 = integer_nthroot(x//j**2, 3)[0]+1
w = squarefreepi(k2-1)
c -= j*(w-l)
l, j = w, isqrt(x//k2**3)
return c+l
return (m:=bisection(f, n, n))//prod(primefactors(m)) # Chai Wah Wu, Sep 14 2024
CROSSREFS
KEYWORD
nonn,look
AUTHOR
Amiram Eldar, Feb 17 2019
STATUS
approved