The Sum Of The First Billion Primes

September 11, 2012

The first ten prime numbers are 2, 3, 5, 7, 11, 13, 17, 19, 23, and 29. Their sum is 129.

Your task is to write a program that calculates the sum of the first billion primes. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

About these ads

Pages: 1 2

9 Responses to “The Sum Of The First Billion Primes”

  1. Kim Walisch said
    //////////////////////////////////////////////////////////////////////
    // sum_primes.cpp
    
    #include 
    #include 
    #include 
    
    uint64_t sum = 0;
    
    void callback(uint64_t prime)
    {
      sum += prime;
    }
    
    int main()
    {
      PrimeSieve ps;
      ps.generatePrimes(2, (uint64_t) 22801763489.0, callback);
      std::cout << "Sum = " << sum << std::endl;
      return 0;
    }
    
    //////////////////////////////////////////////////////////////////////
    

    The above C++ code sums primes using my primesieve library (https://primesieve.googlecode.com). To compile it on a Linux/Unix system use:

    $ svn checkout http://primesieve.googlecode.com/svn/trunk/ primesieve
    $ cd primesieve
    $ make lib
    $ sudo make install
    
    $ g++ -O2 sum_primes.cpp -o sum_primes -lprimesieve
    $ ./sum_primes
    Sum = 11138479445180240497
    

    This takes about 10 seconds on my Intel Core-i5 CPU.

  2. A (inefficient) Python solution:

    import time
    
    def primes():
        yield 2
        yield 3
        stack = [3]
        while True:
            n = stack[-1] + 2
            is_prime = False
            while not is_prime:
                is_prime = True
                for p in stack:
                    if p * p > n:
                        break
                    if n % p == 0:
                        is_prime = False
                        break
                if is_prime:
                    stack.append(n)
                else:
                    n += 2
            yield stack[-1]
    
    def sum_primes(n):
        result = 0
        count = 0
        for p in primes():
            result += p
            count += 1
            if count == n:
                break
        return result
    
    if __name__ == '__main__':
        start = time.clock()
        print(sum_primes(1000000))
        print("elapsed %.2fs" % (time.clock() - start))
    

    On my machine (Intel Core i3) this takes around 87 seconds to complete.

  3. Jeff said

    Java noob’s super inefficient, but extensively documented, Java

    public class AddPrimes {
    
        public static void main(String[] args) {
            
            //Initialize number as 2, the first possible prime number larger than one
            number = 2;
            //Main Loop; loop until 1 billion prime numbers have been discovered
            while(primeCounter < 1000000){
                //Initialize fields for Calculation Loop.
                m = Math.sqrt(number);
                divisor = 2;
                //Calculation Loop
                while(divisor <= m){
                    //Initialize quotient
                    quotient = number /divisor;
                    //Initialize isInteger
                    isInteger = checkIsInteger(quotient);
                    //If the quotient is an integer, add it to the list of factors.  Increase divisor by 1.
                    if(isInteger == true ){
                        factors.add(quotient);
                        divisor++;
                    //If the quotient is not an integer, increase the divisor by 1.
                    }else{
                        divisor++;
                    }
                }
                //Once all quotients are analyzed, if the list of factors is empty, then the number is prime.
                //Add the number to the sum of primes.  Increase primeCounter by 1.  Increase number by 1.
                if(factors.isEmpty() == true){
                    sum += number;
                    primeCounter ++;
                    number++;
                //If the list of factors is not empty, then the number is not prime.  Clear the list of factors.
                //Increase number by 1.
                }else{
                    factors.clear();
                    number++;
                }
            }
            //Display the sum of the total number of primes discovered.
            System.out.println(sum);
        }
        /**
         * Determines whether or not the given double is an integer.
         * @param d - double
         * @return true if d is an integer or false if d is not an integer
         */
        public static boolean checkIsInteger(double d){
            boolean intCheck;
            String doubleText = Double.toString(d);
            if(doubleText.endsWith("0")){
                intCheck = true;
            }else{
                intCheck = false;
            }
            return intCheck;
        }
        
        //Necessary fields
            /**
             * The current number being tested for primality.
             */
            static double number;
            /**
             * The sum of all numbers determined to be prime.
             */ 
            static int sum;
            /**
             * The square root of the current number.
             */ 
            static double m;
            /**
             * The current divisor.
             */
            static int divisor;
            /**
             * The quotient of number / divisor.
             */
            static double quotient;
            /**
             * A counter that records the total number of discovered prime numbers.
             */
            static int primeCounter;
            /**
             * Boolean placeholder for the result of the checkInteger method.
             */
            static boolean isInteger;
            /**
             * List of any factors discovered for number.
             */
            static ArrayList factors = new ArrayList();
    }
    
  4. David said

    Done in SwiftForth with a bitset library I wrote a long time ago. Took 40 minutes or 410K primes per second. (While also watching web video, etc. on the computer.) It uses the compressed 33K file from previous exercise to store primes < 1 million, which is used as the sieving primes, and the segmented sieve from another previous exercise to sieve blocks after 1 million, in increments of 100,000.

    requires bitset

    \ Table lookup for "small" primes, (33K table for primes < 1e6)

    1,000,003 drop Constant FILE_PRIME_MAX

    : load-primes ( str — addr )
    r/o bin open-file throw locals| fd |
    fd file-size throw drop ( — size )
    s" CREATE prime-table" evaluate
    here over allot
    swap fd read-file throw drop
    fd close-file throw ;

    s" ../prime1e6.dat" load-primes

    hex
    create masks[]
    FF c, \ 0
    FE c, FE c, FE c, FE c, FE c, FE c, \ 1 – 6
    FC c, FC c, FC c, FC c, \ 7 – 10
    F8 c, F8 c, \ 11 – 12
    F0 c, F0 c, F0 c, F0 c, \ 13 – 16
    E0 c, E0 c, \ 17 – 18
    C0 c, C0 c, C0 c, C0 c, \ 19 – 22
    80 c, 80 c, 80 c, 80 c, 80 c, 80 c, \ 23 – 28
    00 c, \ 29

    decimal

    create offset_mod_11
    -1 , 1 , 7 , -1 , 11 , 17 , -1 , 29 , 13 , 23 , 19 ,

    : lowbit>offset
    dup negate and 11 mod cells offset_mod_11 + @ ;

    : next-small-prime ( n — n )
    dup 2 < IF drop 2 EXIT THEN
    dup 3 < IF drop 3 EXIT THEN
    dup 5 < IF drop 5 EXIT THEN
    dup FILE_PRIME_MAX >= abort" : precomputed prime table exceeded"

    30 /mod prime-table +
    swap masks[] + c@ \ addr mask
    over c@ and lowbit>offset
    BEGIN dup 0< WHILE
    drop
    1+ dup c@ lowbit>offset
    REPEAT
    swap prime-table – 30 * + ;

    \ segmented sieve

    100,000 drop Constant SIEVE_BLOCK_SIZE
    SIEVE_BLOCK_SIZE 2/ bset sieve

    2variable sieve_start \ used by the generator
    : sieve_end ( — d )
    sieve_start 2@ [ SIEVE_BLOCK_SIZE 1- ] Literal M+ ;

    : d/mod ( d n — d%n d/n )
    locals| n |
    \ FORTH does not have remainder function for double division
    \ so we multiply the quotient by the divisor and subtract…
    \
    2dup 1 n M*/ 2dup 2>r \ get quotient & save a copy
    n 1 M*/ d- drop \ compute remainder
    2r> ; \ restore quotient

    : dmod ( d m — n )
    locals| m |
    m d/mod 2drop \ note: FM/MOD is faster but overflows…
    dup 0< IF m + THEN ;

    : offset ( n — n ) \ first offset into bit array of prime #
    sieve_start 2@ third 1+ M+ D2/ dnegate rot dmod ;

    : seg-sieve ( d — )
    sieve_start 2!
    sieve bits erase
    sieve bset-invert! \ set everything to TRUE

    3
    BEGIN
    dup offset
    BEGIN dup sieve nip ( length ) < WHILE
    dup sieve -bit
    over +
    REPEAT drop
    next-small-prime
    sieve_end third dup M* D< UNTIL
    drop ;

    \ gen-primes — generate primes and call callback for each one.
    \ callback returns 0 when no more primes desired
    \ After ~1e12 primes generated will throw exception…

    variable 'with_prime
    variable complete?

    : sieve_cb ( n — n )
    complete? @ IF
    drop
    ELSE
    sieve_start 2@ rot 2* 1+ M+ 'with_prime @execute not complete? !
    THEN ;

    : gen-primes ( xt — )
    'with_prime ! false complete? !

    2
    BEGIN
    dup s>d 'with_prime @execute not IF drop EXIT THEN
    next-small-prime
    dup 1000000 > UNTIL drop

    1,000,000
    BEGIN
    2dup seg-sieve
    sieve ['] sieve_cb bset-foreach
    complete? @ IF 2drop EXIT THEN
    SIEVE_BLOCK_SIZE M+
    AGAIN ;

    \ sum first N primes

    variable #primes
    2variable sum

    : sum_cb ( d — ? )
    sum 2@ D+ sum 2!
    -1 #primes +!
    #primes @ 0> ;

    : sum_primes ( n — d )
    0, sum 2! #primes !
    ['] sum_cb gen-primes
    sum 2@ ;

    counter 1,000,000,000 drop sum_primes du. timer 11138479445180240497 2401035 ok

  5. David said

    Trying this again, hope it formats better…

    Done in SwiftForth with a bitset library I wrote a long time ago. Took 40 minutes or 410K primes per second. (While also watching web video, etc. on the computer.) It uses the compressed 33K file from previous exercise to store primes < 1 million, which is used as the sieving primes, and the segmented sieve from another previous exercise to sieve blocks after 1 million, in increments of 100,000.

    requires bitset

    \ Table lookup for "small" primes, (33K table for primes < 1e6)

    1,000,003 drop Constant FILE_PRIME_MAX

    : load-primes ( str — addr )
    r/o bin open-file throw locals| fd |
    fd file-size throw drop ( — size )
    s" CREATE prime-table" evaluate
    here over allot
    swap fd read-file throw drop
    fd close-file throw ;

    s" ../prime1e6.dat" load-primes

    hex
    create masks[]
    FF c, \ 0
    FE c, FE c, FE c, FE c, FE c, FE c, \ 1 – 6
    FC c, FC c, FC c, FC c, \ 7 – 10
    F8 c, F8 c, \ 11 – 12
    F0 c, F0 c, F0 c, F0 c, \ 13 – 16
    E0 c, E0 c, \ 17 – 18
    C0 c, C0 c, C0 c, C0 c, \ 19 – 22
    80 c, 80 c, 80 c, 80 c, 80 c, 80 c, \ 23 – 28
    00 c, \ 29

    decimal

    create offset_mod_11
    -1 , 1 , 7 , -1 , 11 , 17 , -1 , 29 , 13 , 23 , 19 ,

    : lowbit>offset
    dup negate and 11 mod cells offset_mod_11 + @ ;

    : next-small-prime ( n — n )
    dup 2 < IF drop 2 EXIT THEN
    dup 3 < IF drop 3 EXIT THEN
    dup 5 < IF drop 5 EXIT THEN
    dup FILE_PRIME_MAX >= abort" : precomputed prime table exceeded"

    30 /mod prime-table +
    swap masks[] + c@ \ addr mask
    over c@ and lowbit>offset
    BEGIN dup 0< WHILE
    drop
    1+ dup c@ lowbit>offset
    REPEAT
    swap prime-table – 30 * + ;

    \ segmented sieve

    100,000 drop Constant SIEVE_BLOCK_SIZE
    SIEVE_BLOCK_SIZE 2/ bset sieve

    2variable sieve_start \ used by the generator
    : sieve_end ( — d )
    sieve_start 2@ [ SIEVE_BLOCK_SIZE 1- ] Literal M+ ;

    : d/mod ( d n — d%n d/n )
    locals| n |
    \ FORTH does not have remainder function for double division
    \ so we multiply the quotient by the divisor and subtract…
    \
    2dup 1 n M*/ 2dup 2>r \ get quotient & save a copy
    n 1 M*/ d- drop \ compute remainder
    2r> ; \ restore quotient

    : dmod ( d m — n )
    locals| m |
    m d/mod 2drop \ note: FM/MOD is faster but overflows…
    dup 0< IF m + THEN ;

    : offset ( n — n ) \ first offset into bit array of prime #
    sieve_start 2@ third 1+ M+ D2/ dnegate rot dmod ;

    : seg-sieve ( d — )
    sieve_start 2!
    sieve bits erase
    sieve bset-invert! \ set everything to TRUE

    3
    BEGIN
    dup offset
    BEGIN dup sieve nip ( length ) < WHILE
    dup sieve -bit
    over +
    REPEAT drop
    next-small-prime
    sieve_end third dup M* D< UNTIL
    drop ;

    \ gen-primes — generate primes and call callback for each one.
    \ callback returns 0 when no more primes desired
    \ After ~1e12 primes generated will throw exception…

    variable 'with_prime
    variable complete?

    : sieve_cb ( n — n )
    complete? @ IF
    drop
    ELSE
    sieve_start 2@ rot 2* 1+ M+ 'with_prime @execute not complete? !
    THEN ;

    : gen-primes ( xt — )
    'with_prime ! false complete? !

    2
    BEGIN
    dup s>d 'with_prime @execute not IF drop EXIT THEN
    next-small-prime
    dup 1000000 > UNTIL drop

    1,000,000
    BEGIN
    2dup seg-sieve
    sieve ['] sieve_cb bset-foreach
    complete? @ IF 2drop EXIT THEN
    SIEVE_BLOCK_SIZE M+
    AGAIN ;

    \ sum first N primes

    variable #primes
    2variable sum

    : sum_cb ( d — ? )
    sum 2@ D+ sum 2!
    -1 #primes +!
    #primes @ 0> ;

    : sum_primes ( n — d )
    0, sum 2! #primes !
    ['] sum_cb gen-primes
    sum 2@ ;

    counter 1,000,000,000 drop sum_primes du. timer 11138479445180240497 2401035 ok

  6. [...] This problem from Programming Praxis came about in the comments to my last post and intrigued me. So today, we are trying to sum the first one billion primes. Summing the first hundred, thousand, even million primes isn’t actually that bad. But it takes a bit more effort when you scale it up to a billion. And why’s that? [...]

  7. JP said

    I know it’s an older, but after your comments on my solution to the Pythagorean Triples hinting at this one, I had to try it. I ended up working on it way more than I probably should have with quals coming up, but what I ended up with is 5 different ways (including all three sieves mentioned on Wikipedia) of summing the first n primes with different data structure variants on the two fastest ones.

    It still ends up taking an hour and a half on the quickest run, but it was really interesting to work out all of the ways to do it. Perhaps I’ll see if I can optimize even more at some point. Although 10 seconds? Likely not. :)

    Linky: The Sum of the First Billion Primes

  8. danaj said

    Five methods in Perl.

    The first is calling Math::Prime::Util::next_prime in a dumb loop. Not very fast since it’s doing trial division and/or MR tests. 400k-500k / second.

    We can speed that up by having MPU precalculate all the primes, obviously at the expense of memory (about 760MB for sieving 22.8B numbers). Using nth_prime_upper() to get an upper bound (0.01% larger than the real value, but obtained in microseconds) we can do this pretty fast. A reasonably fast solution at over 3M primes/s.

    Next is using Math::Prime::Util to do it more efficiently — get primes in segments until we’ve reached our limit. Not bad — 66s on my machine to sum the first 1 billion. Over 15M/s.

    Now we move to pure Perl. First let’s use an unsegmented sieve. This one is a modification of one from c2.com’s wiki, which is a very simple SoE, and much faster than most Perl sieves found on the web (I have a collection of 17 of them — some shockingly bad, some quite decent). It runs at about 260k/s. But the biggest problem is that it isn’t segmented, nor is it storing the numbers efficiently. My 32GB machine started swapping when calculating 10^8.

    Lastly is a segmented sieve in pure Perl using strings. A bit wacky and too verbose, but it can be faster than the obvious solutions depending on memory speed. As a basic unsegmented sieve, it’s about 2x faster than any of the other Perl solutions on my workstation. On my old Atom netbook it’s a bit slower than the above sieve. I wrote it earlier this year so Math::Prime::Util could have something reasonable as a fallback if the XS failed to compile on a platform. Since it’s segmented it will work for our task without going crazy with memory. Over 500k/s.

    None of them compare to primesieve of course (the serial version of Kim’s program takes all of 7 seconds on my machine, the parallel version runs in 1.1s). The callback is spiffy also. I thought about adding something like it to Math::Prime::Util (e.g. prime_foreach 2 100 { …sub… }) and still might, but calling Perl subroutines is far slower than C — it’s over a minute just to make 1B empty subroutine calls so it seems of limited use for this example.

    Timings:

           MPU next     next+precalc   MPU sieve      Perl big    Perl segment
          ------------  ------------  ------------  ------------  ------------
    10^6      1.6s          0.3s          0.1s          2.9s          1.4s
    10^7     20.4s          2.3s          0.6s         37.7s         16.0s
    10^8    252.1s         26.0s          6.1s          n/a         197.6s
    10^9                  304.9s         65.8s          n/a        2548.0s
    

    Code:

    #!/usr/bin/perl
    use warnings;
    use strict;
    use Math::Prime::Util qw/:all -nobigint/;
    use List::Util qw/sum/;
    
    my $n = $ARGV[0] || 1_000_000;
    #print "sum of the first $n primes: ", prime_sum_next($n, 1), "\n";
    #print "sum of the first $n primes: ", prime_sum_sieve($n), "\n";
    #print "sum of the first $n primes: ", prime_sum_ns($n), "\n";
    print "sum of the first $n primes: ", prime_sum($n), "\n";
    
    sub prime_sum_next {
      my ($np, $do_precalc) = @_;
      prime_precalc(nth_prime_upper($np)) if $do_precalc;
      my $sum = 0;
      my $p = 0;
      while ($np-- > 0) {
        $p = next_prime($p);
        $sum += $p;
      }
      $sum;
    }
    
    sub prime_sum_sieve {
      my $np = shift;
      my $sum = 0;
      my $low = 1;
      my $segsize = 5_000_000;
      while ($np > 0) {
        my $chunkref = primes($low, $low+$segsize);
        $low += $segsize + 2;
        $#$chunkref = $np-1 if scalar @$chunkref > $np;
        $sum += int(sum @$chunkref);
        $np -= scalar @$chunkref;
      }
      $sum;
    }
    
    sub prime_sum_ns {   # Basic unsegmented sieve
      my $np = shift;
      return if !defined $np || $np < 0;
      return (0,2,2+3,2+3+5)[$np] if $np <= 3;
      my $max = nth_prime($np);
      my @c;
      for(my $t = 3; $t*$t <= $max; $t += 2) {
         if (!$c[$t]) {
             for(my $s = $t*$t; $s <= $max; $s += $t*2) { $c[$s]++ }
         }
      }
      my $sum = 2;
      for (my $t = 3; $t <= $max; $t += 2) {
        $c[$t] || do {$sum += $t;}
      }
      $sum;
    }
    
    sub prime_sum {   # Using string segmented sieve
      my $np = shift;
      return if !defined $np || $np < 0;
      return (0,2,2+3,2+3+5)[$np] if $np <= 3;
      $np -= 3;
    
      my $sum = 2+3+5;
      my $segsize = 1_000_000;
      my $low = 7;
      while ($np > 0) {
        my $high = $low + $segsize;
        my $sref = _pp_segment_sieve($low, $high);
        my $p = $low - 2;
        foreach my $s (split("0", $$sref, -1)) {
          $p += 2 + 2 * length($s);
          last if $p > $high;
          $sum += $p;
          last if --$np <= 0;
        }
        $low = $high + 2;
      }
      $sum;
    }
    
    sub _pp_segment_sieve {
      my($beg,$end) = @_;
      my $range = int( ($end - $beg) / 2 ) + 1;
      # Prefill with 3 and 5 already marked, and offset to the segment start.
      my $whole = int( ($range+14) / 15);
      my $startp = ($beg % 30) >> 1;
      my $sieve = substr("011010010010110", $startp) . "011010010010110" x $whole;
      # Set 3 and 5 to prime if we're sieving them.
      substr($sieve,0,2) = "00" if $beg == 3;
      substr($sieve,0,1) = "0"  if $beg == 5;
      # Get rid of any extra we added.
      substr($sieve, $range) = '';
    
      # If the end value is below 7^2, then the pre-sieve is all we needed.
      return \$sieve if $end < 49;
    
      my $limit = int(sqrt($end)) + 1;
      # For large value of end, it's a huge win to just walk primes.
      my $primesieveref = _pp_string_sieve($limit);
      my $p = 7-2;
      foreach my $s (split("0", substr($$primesieveref, 3), -1)) {
        $p += 2 + 2 * length($s);
        my $p2 = $p*$p;
        last if $p2 > $end;
        if ($p2 < $beg) {
          $p2 = int($beg / $p) * $p;
          $p2 += $p if $p2 < $beg;
          $p2 += $p if ($p2 % 2) == 0;   # Make sure p2 is odd
        }
        # With large bases and small segments, it's common to find we don't hit
        # the segment at all.  Skip all the setup if we find this now.
        if ($p2 <= $end) {
          # Inner loop marking multiples of p
          # (everything is divided by 2 to keep inner loop simpler)
          my $filter_end = ($end - $beg) >> 1;
          my $filter_p2  = ($p2  - $beg) >> 1;
          while ($filter_p2 <= $filter_end) {
            substr($sieve, $filter_p2, 1) = "1";
            $filter_p2 += $p;
          }
        }
      }
      \$sieve;
    }
    
    sub _pp_string_sieve {
      my($end) = @_;
      return 0 if $end < 2;
      return 1 if $end < 3;
      $end-- if ($end & 1) == 0;
    
      # string with prefill
      my $whole = int( ($end>>1) / 15);
      die "Sieve too large" if $whole > 1_145_324_612;  # ~32 GB string
      my $sieve = '100010010010110' . '011010010010110' x $whole;
      substr($sieve, ($end>>1)+1) = '';
      my $n = 7;
      while ( ($n*$n) <= $end ) {
        my $s = $n*$n;
        my $filter_s   = $s >> 1;
        my $filter_end = $end >> 1;
        while ($filter_s <= $filter_end) {
          substr($sieve, $filter_s, 1) = '1';
          $filter_s += $n;
        }
        do { $n += 2 } while substr($sieve, $n>>1, 1);
      }
      return \$sieve;
    }
    
  9. JP said

    Dang. After some comments over on my original post, I decided to go back and add a the segmented version of the Sieve of Eratosthenes into the mix. I’m glad I did.

    While my previous times were 2 hr 42 min for Eratosthenes and 1 hr 28 min for Atkin, the new segmented sieve blows both out of the water at only 25 min 12 sec. It’s still nowhere near as quick as some people’s, but it’s at least a good example of what being conscious of what memory/caching behavior can do to runtimes.

    Linky: One Billion Primes – Segmented Sieve

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 576 other followers

%d bloggers like this: