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;

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?

  • $\begingroup$ Your first code does not work out of the box, quit the kernel, copy it and check yourself. Also, it would be nice to describe what is an input for desired procedure and expected output. $\endgroup$
    – Kuba
    Jun 9, 2017 at 8:08
  • $\begingroup$ Could you say in words what this code returns? $\endgroup$
    – Greg Hurst
    Jun 9, 2017 at 18:25
  • $\begingroup$ @Kuba I have added the necessary variables. Also there is no input other than changing n manually. $\endgroup$
    – qwr
    Jun 9, 2017 at 19:08
  • $\begingroup$ @ChipHurst all $pqr \le n$, with $p<q<r$ and $p,q,r$ prime $\endgroup$
    – qwr
    Jun 9, 2017 at 19:18

2 Answers 2


The following is about 50 times faster on my machine.

Block[{n = 10^9, CbrtNindex},
      CbrtNindex = PrimePi[CubeRoot[n]];
          PrimePi[n/(Prime[i]*Prime[j])] - j,
          {i, 1, CbrtNindex},
          {j, i + 1, PrimePi[Sqrt[n/Prime[i]]]}]
  • $\begingroup$ The code (alone) does not run for me, but replacing p with Prime[i] does work well. But, if I may ask, what is the purpose of Block here? n and CbrtNindex do not change throughout the sum. $\endgroup$
    – qwr
    Jun 9, 2017 at 19:02
  • $\begingroup$ The code runs about 3x faster than the iterative version on my machine $\endgroup$
    – qwr
    Jun 9, 2017 at 21:06
  • $\begingroup$ The Block ensures that values of n and CbrtNindex in the enclosing context are not altered by the assignments in this code. In short, it ensures that assignments to local variables do not clobber globals. $\endgroup$ Apr 15, 2021 at 17:51

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;

  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.)


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.