Good algorithms to determine whether one number is prime?

Question

Possible Duplicate:
Checking if an int is prime more efficiently

I need to test some very large integer to see whether it is a prime. Could you provide some good algorithms or library routines?

EDIT: C/C++ would be OK.

Thanks.

Answer 1

The Miller-Rabin test is very fast, and it can be both fast and deterministic for certain ranges of numbers.

Answer 2

The easiest option would be use use an existing big integer library. It won't have bugs, and it will provide all the supporting functions.

If you are writing your own implementation (i.e. for an assignment), I'd recommend working from a pseudo code algorithm in a book so that you understand what you are doing.

That being said, one of the simplest methods is to use Jacobi and Legendre, and compare for equality. I just submitted an assignment for RSA encryption. Here is what I did for single precision, however the algorithms are general and work for multiple precision integers too.

typedef uint64_t BigIntT;
typedef int64_t SBigIntT;

// This function calculations the power of b^e mod phi
// As long as 
//      b*b is smaller than max(BigIntT) 
//      b*phi is smaller than max(BigIntT)
// we will not have overflow.
BigIntT calculatePower (BigIntT b, BigIntT e, BigIntT m) {
    BigIntT result = 1;

    while (e != 0) {
        if (e & 1) {
            result = (result * b) % m;
        }

        e = e >> 1;
        b = (b * b) % m;
    }

    return result;
}

// This function implements simple jacobi test.
// We can expect compiler to perform tail-call optimisation.
SBigIntT jacobi (SBigIntT a, SBigIntT b) {
    if (a == 0 || a == 1) {
        return a;
    } else if (a % 2 == 0) {
        if (((b*b - 1) / 8) % 2 == 0) {
            return jacobi(a/2, b);
        } else {
            return -jacobi(a/2, b);
        }
    } else if ((((a-1) * (b-1)) / 4) % 2 == 0) {
        return jacobi(b % a, a);
    } else {
        return -jacobi(b % a, a);
    }
}

// This function implements : http://en.wikipedia.org/wiki/Solovay-Strassen_primality_test
bool testPrime (BigIntT p) {
    int tests = 10;

    if (p == 2) {
        return true;
    }

    while (tests-- > 0) {
        BigIntT a = generateRandomNumber(2, p);

        if (greatestCommonDivisor(a, p) == 1) {
            BigIntT l = calculatePower(a, (p-1)/2, p);
            SBigIntT j = jacobi(a, p);

            // j % p == l
            if ((j == -1) && (l == p-1) || (j == l)) {
                // So far so good...
            } else {
                // p is composite
                return false;
            }
        } else {
            // p is composite
            return false;
        }
    }

    return true;
}