login
A286708
Powerful numbers (A001694) that are not prime powers (A000961).
53
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,1
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.
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.
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
EXAMPLE
-------------------------------
| 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}, ...}.
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 *)
PROG
(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
KEYWORD
nonn
AUTHOR
Ilya Gutkovskiy, May 13 2017
STATUS
approved