`````` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 `````` ``````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より大きい．