読者です 読者をやめる 読者になる 読者になる

0CTF Writeup: oneTimePad2

Japanese Version

Task

oneTimePad1と同じように暗号化スクリプトoneTimePad2.pyとciphertxtが入ったzipファイルがもらえる。
今度は\(GF(2^{128})\)のブロック暗号システムのようだ。

\begin{align}
C_i &= M_i + R_i \\
\begin{bmatrix}
R_{i+1} \\
1
\end{bmatrix} &=
\begin{bmatrix}
A & B\\
0 & 1\\
\end{bmatrix}^{N_i}
\begin{bmatrix}
R_{i} \\
1
\end{bmatrix} \\
N_{i+1} &= N_i^2\\
\end{align}
\(C_1 \dots, C_5 , M_1, \dots, M_4, A, B \)が既知である。\(M_5\)にフラグがあるようだ。
秘密鍵は\(N_1\)なのでこれを解読したい。

解析

まず、行列の累乗を計算する。
\begin{align}
\begin{bmatrix}
A & B\\
0 & 1\\
\end{bmatrix}^{n} &=
\begin{bmatrix}
A^n & A^{n-1} B + A^{n-2} B + \dots + B \\
0 & 1
\end{bmatrix} \\
&=
\begin{bmatrix}
A^n & \frac{(A^{n} - 1)B}{A - 1}\\
0 & 1
\end{bmatrix}
\end{align}
したがって
\[
R_{2} = A^{N_1} R_{1} + \frac{(A^{N_1} - 1)B}{A - 1}
\]
これを\(A^{N_1}\)について解くと
\[
A^{N_1} = (R_2 + B (A - 1)^{-1})(R_1 + B(A-1)^{-1})^{-1}
\]
よって\(A^{N_1}\)の値は計算できる。
\(N_1\)を計算するのは離散対数問題といって一般の体では計算するのは難しいが、少なくとも\(GF(2^{128})\)なら離散対数は現実的な時間で計算できるようだ。

解法

pari/gpで離散対数を計算する処理までを行った。

gistc2a4481d5f18d35e913503664320a016

$ gp -q doit.gp  
x^123 + x^117 + x^116 + x^114 + x^113 + x^112 + x^110 + x^109 + x^108 + x^107 + x^104 + x^99 + x^98 + x^97 + x^95 + x^93 + x^91 + x^89 + x^87 + x^84 + x^83 + x^82 + x^81 + x^80 + x^78 + x^74 + x^69 + x^68 + x^66 + x^62 + x^57 + x^53 + x^52 + x^51 + x^45 + x^43 + x^42 + x^41 + x^38 + x^37 + x^35 + x^33 + x^32 + x^31 + x^30 + x^29 + x^28 + x^25 + x^22 + x^19 + x^18 + x^17 + x^15 + x^14 + x^13 + x^12 + x^10 + x^8 + x^5 + x^2 + 1
76716889654539547639031458229653027958

10秒くらい待つと離散対数が計算できたので元のプログラムを改変して復号した。

gistedf58b1bf3484af727d98a7ebda023ff

$ python exploit.py
['\r\xa8\xe9\xe8J\x99\xd2M\x0fx\x8cqn\xf9\xe9\x9c', '\xc4G\xc3\xcf\x12\xc7\x16 m\xee\x92\xb9\xceY\x1d\xc0', 'r-BF)\x18b\x11 \xec\xe6\x8a\xc6NI:', 'A\xea:p\xdd\x7f\xe2\xb1\xd1\x16\xacH\xf0\x8d\xbf+', '&\xbdc\x83O\xa5\xb4\xcbu\xe3\xc6\rIg`\x92', '\x1b\x91\xdf^^c\x1e\x8e\x9eP\xc9\xd8\x03P$\x9c']
['One-Time Pad is ', 'used here. You w', "on't know that t", 'he flag is flag{']
One-Time Pad is used here. You won't know that the flag is flag{LCG1sN3ver5aFe!!}.��Y���+�]j�
 

flag文字列が自分の解法とあまり一致していない気がするので離散対数をもとめずとも計算できるのかもしれない。

English Version

Task

In the provided zip file, there are an encryption script 'oneTimePad2.py' and 'ciphertext'.
In this time, the crypt system is built on \(GF(2^{128})\)

\begin{align}
C_i &= M_i + R_i \\
\begin{bmatrix}
R_{i+1} \\
1
\end{bmatrix} &=
\begin{bmatrix}
A & B\\
0 & 1\\
\end{bmatrix}^{N_i}
\begin{bmatrix}
R_{i} \\
1
\end{bmatrix} \\
N_{i+1} &= N_i^2\\
\end{align}
We know \(C_1 \dots, C_5, M_1, \dots, M_4, A, B \).
The flag must be in \(M_5\).
We want to detect the value of \(N_1\) in order to break this system.

Analysis

First, let's compute the matrix power.
\begin{align}
\begin{bmatrix}
A & B\\
0 & 1\\
\end{bmatrix}^{n} &=
\begin{bmatrix}
A^n & A^{n-1} B + A^{n-2} B + \dots + B \\
0 & 1
\end{bmatrix} \\
&=
\begin{bmatrix}
A^n & \frac{(A^{n} - 1)B}{A - 1}\\
0 & 1
\end{bmatrix}
\end{align}
Therefore,
\[
R_{2} = A^{N_1} R_{1} + \frac{(A^{N_1} - 1)B}{A - 1}
\]
By solving this equation for \(A^{N_1}\),
\[
A^{N_1} = (R_2 + B (A - 1)^{-1})(R_1 + B(A-1)^{-1})^{-1}
\]
Then, we can conpute the value of \(A^{N_1}\).
Task to compute the value of \(N_1\) from \(A^{N_1}\) is called 'discrete logarithm problem', which is known to be difficult for general fields.
Fortunately, in this field \(GF(2^{128})\), the discrete logarithm is computable in a reasonable time.

Solution

We compute the descrete logarithm with PARI/GP.

gistc2a4481d5f18d35e913503664320a016

