OFFSET
1,1
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
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
EXAMPLE
A006881 begins 6 10 14 15 ... so this sequence begins 36 100 196 225 ...
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 *)
PROG
(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
(Magma) [k^2:k in [1..150]| IsSquarefree(k) and #PrimeDivisors(k) eq 2]; // Marius A. Burtea, Oct 24 2019
(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
CROSSREFS
KEYWORD
easy,nonn
AUTHOR
Alford Arnold, Jul 06 2003
STATUS
approved