%I #16 Sep 15 2024 22:01:19
%S 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,
%T 128,48,17,54,49,19,28,40,72,21,22,50,256,23,96,125,108,45,26,243,56,
%U 80,29,144,30,31,44,162,100,512,33,75,192,34,35,216,63,121,52
%N a(n) = A001694(n)/A007947(A001694(n)), the powerful numbers divided by their squarefree kernel.
%C A permutation of the positive integers.
%H Robert Israel, <a href="/A306458/b306458.txt">Table of n, a(n) for n = 1..10000</a>
%F A064549(a(n)) = A001694(n).
%p N:= 10^4: # to get terms corresponding to powerful numbers <= N
%p rad:= n -> convert(numtheory:-factorset(n), `*`):
%p S:= {1}:
%p p:= 1:
%p do
%p p:= nextprime(p);
%p if p^2 > N then break fi;
%p S:= S union map(t -> seq(t*p^i,i=2..floor(log[p](N/t))),select(`<=`,S,N/p^2));
%p od:
%p map(t -> t/rad(t), sort(convert(S,list))); # _Robert Israel_, Mar 20 2019
%t 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 *)
%o (PARI) apply(x->(x/factorback(factorint(x)[, 1])), select(x->ispowerful(x), vector(1600, k, k))) \\ _Michel Marcus_, Feb 17 2019
%o (Python)
%o from math import isqrt, prod
%o from sympy import mobius, integer_nthroot, primefactors
%o def A306458(n):
%o def squarefreepi(n): return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+1)))
%o def bisection(f,kmin=0,kmax=1):
%o while f(kmax) > kmax: kmax <<= 1
%o while kmax-kmin > 1:
%o kmid = kmax+kmin>>1
%o if f(kmid) <= kmid:
%o kmax = kmid
%o else:
%o kmin = kmid
%o return kmax
%o def f(x):
%o c, l = n+x-squarefreepi(integer_nthroot(x,3)[0]), 0
%o j = isqrt(x)
%o while j>1:
%o k2 = integer_nthroot(x//j**2,3)[0]+1
%o w = squarefreepi(k2-1)
%o c -= j*(w-l)
%o l, j = w, isqrt(x//k2**3)
%o return c+l
%o return (m:=bisection(f,n,n))//prod(primefactors(m)) # _Chai Wah Wu_, Sep 14 2024
%Y Cf. A001694, A007947, A064549.
%K nonn,look
%O 1,2
%A _Amiram Eldar_, Feb 17 2019