$ gp -q doit.gp  
x^123 + x^117 + x^116 + x^114 + x^113 + x^112 + x^110 + x^109 + x^108 + x^107 + x^104 + x^99 + x^98 + x^97 + x^95 + x^93 + x^91 + x^89 + x^87 + x^84 + x^83 + x^82 + x^81 + x^80 + x^78 + x^74 + x^69 + x^68 + x^66 + x^62 + x^57 + x^53 + x^52 + x^51 + x^45 + x^43 + x^42 + x^41 + x^38 + x^37 + x^35 + x^33 + x^32 + x^31 + x^30 + x^29 + x^28 + x^25 + x^22 + x^19 + x^18 + x^17 + x^15 + x^14 + x^13 + x^12 + x^10 + x^8 + x^5 + x^2 + 1
76716889654539547639031458229653027958

It took around ten seconds before we have the discrete logarithm.
Then, we decrypted the cipher text with a modified version of the encryption program.

gistedf58b1bf3484af727d98a7ebda023ff

$ python exploit.py
['\r\xa8\xe9\xe8J\x99\xd2M\x0fx\x8cqn\xf9\xe9\x9c', '\xc4G\xc3\xcf\x12\xc7\x16 m\xee\x92\xb9\xceY\x1d\xc0', 'r-BF)\x18b\x11 \xec\xe6\x8a\xc6NI:', 'A\xea:p\xdd\x7f\xe2\xb1\xd1\x16\xacH\xf0\x8d\xbf+', '&\xbdc\x83O\xa5\xb4\xcbu\xe3\xc6\rIg`\x92', '\x1b\x91\xdf^^c\x1e\x8e\x9eP\xc9\xd8\x03P$\x9c']
['One-Time Pad is ', 'used here. You w', "on't know that t", 'he flag is flag{']
One-Time Pad is used here. You won't know that the flag is flag{LCG1sN3ver5aFe!!}.��Y���+�]j�

Now, we have the flag, but there may be another solution which does not require to compute discrete logarithms, since our solution seems to be different from what the flag suggests.

0CTF Writeup: oneTimePad1

これは0CTFのoneTimePad1という問題のWriteupです。

Japanese Version

Task

zipファイルを開くと暗号化スクリプトoneTimePad.pyと暗号文ciphertextがある。
暗号化の仕組みはブロック暗号で\(GF(2^{256})\)上で次のように計算する。

\begin{align}
C_1 &= M_1 + R_1 \\
C_2 &= M_2 + R_2 \\
C_3 &= M_3 + R_3
\end{align}
ここで\(C_1, C_2, C_3\)が暗号文、 \(M_1, M_2, M_3\)が平文となる。
\(C_1, C_2, C_3, M_2, M_3\)は既知で\(M_1\)を復号するのが問題のようだ。

解析

もちろん、暗号化の肝は\(R_1, R_2, R_3\)という擬似乱数生成器の生成した数列にある。
擬似乱数生成器は、秘密鍵\(K\)と現在の乱数\(R_i\)から次の乱数\(R_{i+1}\)を次のように生成する。
\[
R_{i+1} = (R_{i} + K)^2
\]
今\(GF(2^{256})\)上で考えているので、\(2 = 0, X + Y = X - Y \)が成り立つことに注意すると、
\[
R_{i+1} = R_{i}^2 + K^2
\]
となる。\(C_2, C_3, M_2, M_3\)はわかっているので、\(K^2 = R_2 + R_1^2\)の右辺は計算できる。
これを\(R_1 = R_0^2 + K^2\)に代入すると
\[
R_0 = \sqrt{ R_1 + R_2 + R_1^2 }
\]
である。右辺の平方根は平方剰余というやつでPARI/GPというツールを使うと計算できる。

解法


gist8bcd0ce53f61a87c8742e697917cfae1

$ gp -q doit.gp 
52553756545437879763013576434186236591362916937163148104579268147006553551449
$ python
>>> hex(52553756545437879763013576434186236591362916937163148104579268147006553551449)[2:-1].decode('hex')
't0_B3_r4ndoM_en0Ugh_1s_nec3s5arY'

English Version

Task

There are an encryption script 'oneTimePad.py' and 'ciphertext' in the provided zip file.

The crypto system is described over \(GF(2^{256})\)
\begin{align}
C_1 &= M_1 + R_1 \\
C_2 &= M_2 + R_2 \\
C_3 &= M_3 + R_3
\end{align}
Here, \(C_1, C_2, C_3\) are cipher texts and \(M_1, M_2, M_3\) are plain texts.
We know \(C_1, C_2, C_3, M_2, M_3\), and the task of this problem is to decrypt \(C_1\).

Analysis

Of course, the clue of this crypto system is pseudo random numbers \(R_1, R_2, R_3\).
The random number generator generates the next number \(R_{i+1}\) from Secret key \(K\)
and the current number \(R_i\) in the following way:
\[
R_{i+1} = (R_{i} + K)^2
\]
We note that equations like \(2 = 0, X + Y = X - Y\) hold because the field is \(GF(2^{256})\).
Then, we have,
\[
R_{i+1} = R_{i}^2 + K^2
\]
The rhs of \(K^2 = R_2 + R_1^2\) is computable since we know \(C_2, C_3, M_2, M_3\).
By substituting it to \(R_1 = R_0^2 + K^2\), we have
\[
R_0 = \sqrt{ R_1 + R_2 + R_1^2 }
\]
The square root in its rhs is so called quadratic residue and it is computable with
tools like PARI/GP.

Solution


gist8bcd0ce53f61a87c8742e697917cfae1

$ gp -q doit.gp 
52553756545437879763013576434186236591362916937163148104579268147006553551449
$ python
>>> hex(52553756545437879763013576434186236591362916937163148104579268147006553551449)[2:-1].decode('hex')
't0_B3_r4ndoM_en0Ugh_1s_nec3s5arY'

State Monadと正格性について

今回の話題は、State MonadはStrictなものを使おうという話である。

StateT s m a

HaskellでState Monadを使う際、一番よく使うのは
transformer packageのStateT s m aだと思う。
transformers: Concrete functor and monad transformers
あるいは、そのwrapper libraryとしてmtlを使うかもしれない。
mtl: Monad classes, using functional dependencies
どちらにしろ実体は同じだ。
このStateTにはLazyなバージョンとStrictなバージョンの2種類がある。
二つともデータ型の定義は同じである。

