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 '17 at 8:08
  • $\begingroup$ Could you say in words what this code returns? $\endgroup$ – Chip Hurst Jun 9 '17 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 '17 at 19:08
  • $\begingroup$ @ChipHurst all $pqr \le n$, with $p<q<r$ and $p,q,r$ prime $\endgroup$ – qwr Jun 9 '17 at 19:18

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]]]}]
| improve this answer | |
  • $\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 '17 at 19:02
  • $\begingroup$ The code runs about 3x faster than the iterative version on my machine $\endgroup$ – qwr Jun 9 '17 at 21:06

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

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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