A303606 Powers of composite squarefree numbers that are not squarefree.

36, 100, 196, 216, 225, 441, 484, 676, 900, 1000, 1089, 1156, 1225, 1296, 1444, 1521, 1764, 2116, 2601, 2744, 3025, 3249, 3364, 3375, 3844, 4225, 4356, 4761, 4900, 5476, 5929, 6084, 6724, 7225, 7396, 7569, 7776, 8281, 8649, 8836, 9025, 9261, 10000, 10404, 10648, 11025, 11236

Offset

1,1

Links

Amiram Eldar, Table of n, a(n) for n = 1..10000

Eric Weisstein's World of Mathematics, Prime Power.

Eric Weisstein's World of Mathematics, Squarefree.

Formula

Sum_{n>=1} 1/a(n) = Sum_{n>=1} 1/((A120944(n)-1)*A120944(n)) = Sum_{k>=2} (zeta(k)/zeta(2*k) - P(k) - 1) = 0.07547719891508850482..., where P(k) is the prime zeta function. - Amiram Eldar, Feb 12 2021

Example

196 is in the sequence because 196 = 2^2*7^2.
4900 is in the sequence because 4900 = 2^2*5^2*7^2.

Maple

isA303606 := proc(n) local f, k; f := ifactors(n)[2];
seq(f[k][2], k = 1..nops(f)); nops({%}) = 1 and nops(f) > 1 and %[1] > 1 end:
select(isA303606, [seq(1..12000)]);  # Peter Luschny, Sep 15 2025

Mathematica

Select[Range[12000], Length[Union[FactorInteger[#][[All, 2]]]] == 1 && ! SquareFreeQ[#] && ! PrimePowerQ[#] &]
seq[max_] := Module[{sp = Select[Range[Floor@Sqrt[max]], SquareFreeQ[#] && PrimeNu[#] > 1 &], s = {}}, Do[s = Join[s, sp[[k]]^Range[2, Floor@Log[sp[[k]], max]]], {k, 1, Length[sp]}]; Union@s]; seq[10^4] (* Amiram Eldar, Feb 12 2021 *)

Prog

(Python)
from math import isqrt
from sympy import mobius, primepi, integer_nthroot
def A303606(n):
    def g(x): return int(sum(mobius(k)*(x//k**2) for k in range(1, isqrt(x)+1))-primepi(x))
    def f(x): return n-3+x+(y:=x.bit_length())-sum(g(integer_nthroot(x, k)[0]) for k in range(2, y))
    kmin, kmax = 1, 2
    while f(kmax) >= kmax:
        kmax <<= 1
    while True:
        kmid = kmax+kmin>>1
        if f(kmid) < kmid:
            kmax = kmid
        else:
            kmin = kmid
        if kmax-kmin <= 1:
            break
    return kmax # Chai Wah Wu, Aug 19 2024
(SageMath)
def is_A303606(n: int) -> int:
    if n == 1: return False
    factors = list(factor(n))
    exponents = [f[1] for f in factors]
    return 1 != len(factors) and 1 == len(set(exponents)) and 1 != max(exponents)
print([n for n in range(1, 12000) if is_A303606(n)])  # Peter Luschny, Sep 19 2025

Crossrefs

Intersection of A024619 and A072777.

Intersection of A072774 and A126706.

Intersection of A013929 and A182853.

Cf. A000469, A001597, A005117, A120944, A286708.

Sequence in context: A375144 A395414 A340017 * A390126 A303661 A390956

Adjacent sequences: A303603 A303604 A303605 * A303607 A303608 A303609

Keyword

nonn

Author

Ilya Gutkovskiy, Apr 26 2018

Status

approved