Which is the fastest algorithm to find prime numbers?

Question

Which is the fastest algorithm to find out prime numbers using C++? I have used sieve's algorithm but I still want it to be faster!

Given

p1, p2, p3, p4, p5, p6, p7, p8, p9 ..., pn

How can we find quickly at least a new prime?

Answer 1

I've implemented the sieve of erastosthenes in C++20 using with the punching of the non-prime multiples done multithreaded:

#include <iostream>
#include <cstdlib>
#include <charconv>
#include <cstring>
#include <vector>
#include <algorithm>
#include <cmath>
#include <bit>
#include <fstream>
#include <string_view>
#include <thread>
#include <utility>
#include <new>
#include <span>
#include <array>
#include <cassert>
#include <sstream>
#if defined(_MSC_VER) && (defined(_M_IX86) || defined(_M_X64))
    #include <intrin.h>
#elif (defined(__GNUC__) || defined(__clang__)) && (defined(__i386__) || defined(__x86_64__))
    #include <cpuid.h>
#endif

#if defined(_MSC_VER)
    #pragma warning(disable: 26495) // uninitialized member
#endif

using namespace std;

#if defined(__cpp_lib_hardware_interference_size)
constexpr ptrdiff_t CL_SIZE = hardware_destructive_interference_size;
#else
constexpr ptrdiff_t CL_SIZE = 64;
#endif

size_t getL2Size();
bool smt();

