Miller Rabin Overview
The Miller Rabin primality test is a probabilistic primality test. It can tell you if a number is probably prime. The test is based on something that comes as a consequence of Fermat’s Little Theorem which states
\[a^{p-1} \equiv 1 \pmod p\]Which implies \(a^{p-1} - 1\equiv 0 \pmod p \). So p must divide \( a^{p-1}-1 \) if it is prime.
Since p is an odd prime p-1 is even, so \(p-1 = 2^k * q \) for some odd q. This means we can take the square root of (\( a^{p-1} \)) \( k \) times. So \( (a^{p-1}-1) \) can be written as a difference of squares.
\begin{align}
a^{p-1}-1 &= (a^{\frac{p-1}{2}} + 1)(a^{\frac{p-1}{2}} - 1)
&= (a^{\frac{p-1}{2}} + 1)(a^{\frac{p-1}{4}} + 1)(a^{\frac{p-1}{4}} - 1)
&= (a^{\frac{p-1}{2}} + 1)(a^{\frac{p-1}{4}} + 1) \cdots (a^{\frac{p-1}{2^k}} + 1)(a^{\frac{p-1}{2^k}} - 1)
\end{align}
Thus if p is prime it must divide one of the above factors.
Further Explanation
We can see that the last term \( (a^{\frac{p-1}{2^k}} - 1) \) is congruent to 1 mod n
\[a^{\frac{p-1}{2^k}} - 1 \equiv 0 \pmod p \to a^{\frac{p-1}{2^k}} \equiv 1 \pmod p\]We can also see every other term is congruent to -1
\[a^{\frac{p-1}{2^k}} + 1 \equiv 0 \pmod p \to a^{\frac{p-1}{2^k}} \equiv -1 \pmod p\]And \( -1 \pmod p = p -1 \). So if \( a \equiv -1 \pmod p \) then \( a \pmod p = (p - 1) \). We’ll use this in the actual implementation.
Algorithm
In this context test fails means its probably prime.
- If n is even or \( 1 < \gcd(a,n) < n \) return composite.
- Write \( n-1 = 2^k q \) with q odd.
- Set \( a \equiv a^q \pmod n \).
- If \( a \equiv 1 \pmod n \), return Test Fails.
- Loop \(i = 0,1,2,…,k-1 \)
- If \( a \equiv -1 \pmod n \), return Test Fails.
- Set \( a \equiv a^2 \pmod n \).
- Increment i and loop again at step 5.
- Return composite.
Implementation
The most important part of the code is the strongPrimeTest function.
strongPrimeTest :: Integer -> Integer -> Bool
strongPrimeTest n b | (gcd n b) /= 1 = False
| mpow b q n == 1 = True
| mpow b q n == (n - 1) = True
| otherwise = testSquares (k - 1) n $ mpow b q n
where
k = findk $ n - 1
q = findq $ n - 1
testSquares i n b | i < 0 = False
| (mpow b 2 n) == (n - 1) = True
| otherwise = testSquares (i - 1) n (mpow b 2 n)
This is a fairly direct translation of the algorithm. mpow = powMod is a modular power function
for computing modular exponents effectively. b is the base being tested as a witness for the primality
of n.
The millerRabin routine selects random bases to test against its input n. It also handles some
corner cases, such as is it negative, is it even, etc.
millerRabin :: (RandomGen g) => Integer -> g -> Bool
millerRabin n g | n < 2 = False --2 is first prime and all primes are positive
| n == 2 = True --2 is prime
| even n = False --all other primes are odd
| otherwise = checkRandomBases n 10 [x | x <- (randInt (n - 1) g), odd x]
where
checkRandomBases :: Integer -> Integer -> [Integer] -> Bool
checkRandomBases n iters (r:rs) | iters == 0 = True --all tests passed
| strongPrimeTest n r == False = False
| otherwise = checkRandomBases n (iters - 1) rs
This is the actual implementation in use in a program to generate large primes.
import Prelude
import Math.NumberTheory.Roots
import System.Random
import System.IO
import System.Environment
import Math.NumberTheory.Powers.Modular
-------------------------UTILITIES---------------------------
isqrt = integerSquareRoot
mpow = powMod --modular power function, this one is deprecated, but it was easier to get working
-- mpow b e m = (b^e) `mod` m --modular power --this one is naive it will die for large nums
-------------------------RANDOM NUMBERS-----------------------
randBits :: (RandomGen g) => Integer -> g -> [Integer] --return stream of random numbers of n bits
randBits nbits = randomRs (2^nbits, 2^(nbits+1) - 1)
randInt :: (RandomGen g) => Integer -> g -> [Integer] --return stream of random nums in range
randInt n = randomRs (2, n - 1)
-------------------------MILLER/RABIN-------------------------
--find k such that n = 2^k * q where q is some odd
findk :: Integer -> Integer
findk n = if n==0 then 0 else findkhelp n 0
where
findkhelp n k | odd n = k
| otherwise = findkhelp (div n 2) (k+1)
findq :: Integer -> Integer --find q s.t n = 2^k * q, constraints n > 0
findq n = if even n then findq (div n 2) else n
strongPrimeTest :: Integer -> Integer -> Bool
strongPrimeTest n b | (gcd n b) /= 1 = False
| mpow b q n == 1 = True
| mpow b q n == (n - 1) = True
| otherwise = testSquares (k - 1) n $ mpow b q n --loop part of algorithm
where
k = findk $ n - 1
q = findq $ n - 1
testSquares i n b | i < 0 = False
| (mpow b 2 n) == (n - 1) = True
| otherwise = testSquares (i - 1) n (mpow b 2 n)
millerRabin :: (RandomGen g) => Integer -> g -> Bool
millerRabin n g | n < 2 = False --2 is first prime and all primes are positive
| n == 2 = True --2 is prime
| even n = False --all other primes are odd
| otherwise = checkRandomBases n 10 [x | x <- (randInt (n - 1) g), odd x]
where
checkRandomBases :: Integer -> Integer -> [Integer] -> Bool
checkRandomBases n iters (r:rs) | iters == 0 = True --all tests passed
| strongPrimeTest n r == False = False
| otherwise = checkRandomBases n (iters - 1) rs
--Alias for millerRabin
isprime :: (RandomGen g) => Integer -> g -> Bool
isprime = millerRabin
--Driver
genPrime :: (RandomGen g) => Integer -> g -> Integer
genPrime nBits g | nBits < 2 = error "Smallest prime cannot fit in less than 2 bits"
| otherwise = findPrime (randBits nBits g) g
where findPrime (r:rs) g | isprime r g = r
| otherwise = findPrime rs g
printUsage = do
putStr "Usage: genPrime [nbits]\n"
putStr " Prints a prime number with the specified number of bits\n"
putStr " example: 'genPrime 100' creates and prints a 100 bit prime number\n"
main = do
g <- getStdGen
args <- getArgs
if length args /= 1 then
printUsage
else do
let nBits = read (args !! 0) :: Integer
print ( genPrime nBits g)