Fast deterministic primality test for 32bit numbers

isprime.c

#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>

/*
 * returns a^e % n
 */
static unsigned int
powm(unsigned int a, unsigned int e, unsigned int n)
{
        uint64_t x = a;
        uint64_t y = 1;

        while (e) {
                if ((e & 1)) {
                        y *= x;
                        y %= n;
                }
                e >>= 1;
                x *= x;
                x %= n;
        }

        return (unsigned int)y;
}

/*
 * returns true if a < n,
 * a^d % n != 1, and
 * a^(2^r * d) % n != n - 1 for all r in [0, s - 1]
 */
static int
test_for_composite(unsigned int a, unsigned int s, unsigned int d, unsigned int n)
{
        uint64_t x;

        if (a >= n || (x = powm(a, d, n)) == 1)
                return 0;

        while (s--) {
                if (x == n - 1)
                        return 0;
                x *= x;
                x %= n;
        }

        return 1;
}

/*
 * this function is valid for n < 4,759,123,141,
 * but uint32_t only goes to 4,294,967,295 ;)
 */
static int
is_prime(uint32_t n)
{
        unsigned int d;
        unsigned int s;

        /* if n is 2 return true */
        if (n == 2)
                return 1;

        /* otherwise if n is even return false */
        if (n == 1 || (n & 1) == 0)
                return 0;

        /* n - 1 = 2^s * d */
        for (s = 0, d = n - 1; (d & 1) == 0; s++, d >>= 1);

        if (test_for_composite( 2, s, d, n) ||
            test_for_composite( 7, s, d, n) ||
            test_for_composite(61, s, d, n))
                return 0;

        return 1;
}

int
main(int argc, char *argv[])
{
        uint32_t i;

        for (i = 0; i < 1000000; i++) {
                if (is_prime(i))
                        printf("%u is prime\n", i);
        }

        return EXIT_SUCCESS;
}