Miller-Rabin法で擬素数を探す

Haskell

Miller-Rabbin法による擬素数探索をHaskellで書いてみた。

感想
1. 多倍長計算を気にせずできるInteger型素晴らしい！
2. とりあえず乱択のパラメーターを乱数にしているけど、IOモナド破壊しても大して困らなそうなので破壊したいと思います。

改善点
1. 実際に暗号などの目的で素数探索を使うとなると、HaskellのStdGenの乱数は暗号的にはちょっと信用ならないので乱数自体は改善する必要がある。Miller-Rabinの方のテスターは別にいいと思うのですが、素数探索の開始点は十分品質の良い暗号的乱数を選ぶべし。
2. 速度を稼ぐのであれば、事前にエラトステネスのふるいを使ってある程度ふるい落としたほうが早いらしい。(そりゃそうか)
3. 並列化を睨むのであれば逐次処理であるfind_next_prime処理を並列に投機実行できるようにするのが望ましいが、この部分をもうちょっと抽象的に書けないか。
4. Miller Rabinテスト自体は乱択なのでIOを含むが、大したIOじゃないので破壊したい！

疑問点
1. 一旦余再帰で書いた関数を末尾再帰型にリファクタしたんですけど、これ、普通なんですかね。パフォーマンスとか比べてないけどやっぱりやるべきなんだろうか。
2. Haskellコードを美しく書くコツってなんですかね。。。。今時点での実装能力ではレイアウトにどうしても気が行ってしまってコードの読みやすさなどを考慮できていない気がします。

import System.Random

modexp m p n = modexp_sub m p n m 1
  where
    modexp_sub m p n c r
      | p == 0 = r
      | otherwise =
         let res = if ( p `mod` 2 == 1 ) then (c * r) `mod` n else r in
         modexp_sub m (p `div`2) n ((c * c)`mod` n) res

find_s_k r = find_s_k_sub (r-1) 1 1
  where
    find_s_k_sub r ss s
      | (r`div` ss) `mod` 2 == 1 = (s,r `div` ss)
      | otherwise = find_s_k_sub r (2*ss) (1+ s)

millerrabbintest r t =
  let (s,k) = find_s_k r in
  let rabbintest_sub j b' =
        if (b' == r-1)||(b' == 1)
        then True
        else
            if j < s
            then rabbintest_sub (j+1)((b'*b')`mod` r)
            else False
  in
  let rabbintest a k r = rabbintest_sub 0 $ modexp a k r in
  let millerrabbin_sub r t i =
        if (i > t )
        then
          do
            putChar '\n'
            return True
        else
          do
            a <- randomRIO ( 1, r-1)
            ires <- return $ rabbintest a k r
            putChar '+'
            case ires of
              False -> return False
              True -> millerrabbin_sub r t (i+1)
  in
  millerrabbin_sub r t 1

find_next_prime n =
  do
    res <- millerrabbintest n 10
      case res of
      True -> return n
      False -> find_next_prime (n + 1)

relative_prime p q
  | p == 0 || q == 0 = False
  | p == 1 || q == 1 = True
  | p > q = relative_prime q (p `mod` q)
  | otherwise = relative_prime p (q `mod` p)

find_prime keylen =
  let randomv = let min = 2 ^ (keylen - 1)
                    max = 2 * min
                in
                randomRIO ( min,max)
  in
  do
    v <- randomv
    find_next_prime v

書いてから下記のサイトにmiller rabinテストの実装があることに気づいた。。。ただ、他のサイトでもそうなのだけど、乱択の際にIOとかが混入しちゃうのを嫌って乱択のところを端折るのがセオリーみたいです。

http://www.haskell.org/haskellwiki/Testing_primality#Miller-Rabin_Primality_Test

　　　　　　　 　∧＿∧　　　┌────────────
　　　　　　 ◯（ ´∀｀ ）◯ ＜ 僕は、unsafePerformIOちゃん！
　　　　　　　 ＼　　　 ／　　└────────────
　　　　　　　_／ ＿＿ ＼_
　　　　　　（＿／　　 ＼＿）
　　　　　　　　 　　ｌｌｌ