%I #26 Aug 30 2024 18:03:33
%S 900,1764,4356,4900,6084,10404,11025,12100,12996,16900,19044,23716,
%T 27225,28900,30276,33124,34596,36100,38025,49284,52900,53361,56644,
%U 60516,65025,66564,70756,74529,79524,81225,81796,84100,96100,101124
%N a(n) = A007304(n)^2.
%C Numbers that are the product of exactly 3 distinct squares of primes (p^2*q^2*r^2).
%H T. D. Noe, <a href="/A162143/b162143.txt">Table of n, a(n) for n = 1..1000</a>
%H <a href="/index/Pri#prime_signature">Index to sequences related to prime signature</a>
%F A050326(a(n)) = 8. - _Reinhard Zumkeller_, May 03 2013
%F Sum_{n>=1} 1/a(n) = (P(2)^3 + 2*P(6) - 3*P(2)*P(4))/6 = (A085548^3 + 2*A085966 - 3*A085548*A085964)/6 = 0.0036962441..., where P is the prime zeta function. - _Amiram Eldar_, Oct 30 2020
%e 900 = 2^2*3^2*5^2, 1764 = 2^2*3^2*7^2, 4356 = 2^2*3^2*11^2, ..
%t fQ[n_]:=Last/@FactorInteger[n]=={2,2,2}; Select[Range[100000], f]
%o (Python)
%o from math import isqrt
%o from sympy import primepi, primerange, integer_nthroot
%o def A162143(n):
%o def f(x): return int(n+x-sum(primepi(x//(k*m))-b for a,k in enumerate(primerange(integer_nthroot(x,3)[0]+1),1) for b,m in enumerate(primerange(k+1,isqrt(x//k)+1),a+1)))
%o def bisection(f,kmin=0,kmax=1):
%o while f(kmax) > kmax: kmax <<= 1
%o while kmax-kmin > 1:
%o kmid = kmax+kmin>>1
%o if f(kmid) <= kmid:
%o kmax = kmid
%o else:
%o kmin = kmid
%o return kmax
%o return bisection(f)**2 # _Chai Wah Wu_, Aug 29 2024
%Y Cf. A006881, A050326, A162142.
%Y Cf. A085548, A085964, A085966.
%K nonn
%O 1,1
%A _Vladimir Joseph Stephan Orlovsky_, Jun 25 2009
%E Edited by _N. J. A. Sloane_, Jun 27 2009