Most efficient code for the first 10000 prime numbers?

Viewed 70519

I want to print the first 10000 prime numbers. Can anyone give me the most efficient code for this? Clarifications:

  1. It does not matter if your code is inefficient for n >10000.
  2. The size of the code does not matter.
  3. You cannot just hard code the values in any manner.
30 Answers

The Sieve of Atkin is probably what you're looking for, its upper bound running time is O(N/log log N).

If you only run the numbers 1 more and 1 less than the multiples of 6, it could be even faster, as all prime numbers above 3 are 1 away from some multiple of six. Resource for my statement

I recommend a sieve, either the Sieve of Eratosthenes or the Sieve of Atkin.

The sieve or Eratosthenes is probably the most intuitive method of finding a list of primes. Basically you:

  1. Write down a list of numbers from 2 to whatever limit you want, let's say 1000.
  2. Take the first number that isn't crossed off (for the first iteration this is 2) and cross off all multiples of that number from the list.
  3. Repeat step 2 until you reach the end of the list. All the numbers that aren't crossed off are prime.

Obviously there are quite a few optimizations that can be done to make this algorithm work faster, but this is the basic idea.

The sieve of Atkin uses a similar approach, but unfortunately I don't know enough about it to explain it to you. But I do know that the algorithm I linked takes 8 seconds to figure out all the primes up to 1000000000 on an ancient Pentium II-350

Sieve of Eratosthenes Source Code: http://web.archive.org/web/20140705111241/http://primes.utm.edu/links/programs/sieves/Eratosthenes/C_source_code/

Sieve of Atkin Source Code: http://cr.yp.to/primegen.html

GateKiller, how about adding a break to that if in the foreach loop? That would speed up things a lot because if like 6 is divisible by 2 you don't need to check with 3 and 5. (I'd vote your solution up anyway if I had enough reputation :-) ...)

ArrayList primeNumbers = new ArrayList();

for(int i = 2; primeNumbers.Count < 10000; i++) {
    bool divisible = false;

    foreach(int number in primeNumbers) {
        if(i % number == 0) {
            divisible = true;
            break;
        }
    }

    if(divisible == false) {
        primeNumbers.Add(i);
        Console.Write(i + " ");
    }
}

@Matt: log(log(10000)) is ~2

From the wikipedia article (which you cited) Sieve of Atkin:

This sieve computes primes up to N using O(N/log log N) operations with only N1/2+o(1) bits of memory. That is a little better than the sieve of Eratosthenes which uses O(N) operations and O(N1/2(log log N)/log N) bits of memory (A.O.L. Atkin, D.J. Bernstein, 2004). These asymptotic computational complexities include simple optimizations, such as wheel factorization, and splitting the computation to smaller blocks.

Given asymptotic computational complexities along O(N) (for Eratosthenes) and O(N/log(log(N))) (for Atkin) we can't say (for small N=10_000) which algorithm if implemented will be faster.

Achim Flammenkamp wrote in The Sieve of Eratosthenes:

cited by:

@num1

For intervals larger about 10^9, surely for those > 10^10, the Sieve of Eratosthenes is outperformed by the Sieve of Atkins and Bernstein which uses irreducible binary quadratic forms. See their paper for background informations as well as paragraph 5 of W. Galway's Ph.D. thesis.

Therefore for 10_000 Sieve of Eratosthenes can be faster then Sieve of Atkin.

To answer OP the code is prime_sieve.c (cited by num1)

Using GMP, one could write the following:

#include <stdio.h>
#include <gmp.h>

int main() {
  mpz_t prime;
  mpz_init(prime);
  mpz_set_ui(prime, 1);
  int i;
  char* num = malloc(4000);
  for(i=0; i<10000; i++) {
    mpz_nextprime(prime, prime);
    printf("%s, ", mpz_get_str(NULL,10,prime));
  }
}

On my 2.33GHz Macbook Pro, it executes as follows:

time ./a.out > /dev/null

real    0m0.033s
user    0m0.029s
sys    0m0.003s

Calculating 1,000,000 primes on the same laptop:

time ./a.out > /dev/null

real    0m14.824s
user    0m14.606s
sys     0m0.086s

GMP is highly optimized for this sort of thing. Unless you really want to understand the algorithms by writing your own, you'd be advised to use libGMP under C.

Not efficient at all, but you can use a regular expression to test for prime numbers.

/^1?$|^(11+?)\1+$/

This tests if, for a string consisting of k “1”s, k is not prime (i.e. whether the string consists of one “1” or any number of “1”s that can be expressed as an n-ary product).

I have adapted code found on the CodeProject to create the following:

ArrayList primeNumbers = new ArrayList();

for(int i = 2; primeNumbers.Count < 10000; i++) {
    bool divisible = false;

    foreach(int number in primeNumbers) {
        if(i % number == 0) {
            divisible = true;
        }
    }

    if(divisible == false) {
        primeNumbers.Add(i);
        Console.Write(i + " ");
    }
}

