A286708 Powerful numbers (A001694) that are not prime powers (A000961).

36, 72, 100, 108, 144, 196, 200, 216, 225, 288, 324, 392, 400, 432, 441, 484, 500, 576, 648, 675, 676, 784, 800, 864, 900, 968, 972, 1000, 1089, 1125, 1152, 1156, 1225, 1296, 1323, 1352, 1372, 1444, 1521, 1568, 1600, 1728, 1764, 1800, 1936, 1944, 2000, 2025, 2116, 2304, 2312, 2500, 2592, 2601, 2700, 2704, 2744

Offset: 1

Author

Ilya Gutkovskiy, May 13 2017

Keywords

• nonn

Comments

If a prime p divides a(n) then p^2 must also divide a(n) and number of distinct primes dividing a(n) > 1.

Intersection of A001694 and A024619.

Examples

			-------------------------------
| n | a(n) | prime            |
|   |      | factorization    |
|------------------------------
| 1 | 36   | {{2, 2}, {3, 2}} |
| 2 | 72   | {{2, 3}, {3, 2}} |
| 3 | 100  | {{2, 2}, {5, 2}} |
| 4 | 108  | {{2, 2}, {3, 3}} |
| 5 | 144  | {{2, 4}, {3, 2}} |
| 6 | 196  | {{2, 2}, {7, 2}} |
| 7 | 200  | {{2, 3}, {5, 2}} |
| 8 | 216  | {{2, 3}, {3, 3}} |
| 9 | 225  | {{3, 2}, {5, 2}} |
-------------------------------
a(n) = p_1^e_1*p_2^e_2*... : {{p_1, e_1}, {p_2, e_2}, ...}.
		

Links

• Michael De Vlieger, Table of n, a(n) for n = 1..10000 (first 5997 terms from Robert Israel)
• Eric Weisstein's World of Mathematics, Prime Power.
• Eric Weisstein's World of Mathematics, Powerful Number.
• Index entries for sequences related to powerful numbers

Crossrefs

Cf. A000961, A001221, A001694, A024619, A082695, A136141, A131605.

Programs

Maple

N:= 10000:
S:= {1}: P:= {1}:
p:= 1:
do
  p:= nextprime(p);
  if p^2 > N then break fi;
  S:= map(s -> (s, seq(s*p^k, k = 2 .. floor(log[p](N/s)))), S);
  P:= P union {seq(p^k, k=2..floor(log[p](N)))}:
od:
sort(convert(S minus P, list)); # Robert Israel, May 14 2017

Mathematica

Select[Range@2750, Min@FactorInteger[#][[All, 2]] > 1 && ! PrimePowerQ[#] &]
(* Second program *)
nn = 2^25; Select[Rest@ Union@ Flatten@ Table[a^2*b^3, {b, nn^(1/3)}, {a, Sqrt[nn/b^3]}], ! PrimePowerQ[#] &] (* Michael De Vlieger, Jun 22 2022 *)

Python

from sympy import primefactors, factorint
print([n for n in range(4,2745) if len(primefactors(n)) > 1 and min(list(factorint(n).values())) > 1]) # Karl-Heinz Hofmann, Feb 07 2023

Python

from math import isqrt
from sympy import integer_nthroot, primepi, mobius
def A286708(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, 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)
        c -= squarefreepi(integer_nthroot(x,3)[0])-l
        return c+1+sum(primepi(integer_nthroot(x, k)[0]) for k in range(2, x.bit_length()))
    return bisection(f,n,n) # Chai Wah Wu, Sep 10 2024

Formula

Sum_{n>=1} 1/a(n) = zeta(2)*zeta(3)/zeta(6) - Sum_{p prime} 1/(p*(p-1)) - 1 = A082695 - A136141 - 1 = 0.17043976777096407719... - Amiram Eldar, Feb 12 2021