newtype StateT s m a = StateT { runStateT :: s -> m (a, s ) }

違いはMonad (StateT s m)のインスタンス宣言である。

-- Strict
instance Monad m => Monad (StateT s m) where
    return x = StateT (\s -> return (x, s))
    StateT mv >>= f = StateT $ \s -> do
        (v,s1) <- mv s
        runStateT (f v) s1
-- Lazy
instance Monad m => Monad (StateT s m) where
    return x = StateT (\s -> return (x, s))
    StateT mv >>= f = StateT $ \s -> do
        ~(v,s1) <- mv s
        runStateT (f v) s1

Lazyの方は~(チルダ)パターンという見慣れないものが使われているが、これは概ね以下と同じだ。

    p <- mv s
    let (v,s1) = p

違いは、mv sの結果をバインドするときに、タプルのWHNFまで簡約するか、あるいは遅延させるかという点である。

正格性

ここでState MonadのExampleとして、状態を使って1からnまでの和を求めるプログラムを書いてみよう。

-- Main.hs
import Data.Functor.Identity
import Control.Monad.Trans.State.Lazy (StateT(..), get, put)
main :: IO ()
main = do
  n <- readLn
  let (_,s) = runIdentity (runStateT (sumState n) 0)
  print s

sumState :: Monad m => Int -> StateT Int m ()
sumState n = mapM_ (\i -> modify (+i)) [1..n]

modify :: Monad m => (s -> s) -> StateT s m ()
modify f = do
    v <- get
    put (f v)

しかし、このプログラムをコンパイルして実行してみると非常に遅いことがわかる。

$ ghc -O Main.hs
$ echo 10000000 | time ./Main
50000005000000
        2.27 real         1.41 user         0.72 sys

これはなぜかというと、modify fが正格でないため、Stateに(...((0 + 1) + 2) + ... + n)というthunkが積まれてしまうからである。
この現象はControl.Monad.Trans.State.LazyとControl.Monad.Trans.State.Strictの両方で発生する。
状態を正格に評価するために、modifyの定義を少し変更する。

modify f = do
    v <- get
    put $! f v

($!)は正格な関数適用演算子である。こうすることで、modify (+i)の部分が実行されるたびにStateを正格評価するように思える。
実際、この関数はControl.Monad.Trans.State.(Strict | Lazy).modfy'としてライブラリにも定義されている。

実際、Control.Monad.Trans.State.Strict.StateTではこの変更で期待通りに動く。

$ echo 10000000 | time ./Main
50000005000000
        0.02 real         0.01 user         0.00 sys

しかし、残念ながら、Control.Monad.Trans.State.Lazy.StateTでは上手くいかない。

echo 10000000 | time ./Main
50000005000000
        1.94 real         1.27 user         0.50 sys

また、奇妙なことだが、LazyなStateTでも、baseモナドをIdentityからIOに変えたりすると状態が正格に計算される。

import Control.Monad.Trans.State.Lazy (StateT(..), get, put)
main :: IO ()
main = do
  n <- readLn
  (_,s) <- runStateT (sumState n) 0
  print s
$ echo 10000000 | time ./Main 
50000005000000
        0.02 real         0.01 user         0.00 sys

謎解き

これはなぜだろうか、この謎を解くためにはmodify fがどのようにコンパイルされるかをみていく必要がある。

 modify f 
== get >>= (\v -> put $! (f v)) -- unfold "modify"
== StateT (\s -> return (s,s)) >>= (\v -> put $! (f v))  -- unfold "get"
== StateT (\s -> do  -- unfold (>>=)
            ~(a,s1) <- return (s,s)
            runStateT ((\v -> put $! (f v)) a) s1)
== StateT (\s -> do  -- beta reduction
            ~(a,s1) <- return (s,s)
            runStateT (put $! (f a)) s1)
== StateT (\s -> do -- unfold ($!) and "put"
            ~(a,s1) <- return (s,s)
            runStateT (let !v = f a in StateT (\_ -> return ((),v))) s1)
== StateT (\s -> do -- simplify
            ~(a,s1) <- return (s,s)
            (let !v = f a in return ((),v)))
== StateT (\s -> return (s,s) >>= (\p -> -- desugar do notation
                 let (a, s1) = p in
                 let !v = f a in
                 return ((),v)))
== StateT (\s -> let (a, s1) = (s,s) in -- apply monad law (left identity)
                 let !v = f a in
                 return ((),v)))
== StateT (\s -> let !v = f s in return ((),v)) -- simplify

今の所、特に問題のある部分はない。
しかし、baseモナドによっては問題が明らかになる。
その前にIdentity Monadの定義を復習しておく。

newtype Identity a = Identity { runIdentity :: a }
instance Monad Identity where
  return x = Identity x
  m >>= f = f (runIdentity a)

さて、次の項をIdentityモナド上で簡約してみよう。

(modify f) >> action :: StateT s Identity a
== StateT (\s -> let !v = f s in return ((),v)) >> action 
== StateT (\s -> let !v = f s in return ((),v)) >>= (\_ -> action) -- unfold (>>)
== StateT (\s -> do  -- unfold (>>=)
            ~(_,s1) <- let !v = f s in return ((),v)
            runStateT action s1)
== StateT (\s -> -- desugar do notation
           (let !v = f s in return ((),v)) >>= (\p ->
           let (_,s1) = p in
           runStateT action s1))
== StateT (\s -> -- unfold (>>=)
           (\p ->
           let (_,s1) = p in
           runStateT action s1) (let !v = f s in return ((),v)))
== StateT (\s -> -- beta reduction
           let (_,s1) = let !v = f s in return ((),v) in
           runStateT action s1)

