A376720 Product of numbers m that are neither squarefree nor prime powers and rad(m), where rad = A007947.

72, 108, 200, 144, 392, 216, 400, 968, 675, 288, 500, 1352, 324, 784, 1800, 1323, 2312, 432, 1125, 2888, 800, 3528, 1936, 2700, 4232, 576, 1372, 3267, 1000, 2704, 648, 1568, 6728, 4563, 3600, 7688, 5292, 8712, 2025, 4624, 9800, 864, 3087, 10952, 4500, 5776, 7803

Offset: 1

Author

Michael De Vlieger, Oct 05 2024

Keywords

• nonn
• easy

Comments

Term a(n) = k is powerful but not a prime power (i.e., in A286708) such that k/rad(k) is not squarefree, where rad = A007947 and k/rad(k) = A003557(k).

Permutation of A372404.

Examples

			Let b(n) = A126706(n).
Table of b(n) and a(n) for n <= 12:
   n   b(n)                a(n)
  -----------------------------------------
   1    12 = 2^2 * 3        72 = 2^3 * 3^2
   2    18 = 2 * 3^2       108 = 2^2 * 3^3
   3    20 = 2^2 * 5       200 = 2^3 * 5^2
   4    24 = 2^3 * 3       144 = 2^4 * 3^2
   5    28 = 2^2 * 7       392 = 2^3 * 7^2
   6    36 = 2^2 * 3^2     216 = 2^3 * 3^3
   7    40 = 2^3 * 5       400 = 2^4 * 5^2
   8    44 = 2^2 * 11      968 = 2^3 * 11^2
   9    45 = 3^2 * 5       675 = 3^3 * 5^2
  10    48 = 2^4 * 3       288 = 2^5 * 3^2
  11    50 = 2 * 5^2       500 = 2^2 * 5^3
  12    52 = 2^2 * 13     1352 = 2^3 * 13^2
		

Links

• Michael De Vlieger, Table of n, a(n) for n = 1..65536
• Michael De Vlieger, Log log scatterplot of a(n), n = 1..2^16.

Crossrefs

Cf. A003557, A007947, A126706, A286708, A372404.

Programs

Mathematica

Map[#*Times @@ FactorInteger[#][[All, 1]] &, Select[Range[12, 160], Nor[PrimePowerQ[#], SquareFreeQ[#]] &] ]

Python

from math import prod, isqrt
from sympy import primepi, integer_nthroot, mobius, primefactors
def A376720(n):
    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): return int(n+sum(primepi(integer_nthroot(x,k)[0]) for k in range(2,x.bit_length()))+sum(mobius(k)*(x//k**2) for k in range(1, isqrt(x)+1)))
    return (m:=bisection(f,n,n))*prod(primefactors(m)) # Chai Wah Wu, Oct 05 2024

Formula

a(n) = m * rad(m) for m in A126706.