OFFSET
1,2
COMMENTS
A permutation of the positive integers.
LINKS
Robert Israel, Table of n, a(n) for n = 1..10000
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