Two Million: The Sum of Primes Below It vs. The Sum of the First That Many Primes

两百万前素数之和与前两百万素数之和 苏剑林 · 2014-06-10 · Mathematics prime number number theory Python

The title mentions two rather fun programming problems. If the wording feels dizzying, let me spell it out more clearly:

The sum of primes below two million
refers to
the sum
of all primes not exceeding two million;
the sum of the first two million primes
refers to
the sum
of the first two million primes.

I first saw this problem on Zimou's blog. The first problem is Project Euler Problem 10, and the second one is something Zimou and I explored together for fun. For Zimou's research and code, you can go study his blog. Here I'll share my own thoughts. more

Both problems are fundamentally about constructing a table of primes — the first one is a bit simpler, the second a bit more involved and computationally heavier. Building a prime table and testing primality both rely on the basic "Sieve of Eratosthenes," i.e., dividing $n$ by primes from 2 up to $\sqrt{n}$. For convenience of implementation, people usually just divide by all integers from 2 to $\sqrt{n}$. But the numbers involved in the problems discussed here are quite large, so improving efficiency is essential. Hence we need to implement a version that only divides $n$ by primes from 2 to $\sqrt{n}$.

It's worth noting that many readers have a misconception about the sieve method. Typically, people test the primality of $n$ by dividing $n$ by primes from 2 to $\sqrt{n}$, then apply this test to every number from 1 to 2 million and sum up all the primes found. This works, but it's hugely inefficient! The original Sieve of Eratosthenes doesn't work by division at all — it's described as follows:

Sieve of Eratosthenes
To find all primes not exceeding n, first find all primes $2,3,5,\dots,p_m$ between 2 and $\sqrt{n}$, then successively strike out, from 1 to n, the multiples of 2 (except 2 itself), the multiples of 3 (except 3 itself), ..., the multiples of $p_m$ (except $p_m$ itself).

In other words, building a prime table with the Sieve of Eratosthenes requires only multiplication and deletion operations! Multiplication is far more efficient than division, and this algorithm also performs many fewer operations overall than testing each number one by one.

The implementation idea is: first build an array containing all integers from 1 to $n$, then successively set to 0 every entry that's a multiple of 2, every entry that's a multiple of 3, ..., and every entry that's a multiple of $p_m$. But how do we get the primes from 1 to $\sqrt{n}$ in the first place? Simple — we just pull them from the array itself as we go, judging and collecting at the same time. The rate at which primes are generated grows faster than the rate at which we need them. (What this means is: for example, within 20, we only need to strike out multiples of 2 and 3 — excluding themselves — to obtain all the primes 2, 3, 5, 7, ..., 19; the largest prime reached is 19, which is already enough to test primality for all numbers up to 400. To get each prime $p_i$ between 2 and $\sqrt{n}$, we simply take the nonzero entries near the front of the array, since composite numbers have already been zeroed out. With a little thought, the reader will see that this process advances progressively.)

At this point the basic programming approach is clear. However, we also need to account for the fact that both problems already involve fairly large numbers — especially the second one, since the first two million primes extend to roughly 35 million. This brings up the question of tooling. I initially wanted to use C++ for its efficiency, but C++ arrays have size limits, and integer types have limits too, making overflow easy; so I considered using the GMP big-number library with C++, but my C++ skills were too weak and I could never get it working properly (apparently the GMP authors have never been willing to support VC or VS, and all the workarounds online are community patches; I also considered writing my own simple big-integer class, but didn't have the energy to fuss over it). I also thought about MATLAB, but MATLAB doesn't support high-precision big-integer arithmetic — it only gives approximate values, which isn't quite suitable. I then considered Mathematica, which is itself a symbolic computation tool and supports arbitrarily large integers, but I found that Mathematica is also quite slow when handling arrays with millions of elements. In the end I settled on Python.

I had tried Python a long time ago but set it aside for ages. Today I discovered that Python natively supports big-integer arithmetic! Although it's not as fast as C++, it handles the programming tasks in this post comfortably, and after some optimization the runtime was significantly reduced. I've attached the code at the end of the post for anyone interested.

Test results:

In Python 3.3

Sum of primes below two million:
142913828922
time: 2.4048174478605646
Sum of the first two million primes:
31381137530481
time: 46.75734807838953

Apparently Python 2.7 runs faster, but in my own testing on my machine, I didn't notice a clear difference. Computing the sum of the first two million primes took about 50 seconds on average, which is a bit faster than Zimou's algorithm — probably because I constructed the prime table directly rather than testing numbers one by one. What we have in common is that both Zimou and I only test divisibility by primes between 2 and $\sqrt{n}$, rather than by all integers.

Looking forward to seeing an even more efficient Python algorithm~

Sum of primes up to 2 million

#200万以内的素数之和
import time
start=time.clock()
import math

n=2000000 #定义上限
prime=[i for i in range(1,n+1)] #定义整数表

r=int(math.sqrt(n))

#下面是用删除式的方法把整数表中的合数删除掉
for j in range(2,r+1):
    if prime[j-1] != 0:
        s=j*j
        while s <= n:
            prime[s-1]=0
            s=s+j

print(sum(prime)-1) #求和        
end=time.clock()

print("time:",end-start)

Sum of the first 2 million primes

#前200万个素数之和
import time
start=time.clock()
import math
n=35000000 
#定义上限，两百万个素数大约是前3500万的素数
#这是根据公式pi(n)约等于n/ln(n)得到的。
prime=[i for i in range(1,n+1)]    #定义整数表

r=int(math.sqrt(n))

#下面是用删除式的方法把整数表中的合数删除掉
for j in range(2,r+1):
    if prime[j-1] != 0:
        s=j*j
        while s <= n:
            prime[s-1]=0
            s=s+j

#下面开始确定素数个数到200万个。
prime.sort() #从小到大重新排列（让0位于前面）
z=prime.count(0) #统计0的个数
print(sum(prime[z+1:z+1+2000000])) #计算前两百万个素数之和。
end=time.clock()

print("time:",end-start)