Efficiently get primes in range

Question

I first had an iterative method, which ran sufficiently fast:

n = 10^12;
primes = Prime@Range@PrimePi@Sqrt[n];
primesCbrt = Prime@Range@PrimePi@CubeRoot[n];
count = 0;
For[pi = 1, pi <= Length[primesCbrt], pi++,
  p = primesCbrt[[pi]];
  For[qi = pi + 1, primes[[qi]] <= Sqrt[n/p], qi++,
   count += PrimePi[n/(p*primes[[qi]])] - qi;
   ];
  ];
count

For n=10^3, count=135, for n=10^6, count=206964, and for n=10^12, count=190614467420.

It is much cleaner to rewrite it in summation notation, but now it won't finish at all:

Sum[PrimePi[n/(p q)] - PrimePi[q], {p, primesCbrt}, {q, Prime[Range[PrimePi[p] + 1, PrimePi[Sqrt[n/p]]]]}]

This is directly $$\sum_{p < \sqrt[3]{n}} \sum_{p < q < \sqrt{n/p}} \pi(\frac{n}{pq}) - \pi(q)$$

I think the slowest part is getting the range of primes with Prime[Range[PrimePi[x]]], which if I understand correctly, has to calculate PrimePi and Prime, which seems much slower than using a pre-generated list of primes and iterating through, taking those that fall in the range. How can I generate primes in a range without using iterative functions?

Answer 1

The following is about 50 times faster on my machine.

Block[{n = 10^9, CbrtNindex},
      CbrtNindex = PrimePi[CubeRoot[n]];
      Sum[
          PrimePi[n/(Prime[i]*Prime[j])] - j,
          {i, 1, CbrtNindex},
          {j, i + 1, PrimePi[Sqrt[n/Prime[i]]]}]
]

Answer 2

On my system (v10.1 under Windows) this is a modest improvement on Kenny Colnago's code. I added vectorization and use of machine precision Reals at an intermediate step rather than Rational values.

n = 1*^10;

Sum[
  Tr[ PrimePi[n/(N@Prime[i]*Prime[j])] - j ],
  {i, PrimePi @ CubeRoot @ n},
  {j, { Range[ i + 1, PrimePi @ Sqrt[n/N@Prime[i]] ] }}
] // RepeatedTiming

{1.1, 1997171674}

(I get {1.34, 1997171674} for Kenny Colnago's code.)