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

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?
• FWIW, Total[Prime /@ Range[PrimePi[10^9]]] evaluates to 24739512092254535 — i.e. the same as your trial 2 result. The computation took a few minutes, and the peak kernel memory usage was about 4.7GB. – Stephen Luttrell Apr 19 '15 at 14:01
• Approximation to 5-7 digits on the OP's values for n: sumP = #*RiemannR[#] - (NIntegrate[RiemannR[x], {x, 2, #}]) &. It takes about 0.15 sec. on all (those) inputs. – Michael E2 Apr 20 '15 at 0:06
• Dynamic programming (DP) is, "a method for solving a complex problem by breaking it down into a collection of simpler subproblems." As an example, solving the problem of making $n$ cents with the fewest number of coins has a DP algorithm where we find the number of ways to make each amount $\leq n$ using only pennies, then with pennies and nickels, then wtih dimes, etc. Each solution is a simple function of two of the previously calculated solutions. (I learned about DP while competing in USACO, a good way to build your algorithm skills.) – 2012rcampion Apr 20 '15 at 2:21
• @2012rcampion,OK, thanks a lot :) – xyz Apr 20 '15 at 5:05

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.

• You'll get a nice speed bump changing TakeWhile... to Pick[v, UnitStep@Subtract[v, p^2], 1]. +1 on translating that. – ciao Apr 20 '15 at 4:45

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:

• +1 This is why you are the Wiz. Man, my answer is huge compared to yours and not nearly as good. LOL! – Edmund Apr 19 '15 at 16:50

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]

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

• @ShutaoTang Please see my revision. I just edited the code to be much faster. – Chip Hurst Apr 22 '15 at 14:37

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[] + 1, next + running[]}], {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.