int main( int argc, char **argv )
{
    try
    {
        if( argc < 2 )
            return EXIT_FAILURE;
        auto parse = []( char const *str, auto &value )
        {
            bool hex = str[0] == '0' && tolower( str[1] ) == 'x';
            str += 2 * hex;
            auto [ptr, ec] = from_chars( str, str + strlen( str ), value, !hex ? 10 : 16 );
            return ec == errc() && !*ptr;
        };
        size_t end;
        if( !parse( argv[1], end ) )
            return EXIT_FAILURE;
        if( end < 2 )
            return EXIT_FAILURE;
        if( (ptrdiff_t)end++ < 0 )
            throw bad_alloc();
        unsigned
            hc = jthread::hardware_concurrency(),
            nThreads;
        if( argc < 4 || !parse( argv[3], nThreads ) )
            nThreads = hc;
        if( !nThreads )
            return EXIT_FAILURE;
        using word_t = size_t;
        constexpr size_t
            BITS_PER_CL = CL_SIZE * 8,
            BITS = sizeof(word_t) * 8;
        size_t nBits = end / 2;
        union alignas(CL_SIZE) ndi_words_cl { word_t words[CL_SIZE / sizeof(word_t)]; ndi_words_cl() {} };
        vector<ndi_words_cl> sieveCachelines( (nBits + BITS_PER_CL - 1) / BITS_PER_CL );
        span<word_t> sieve( &sieveCachelines[0].words[0], (nBits + BITS - 1) / BITS );
        assert(to_address( sieve.end() ) <= &to_address( sieveCachelines.end() )->words[0]);
        fill( sieve.begin(), sieve.end(), -1 );
        auto punch = [&]( size_t start, size_t end, size_t prime, auto )
        {
            size_t
                bit = start / 2,
                bitEnd = end / 2;
            if( bit >= bitEnd )
                return;
            if( prime >= BITS ) [[likely]]
                do [[likely]]
                    sieve[bit / BITS] &= rotl( (word_t)-2, bit % BITS );
                while( (bit += prime) < bitEnd );
            else
            {
                auto pSieve = sieve.begin() + bit / BITS;
                do [[likely]]
                {
                    word_t
                        word = *pSieve,
                        mask = rotl( (word_t)-2, bit % BITS ),
                        oldMask;
                    do
                        word &= mask,
                        oldMask = mask,
                        mask = rotl( mask, prime % BITS ),
                        bit += prime;
                    while( mask < oldMask );
                    *pSieve++ = word;
                } while( bit < bitEnd );
            }
        };
        size_t sqrt = (ptrdiff_t)ceil( ::sqrt( (ptrdiff_t)end ) );
        for( size_t bSqrt = sqrt / 2, bit = 1; bit < bSqrt; ++bit ) [[likely]]
        {
            auto pSieve = sieve.begin () + bit / BITS;
            do [[likely]]
                {
                    if( word_t w = *pSieve >> bit % BITS; w ) [[likely]]
                    {
                        bit += countr_zero ( w );
                        break;
                    }
                    ++pSieve;
                    bit = bit + BITS & -(ptrdiff_t)BITS;
                } while( bit < bSqrt );
            if( bit >= bSqrt ) [[unlikely]]
                break;
            size_t prime = 2 * bit + 1;
            punch ( prime *prime, sqrt, prime, false_type () );
        }
        auto scan = [&]( size_t start, size_t end, auto fn )
        {
            size_t
                bit = start / 2,
                bitEnd = end / 2;
            if( bit >= bitEnd )
                return;
            auto pSieve = sieve.begin() + bit / BITS;
            word_t word = *pSieve++ >> bit % BITS;
            for( ; ; )
            {
                size_t nextWordBit = bit + BITS & -(ptrdiff_t)BITS;
                for( unsigned gap; word; word >>= gap, word >>= 1, ++bit )
                {
                    gap = countr_zero( word );
                    bit += gap;
                    if( bit >= bitEnd ) [[unlikely]]
                        return;
                    fn( bit * 2 + 1, true_type() );
                }
                bit = nextWordBit;
                if( bit >= bitEnd ) [[unlikely]]
                    return;
                word = *pSieve++;
            }
        };
        vector<jthread> threads;
        threads.reserve( nThreads );
        static auto dbl = []( ptrdiff_t x ) { return (double)x; };
        auto initiate = [&]( size_t start, auto fn )
        {
            double threadRange = dbl( end - start ) / nThreads;
            ptrdiff_t t = nThreads;
            size_t trEnd = end;
            size_t prevStart, prevEnd;
            bool prevValid = false;
            do [[likely]]
            {
                size_t trStart = start + (ptrdiff_t)(--t * threadRange + 0.5);
                trStart &= -(2 * CL_SIZE * 8);
                trStart = trStart >= start ? trStart : start;
                if( trStart < trEnd )
                {
                    if( prevValid )
                        threads.emplace_back( fn, prevStart, prevEnd );
                    prevStart = trStart;
                    prevEnd = trEnd;
                    prevValid = true;
                }
                trEnd = trStart;
            } while( t );
            if( prevValid )
                fn( prevStart, prevEnd );
            threads.resize( 0 );
        };
        double maxCacheRange = dbl( getL2Size() * 8 * 2 ) / 3 * (smt() && nThreads > hc / 2 ? 1 : 2);
        initiate( sqrt,
            [&]( size_t rStart, size_t rEnd )
            {
                double thisThreadRange = dbl( rEnd - rStart );
                ptrdiff_t nCachePartitions = (ptrdiff_t)ceil( thisThreadRange / maxCacheRange );
                double cacheRange = thisThreadRange / dbl( nCachePartitions );
                for( size_t p = nCachePartitions, cacheEnd = rEnd, cacheStart; p--; cacheEnd = cacheStart ) [[likely]]
                {
                    cacheStart = rStart + (ptrdiff_t)((double)(ptrdiff_t)p * cacheRange + 0.5);
                    cacheStart &= -(2 * CL_SIZE * 8);
                    cacheStart = cacheStart >= sqrt ? cacheStart : sqrt;
                    scan( 3, sqrt,
                        [&]( size_t prime, auto )
                        {
                            size_t start = (cacheStart + prime - 1) / prime * prime;
                            start += start % 2 ? 0 : prime;
                            punch( start, cacheEnd, prime, true_type() );
                        } );
                }
            } );
        if( bool count = false; argc >= 3 && (!*argv[2] || (count = strcmp( argv[2], "*" ) == 0)) )
        {
            if( !count )
                return EXIT_SUCCESS;
            atomic<size_t> aNPrimes( 1 );
            initiate( 3,
                [&]( size_t rStart, size_t rEnd )
                {
                    size_t nPrimes = 0;
                    scan( rStart, rEnd, [&]( ... ) { ++nPrimes; } );
                    aNPrimes.fetch_add( nPrimes, memory_order::relaxed );
                } );
            cout << aNPrimes.load( memory_order::relaxed ) << " primes" << endl;
            return EXIT_SUCCESS;
        }
        ofstream ofs;
        ofs.exceptions( ofstream::failbit | ofstream::badbit );
        ofs.open( argc >= 3 ? argv[2] : "primes.txt", ofstream::binary | ofstream::trunc );
        constexpr size_t
            BUF_SIZE = 0x100000,
            AHEAD = 32;
        union ndi_char { char c; ndi_char() {} };
        vector<ndi_char> rawPrintBuf( BUF_SIZE + AHEAD - 1 );
        span printBuf( &to_address( rawPrintBuf.begin() )->c, &to_address( rawPrintBuf.end() )->c );
        auto
            bufBegin = printBuf.begin(),
            bufEnd = bufBegin,
            bufFlush = bufBegin + BUF_SIZE;
        auto print = [&]() { ofs << string_view( bufBegin, bufEnd ); };
        auto printPrime = [&]( size_t prime, auto )
        {
            auto [ptr, ec] = to_chars( &*bufEnd, &bufEnd[AHEAD - 1], prime );
            if( (bool)ec ) [[unlikely]]
                throw system_error( (int)ec, generic_category(), "converson failed" );
            bufEnd = ptr - &*bufBegin + bufBegin; // pointer to iterator - NOP
            *bufEnd++ = '\n';
            if( bufEnd >= bufFlush ) [[unlikely]]
                print(), bufEnd = bufBegin;
        };
        printPrime( 2, false_type() );
        scan( 3, end, printPrime );
        print();
    }
    catch( exception const &exc )
    {
        cout << exc.what() << endl;
    }
}

