Miller-Rabin法で擬素数を探す
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ちゃん! \ / └──────────── _/ __ \_ (_/ \_) lll
Euler 4 を Haskell で高速化してみたが、、、
下記のような記事があって、
HaskellとOCamlの速度を比較してみた - let start programming = fun life -> life * 100;;
キャミバ様とガリグ先生がOCamlでコードを戦わせていたようだったので、こりゃ勝てないと思ってHaskellで高速化をしてみた。
結果、ぼくのコードは10倍に膨れ上がってしまった。
改善の趣旨としては、この問題は最大値を求める問題なので、探索空間をおおよそ大きい方からのリストにして、最初に見つかった解に対して、明らかに小さくなってしまう部分を切り捨てる。という方針。
is_circle n = let is_circle_sub m n r = if m == 0 then r == n else is_circle_sub (m `div` 10) n (r * 10 + m `mod` 10) in is_circle_sub n n 0 declist a b limit tail = if b <= a then (a,b):(declist (a-1) (b+1) limit tail ) else tail -- decendant list of pair (a,b) on value a + b declist_rec a b limit = declist a b limit $ if(a + b > 2) then if (a+b-1) > limit then declist_rec limit (a + b -1 -limit) limit else declist_rec (a+b-2) 1 limit else [] mult_list = map (\ (x,y) -> x * y) $ declist_rec 999 999 999 index_mult_list = zip [1..] mult_list main = let (index,first_sol) = head $ filter (is_circle.snd) index_mult_list root = floor $ sqrt $ fromInteger first_sol search_limit = (\x->x * x) $ (\x -> (999 - x +1)) $ root in do print $ maximum $ filter is_circle $ drop (index-1) $ take search_limit mult_list
これがhideshi_oさんのオリジナルのコード
--Haskell import Data.List main = do print $ last $ filter (\x -> (reverse $ show x) == show x) $ sort ([x * y | x <- [100..999], y <- [100..999]])
さて、実行結果は
オリジナルコード:
[seg@seg-PC ~/haskell-sample]$ time ./Euler0 906609 0.000u 0.046s 0:04.09 0.9% 0+0k 0+0io 1355pf+0w
今回アルゴリズム改善したコード:
[seg@seg-PC ~/haskell-sample]$ time ./Euler 906609 0.000u 0.031s 0:00.04 75.0% 0+0k 0+0io 1356pf+0w
爆速になった。100万要素リストのトラバースをしないのだから当然と言っては当然。実際に探索した空間は2300くらい。500倍の差がつかないのは他に余計な計算をたくさんしているせいか。あ、takeを使わないでリストをたぐるだけの関数書けば速くなりそうな気がしますね。
ただ、元のコードにあった簡潔さはもうみじんもないし、コード行数も実質15倍。トレードオフといえばそれまでだけどなんだかなぁ、という気分にもなる。
元の趣旨であったと思われる「簡潔で見やすいコードの速度を比べてみる」というところからは明らかに外れたコード修正でした。
for による組み合わせ回路
昨日の回路を少し改善してみた。
functionとforを組み合わせることで超シンプルな回路として書けた。
wireで実装するよりシミュレーション速度も速い。(仕方ないか。。。)
合成後の速度は分からないが、、、
なぜ自分がわざわざこんなことをやっているかというと、JavaRockがスゲー!って思ったからです。
JavaRock
module modexp1 ( e ,n,m, r); input[511:0] e,n,m; output[511:0] r; assign r = modexp(e,n,m); function [511:0] modexp; input [511:0] e,n,m; integer i; reg [1023:0] m1; reg [1023:0] r1; begin for(i=0;i<512;i=i+1) begin if ( i == 0) begin m1 = m; r1 = (e[0]? m : 1) % n; end else begin m1 = twice_mod(m1,n); r1 = mult_non_zero(r1, e[i],m1) % n; end end modexp = r1; end endfunction function [1023:0] mult; input [511:0] m,n; mult = m * n; endfunction function [1023:0] mult_non_zero; input [511:0] r; input [1:0] b; input [511:0] m; mult_non_zero = mult(r , (b == 1'b1 ? m : 1) ); endfunction function [511:0] twice_mod; input [511:0] m, n; twice_mod = mult(m,m) % n; endfunction endmodule
組み合わせ回路に挑戦
自作回路を作ってみた。
使用したのはicarus-verilog。
このicarusシミュレータはWindowsでも動作するのだが、ためしてみるとシミュレーションにめっちゃ時間が掛かる。
実装してみたのは冪剰余演算。512bit整数で3回ほど回すシミュレーションをしてみたのだが、、、
0.000u 0.062s 3:06.88 0.0% 0+0k 0+0io 1313pf+0w
つらい。。。
ソースコード:
https://github.com/dec9ue/verilog-study/blob/master/modexp1.v
はてな記法
いままで使ってなかった。。。。
ちょっとはてな固有でやだけど楽にハイライトできるから使ってみるかな。。。
-- Simple example for Bang-Patterns leak import Control.Monad import Data.Foldable ( foldr' ) import System.IO main :: IO () main = do result <- takeStat (return 100) 10000000 print $ "result is " ++ show result takeStat :: IO Int -> Int -> IO Int takeStat io count = foldr' (liftM2 (+)) (return 0) $ replicate count io
でも未だGist埋め込みと迷っている。。。。
世界一難しい数独@Alloy(cont'd)
下記のような記事を見つけた。
なんと以前の記事を見てdisktnkさんという方が「遅すぎる」と検証してくださったようだ。。。
Lets's Compare!
ぼくの実行結果
Executing "Run run$1 for 13" Solver=sat4j Bitwidth=0 MaxSeq=0 SkolemDepth=1 Symmetry=20 77305 vars. 2201 primary vars. 146903 clauses. 3635ms. . found. Predicate is consistent. 82243ms.
disktnkさんの実行結果
Executing "Run game" Solver=sat4j Bitwidth=0 MaxSeq=0 SkolemDepth=1 Symmetry=20 2725 vars. 729 primary vars. 3303 clauses. 508ms. . found. Predicate is consistent. 2681ms.
なんだこの複雑度の差は、、、、述語数が30~40倍くらいちがう。僕の実行環境がNetbookでバカ遅いという状況をはるかに超えている。。。たしかに余計な述語が多かったんだよな。
精進します。
入門 Verilog
長らく堪えていたが、耐え切れなくなってポチってしまった。夏休みの積ん読リスト入り。