C - determine if a number is prime

Viewed 202658

I am trying to come up with a method that takes an integer and returns a boolean to say if the number is prime or not and I don't know much C; would anyone care to give me some pointers?

Basically, I would do this in C# like this:

static bool IsPrime(int number)
{
    for (int i = 2; i < number; i++)
    {
        if (number % i == 0 && i != number)
            return false;
    }
    return true;
}
12 Answers

I would just add that no even number (bar 2) can be a prime number. This results in another condition prior to for loop. So the end code should look like this:

int IsPrime(unsigned int number) {
    if (number <= 1) return 0; // zero and one are not prime
    if ((number > 2) && ((number % 2) == 0)) return 0; //no even number is prime number (bar 2)
    unsigned int i;
    for (i=2; i*i<=number; i++) {
        if (number % i == 0) return 0;
    }
    return 1;
}

Avoid overflow bug

unsigned i, number;
...
for (i=2; i*i<=number; i++) {  // Buggy
for (i=2; i*i<=number; i += 2) {  // Buggy
// or
for (i=5; i*i<=number; i+=6) { // Buggy

These forms are incorrect when number is a prime and i*i is near the max value of the type.

Problem exists with all integer types, signed, unsigned and wider.

Example:

Let UINT_MAX_SQRT as the floor of the square root of the maximum integer value. E.g. 65535 when unsigned is 32-bit.

With for (i=2; i*i<=number; i++), this 10-year old failure occurs because when UINT_MAX_SQRT*UINT_MAX_SQRT <= number and number is a prime, the next iteration results in a multiplication overflow. Had the type been a signed type, the overflow is UB. With unsigned types, this itself is not UB, yet logic has broken down. Interations continue until a truncated product exceeds number. An incorrect result may occur. With 32-bit unsigned, try 4,294,967,291‬ which is a prime.

If some_integer_type_MAX been a Mersenne Prime, i*i<=number is never true.


To avoid this bug, consider that number%i, number/i is efficient on many compilers as the calculations of the quotient and remainder are done together, thus incurring no extra cost to do both vs. just 1.

A simple full-range solution:

bool IsPrime(unsigned number) {
    for(unsigned i = 2; i <= number/i; i++){
        if(number % i == 0){
            return false;
        }
    }
    return number >= 2;
}

To check if a number is prime i'm using a Miller/Rabin Algorithm.

#include <stdlib.h>

typedef size_t positive_number; // also try __uint128_t

static inline positive_number multiplication_modulo(positive_number lhs, positive_number rhs, positive_number mod) {
    positive_number res = 0; // we avoid overflow in modular multiplication
    for (lhs %= mod, rhs %= mod; rhs; (rhs & 1) ? (res = (res + lhs) % mod) : 0, lhs = (lhs << 1) % mod, rhs >>= 1);
    return res; // <= (lhs * rhs) % mod
}

static int is_prime(positive_number n, int k) {
    positive_number a = 0, b, c, d, e, f, g; int h, i;
    if ((n == 1) == (n & 1)) return n == 2;
    if (n < 51529) // fast constexpr check for small primes (removable)
        return (n & 1) & ((n < 6) * 42 + 0x208A2882) >> n % 30 && (n < 49 || (n % 7 && n % 11 && n % 13 && n % 17 && n % 19 && n % 23 && n % 29 && n % 31 && n % 37 && (n < 1369 || (n % 41 && n % 43 && n % 47 && n % 53 && n % 59 && n % 61 && n % 67 && n % 71 && n % 73 && ( n < 6241 || (n % 79 && n % 83 && n % 89 && n % 97 && n % 101 && n % 103 && n % 107 && n % 109 && n % 113 && ( n < 16129 || (n % 127 && n % 131 && n % 137 && n % 139 && n % 149 && n % 151 && n % 157 && n % 163 && n % 167 && ( n < 29929 || (n % 173 && n % 179 && n % 181 && n % 191 && n % 193 && n % 197 && n % 199 && n % 211 && n % 223))))))))));
    for (b = c = n - 1, h = 0; !(b & 1); b >>= 1, ++h);
    for (; k--;) {
        for (g = 0; g < sizeof(positive_number); ((char*)&a)[g++] = rand()); // random number
        do for (d = e = 1 + a % c, f = n; (d %= f) && (f %= d););
        while (d > 1 && f > 1);
        for (d = f = 1; f <= b; f <<= 1);
        for (; f >>= 1; d = multiplication_modulo(d, d, n), f & b && (d = multiplication_modulo(e, d, n)));
        if (d == 1) continue;
        for (i = h; i-- && d != c; d = multiplication_modulo(d, d, n));
        if (d != c) return 0;
    }
    return 1;
}

The test is to print the first primes :

#include <stdio.h>
int main() {
    // C Fast Iterative Algorithm
    // The First 10,000 Primes
    for (int i = 0 ; i < 104730 ; ++i)
        if (is_prime(i, 20))
            printf("%d %c", i, (i+1) % 10 ? ' ' : '\n');

    if (is_prime(9223372036854775783UL, 12))
        if (is_prime(9223372036854775643UL, 12))
            if (!is_prime(3037000493ULL * 3037000453ULL, 12))
                printf("Done.\n");
}

You can put it into a primes.c file then compile + execute :

gcc -O3 -std=c99 -Wall -pedantic primes.c ; ./a.out ;

This variant of the Fermat test has a polynomial runtime in log(n).
The __uint128_t type may be available with 128-bit integers GCC extension.

Related