A306458 a(n) = A001694(n)/A007947(A001694(n)), the powerful numbers divided by their squarefree kernel.

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

Author

Amiram Eldar, Feb 17 2019

Keywords

• nonn
• look

Comments

A permutation of the positive integers.

Links

• Robert Israel, Table of n, a(n) for n = 1..10000

Crossrefs

Cf. A001694, A007947, A064549.

Programs

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

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

Formula

A064549(a(n)) = A001694(n).