Fast primality testing for very large primes

Question

I'm working on a project that requires me to find whether extremely large numbers are prime numbers or not. Of course, I've read how to find prime numbers and have come up with a very simple brute force method:

def is_prime_brute_force(p):
    if p == 2 or p == 3:
        return true
    if p == 1 or p % 2 == 0 or any(p % i == 0 for i in range(3, floor_sqrt(p), 2)):
        return false
    return true

I've also investigated such probabilistic methods as the Miller-Rabin Primality Test and Fermat's little theorem (see here for Rosetta code's implementation of the former).

Though the probabilistic options are an order of magnitude faster than brute force, they're still very slow for very large inputs of n (for example, the known prime 10**9999 + 33603).

I came across an interesting observation (of course I'm not the first to come across such an observation) that all primes fit the equation $p = 6k \pm 1$. In Python, such a function looks like this

def is_prime_eq(p):
    if p == 2 or p == 3:
        return True
    if p == 0 or p == 1:
        return False

    # The same as `return (p % 6 == 1) or (p % 6 == 5)`
    prime_test = lambda p, a, m : (p % a == m) or (p % a == (a-m))
    return prime_test(p, 6, 1)

The above is guaranteed to return true if p is a prime, but a true result does not mean the p is a prime. An easy example is 25 ($25 \equiv 1\ (mod\ 6)$, but clearly $25 = 5^2$).

I'm wondering if there's some more general way to apply this interesting property of primes, perhaps with different values of a to improve the speed of my is_prime function.

Answer 1

First off, note that your "brute force" algorithm has some errors. It should be:

def is_prime_brute_force(p):
    if p == 2 or p == 3:
        return true
    if p == 1 or p % 2 == 0 or any(p % i == 0 for i in range(3, floor_sqrt(p), 2)):
        return false
    return true

In relation to this algorithm, your proposed "faster" algorithm is equivalent to

def is_prime_brute_force(p):
    if p == 2 or p == 3:
        return true
    if p == 1 or p % 2 == 0 or p % 3 == 0:
        return false
    return true

Hopefully you see why this is not terribly helpful. Any composite which is a product of primes $\geq 5$ will evaluate as a prime. Usually we use probabilistic primality tests (e.g. Miller-Rabin) for numbers whose prime divisors are all sufficiently large, so ignoring all prime divisors greater than $3$ makes it fairly useless. It's for this reason I facetiously proposed

def is_prime_brute_force(p):
    return true

in the comments as a much faster algorithm that doesn't lose much accuracy.

---

Primality tests are by their nature rather costly on current hardware. The best you can do is to try to optimize for some given assumptions on the input.

Answer 2

If we do not demand a rigorous primality test, the Miller-Rabin test is one of the best tests currently known. Very reliable is a shortcut called BPSW-test. Complete trial division is of course infeasible for huge numbers, but it is reasonable to rule out small factors before the test.

If the number is huge, a proof of the primality will only be feasible, if the number has a special form.

But even if we only apply a single weak Fermat test, the complexity is not better than $O(\ln(n))$. Apart from applying trial division upto a reasonable limit before the primality test, you cannot do any better with the currently known methods, as Brian already pointed out.