さて、簡約結果を見てみると、 let !v = f s in return ((),v)の部分がさらにlet (_,s1) = ... in ...でくるまれているため、
評価が遅延されていることがわかると思う。この問題はbaseモナドがIOの場合は発生しない。
なぜなら、IO a は RealWorld -> (# RealWorld, a #)という型として評価されるわけだが、(# RealWolrd, a #)というタプルはthunkになり得ないからだ。

ちなみに、Control.Monad.Trans.State.Strictの場合は、以下のようになる。

(modify f) >> action :: StateT s Identity a
== StateT (\s -> let !v = f s in return ((),v)) >> action
== StateT (\s -> let !v = f s in return ((),v)) >>= (\_ -> action)
== StateT (\s -> do
            (_,s1) <- let !v = f s in return ((),v)
            runStateT action s1)
== StateT (\s -> 
           (let !v = f s in return ((),v)) >>= (\(_,s1) ->
           runStateT action s1))
== StateT (\s -> 
           (let !v = f s in Identity ((),v)) >>= (\(_,s1) ->
           runStateT action s1))
== StateT (\s ->
           (\(_,s1) ->
           runStateT action s1) (let !v = f s in Identity ((),v)))
== StateT (\s ->
           case (let !v = f s in Identity ((),v)) of
             (_,s1) -> runStateT action s1)
== StateT (\s ->
           case f s of
             v -> case Identity ((),v) of
               (_,s1) -> runStateT action s1)
== StateT (\s ->
           case f s of
             v -> runStateT action v)

従ってf sの評価はactionの評価より前に行われる。

結論

LazyなStateモナドでmodify'を使うべきではない。使う必要が出てきた場合にはStrictなStateモナドを使おう。

普通はStrictなものを使っておけば良いと思う。LazyなStateモナドが役に立つ場面は本当にあるのだろうか。甚だ疑問である。

僕の考えたさいきょうの抽象構文木データ型

あらすじ

プログラミング言語処理系を作成しようとすると避けては通れないのが、構文木データ型の設計である。
言語処理系では構文解析、アルファ変換、脱糖、正規化など構文木を変換するパスがいくつか存在して、それらの構文ごとにデータ型を設計しなければならない。今回はHaskellの依存型機能を使って、拡張のしやすい構文木データ型の構築方法を示す。

はじめに

今回の記事では単純な型なし関数型言語インタープリタを実装していく。コード全体はgithubに上げている。
github.com

今回のお題

今回は次のような構文の関数型言語を題材にする。評価戦略は値呼びである。

<expr> ::= <id> | <int_literal> | <bool_literal>
        |  <prefix_op> <expr> 
        |  <expr> <infix_op> <expr>
        | "let" <id> "=" <expr> "in" <expr>
        | "fun" <id> "->" <expr>
        | <expr> <expr>
        | "if" <expr> "then" <expr> "else" <expr>
        | "if" <expr> "then" <expr>
        | <expr> "&&" <expr> | <expr> "||" <expr>
<id> ::= [a-zA-Z][a-zA-Z0-9'_]*
<int_literal> ::= [0-9]+
<bool_iteral> ::= "true" | "false"
<prefix_op> ::= "+" | "-"
<infix_op> ::= "*" | "/" | "+" | "-" | "<" | ">" | "=" | "<>" | "<=" | ">="

素朴な実装

このような構文をHaskellのデータ型で表現する際、最も素朴なのは代数的データ型を用いるものだろう。
例えばこんな感じ。

data Expr = Id Id | IntL Int | BoolL Bool 
          | PrefixOp POp Expr | InfixOp IOp Expr Expr
          | Let Id Expr Expr 
          | Fun Id Expr
          | App Expr Expr
          | IfThenElse Expr Expr Expr
          | IfThen Expr Expr
          | Cond CondOp Expr Expr
data CondOp = And| Or
type Id = String
data POp = Plus | Minus
data IOp = Mul | Div | Add | Sub | ...

ただし、この定義には二つの問題点がある。

  • どの構文ノードにも共通して存在するデータを埋め込むのが面倒、例えば、構文木ノードにユニークなラベルを貼ったり、ソースコード上での位置を表すデータなどのメタデータを埋め込みたいとき、全てのコンストラクタに書かなければならない。
  • 拡張性に乏しい。例えば"e0 && e1"は"if e0 then e1 else false"の糖衣構文であるので、脱糖後のデータ構造にはそのような構文は許されないが、脱糖後のデータ型でそのことを表現するためにはほぼ同じ代数的データ型を定義する必要がある。コンストラクタの名前を考えるのが面倒だし、名前空間が汚くなる。

僕の考えたさいきょうの抽象構文木データ型

ここで本記事の目的であるデータ型を紹介しよう。
まず、DataKinds拡張を使って構文木のコンストラクタのラベルを表すConstカインドと
そのシングルトン型SConstをGADTとして定義する。

{-# Language DataKinds, GADTs #-}
module Language.ToyML.Syntax.Base where
import Data.Void(Void)
import Data.Kind(Constraint)

data Const = Var | Literal | Infix | Prefix 
           | App | Abs | Let  
           | IfT | IfTE
           | Cond

data SConst (c :: Const) where
    SVar     :: SConst 'Var
    SLiteral :: SConst 'Literal
    SInfix   :: SConst 'Infix
    SPrefix  :: SConst 'Prefix
    SApp     :: SConst 'App
    SAbs     :: SConst 'Abs
    SLet     :: SConst 'Let
    SIfT     :: SConst 'IfT
    SIfTE    :: SConst 'IfTE
    SCond    :: SConst 'Cond

そして構文木のデータ型はこれだ。

data Exp where
    Exp :: WellFormed c Exp arg => SConst c -> arg -> Exp

つまり、構文木はExp op argという二つ組でopはSConst c型、つまりコンストラクタの種類を表している、
一方argの型はこの定義からはわからないが、WellFormed c Exp argという制約がargの型を定めてくれる。
カギはWellFormed c e argという制約にあるわけだが、これはtype familyという仕組みで次のように定義している。

type family WellFormed (c :: Const) e arg :: Constraint where
    WellFormed 'Var e arg     = arg ~ Ident e
    WellFormed 'Literal e arg = arg ~ Literal e
    WellFormed 'Infix e arg   = arg ~ (InfixOp e,e,e)
    WellFormed 'Prefix e arg  = arg ~ (PrefixOp e, e)
    WellFormed 'App e arg     = arg ~ (e,e)
    WellFormed 'Abs e arg     = arg ~ (Ident e, e)
    WellFormed 'Let e arg     = arg ~ (Ident e, e, e)
    WellFormed 'IfT e arg     = arg ~ (e, e)
    WellFormed 'IfTE e arg    = arg ~ (e, e, e)
    WellFormed 'Cond e arg    = arg ~ (CondOp, e, e)

type family Ident e
type family Literal e
type family InfixOp e
type family PrefixOp e

この定義によると、例えばExp SApp argの場合、argは(Exp, Exp)型となり、またExp SVar argの場合はargはIdent Exp型になる。
Ident Exp型がどのような型になるかもtype familyによって定めており、具体的にはStringになるように次のように定めている。

type instance Ident Exp = String
type instance Literal Exp = Lit
type instance InfixOp Exp = Op
type instance PrefixOp Exp = Op
data Lit = CInt Integer | CBool Bool
data Op = Op { headChar :: Char, opName :: String }

この定義に従うと構文解析後のデータとして"正しい"データ構造を表現できている。例えば次のような例があり、間違ったデータ構造はコンパイルエラーとなる。

exp1 :: Exp
exp1 = Exp SVar "x"

exp2 :: Exp
exp2 = Exp SAbs ("x", Exp SApp (Exp SVar "x", Exp SVar "x")) -- fun x -> x x

exp3 :: Op -> Exp
exp3 op = Exp SPrefix (op, Exp SVar "x") 


-- これはコンパイルエラー(前置演算子なのに二つ子のExpを持つため)
-- exp3bug :: Op -> Exp
-- exp3bug op = Exp SPrefix (op, Exp SVar "x", Exp SVar "x")

このパターンで構文木データ型を定義すると共通するデータを埋め込むのも容易である。
例えば、ソースコード上の位置を埋め込みたければ

data Exp where
    Exp :: WellFormed c Exp arg => SConst c -> arg -> SourcePos -> Exp

と変えるだけで良い。

次に拡張も簡単であることを示すために脱糖後のデータ型をこのパターンを用いて定義してみよう。

data Exp where
    Exp :: (WellFormed c Exp arg, Desugared c Exp arg) => 
           SConst c -> arg -> SourcePos -> Exp

type family Desugared (c :: Const) e arg :: Constraint where
    Desugared 'Infix e arg = Impossible 'Infix
    Desugared 'Prefix e arg = Impossible 'Prefix
    Desugared 'IfT e arg = Impossible 'IfT
    Desugared 'Cond e arg = Impossible 'Cond
    Desugared c e arg = ()

これだけである。
先ほど定義したWellFormed c e argという制約に加えて脱糖後にはInfix, Prefix, IfT, Condのコンストラクタが使えないことをDesugared c e argという制約で表現している。
ここでImpossibleという制約があるわけだが、この制約はこのように定義する。

import Data.Void
class Impossible (c :: Const) where
    impossible :: SConst c -> Void

つまりImpossibleは型クラスであり、そのインスタンスは存在しない。したがって、Impossible制約を持つコンストラクタを使うことはできないようになっている。逆にパターンマッチでこのパターンが来たらData.Void.absurd :: Void -> aという最強関数で安全に型を合わせてしまうことができる。

終わりに

今回紹介したパターンを使えばこのように柔軟な型コンストラクタをHaskellでも実装することができる。パターンマッチの部分もなかなか面白いので興味のある人はレポジトリのコードを読んでみてほしい。
github.com

OS Xにおける共有ライブラリについてのメモ

最近Z3のインストールに複数の意味でハマっていて、その過程で動的ライブラリに対する理解が深まったのでメモしておく。
autotaker.hatenablog.com

動的ライブラリとは

動的ライブラリは静的ライブラリと異なり、実行時にリンクされる。
今まで誤解していたのだが、ライブラリを使用する側からはライブラリが静的か動的かは関係がない。
コンパイラソースコードをオブジェクトファイルにする段階ではライブラリのヘッダファイルさえあればよくて、リンカが実行形式ファイルを出力する段階で初めてシンボルを静的に解決するか動的に解決するかが決まる。

覚えておくと便利なコマンド

環境はOS Xを想定している。

  • 実行形式ファイルあるいは動的ライブラリからそれが参照する動的ライブラリの一覧を取得する。
otool -L ファイル
  • 実行形式ファイルあるいは動的ライブラリからload commandsの一覧を取得する。
otool -l ファイル
  • 動的ライブラリがexportしている関数の一覧を取得する。
gobjdump -T ファイル

動的リンクの仕組み

ところで、動的リンカはどのようにして実行時に動的ライブラリの実体を見つけてくるのだろうか?

例えば、手元にある適当な動的ライブラリが使用する動的ライブラリを見てみよう。

/Users/autotaker/.cabal/lib/x86_64-osx-ghc-7.10.2/z3-4.1.0-CQEXZ6rXqZX8ei09CJh4x5/libHSz3-4.1.0-CQEXZ6rXqZX8ei09CJh4x5-ghc7.10.2.dylib:
	@rpath/libHSz3-4.1.0-CQEXZ6rXqZX8ei09CJh4x5-ghc7.10.2.dylib (compatibility version 0.0.0, current version 0.0.0)
	libz3.dylib (compatibility version 0.0.0, current version 0.0.0)
	@rpath/libHSmtl-2.2.1-2BzSpTumj4ZLLPrpLUWDdr-ghc7.10.2.dylib (compatibility version 0.0.0, current version 0.0.0)
	@rpath/libHStransformers-0.4.2.0-3eG64VdP2vzGjP6wJiCp5X-ghc7.10.2.dylib (compatibility version 0.0.0, current version 0.0.0)
	@rpath/libHScontainers-0.5.6.2-LKCPrTJwOTOLk4OU37YmeN-ghc7.10.2.dylib (compatibility version 0.0.0, current version 0.0.0)
	@rpath/libHSdeepseq-1.4.1.1-LbCWUlehDDeLxurARKDH5o-ghc7.10.2.dylib (compatibility version 0.0.0, current version 0.0.0)
	@rpath/libHSarray-0.5.1.0-E0sTtauuKsGDLZoT7lTbgZ-ghc7.10.2.dylib (compatibility version 0.0.0, current version 0.0.0)
	@rpath/libHSbase-4.8.1.0-GDytRqRVSUX7zckgKqJjgw-ghc7.10.2.dylib (compatibility version 0.0.0, current version 0.0.0)
	@rpath/libHSinteger-gmp-1.0.0.0-2aU3IZNMF9a7mQ0OzsZ0dS-ghc7.10.2.dylib (compatibility version 0.0.0, current version 0.0.0)
	@rpath/libHSghc-prim-0.4.0.0-8TmvWUcS1U1IKHT0levwg3-ghc7.10.2.dylib (compatibility version 0.0.0, current version 0.0.0)
	/usr/lib/libiconv.2.dylib (compatibility version 7.0.0, current version 7.0.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 1226.10.1)

このライブラリではライブラリの指定の仕方に以下の3種類が用いられている。(一般にはもう何種類かあるらしいが説明は割愛)

  • 絶対パス(例:/usr/lib/libiconv.2.dylib )
  • 相対パス(例:libz3.dylib)
  • @rpathで始まるパス(例:@rpath/libHSz3-4.1.0-CQEXZ6rXqZX8ei09CJh4x5-ghc7.10.2.dylib)

この書き方にしたがって実行時に動的ライブラリを探してくる仕組みが異なる。

絶対パスの場合

この場合は特に難しいことはなくそのパスで指定された動的ライブラリをリンクする。

相対パスの場合

この場合はまず、DYLD_LIBRARY_PATHという環境変数に指定されたパスからライブラリを探索する。
そこで見つからなかった場合は、DYLD_FALLBACK_LIBRARY_PATHを探索する。DYLD_FALLBACK_PATHのデフォルト値は$(HOME)/lib:/usr/local/lib/:/lib/:/usr/libとなっている。

@rpathの場合

実行ファイルおよび動的ライブラリにはrpathと呼ばれる、実行時に動的ライブラリを探索するパスのリストが埋め込まれている。動的リンカはそのリストに含まれるパスから動的ライブラリを探索する。

rpathの追加

ldコマンドに-rpath というオプションを渡すと、実行ファイルにそのパスがrpathとして埋め込まれる。
clangでコンパイルする場合は-Wl,-rpath -Wl,とする。

install name

そもそも上記の3種類の指定はどのように決定されるのか?
その答えがinstall nameと呼ばれるものである。動的ライブラリを生成するときにldに-install_name というオプションを渡すことで指定できる。あるいは、install_name_toolというコマンドで生成された動的ライブラリのinstall nameを後から変更することができる。
ldが実行ファイルを生成するとき、動的ライブラリを見つけてきて、そのライブラリのinstall nameを実行ファイルの動的ライブラリのパスとして埋め込む。

以上調べたことをまとめて、実際にいろいろ試せるように共有ライブラリのビルドを行うレポジトリをgithubに公開しておいた。
github.com

OS Xにおける動的リンクのまとめは以下のサイトを参考にした。
mikeash.com: Friday Q&A 2009-11-06: Linking and Install Names

MacにZ3をインストールした。

新しいMacを手に入れたので環境構築を行っている。
その過程で、Z3のインストールにハマったので忘備録を書いておく。

目標

Z3はMicrosoftが開発しているSMTソルバで、様々な言語のバインディングがある。
公式でサポートしているのはC/C++, Java, Python, OCaml等でHaskellでは非公式のラッパーライブラリを用いる。
github.com
z3: Bindings for the Z3 Theorem Prover

今回はOCamlHaskell両方のバインディングをインストールしたい。

トラブル

まず、z3をmlバインディングつきでインストールする。

$ git clone https://github.com/Z3Prover/z3.git
$ cd z3
$ scripts/mk_make.py --ml # --mlはmlバインディングをインストールするフラグ
$ make -C build
$ make -C build install

これでlibz3.dylibが/usr/local/lib以下にインストールされ、mlバインディングのライブラリは~/.opam/`opam switch show`/lib/Z3以下にインストールされる。

次にZ3のHaskellバインディングライブラリをインストールしようとして失敗した。

$ cabal install z3 -f examples
=> ERROR!

これをデバッグしていくと、libz3.dylibに_Z3_get_error_msg_exという関数がないことが原因だとわかった。

$ gobjdump -T /usr/local/lib/libz3.dylib | grep _get_error
0000000000028af0 g       1e SECT   01 0000 [.text] __Z22exec_Z3_get_error_codeR11z3_replayer
0000000000028b40 g       1e SECT   01 0000 [.text] __Z21exec_Z3_get_error_msgR11z3_replayer
0000000000069880 g       1e SECT   01 0000 [.text] __Z21log_Z3_get_error_codeP11_Z3_context
0000000000069a10 g       1e SECT   01 0000 [.text] __Z20log_Z3_get_error_msgP11_Z3_context13Z3_error_code
0000000000034060 g       0f SECT   01 0000 [.text] _Z3_get_error_code
00000000000340f0 g       0f SECT   01 0000 [.text] _Z3_get_error_msg
$

原因究明

Z3のレポジトリを読んでいつこの関数が消えたのかということを調べたら、去年の11月ごろに以下のコミット周辺で消されていることが判明した。
github.com
このコミット以前では次のようなAPIになっていた。

Z3_API char const * Z3_get_error_msg_ex(Z3_context c, Z3_error_code err);
Z3_API char const * Z3_get_error_msg(Z3_error_code err);

実装を読むとZ3_get_error_msgはZ3_get_error_msg_ex(0, err)と同じ振る舞いをするようだ。
ずいぶん前からZ3_get_error_msg自体はdeprecatedになっており、代わりにZ3_get_error_msg_exを使えということになっていた。

コミットの以後はつぎのようなAPIになっている。

Z3_API char const * Z3_get_error_msg(Z3_context c, Z3_error_code err);

つまり、長らく負の遺産となっていたZ3_get_error_msgをまともなAPIに直したということだろう。
しかし、このとき同時にZ3_get_error_msg_exも消し去ってしまうのは後方互換性の考え方からしてちょっと配慮が足りないのではないだろうか。当分の間はZ3_get_error_msg_exはdeprecatedのまま残しておいてwrapperライブラリ開発者が対応する時間を与えるべきだったと思う。

解決策

z3-haskellの開発者にプルリクを送ろうかとも考えたが、z3-haskellmercurialで管理されていたのでしばらくプルリクの送り方を勉強してからにする予定である。
当分は問題の少なそうなコミットを適当に探した結果、c2108f74f1e5c5ab96c6e4cc189ccbc55dca339aをチェックアウトしてからビルドすれば問題なさそうである。
そもそもz3のバージョンが変わらないのにAPIの仕様が変わっているのがすごくよくないと思うので、z3にはstableなブランチを是非とも作って欲しいものだ。

5/30 21:30 追記

そもそも、masterブランチで開発されているのはz3-4.4.2でありRelease前のバージョンだと気付いた。stableなバージョンも普通にgithubのreleasesにあった。
Releases · Z3Prover/z3 · GitHub
master HEADは極めてunstableでないとわかったので大人しく最新releaseであるz3-4.4.1をインストールするのが無難だろう。

$ git clone https://github.com/Z3Prover/z3.git
$ cd z3
$ git checkout refs/tags/z3-4.4.1
$ python scripts/make_mk.py --ml --prefix=/usr/local
$ make -C build
$ make -C build install

vectorを使ったData.List.sortより4倍速いsortアルゴリズムの実装

Data.List.sortがあまりに遅くてつらいので、vectorを使って書いてチューニングしたら約4倍速くなりましたという話をします。その過程でvectorのmonadic indexingとは何かという話をします。

仕様としてData.Listと互換性を持たせるため、次のようなインターフェースにしました。

{-# INLINE sort #-}
sort :: Ord a => [a] -> [a]
sort = sortBy compare
sortBy :: (a -> a -> Ordering) -> [a] -> [a]
sortBy = ...

次に実装を見てきます。

import qualified Data.Vector as V
...
sortBy cmp = V.toList . mergeSortAux . V.fromList
    where 
    mergeSortAux l = ...

入力はリストで与えられるので、まず最初にfromList :: [a] -> Vector aを使ってData.Vector.Vectorに変換し、mergeSortAux :: Vector a -> Vector aを使ってソートした後に、toList :: Vector a -> [a]を用いてリストに戻します。Unboxed Vectorを使うとさらに速くなりますが、ソートできるデータが制限されてしまうので今回はBoxed Vectorを使います。

mergeSortAuxの実装はただのマージソートです。今後ところどころunsafeなんとかを使っていますが、これは今回のアルゴリズムでは境界値外アクセスすることはないので配列の境界チェックを省くために行っているだけです。怖くないよ。

import qualified Data.Vector as V
...
sortBy cmp = ...
    where 
    mergeSortAux l
        | n <= 1 = l
        | otherwise = merge (n1, l1') (n2, l2')
             where
                n  = V.length l
                n1 = div n 2
                n2 = n - n1
                l1 = V.unsafeSlice 0  n1 l
                l2 = V.unsafeSlice n1 n2 l
                l1' = mergeSortAux l1
                l2' = mergeSortAux l2
    merge = ... 

さて、最後にmergeの実装を見てきましょう。

{-# LANGUAGE BangPatterns #-}
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as MV

sortBy cmp = ...
    where
    ....
    merge (!n1,!as) (!n2,!bs) = V.create $ do
        res <- MV.unsafeNew (n1+n2)
        let go i1 i2 = ...
        go 0 0

V.create :: (forall s. ST s (MVector s a)) -> Vector aを使ってMVector s aをVectorに変換します。
do構文の中ではマージ後の配列を格納するres :: MVector s aを確保し、go関数で先頭から順番に値を埋めていく感じです。

        let go i1 i2 | i1 >= n1 && i2 >= n2 = return res
                     | i1 >= n1 = tick2 i1 i2
                     | i2 >= n2 = tick1 i1 i2
                     | otherwise = do
                         v1 <- V.unsafeIndexM as i1
                         v2 <- V.unsafeIndexM bs i2
                         case cmp v1 v2 of
                            GT  -> do
                                 MV.unsafeWrite res (i1 + i2) v2
                                 go i1 (i2 + 1)
                            _   -> do
                                 MV.unsafeWrite res (i1 + i2) v1
                                 go (i1 + 1) i2
            tick1 i1 i2 = do
                V.unsafeIndexM as i1 >>= MV.unsafeWrite res (i1 + i2)
                go (i1 + 1) i2
            tick2 i1 i2 = do
                V.unsafeIndexM bs i2 >>= MV.unsafeWrite res (i1 + i2)
                go i1 (i2 + 1)

goの中でV.unsafeIndexMを使っていますが、これの意味について説明します。

indexMについて

例として以下のようなIndexTest.hsを考えます。

gist4eb6a528eed4ef665e9bcf1d91441094

$ ghc -O2 IndexTest.hs
$ ./IndexTest
test1: safe!
test2: safe!
IndexTest: Prelude.undefined

実行するとtest3でundefinedエラーが発生します。

コンパイルの中間コードをdumpして何が起こっているのかを読みます。

$ touch IndexTest.hs
$ ghc -ddump-simpl -O2 IndexTest.hs > IndexTest.log

ghc -ddump-simpl -O2 IndexTest.hs (GHC 7.10.2)

test1のコンパイル結果を見ましょう。

a1_r5Mg
  :: V.Vector Int
     -> GHC.Prim.State# GHC.Prim.RealWorld
     -> (# GHC.Prim.State# GHC.Prim.RealWorld, () #)
[GblId, Arity=2, Str=DmdType <S,1*H><L,U>]
a1_r5Mg =
  \ (v_a2pg :: V.Vector Int)
    (eta_B1 [OS=OneShot] :: GHC.Prim.State# GHC.Prim.RealWorld) ->
    case v_a2pg
    of _ [Occ=Dead]
    { Data.Vector.Vector ipv_s2CZ ipv1_s2Db ipv2_s2Dc ->
    ((check_r2i0
        lvl2_r5Mf
        (case GHC.Prim.indexArray# @ Int ipv2_s2Dc ipv_s2CZ
         of _ [Occ=Dead] { (# ipv3_a5vA #) ->
         ipv3_a5vA
         }))
     `cast` ...)
      eta_B1
    }

test1 [InlPrag=NOINLINE] :: V.Vector Int -> IO ()
[GblId, Arity=2, Str=DmdType <S,1*H><L,U>]
test1 =
  a1_r5Mg
  `cast` ...

ごちゃごちゃしてますが重要なことはcheck関数呼び出し時に

(check_r2i0
        lvl2_r5Mf
        (case GHC.Prim.indexArray# @ Int ipv2_s2Dc ipv_s2CZ
         of _ [Occ=Dead] { (# ipv3_a5vA #) ->
         ipv3_a5vA
         }))

と(case...)の部分が遅延評価されているということです。
この部分は元のソースコードでは(V.unsafeIndex v 0)の部分に対応します。
このため、余計なオーバーヘッドが生じてしまいます。

次にtest2のコンパイル結果を見ます。

a3_r5Mk
  :: V.Vector Int
     -> GHC.Prim.State# GHC.Prim.RealWorld
     -> (# GHC.Prim.State# GHC.Prim.RealWorld, () #)
[GblId, Arity=2, Str=DmdType <S,1*U(U,A,U)><L,U>]
a3_r5Mk =
  \ (v_a2ph :: V.Vector Int)
    (eta_B1 [OS=OneShot] :: GHC.Prim.State# GHC.Prim.RealWorld) ->
    case v_a2ph
    of _ [Occ=Dead]
    { Data.Vector.Vector ipv_s2Dk ipv1_s2Dl ipv2_s2Dm ->
    case GHC.Prim.indexArray# @ Int ipv2_s2Dm ipv_s2Dk
    of _ [Occ=Dead] { (# ipv3_a5vA #) ->
    ((check_r2i0 lvl4_r5Mj ipv3_a5vA)
     `cast` ...)
      eta_B1
    }
    }

test2 [InlPrag=NOINLINE] :: V.Vector Int -> IO ()
[GblId, Arity=2, Str=DmdType <S,1*U(U,A,U)><L,U>]
test2 =
  a3_r5Mk
  `cast` ...

今度はcheckの呼び出しの前にindexArray#が呼び出されています。
しかし、配列の要素自体は評価しないのでundefinedエラーは発生しません。

最後にtest3のコンパイル結果を見ます。

a2_r5Mi
  :: V.Vector Int
     -> GHC.Prim.State# GHC.Prim.RealWorld
     -> (# GHC.Prim.State# GHC.Prim.RealWorld, () #)
[GblId, Arity=2, Str=DmdType <S,1*U(U,A,U)><L,U>]
a2_r5Mi =
  \ (v_a2pj :: V.Vector Int)
    (eta_B1 [OS=OneShot] :: GHC.Prim.State# GHC.Prim.RealWorld) ->
    case v_a2pj
    of _ [Occ=Dead]
    { Data.Vector.Vector ipv_s2Df ipv1_s2Dg ipv2_s2Dh ->
    case GHC.Prim.indexArray# @ Int ipv2_s2Dh ipv_s2Df
    of _ [Occ=Dead] { (# ipv3_a5vA #) ->
    case ipv3_a5vA of vx_a2zI { GHC.Types.I# ipv4_s3c7 ->
    ((check_r2i0 lvl3_r5Mh vx_a2zI)
     `cast` ...)
      eta_B1
    }
    }
    }

test3 [InlPrag=NOINLINE] :: V.Vector Int -> IO ()
[GblId, Arity=2, Str=DmdType <S,1*U(U,A,U)><L,U>]
test3 =
  a2_r5Mi
  `cast` ...

この場合もcheckの呼び出し前にindexArray#を呼び出していて余計なオーバーヘッドは生じないのですが、さらに配列の要素も評価してしまうのでこの例ではundefinedエラーとなってしまいます。

まとめると

  • check "test1" (V.unsafeIndex v 0)だと配列アクセスのためのサンクが作って関数呼び出し。
  • check "test3" $! (V.unsafeIndex v 0)だと配列アクセスと要素の評価の後、関数呼び出し。
  • V.unsafeIndexM v 0 >>= check "test2"だと配列アクセスの後、関数呼び出し。

ちなみにindexM系がモナド内でしか使えないのはある種のHackなのでモナドの種類は関係ないはずです。

実験

今回書いたソートのコードはgistにまとめています。
gist1b847adefb53e04eedd1b7adabfcc330

次のようなテストプログラムを書いて実験を行いました。

gistecc5317d3f07cc2890f68da58bdec379

入力として30万要素のランダムな整数列を与えてみます。

$ ./Main < test300000.in > test300000.out
[2016-04-04 04:38:12.623993 UTC] Parsing: begin
[2016-04-04 04:38:12.709319 UTC] Parsing: end
[2016-04-04 04:38:12.709578 UTC] Parsing: 0.085326s
[2016-04-04 04:38:12.709794 UTC] Vector: begin
[2016-04-04 04:38:12.976383 UTC] Vector: end
[2016-04-04 04:38:12.976596 UTC] Vector: Sorting: 0.266589s
[2016-04-04 04:38:12.976818 UTC] List: begin
[2016-04-04 04:38:14.018594 UTC] List: end
[2016-04-04 04:38:14.018809 UTC] List: Sorting: 1.041776s
[2016-04-04 04:38:14.019109 UTC] result validation begin
[2016-04-04 04:38:14.022344 UTC] Correct!
[2016-04-04 04:38:14.126167 UTC] output done
[2016-04-04 04:38:14.126446 UTC] Elapsed Time: 1.502174s
[2016-04-04 04:38:14.126715 UTC] Speedup (Data.List.sort / VecSort.sort) = 391%

Data.List.sortは1.04秒かかっていたのに対してVecSort.sortは0.267秒でおよそ4倍高速になりました。もちろん、ソートされたリストの最初の幾つかの要素が欲しい場合などではData.List.sortの方が速いこともありますが、全体をソートしたい場合にはたくさんのメモリを消費するリストを用いたソートよりも、省メモリなVectorを用いたソートを用いたほうがよいと思います。