A085986 Squares of the squarefree semiprimes (p^2*q^2).

36, 100, 196, 225, 441, 484, 676, 1089, 1156, 1225, 1444, 1521, 2116, 2601, 3025, 3249, 3364, 3844, 4225, 4761, 5476, 5929, 6724, 7225, 7396, 7569, 8281, 8649, 8836, 9025, 11236, 12321, 13225, 13924, 14161, 14884, 15129, 16641, 17689, 17956, 19881

Offset: 1

Author

Alford Arnold, Jul 06 2003

Keywords

• easy
• nonn

Comments

This sequence is a member of a family of sequences directly related to A025487. First terms and known sequences are listed below: 1, A000007; 2, A000040; 4, A001248; 6, A006881; 8, A030078; 12, A054753; 16, A030514; 24, A065036; 30, A007304; 32, A050997; 36, this sequence; 48, ?; 60, ?; 64, ?; ....

Subsequence of A077448. The numbers in A077448 but not in here are 1, the squares of A046386, the squares of A067885, etc. - R. J. Mathar, Sep 12 2008

a(4)-a(3)=29 and a(3)+a(4)=421 are both prime. There are no other cases where the sum and difference of two members of this sequence are both prime. - Robert Israel and J. M. Bergot, Oct 25 2019

Examples

			A006881 begins 6 10 14 15 ... so this sequence begins 36 100 196 225 ...
		

Links

• T. D. Noe, Table of n, a(n) for n = 1..1000
• Index to sequences related to prime signature

Crossrefs

Subsequence of A036785 and of A077448.

Subsequence of A062503.

Cf. A025487.

Cf. A085548, A085964.

Programs

Magma

[k^2:k in [1..150]| IsSquarefree(k) and #PrimeDivisors(k) eq 2]; // Marius A. Burtea, Oct 24 2019

Mathematica

f[n_]:=Sort[Last/@FactorInteger[n]]=={2,2}; Select[Range[20000], f] (* Vladimir Joseph Stephan Orlovsky, Aug 14 2009 *)
Select[Range[200],PrimeOmega[#]==2&&SquareFreeQ[#]&]^2 (* Harvey P. Dale, Mar 07 2013 *)

PARI

list(lim)=my(v=List(), x=sqrtint(lim\=1), t); forprime(p=2, x\2, t=p; forprime(q=2, min(x\t,p-1), listput(v, (t*q)^2))); Set(v) \\ Charles R Greathouse IV, Sep 22 2015

PARI

is(n)=factor(n)[,2]==[2,2]~ \\ Charles R Greathouse IV, Oct 19 2015

Python

from math import isqrt
from sympy import primepi, primerange
def A085986(n):
    def f(x): return int(n+x+(t:=primepi(s:=isqrt(x)))+(t*(t-1)>>1)-sum(primepi(x//k) for k in primerange(1, s+1)))
    m, k = n, f(n)
    while m != k:
        m, k = k, f(k)
    return m**2 # Chai Wah Wu, Aug 18 2024

Formula

a(n) = A006881(n)^2.

Sum_{n>=1} 1/a(n) = (P(2)^2 - P(4))/2 = (A085548^2 - A085964)/2 = 0.063767..., where P is the prime zeta function. - Amiram Eldar, Jul 06 2020