在 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)
参考链接
不是最优, 仅供参考.
博文