#if (defined(__GNUC__) || defined(__clang__)) && (defined(__i386__) || defined(__x86_64__)) || defined(_MSC_VER) && (defined(_M_IX86) || defined(_M_X64))
array<unsigned, 4> cpuId( unsigned code )
{
    int regs[4];
#if defined(_MSC_VER)
    __cpuid( regs, code );
#elif defined(__GNUC__) || defined(__clang__)
    __cpuid(code, regs[0], regs[1], regs[2], regs[3]);
#endif
    using u = unsigned;
    return array<u, 4> { (u)regs[0], (u)regs[1], (u)regs[2], (u)regs[3] };
}

bool smt()
{
    if( cpuId( 0 )[0] < 1 )
        return false;
    return cpuId( 1 )[3] >> 28 & 1;
}

size_t getL2Size()
{
    if( cpuId( 0x80000000u )[0] < 0x80000006u )
        return 512ull * 1024;
    return (cpuId( 0x80000006u )[2] >> 16) * 1024;
}
#else
size_t getL2Size()
{
    return 512ull * 1024;
}

bool smt()
{
    return false;
}
#endif

Almost all of the CPU time is spent in the multi-threaded code punching the non-prime multiples. Within the punching threads the bit range is partitioned to fit into the L2-cache. If the number of threads is beyond the half number of logical cores on a SMT-system the thread sub-partitions are half the size. Ive just added storing only the odd primes, i.e no multiples of 2, thereby having about +87% more performance. On my Ryzen 9 7950X Zen4 16-core CPU producing all primes up to 2 ^ 32, that's almost 210E6 primes, takes about 0.15 seconds on all cores with suppressed file output; on one core the algorithm takes about 1.73 seconds. The program's third parameter is the number of threads. With 16 threads on my 16 core Zen4-machine the code I get a scaling of 70% over one thread. Occupying further SMT siblings brings almost no further speedup as the code depends on the throughput of the L2-cache. The file output is done with to_chars and my own buffering to speed up the I/O. Producing a 2.2GB file from the mentioned range takes about an additional two seconds on my computer (64GB, PCIe 4.0 SSD).

Answer 2

