Efficient way to sum all the primes below $N$ million in Mathematica

Question

Today, I saw a problem described as below:

Sum all primes below $N$ million

At the beginning, I feel it is rather simple to calculate it in Mathematica enviroment. unfortuately, I failed >__< (Here, I assume $N=10$)

Trail 1

Total@(Prime /@ Range[PrimePi[10^9]])

General::nomem: The current computation was aborted because there was insufficient memory available to complete the computation.

Trail 2(Clumsy method)

PrimePi[10^9](*50847534*)

step0 = Total@(Prime /@ Range[1, 10^7]);
step1 = Total@(Prime /@ Range[10^7, 2*10^7]);
step2 = Total@(Prime /@ Range[2*10^7, 3*10^7]);
step3 = Total@(Prime /@ Range[3*10^7, 4*10^7]);
step4 = Total@(Prime /@ Range[4*10^7, 50847534]);

res=
 step0 + step1 + step2 + step3 + step4 - Total@(Prime /@ {10^7, 2 10^7, 3 10^7, 4 10^7})

24739512092254535

Then I discoved a most efficient algorithm by Google

It is derived from similar algorithms for counting primes. The advantage is that there is no need to find all the primes to find their sum.

Let $S(v,p)$ be the sum of integers in the range $2$ to $v$ that remain after sieving with all primes smaller or equal than $p$. That is $S(v,p)$ is the sum of integers up to $v$ that are either prime or the product of primes larger than $p$.

$S(v, p)$ is equal to $S(v, p-1)$ if $p$ is not prime or $v$ is smaller than $p^2$. Otherwise ($p$ prime, $p^2 \leq v)\quad S(v,p)$ can be computed from $S(v,p-1)$ by finding the sum of integers that are removed while sieving with $p$. An integer is removed in this step if it is the product of $p$ with another integer that has no divisor smaller than $p$. This can be expressed as

$S(v,p) = S(v, p-1) - p (S(\lfloor v/p \rfloor, p-1) - S(p-1,p-1))$

Here,we need to compute $S(n,\sqrt n)$

Dynamic programming can be used to implement this.

Performance

(A implementation can be found here)

Question

• Is it possible to implement this algorithm in Mathematica by Dynamic programming?
• It is my first time to hear Dynamic programming, so I'd like to know what's is Dynamic programming?

Answer 1

Here's a direct translation of the code you linked:

Module[{n = 10^9, r, v, p, sp},
  r = Floor@Sqrt@n;
  v = Table[Floor@(n/i), {i, 1, r}] ~Join~ Range[-1 + Floor@(n/r), 1, -1];
  ClearAll[s];
  Scan[(s[#] = (# (# + 1))/2 - 1) &, v];
  For[p = 2, p <= r, ++p,
    If[s[p] > s[p - 1],
     Scan[(s[#] -= p (s[Floor[#/p]] - s[p - 1])) &, TakeWhile[v, # >= p^2 &]];
     ]
    ]~Monitor~p;
  s[n]
  ] // AbsoluteTiming

{22.134213, 24739512092254535}

For the "dynamic programming" bit, we are using a function s to store values that we have calculated. You could try a ?s afterwards to see what's there.

Answer 2

No time for a full answer but with four cores in ParallelSum:

ParallelSum[Prime@i, {i, PrimePi[10^9]}] // AbsoluteTiming

MaxMemoryUsed[]

{55.366167, 24739512092254535}

108337144

Pretty fast and very little RAM required.

---

A modest refactoring of wxffles's code in a more native style including rasher's suggestion:

Module[{n = 10^9, r, v, s},
  r = Floor@Sqrt@n;
  v = Quotient[n, Range[r]] ~Join~ Range[Quotient[n, r] - 1, 1, -1];
  s = AssociationThread[v -> (v (v + 1))/2 - 1];
  Do[
    With[{sp1 = s[p - 1]},
     If[s[p] > sp1, 
      Scan[(s[#] -= p Subtract[s[# ~Quotient~ p], sp1]) &, 
        Pick[v, UnitStep @ Subtract[v, p^2], 1]];]],
    {p, 2, r}
  ] ~Monitor~ p;
  s[n]
] // AbsoluteTiming

{13.229757, 24739512092254535}

Explicit Subtract is used for performance; see:

• Why are numeric division and subtraction not handled better in Mathematica?

Answer 3

Dan Fortunato and I made use of Compile for this one.

index = Compile[{{p, _Integer}, {r, _Integer}, {n, _Integer}},
  If[p <= r,
    r + Quotient[n, r] - p,
    Quotient[n, p]
  ],
  CompilationTarget -> "C",
  Parallelization -> True,
  RuntimeOptions -> "Speed"
];

PrimeSum = Compile[{{n, _Integer}},
  Module[{r = Floor[Sqrt[n]], V, S, p = 2, ind = 0, i = 1, v = 0},
    V = Quotient[n, Range[r]] ~Join~ Range[Quotient[n, r] - 1, 1, -1];
    S = Quotient[#(#+1), 2] - 1 & /@ V;
    For[p = 2, p <= r, ++p, 
      If[S[[index[p, r, n]]] > S[[index[p-1, r, n]]],
        i = 1;
        While[(v = V[[i++]]) >= p^2,
          ind = index[v, r, n];
          S[[ind]] -= p (S[[index[Quotient[v, p], r, n]]] - S[[index[p-1, r, n]]])
        ];
      ]
    ];
    First[S]
  ],
  CompilationTarget -> "C",
  Parallelization -> True,
  RuntimeOptions -> "Speed",
  CompilationOptions -> {"InlineCompiledFunctions" -> True, "InlineExternalDefinitions" -> True}
];

Here are some timings

Column[AbsoluteTiming[PrimeSum[10^#]]& /@ Range[11]]

{0.000075, 17}
{6.*10^-6, 1060}
{0.000017, 76127}
{0.000069, 5736396}
{0.000319, 454396537}
{0.001524, 37550402023}
{0.007245, 3203324994356}
{0.036089, 279209790387276}
{0.176848, 24739512092254535}
{0.897198, 2220822432581729238}
{4.595891, -1447107067060386762}

Notice the final answer is negative. This is because the answer is larger than the largest machine integer (2^63 - 1), which we can verify here.

Edit: Note the correct answer for PrimeSum[10^11] is 201467077743744681014.

Answer 4

I have two efficiency increase approaches using existing functions.

1. Distribute the load evenly over the available kernels in parallel. This has gains in speed but memory usage will be about the same.
2. Distribute the load evenly over time. This significantly reduces memory usage but increases time.

Capture the number of primes we want to sum.

nn = PrimePi[10^9];

For option 1 I take the told number of primes and divide them evenly over the number of kernels. I first LaunchKernels to be certain all are ready and $KernelCount is correct. If you are going for option 2 then you may replace $KernelCount with any positive integer of your choice.

LaunchKernels[];
intPart = First@IntegerPartitions[nn, {$KernelCount}, {Floor[nn/$KernelCount],
  Ceiling[nn/$KernelCount]}];
intPart = Rest@FoldList[
 Function[{running, next}, {running[[2]] + 1, next + running[[2]]}], {0, 0}, intPart];

I restrict IntegerPartitions to only select partitions of length equal to the number of kernels (or parts specified if option 2) available (I have 4 kernels) and fix the usable values to create evenly sized partitions. I then use FoldList to get the partition intervals.

For option 1

DistributeDefinitions[intPart];
Plus @@ Parallelize[Total@Prime@Range[Sequence @@ #] & /@ intPart, 
 Method -> "CoarsestGrained"] // AbsoluteTiming
(* {43.8049, 24739512092254535} *)

Before I execute in parallel intPart has to be shared across the kernels with DistributeDefinitions. Then I can execute with Parallelize and instruct that it should be broken up into as many pieces as I have kernels (and I have ensured that I have equal pieces to number of kernels).

For option 2

Plus @@ (Total@Prime@Range[Sequence @@ #] & /@ intPart) // AbsoluteTiming
(* {112.644, 24739512092254535} *)

I have the range broken into evenly sized intervals and now just map the Total over each. This executes each item of the map in serial and so uses less memory overall.

The original bulk version

Total[Prime /@ Range[PrimePi[nn]]] // AbsoluteTiming
(* {113.591, 24739512092254535} *)

Summary

It takes 44 seconds to execute in parallel on my machine whereas the two non-parallelised versions take about 113 seconds. Each of the 4 kernels in option 1 used about 1.1GB of memory at the same time; so 4.4GB at once. However, option 2 used about 1.1GB (1/4th of option 1 while using 4 intervals over the range) for its entire execution but took 1.5 times as long to complete. Now the original bulk version uses the same amount of memory as option 1 but takes as long as option 2.

What is interesting is that option 2 and the bulk version take about the same time because each term is independent and the operator (Plus) is commutative. So with a little partitioning of the range you can get your result with your current memory restrictions.

Hope this helps.