质数/素数 Haskell 求解

在 ProjectEuler 经常遇到质数求解的问题, 需要一个够强大的算法. 尝试自己写了个, 用了 Memoization:

primeList :: [Integer]
primeList
= 2:3:5:[p | p <- [7,9..]
      , all ((0 /=) . mod p) $ takeWhile (<= ceiling $ sqrt $ fromInteger p) primeList ]

在 Hugs 下测试:

Main> primeList !! 1000
7927 :: Integer
(3341273 reductions, 4756123 cells, 5 garbage collections)

貌似不错, 但是很明显, 代码适用范围有限, 例如 primeList !! 10000 在半天后终于出结果. ssword 提供了一个更短的算法, 目测貌似没我的效率高:

primes :: [Integer]
primes = sieve [2..]
    where sieve (p : xs) = p : sieve [x | x <- xs, x `mod` p > 0]

代码相当巧妙, 很不错的筛选法, 但是效率实测... 估计是时间花在 [] 上和不断的回朔上, Lazy 嘛~ 发现用到下个值只好回朔到 2, 3, 5, 7, 11, 13....

Main> primes !! 1000
7927 :: Integer
(35111063 reductions, 51142767 cells, 59 garbage collections)

以上算法都只是处理小质数数列, 遇到大质数.... 全部死掉. PS, stack overflow. 熟悉的错误. 我找到了这个 http://en.wikipedia.org/wiki/Miller-Rabin_test, 感觉不错, 下面的 Haskell 算法是自己写的, 很丑陋.... 根据 Python 代码改的. Miller Rabin 筛选法.

-- Miller-Rabin methold By Andelf
-- true for p < 10^16 except `46 856 248 255 981’, so you can add 13
-- optimize: order of and, or, any, all/ all judge to odd/even
isPrime :: Integer -> Bool
isPrime p
    | p == 2          = True 
    | even p || p < 2 = False
    | otherwise    
        = all (isPrimeTest p) (takeWhile (< p) [2, 3, 7, 61, 24251])
    where 
      isPrimeTest :: Integer -> Integer -> Bool   
      isPrimeTest n a 
          = odd d || t == n - 1 
          where 
            (t, d) = loop (pow a (n-1) n) (shiftTillOdd (n - 1))
            loop t d = if t /= 1 && d /= n - 1 && t /= n - 1 
                       then loop (t * t `mod` n) (d + d) 
                       else (t, d)
            shiftTillOdd = until odd $ flip div 2
            pow :: Integer -> Integer -> Integer -> Integer 
            pow a d n
                | even d            = pow (a * a `mod` n) (d `div` 2) n `mod` n 
                | d == 1 || a == 1  = a
                | otherwise         = pow (a * a `mod` n) (d `div` 2) n * a `mod` n

本来想用 AKS 的, 后来发现在多项式判断那块时间太长, 所以还是采用了伪素数算法, 你所看到的第 2 行注释, p < 10^16 except `46 856 248 255 981’ 是没问题的. 代码经过简单优化, guard 判断的分支顺序可能很诡异, mod, even, odd, 等函数的使用可能也很诡异, 但请相信这个是在 Hugs 下测试过的, 目前是最优. div/(+) 就是比 Data.Bits 快很多. 很诡异. ghc 同学请自己权衡. 再次动手可能在 pow 函数上. [2, 3, 7, 61, 24251] 的选择请参考链接.

Main> head $ filter isPrime [10 ^ 100..]
100000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000267 :: Integer
(8965148 reductions, 25399307 cells, 29 garbage collections)

参考链接

• http://en.wikipedia.org/wiki/Miller-Rabin_test
• http://en.literateprograms.org/Miller-Rabin_primality_test_(Python) 很多语言的实现
• http://primes.utm.edu/prove/prove2_3.html

不是最优, 仅供参考.