Rabin-Miller is a standard probabilistic primality test. You run it K times and the input number is either definitely composite, or it is probably prime with probability of error 4^-K (a few hundred iterations and it's almost certainly telling you the truth).

There is a non-probabilistic (deterministic) variant of Rabin Miller.

The Great Internet Mersenne Prime Search (GIMPS) which has found the world's record for largest proven prime (2^74,207,281 - 1 as of June 2017), uses several algorithms, but these are primes in special forms. However the GIMPS page above does include some general deterministic primality tests. They appear to indicate that which algorithm is "fastest" depends upon the size of the number to be tested. If your number fits in 64 bits then you probably shouldn't use a method intended to work on primes of several million digits.

Answer 3

I wrote it today in C, compiled with tcc, figured out during preparation of competitive exams several years back. I don't know if anyone has written it alredy. It's really fast (but you should decide whether it is fast or not). It took one or two minutes to find out about 100,004 prime numbers between 10 and 10,000,000 on i7 processor with average 32% CPU use.

As you know, only those numbers can be prime which have as their last digit either 1, 3, 7, or 9; and to check if that number is prime or not, you have to divide that number by previously found prime numbers only. So first take a group of four numbers = {1,3,7,9}, test it by dividing by known prime numbers, if reminder is non-zero then the number is prime; add it to prime number array. Then add 10 to the group so it becomes {11,13,17,19} and repeat the process.

#include <stdio.h>
int main() {    
    int nums[4]={1,3,7,9};
    int primes[100000];
    primes[0]=2;
    primes[1]=3;
    primes[2]=5;
    primes[3]=7;
    int found = 4;
    int got = 1;
    int m=0;
    int upto = 1000000;
    for(int i=0;i<upto;i++){
        //printf("iteration number: %d\n",i);
        for(int j=0;j<4;j++){
            m = nums[j]+10;
            //printf("m = %d\n",m);
            nums[j] = m;
            got = 1;
            for(int k=0;k<found;k++){
                //printf("testing with %d\n",primes[k]);
                if(m%primes[k]==0){
                    got = 0;
                    //printf("%d failed for %d\n",m,primes[k]);
                    break;
                }
            }
            if(got==1){
                //printf("got new prime: %d\n",m);
                primes[found]= m;
                found++;
            }
        }
    }
    printf("found total %d prime numbers between 1 and %d",found,upto*10);
    return 0;
}

Answer 4

This is the fastest algorithm to find all prime numbers from 1 to n (on my PC it found all from 1 to 1000000 in just 0.004 seconds).

#include <iostream>
#include <fstream>

using namespace std;

double FindPrime(bool* array, int size){
    clock_t start;
    double runtime;
    for (int i = 2; i < size; i++)
        array[i] = true;
    start = clock();
    for (int i = 2; i <= size; i++)
        if (array[i])
            for (int j = 2 * i; j < size; j += i)
                array[j] = false;
    runtime = (double)(clock() - start) / CLOCKS_PER_SEC;
    return runtime;
}


int main() {
    ofstream fout("prime.txt");
    int n = 0;
    cout << "Enter the upper limit of prime numbers searching algorithm:";
    cin >> n;
    bool* array = new bool[n + 1];
    double duration = FindPrime(array, n + 1);
    printf("\n%f seconds.\n", duration);
    for (int i = 2; i <= n; i++)
        if (array[i])
            fout << i << endl;
    fout.close();
    
    return 0;
}

Answer 5

Fermat's formula is the simplest one to use but may have memory efficiency problems over the square root method and AKS method due to the usage of power calculations.

Python implementation:

def prime(n):
    # you can replace 2 with any coprime number
    if (n < 2):
        return False
    if 2**(n-1)%(n) == 1:
        return True
    return False

C:

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

bool prime(int n) {
    if (n < 2) return false;
    if ((int)pow(2, n-1) % n == 1) return true;
    return false;
}

int main() {
    int n = 5; // Example number to check
    if (prime(n)) {
        printf("%d is prime.\n", n);
    } else {
        printf("%d is not prime.\n", n);
    }
    return 0;
}

Answer 6

Until now, I believe that the fastest prime number testing algorithm is Strong Probable Prime (SPRP). I am quoting from Nvidia CUDA forums:

One of the more practical niche problems in number theory has to do with identification of prime numbers. Given N, how can you efficiently determine if it is prime or not? This is not just a thoeretical problem, it may be a real one needed in code, perhaps when you need to dynamically find a prime hash table size within certain ranges. If N is something on the order of 2^30, do you really want to do 30000 division tests to search for any factors? Obviously not.

The common practical solution to this problem is a simple test called an Euler probable prime test, and a more powerful generalization called a Strong Probable Prime (SPRP). This is a test that for an integer N can probabilistically classify it as prime or not, and repeated tests can increase the correctness probability. The slow part of the test itself mostly involves computing a value similar to A^(N-1) modulo N. Anyone implementing RSA public-key encryption variants has used this algorithm. It's useful both for huge integers (like 512 bits) as well as normal 32 or 64 bit ints.

The test can be changed from a probabilistic rejection into a definitive proof of primality by precomputing certain test input parameters which are known to always succeed for ranges of N. Unfortunately the discovery of these "best known tests" is effectively a search of a huge (in fact infinite) domain. In 1980, a first list of useful tests was created by Carl Pomerance (famous for being the one to factor RSA-129 with his Quadratic Seive algorithm.) Later Jaeschke improved the results significantly in 1993. In 2004, Zhang and Tang improved the theory and limits of the search domain. Greathouse and Livingstone have released the most modern results until now on the web, at http://math.crg4.com/primes.html, the best results of a huge search domain.

See here for more info: http://primes.utm.edu/prove/prove2_3.html and http://forums.nvidia.com/index.php?showtopic=70483

If you just need a way to generate very big prime numbers and don't care to generate all prime numbers < an integer n, you can use Lucas-Lehmer test to verify Mersenne prime numbers. A Mersenne prime number is in the form of 2^p -1. I think that Lucas-Lehmer test is the fastest algorithm discovered for Mersenne prime numbers.

And if you not only want to use the fastest algorithm but also the fastest hardware, try to implement it using Nvidia CUDA, write a kernel for CUDA and run it on GPU.

You can even earn some money if you discover large enough prime numbers, EFF is giving prizes from $50K to $250K: https://www.eff.org/awards/coop

Answer 7

Another Python implementation that is more straightforward and slightly faster* than the one in the answer of Death Mask Salesman:

import numpy as np

def prime_numbers(limit: int) -> list[int]:
    """Provide a list of all prime numbers <= the limit."""
    is_prime = np.full((limit + 1, ), True)
    is_prime[0:2] = False
    for n in range(2, limit + 1):
        if is_prime[n]:
            is_prime[n**2::n] = False
    return list(np.where(is_prime)[0])

You could further optimize by, e.g., excluding the 2, or by hard-coding even more prime numbers, but I wanted to keep it simple.

---

*exemplary runtime comparison (note: I used the optimized form of the other implementation, see my comment):

Answer 8

solution to find factors:

def divisors(integer):
    result = set()
    i = 2
    j = integer/2
    while(i <= j):
        if integer % i == 0:
            result.add(i)
            #it dont need to 
            result.add(integer//i)
        i += 1
        j = integer//i
    if len(result) > 0:
        return f"not  prime {sorted(result)}"
    else:
        return f"{integer} is prime"

--- Tests ---- import time

start_time = time.time()
print(divisors(180180180180))
print("--- %s seconds ---" % (time.time() - start_time))

--- 0.06314539909362793 seconds ---

start_time = time.time()
print(divs(180180180180180))
print("--- %s seconds ---" % (time.time() - start_time))

--- 1.5997519493103027 seconds ---

start_time = time.time()
print(divisors(1827))
print("--- %s seconds ---" % (time.time() - start_time))

--- 0.0 seconds ---

start_time = time.time()
print(divisors(104729))
print("--- %s seconds ---" % (time.time() - start_time))

--- 0.0 seconds ---

with this code:

def divs(integer):
    result = set()
    i = 2
    j = integer / 2
    loops = 0
    while (i <= j):
        if integer % i == 0:
            print(f"loops:{loops}")
            return f"{integer} is not a prime"
        i += 1
        j = integer // i
        loops += 1
    print(f"loops:{loops}")
    
    return f"{integer} is prime"

--- Tests ---

start_time = time.time()
print(divs(180180180180180180180180))
print("--- %s seconds ---" % (time.time() - start_time))

--- 0.0 seconds ---

Answer 9

This is an implementation of the Sieve of Eratosthenes in Python I've been toying with.

def eratosthenes(maximum: int) -> list[int | None]:
    """
    Find all the prime numbers between 2 and `maximum`.

    Args:
        maximum: The maximum number to check.

    Returns:
        A list of primes between 2 and `maximum`.
    """

    if maximum < 2:
        return []

    # Discard even numbers by default.
    sequence = dict.fromkeys(range(3, maximum+1, 2), True)

    for num, is_prime in sequence.items():
        # Already filtered, let's skip it.
        if not is_prime:
            continue

        # Avoid marking the same number twice.
        for num2 in range(num ** 2, maximum+1, num):
            # Here, `num2` might contain an even number - skip it.
            if num2 in sequence:
                sequence[num2] = False

    # Re-add 2 as prime and filter out the composite numbers.
    return [2] + [num for num, is_prime in sequence.items() if is_prime]

The code appears to take roughly 16s for 10000000 numbers on a humble Samsung Galaxy A40.

Suggestions are welcome!

Answer 10

I found this solution pretty fast but it comes with consequences, So this is called Fermat's Little Theorem. If we take any number p and put that in (1^p)-1 or (2^p)-2...(n^p)-n likewise and the number we get is divisible by p then it's a prime number. Talking about consequences, it's not 100% right solution. There are some numbers like 341(not prime) it will pass the test with (2^341)-2 but fails on (3^341)-3, so it's called a composite number. We can have two or more checks to make sure they pass all of them. There is one more kind of number which are not prime but also pass all the test case:( 561, 1729 Ramanujan taxi no etc.

The good thing: With (2^p)-2 in first 25 billion numbers only 2183 fails this case.

#include <iostream>
#include <math.h>
using namespace std;

int isPrime(int p)
{
    int tc = pow(2, p) - 2;
    if (tc % p == 0)
    {
        cout << p << "is Prime ";
    }
    else
    {
        cout << p << "is Not Prime";
    }
    return 0;
}

int main()
{
    int p;
    cin >> p;
    isPrime(p);
    return 0;
} 

Answer 11

I recently wrote this code to find the sum of numbers. It can be easily modified to find if a number is prime or not instead. Benchmarks are on top of the code.

// built on core-i2 e8400
// Benchmark from PowerShell
// Measure-Command { ExeName.exe }
// Days              : 0
// Hours             : 0
// Minutes           : 0
// Seconds           : 23
// Milliseconds      : 516
// Ticks             : 235162598
// TotalDays         : 0.00027217893287037
// TotalHours        : 0.00653229438888889
// TotalMinutes      : 0.391937663333333
// TotalSeconds      : 23.5162598
// TotalMilliseconds : 23516.2598
// built with latest MSVC
// cl /EHsc /std:c++latest main.cpp /O2 /fp:fast /Qpar

#include <cmath>
#include <iostream>
#include <vector>

inline auto prime = [](std::uint64_t I, std::vector<std::uint64_t> &cache) -> std::uint64_t {
    std::uint64_t root{static_cast<std::uint64_t>(std::sqrtl(I))};
    for (std::size_t i{}; cache[i] <= root; ++i)
        if (I % cache[i] == 0)
            return 0;

    cache.push_back(I);
    return I;
};

inline auto prime_sum = [](std::uint64_t S) -> std::uint64_t {
    std::uint64_t R{5};
    std::vector<std::uint64_t> cache;
    cache.reserve(S / 16);
    cache.push_back(3);

    for (std::uint64_t I{5}; I <= S; I += 8)
    {
        std::uint64_t U{I % 3};
        if (U != 0)
            R += prime(I, cache);
        if (U != 1)
            R += prime(I + 2, cache);
        if (U != 2)
            R += prime(I + 4, cache);
        R += prime(I + 6, cache);
    }
    return R;
};

int main()
{
    std::cout << prime_sum(63210123);
}

Answer 12

There is a 100% mathematical test that will check if a number P is prime or composite, called AKS Primality Test.

The concept is simple: given a number P, if all the coefficients of (x-1)^P - (x^P-1) are divisible by P, then P is a prime number, otherwise it is a composite number.

For instance, given P = 3, would give the polynomial:

   (x-1)^3 - (x^3 - 1)
 = x^3 + 3x^2 - 3x - 1 - (x^3 - 1)
 = 3x^2 - 3x

And the coefficients are both divisible by 3, therefore the number is prime.

And example where P = 4, which is NOT a prime would yield:

   (x-1)^4 - (x^4-1)
 = x^4 - 4x^3 + 6x^2 - 4x + 1 - (x^4 - 1)
 = -4x^3 + 6x^2 - 4x

And here we can see that the coefficients 6 is not divisible by 4, therefore it is NOT prime.

The polynomial (x-1)^P will P+1 terms and can be found using combination. So, this test will run in O(n) runtime, so I don't know how useful this would be since you can simply iterate over i from 0 to p and test for the remainder.

Answer 13

If it has to be really fast you can include a list of primes:
http://www.bigprimes.net/archive/prime/

If you just have to know if a certain number is a prime number, there are various prime tests listed on wikipedia. They are probably the fastest method to determine if large numbers are primes, especially because they can tell you if a number is not a prime.

Answer 14

I always use this method for calculating primes numbers following with the sieve algorithm.

void primelist()
 {
   for(int i = 4; i < pr; i += 2) mark[ i ] = false;
   for(int i = 3; i < pr; i += 2) mark[ i ] = true; mark[ 2 ] = true;
   for(int i = 3, sq = sqrt( pr ); i < sq; i += 2)
       if(mark[ i ])
          for(int j = i << 1; j < pr; j += i) mark[ j ] = false;
  prime[ 0 ] = 2; ind = 1;
  for(int i = 3; i < pr; i += 2)
    if(mark[ i ]) ind++; printf("%d\n", ind);
 }

Answer 15

I know it's somewhat later, but this could be useful to people arriving here from searches. Anyway, here's some JavaScript that relies on the fact that only prime factors need to be tested, so the earlier primes generated by the code are re-used as test factors for later ones. Of course, all even and mod 5 values are filtered out first. The result will be in the array P, and this code can crunch 10 million primes in under 1.5 seconds on an i7 PC (or 100 million in about 20). Rewritten in C it should be very fast.

var P = [1, 2], j, k, l = 3

for (k = 3 ; k < 10000000 ; k += 2)
{
  loop: if (++l < 5)
  {
    for (j = 2 ; P[j] <= Math.sqrt(k) ; ++j)
      if (k % P[j] == 0) break loop

    P[P.length] = k
  }
  else l = 0
}

Answer 16

#include<iostream>
using namespace std;

void main()
{
    int num,i,j,prime;
    cout<<"Enter the upper limit :";
    cin>>num;

    cout<<"Prime numbers till "<<num<<" are :2, ";

    for(i=3;i<=num;i++)
    {
        prime=1;
        for(j=2;j<i;j++)
        {
            if(i%j==0)
            {
                prime=0;
                break;
            }
        }

        if(prime==1)
            cout<<i<<", ";

    }
}

Answer 17

#include <iostream>

using namespace std;

int set [1000000];

int main (){

    for (int i=0; i<1000000; i++){
        set [i] = 0;
    }
    int set_size= 1000;
    set [set_size];
    set [0] = 2;
    set [1] = 3;
    int Ps = 0;
    int last = 2;

    cout << 2 << " " << 3 << " ";

    for (int n=1; n<10000; n++){
        int t = 0;
        Ps = (n%2)+1+(3*n);
        for (int i=0; i==i; i++){
            if (set [i] == 0) break;
            if (Ps%set[i]==0){
                t=1;
                break;
            }
        }
        if (t==0){
            cout << Ps << " ";
            set [last] = Ps;
            last++;
        }
    }
    //cout << last << endl;


    cout << endl;

    system ("pause");
    return 0;
}

Answer 18

#include<stdio.h>
main()
{
    long long unsigned x,y,b,z,e,r,c;
    scanf("%llu",&x);
    if(x<2)return 0;
    scanf("%llu",&y);
    if(y<x)return 0;
    if(x==2)printf("|2");
    if(x%2==0)x+=1;
    if(y%2==0)y-=1;
    for(b=x;b<=y;b+=2)
    {
        z=b;e=0;
        for(c=2;c*c<=z;c++)
        {
            if(z%c==0)e++;
            if(e>0)z=3;
        }
        if(e==0)
        {
            printf("|%llu",z);
            r+=1;
        }
    }
    printf("|\n%llu outputs...\n",r);
    scanf("%llu",&r);
}    

Answer 19

It depends on your application. There are some considerations:

• Do you need just the information whether a few numbers are prime, do you need all prime numbers up to a certain limit, or do you need (potentially) all prime numbers?
• How big are the numbers you have to deal with?

The Miller-Rabin and analogue tests are only faster than a sieve for numbers over a certain size (somewhere around a few million, I believe). Below that, using a trial division (if you just have a few numbers) or a sieve is faster.

Answer 20

Is your problem to decide whether a particular number is prime? Then you need a primality test (easy). Or do you need all primes up to a given number? In that case prime sieves are good (easy, but require memory). Or do you need the prime factors of a number? This would require factorization (difficult for large numbers if you really want the most efficient methods). How large are the numbers you are looking at? 16 bits? 32 bits? bigger?

One clever and efficient way is to pre-compute tables of primes and keep them in a file using a bit-level encoding. The file is considered one long bit vector whereas bit n represents integer n. If n is prime, its bit is set to one and to zero otherwise. Lookup is very fast (you compute the byte offset and a bit mask) and does not require loading the file in memory.

Answer 21

A very fast implementation of the Sieve of Atkin is Dan Bernstein's primegen. This sieve is more efficient than the Sieve of Eratosthenes. His page has some benchmark information.