GHCソースコードリーディング勉強会 第1回 行ってきた

第0回に続いて、
GHCソースコードリーディング勉強会 第1回 #readghc に行ってきた。

というか発表しました。

発表資料：
GCをみればRTSが見えてくる、かも。。。

まー、ちょっと読み込み不足なところもあって、メンバー見たらかなりこわい人達だったので、これは殺される！って思ったんだけど、予想通りやさしく殺していただきました。

いちおー、下記のような質問がこたえられなかったんで、これはおいおい調べておきたいと思います。

1. nurseryとgeneration 0は同じ？(いちおう、Commentaryのここによると、nurseryからgen 0へコピーする、といった記述があって、まるでnurseryという明示的な領域があるように書いてあるが、、、)
2. Promotionの前にAgingをやっているはずなのだが、相当するコードが見当たらない。

自分のほか、3名の方が発表されていました。

@nothingcosmosさん

飛び込みLTスゲー。GHCのLLVMへの対応について。LLVMに対応するためにあの手この手使ってるサマがよくわかりました。。。というか、RTSの型超越的自由奔放実装はやはり良くないのでは？と思ってしまいました。

@khibinoさん

Core最適化処理を追加する！GHCのterm re-writingの処理を実装されてて面白かった。ツリー変換ぐらいまでならHaskellでサクッと(とは行かないかもだが)書けちゃうのはGHCの一つの魅力かなと思いました。

@master_qさん

gdbデバッグとCommentary GeneratedCode (I know kung fu)の解説。20分で話しきれるわけもなく速すぎてわからんかったわ。

最後に、当初の時間余りの見込みをはるかに通り越してハーフウェイギリギリまでやっちゃってすいませんでした。みんなmaster_q師匠のkung-fu聞きたかったに違いない。。。

STG祭りとCapabilityまわりの挙動について今後期待します。

Low Layer楽しいよ！

• • • • -

懇親会ではsherarcyさんに並列処理(Par Monad)教えてもらったんで、近々にSMPコードを書きたい。

• • • • -

資料を追記(2012/10/01)
議事録： HaskellJP wiki - Workshop/ReadGHC/1
http://wiki.haskell.jp/Workshop/ReadGHC/1

readghc 第一回 - Togetter
http://togetter.com/li/381712

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ちゃん！
　　　　　　　 ＼　　　 ／　　└────────────
　　　　　　　_／ ＿＿ ＼_
　　　　　　（＿／　　 ＼＿）
　　　　　　　　 　　ｌｌｌ

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)

下記のような記事を見つけた。

  星11個の数独の問題をAlloyで

なんと以前の記事を見て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でバカ遅いという状況をはるかに超えている。。。たしかに余計な述語が多かったんだよな。

精進します。