A131605 Perfect powers of nonprimes (m^k where m is a nonprime positive integer and k >= 2).
1, 36, 100, 144, 196, 216, 225, 324, 400, 441, 484, 576, 676, 784, 900, 1000, 1089, 1156, 1225, 1296, 1444, 1521, 1600, 1728, 1764, 1936, 2025, 2116, 2304, 2500, 2601, 2704, 2744, 2916, 3025, 3136, 3249, 3364, 3375, 3600, 3844, 3969, 4225, 4356, 4624
Offset: 1
Keywords
Links
- Michael De Vlieger, Table of n, a(n) for n = 1..10000 (Terms a(n), n = 1..1323 from Klaus Brockhaus, and n = 1324..8649 from Daniel Forgues.)
Crossrefs
Programs
-
Mathematica
With[{nn = 2^20}, {1}~Join~Select[Union@ Flatten@ Table[a^2*b^3, {b, Surd[nn, 3]}, {a, Sqrt[nn/b^3]}], And[Length[#2] > 1, GCD @@ #2 > 1] & @@ {#, FactorInteger[#][[;; , -1]]} &] ] (* Michael De Vlieger, Aug 11 2025 *)
-
PARI
isok(n) = if (n == 1, return (1), return (ispower(n, ,&np) && (! isprime(np)))); \\ Michel Marcus, Jun 12 2013
-
Python
from sympy import mobius, integer_nthroot, primepi def A131605(n): def f(x): return int(n-2+x+sum(mobius(k)*((a:=integer_nthroot(x,k)[0])-1)+primepi(a) for k in range(2,x.bit_length()))) 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 14 2024
Comments