Почему мой алгоритм Миллера Рабина не работает (Haskell)?

Question

Я реализовал тест Миллера Рабина в Хаскелле. Я попытался строго следовать псевдокоду, как это указано, например, в записи wikipedia для теста Миллера Рабина. Теперь я нашел онлайн, что для определенных выборов свидетелей тест является детерминированным до определенных заданных границ. Меня интересуют простые слова под 2^64, и поэтому я нашел достаточные оценки в этом посте What witnesses do i need for Rabin-Miller test for numbers up to 10¹⁸?. Тем не менее, этот код, похоже, работает для большинства небольших простых чисел, которые я тестировал, но не для некоторых более крупных. Например, я попробовал десятизначную цифру 5915587277, и тест возвращает false. Я считаю, что моя реализация правильная, но, надеюсь, кто-то может заметить, где я ошибся, и неправильно понял что-то о МР-тесте. Заранее благодарю за любую помощь. Кроме того, извините за грязный код.

isPrime :: Int -> Bool 
isPrime n = millerRabinTest n (factorizeN (n-1)) 

{- factorizeN finds a number s and odd number d such that n -1 = (2^s)d by 
succesively dividing n by two if it is even. -} 
factorizeN :: Int -> (Int, Int) 
factorizeN n = fN n 0 
    where 
    fN n s | even n = fN (n `div` 2) (s + 1) 
      | otherwise = (n,s) 

{- this is the main function. it takes w values from a set of witnesses 
and checks if n passes the test. If it doesn't, n is not prime, if it does 
for all w, it is probably prime. -} 
millerRabinTest :: Int -> (Int,Int) -> Bool 
millerRabinTest n (d,s) = and [test n (expmod w d n) s | w <- onesToCheck] 

{- this is the test that is used in the millerRabinTest function. it sees if 
w^d = 1 mod n or n-1 mod n, if not it multiplies by two 
and checks again for a total of s-1 times. If it is never true then the number 
is not prime -} 
test :: Int -> Int -> Int -> Bool 
test n w s | w `elem` [1,n-1] = True 
      | otherwise  = or [ (expmod w (2^k) n) `elem` [1,n-1]| k <- [1..s]] 

{- set of witnesses that should make the Miller Rabin test deterministic if 
n < 2^64. -} 
onesToCheck :: [Int] 
onesToCheck = [2,325,9375,28178,450775,9780504,1795265022] 

{- function that calculates a^e mod n. -} 
expmod :: Int -> Int -> Int -> Int 
expmod a e n | e == 1   = a `mod` n 
       | (e `mod` 2) == 0 = (expmod ((a*a) `mod` n) (e `div` 2) n) 
       | otherwise  = (a*(expmod ((a*a) `mod` n) (e `div` 2) n)) `mod` n 

Answer 1

Возможно, ваш Int переполнен в expmod при вычислении a*a. Int - это целое число размером от машины, не более 64 бит. Вы должны заменить некоторые из вхождений Int в свою программу с помощью Integer, цельного типа произвольной точности.