Testing this on my ASP.NET Server took the rountine about 1 minute to run.

Sieve of Eratosthenes is the way to go, because of it's simplicity and speed. My implementation in C

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>

int main(void)
{
    unsigned int lim, i, j;

    printf("Find primes upto: ");
    scanf("%d", &lim);
    lim += 1;
    bool *primes = calloc(lim, sizeof(bool));

    unsigned int sqrtlim = sqrt(lim);
    for (i = 2; i <= sqrtlim; i++)
        if (!primes[i])
            for (j = i * i; j < lim; j += i)
                primes[j] = true;

    printf("\nListing prime numbers between 2 and %d:\n\n", lim - 1);
    for (i = 2; i < lim; i++)
        if (!primes[i])
            printf("%d\n", i);

    return 0;
}

CPU Time to find primes (on Pentium Dual Core E2140 1.6 GHz, using single core)

~ 4s for lim = 100,000,000

Adapting and following on from GateKiller, here's the final version that I've used.

    public IEnumerable<long> PrimeNumbers(long number)
    {
        List<long> primes = new List<long>();
        for (int i = 2; primes.Count < number; i++)
        {
            bool divisible = false;

            foreach (int num in primes)
            {
                if (i % num == 0)
                    divisible = true;

                if (num > Math.Sqrt(i))
                    break;
            }

            if (divisible == false)
                primes.Add(i);
        }
        return primes;
    }

It's basically the same, but I've added the "break on Sqrt" suggestion and changed some of the variables around to make it fit better for me. (I was working on Euler and needed the 10001th prime)

The Sieve seems to be the wrong answer. The sieve gives you the primes up to a number N, not the first N primes. Run @Imran or @Andrew Szeto, and you get the primes up to N.

The sieve might still be usable if you keep trying sieves for increasingly larger numbers until you hit a certain size of your result set, and use some caching of numbers already obtained, but I believe it would still be no faster than a solution like @Pat's.

In Python

import gmpy
p=1
for i in range(10000):
    p=gmpy.next_prime(p)
    print p 

I have written this using python, as I just started learning it, and it works perfectly fine. The 10,000th prime generate by this code as same as mentioned in http://primes.utm.edu/lists/small/10000.txt. To check if n is prime or not, divide n by the numbers from 2 to sqrt(n). If any of this range of number perfectly divides n then it's not prime.

import math
print ("You want prime till which number??")
a = input()
a = int(a)
x = 0
x = int(x)
count = 1
print("2 is prime number")
for c in range(3,a+1):
    b = math.sqrt(c)
    b = int(b)
    x = 0
    for b in range(2,b+1):
        e  = c % b
        e = int(e)
        if (e == 0):
            x = x+1
    if (x == 0):
        print("%d is prime number" % c)
        count = count + 1
print("Total number of prime till %d is %d" % (a,count))

Here the code that I made :


enter code here
#include <cmath>
#include <cstdio>
#include <vector>
#include <iostream>
#include <algorithm>
using namespace std;


int main() {
/* Enter your code here. Read input from STDIN. Print output to STDOUT*/   

unsigned long int n;

int prime(unsigned long int);

scanf("%ld",&n);

unsigned long int val;

for(unsigned long int i=0;i<n;i++)
{
    int flag=0;

    scanf("%ld",&val);

    flag=prime(val);

    if(flag==1)
        printf("yes\n");

    else
        printf("no\n");
}

return 0;

}

int prime(unsigned long int n)
{

if(n==2) return 1;

else if (n == 1||n%2==0)  return 0;

for (unsigned long int i=3; i<=sqrt(n); i+=2)
    if (n%i == 0)
        return 0;

return 1;
}

Using the Array.prototype.find() method in Javascript. 2214.486 ms

function isPrime (number) {

  function prime(element) {
    let start = 2;
    while (start <= Math.sqrt(element)) {
      if (element % start++ < 1) {
        return false;
      }
    }
    return element > 1;
  }

  return [number].find(prime)

}

function logPrimes (n) {

  let count = 0
  let nth = n

  let i = 0
  while (count < nth) {
    if (isPrime(i)) {
      count++
      console.log('i', i) //NOTE: If this line is ommited time to find 10,000th prime is 121.157ms
      if (count === nth) {
        console.log('while i', i)
        console.log('count', count)
      }
    }
    i++
  }

}

console.time(logPrimes)

logPrimes(10000)

console.timeEnd(logPrimes) // 2214.486ms

I can give you some tips, you have to implement it.

  1. For each number, get the half of that number. E.g. for checking 21, only obtain the remainder by dividing it from range 2-10.
  2. If its an odd number, only divide by odd number, and vice versa. Such as for 21, divide with 3, 5, 7, 9 only.

Most efficient method I got up to so far.

