A153158 a(n) = A007916(n)^2.

4, 9, 25, 36, 49, 100, 121, 144, 169, 196, 225, 289, 324, 361, 400, 441, 484, 529, 576, 676, 784, 841, 900, 961, 1089, 1156, 1225, 1369, 1444, 1521, 1600, 1681, 1764, 1849, 1936, 2025, 2116, 2209, 2304, 2500, 2601, 2704, 2809, 2916, 3025, 3136, 3249, 3364, 3481

Offset: 1

Author

Vladimir Joseph Stephan Orlovsky, Dec 19 2008

Keywords

• nonn

Comments

A378168(n) is the number of terms <= 10^n. - Chai Wah Wu, Nov 21 2024

Examples

			2^2 = 4, 3^2 = 9, 4^2 = 16 = 2^4 is not in the sequence, 5^2 = 25, 6^2 = 36, ...
		

Links

• Reinhard Zumkeller, Table of n, a(n) for n = 1..1000

Crossrefs

Cf. A007916, A153147, A153157, A153159, A153160, A062503, A076467, A378168, A378287.

Programs

Haskell

a153158 n = a153158_list !! (n-1)
a153158_list = filter ((== 2) . foldl1 gcd . a124010_row) [2..]
-- Reinhard Zumkeller, Apr 13 2012

Maple

q:= n-> is(igcd(seq(i[2], i=ifactors(n)[2]))=2):
select(q, [i^2$i=2..60])[];  # Alois P. Heinz, Nov 26 2024

Mathematica

Select[Range[2,100],GCD@@Last/@FactorInteger@#==1&]^2

Python

from sympy import mobius, integer_nthroot
def A153158(n):
    def f(x): return int(n+1-sum(mobius(k)*(integer_nthroot(x,k)[0]-1) for k in range(2,x.bit_length())))
    m, k = n, f(n)
    while m != k:
        m, k = k, f(k)
    return m**2 # Chai Wah Wu, Aug 13 2024

Formula

GCD(exponents in prime factorization of a(n)) = 2, cf. A124010. - Reinhard Zumkeller, Apr 13 2012

Sum_{n>=1} 1/a(n) = zeta(2) - 1 - Sum_{k>=2} mu(k)*(1 - zeta(2*k)) = 0.5444587396... - Amiram Eldar, Jul 02 2022

Intersection of A000290 and A378287. Squares that are not of the form m^k for some k>=3. - Chai Wah Wu, Nov 21 2024

Extensions

Edited by Ray Chandler, Dec 22 2008