(Edit by Excale, rajouté l'explication qui était dans le mail):
The points said something about an explaination, and I wasn't sure where to put it, so I just put it here:

It uses Fermat primality test with a = 2.
probablyPrime(p):
return modPow(2, p - 1, p) == 1

To compute modPow, I use Exponentiation by squaring.
modPow(base, exp, mod):
result = base
for each bit in base:
if bit is set:
result = base * result % mod
base = base * base % mod

To speed up this computation, I use Montgomery multiplication.
montgomery(first, second, mod):
result = 0
for each bit in first:
if bit is set:
a += second
if a is odd:
a += mod

In order to use Montgomery multiplication, I must put the operands into Montgomery form before and reverse the process afterwards.  If x' is the Montgomery form of x, then x' = x * r (mod m).  I chose r to be 2^32 since it just has to be a power of two larger than any operand.  Since montgomery(a, b, m) = (a * r) * (b * r) / r = a * b * r (mod m), I can reuse montgomery to compute the Montgomery form of the input.  Lastly, I can avoid inverting it by...
probablyPrime2(p):
base = 2
exp = p - 1
mod = p
result = base * 2^32 % m
base' = montgomery(base, r, mod)
for each bit in base:
if bit is set:
result = montgomery(base, result, mod)
base = montgomery(base, base, mod)
return result == base * 2^32 % m