Since you want first 10000 primes only, rather than coding complex algorithm I'll suggest the following

boolean isPrime(int n){
//even but is prime
    if(n==2)
        return true;
//even numbers filtered already 
    if(n==0 || n==1 || n%2==0)
        return false;

// loop for checking only odd factors
// i*i <= n  (same as i<=sqrt(n), avoiding floating point calculations)
    for(int i=3 ; i*i <=n ; i+=2){
    // if any odd factor divides n then its not a prime!
        if(n%i==0)
            return false;
    }
// its prime now
    return true;
}

now call is prime as you need it

for(int i=1 ; i<=1000 ; i++){
    if(isPrime(i)){
       //do something
    }
}

This is an old question, but there's something here everyone's missing...

For primes this small, trial division isn't that slow... there are only 25 primes under 100. With so few primes to test, and such small primes, we can pull out a neat trick!

If a is coprime to b, then gcd a b = 1. Coprime. Fun word. Means it doesn't share any prime factors. We can thus test for divisibility by several primes with one GCD call. How many? Well, the product of the first 15 primes is less than 2^64. And the product of the next 10 is also less than 2^64. That's all 25 that we need. But is it worth it?

Let's see:

check x = null $ filter ((==0) . (x `mod`)) $ [<primes up to 101>]
Prelude> length $ filter check [101,103..85600]
>>> 9975
(0.30 secs, 125,865,152 bytes

a = 16294579238595022365 :: Word64
b = 14290787196698157718
pre = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97]
primes = (pre ++) $ filter ((==1) . gcd a) $ filter ((==1) . gcd b) [99,101..85600]
main = print $ length primes

Prelude> main
>>> 10000
(0.05 secs, 36,387,520 bytes)

A 6 fold improvement there.

(length is to force the list to be computed. By default Haskell prints things 1 Unicode character at a time and so actually printing the list will either dominate the time or dominate the amount of actual code used.)

Of course, this is running in GHCi - a repl running interpreted code - on an old laptop and it is not interpreting any of these numbers as int64s or even BigInts, nor will it even if you ask it to (well, you can force it, but it's ugly and doesn't really help). It is interpreting every single number there as generalized Integer-like things that can be specialized to some particular type via dictionary lookup, and it is traversing a linked list (which is not fused away here as it's not compiled) 3 times. Interestingly, hand fusing the two filters actually slows it down in the REPL.

Let's compile it:

...\Haskell\8.6\Testbed>Primes.exe +RTS -s
10000
606,280 bytes allocated in the heap
Total   time    0.000s  (  0.004s elapsed)

Using the RTS report because Windows. Some lines trimmed because they aren't relevant - they were other GC data, or measurements of only part of the execution, and together add up to 0.004s (or less). It's also not constant folding, because Haskell doesn't actually do much of that. If we constant fold ourselves (main = print 10000), we get dramatically lower allocation:

...Haskell\8.6\Testbed>Primes.exe +RTS -s
10000
47,688 bytes allocated in the heap
Total   time    0.000s  (  0.001s elapsed)

Literally just enough to load the runtime, then discover there's nothing to do but print a number and exit. Let's add wheel factorization:

wheel = scanl (+) 7 $ cycle [4, 2, 4, 2, 4, 6, 2, 6]
primes = (pre ++) $ filter ((==1) . gcd a) $ filter ((==1) . gcd b) $ takeWhile (<85600) wheel

Total   time    0.000s  (  0.003s elapsed)

Cut down approximately 1/3rd relative to our reference of main = print 10000, but there's definitely room for more optimization. It actually stopped to perform a GC in there for example, while with tweaking there shouldn't be any heap use. For some reason, compiling for profiling here actually cuts the runtime down to 2 milliseconds:

Tue Nov 12 21:13 2019 Time and Allocation Profiling Report  (Final)

   Primes.exe +RTS -p -RTS

total time  =        0.00 secs   (2 ticks @ 1000 us, 1 processor)
total alloc =     967,120 bytes  (excludes profiling overheads)

I'm going to leave this as is for now, I'm pretty sure random jitter is starting to dominate.

def compute_primes(bound):
"""
Return a list of the prime numbers in range(2, bound)
Implement the Sieve of Eratosthenes
https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
"""
primeNumber = [True for i in range(bound + 1)]
start_prime_number = 2
primes = []
while start_prime_number * start_prime_number <=bound:
    # If primeNumber[start_prime_number] is not changed, then it is a prime
    if primeNumber[start_prime_number]:
        # Update all multiples of start_prime_number
        for i in range(start_prime_number * start_prime_number, bound + 1, start_prime_number):
            primeNumber[i] = False
    start_prime_number += 1

# Print all prime numbers
for start_prime_number in range(2, bound + 1):
    if primeNumber[start_prime_number]:
        primes.append(start_prime_number)

return primes

print(len(compute_primes(200)))

print(len(compute_primes(2000)))

Related