%I #15 Aug 30 2024 02:53:17

%S 8,12,20,27,28,44,45,52,63,68,76,92,99,116,117,124,125,148,153,164,

%T 171,172,175,188,207,212,236,244,261,268,275,279,284,292,316,325,332,

%U 333,343,356,369,387,388,404,412,423,425,428,436,452,475,477,508,524,531,539,548,549,556,575,596,603

%N Numbers of the form p^2*q where p and q are primes and p <= q.

%H Robert Israel, <a href="/A337877/b337877.txt">Table of n, a(n) for n = 1..10000</a>

%e a(3) = 20 is a term because 20=2^2*5 with 2 <= 5.

%p N:= 3000: # for terms <= N

%p P:= select(isprime, [2,seq(i,i=3..N/2,2)]): nP:= nops(P):

%p R:= NULL:

%p for i from 1 to nP do

%p p2:= P[i]^2;

%p for j from i to nP do

%p x:= p2*P[j];

%p if x > N then break fi;

%p R:= R, x

%p od od:

%p sort([R]);

%o (Python)

%o from sympy import primepi, primerange, integer_nthroot

%o def A337877(n):

%o def f(x): return int(n+x-sum(primepi(x//k**2)-a for a,k in enumerate(primerange(integer_nthroot(x,3)[0]+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) # _Chai Wah Wu_, Aug 29 2024

%Y Contained in A337806.

%K nonn

%O 1,1

%A _Robert Israel_, Sep 27 2020