Solovay-Strassen: 確率的素数判定アルゴリズム

実装したので，メモ．
実装といっても，単純なアルゴリズムなので，コードは至極シンプル．

import System.Random (randomRIO)

solovayStrassen :: Integer -> IO Bool
solovayStrassen n = return.isProbablyPrime n =<< randomRIO (1,n-1)

isProbablyPrime :: (Integral a) => a -> a -> Bool
isProbablyPrime n a
    | gcd a n /= 1 = False
    | otherwise    = powMod a (div (n-1) 2) n == mod (jacobi a n) n

jacobi :: Integral a => a -> a -> a
jacobi a n | a < 0 && mod n 4 == 1 = jacobi (-a) n
           | a < 0 && mod n 4 == 3 = - jacobi (-a) n
           | a < 2 = a
           | even a && (mod n 8 == 1 || mod n 8 == 7) = jacobi (div a 2) n
           | even a && (mod n 8 == 3 || mod n 8 == 5) = - jacobi (div a 2) n
           | mod a 4 == 3 && mod n 4 == 3 = - jacobi (mod n a) a
           | otherwise = jacobi (mod n a) a

-- a^n (mod p)
powMod :: (Integral a, Integral b) => a -> b -> a -> a
powMod a n p | n == 0 = 1
             | even n = powMod (mod (a*a) p) (div n 2) p
             | otherwise = (a * powMod a (n-1) p)`mod` p

main :: IO ()
main = print =<< solovayStrassen (2 ^ m26 - 1)
    where m26 = 23209 -- 26th Mersenne prime
          m27 = 44497 -- 27th Mersenne prime

Haskellが整数乗算にFFTを使用していれば，計算量はO( (log n)^2 log log n)のはず．
ちなみに，素数の入力に対しては常にTrue，合成数に対してはFalseとなる確率が1/2より大きい．