オイラー関数の公式の確率的意味

最近、オイラー関数（Euler's totient function）について少し復習してみた。新しい発見があったのでメモしておく。

オイラー関数とは、

と素因数分解できるある自然数nに対して、

と表され、これはn未満の自然数でnと互いに素なものの数を表しています。

　上記の公式の証明は、以下3つの補題から簡単にできます。

1番目の補題は、pが素数なことから自明。2番目の補題は、p^k未満の数でpを素因数にもつものを数えれば導出できます。3番目の補題は、aと因数を共有しないもの、bと因数を共有しないものを、それぞれaとbの法を元に列挙して、中国の剰余定理を使えば証明できます。

これら3つの補題からオイラー関数の公式は導出可能ですが、詳細は省略。

面白いと思ったのは、オイラー関数を

(n未満のnと互いに素な自然数の個数)

= (1-nまでの自然数の個数)

* (1-nまで自然数のp1で割れない確率)

* (1-nまで自然数のp2で割れない確率)

.....

* (1-nまで自然数のpkで割れない確率)

と捉えることができるということです。確かに言われてみればその通りですが、この発想は面白いと思いました。

C++でハッシュのマップを実装する

C++のset/mapは、二分木で実装されている。要素数が大きくなると重くなってしまうので（平衡木の回転の処理が重たい？）、ハッシュを使ってmapみたいなものを実装してみた（挿入、参照のみ）。

template <class K, class V>
class hashMap {
    static const unsigned int DEFAULT_SIZE = 1000000;
    vector<pair<K, V> > *contents;
    unsigned int _size;
    unsigned int (*_hashValue) (K, unsigned int);

public:
    hashMap(unsigned int (*func) (K, unsigned int) , unsigned int size = DEFAULT_SIZE) {
        _hashValue = func;
        _size = size;
        contents = new vector<pair<K, V> > [_size];
    }

    ~hashMap() {
        delete [] contents;
    }

    bool hashElement(K s) {
        unsigned int h = _hashValue(s, _size);

        for (int i = 0; i < (int)contents[h].size(); i++) {
            if (contents[h][i].first == s)
                return true;
        }

        return false;
    }

    void insert(K x, V n) {
        if (hashElement(x))
            return;

        unsigned int h = _hashValue(x, _size);
        contents[h].push_back(make_pair(x, n));
    }

    pair<K, V> *getElement(K x) {
        unsigned int h = _hashValue(x, _size);
        for (int i = 0; i < (int)contents[h].size(); i++)
            if (contents[h][i].first == x)
                return &contents[h][i];

        return NULL;
    }

    void printHashBalance() {
        int cnt = 0;
        double ret = 0.0;
        int worst = 0;

        for (int i = 0; i < (int)_size; i++) {
            if (contents[i].size()) {
                ++cnt;
                ret += contents[i].size();
                worst = max<int>(worst, contents[i].size());
            }
        }

        if (!cnt)
            cerr << "Hash is empty" << endl;
        else {
            cerr << "collisions occurred(average): " << (ret/cnt) << endl;
            cerr << "collisions occurred(worst): " << worst << endl;
        }
    }
};

これで、参照/挿入はO(1)でできるはず。衝突がどれくらい起きているかを調べるためにprintHashBalance()という関数を定義した。コンストラクタには、使用するハッシュ関数のポインタと、ハッシュ値のサイズを渡す。

適当にテストコードを書いて、標準関数のmapと比較する。WORDS_NUM個のランダムな文字列を、それぞれのコンテナにぶちこんで、最頻出の文字列を求める。

#define WORDS_NUM 10000

/*
* generate random strings of length between 3 and 8 inclusive
*/
vector<string> randomGen() {
    vector<string> ret;

    for (int i = 0; i < WORDS_NUM; i++) {
        int len = 3 + rand() % 6;

        char str[len+1];
        for (int j = 0; j < len; j++)
            str[j] = 97 + rand() % 26;
        str[len] = '\0';

        ret.push_back(str);
    }

    return ret;
}
/*
* timer
*/
class Timer {
    double start_t;
public:
    void start() {
        cout << "Timer Start." << endl;
        start_t = gettimeofday_sec();
    }

    void end() {
        cout << "Timer stopped." << endl;
        cout << "Time Elapsed: " << gettimeofday_sec() - start_t << endl;
    }
};

//#define _MAP
#define _HASH

unsigned int stringHash(string str, unsigned int sz) {
    unsigned int ret = 0;

    for (int i = 0; i < (int)str.size(); i++)
        ret = 26 * ret + str[i] - 97;

    return ret % sz;
}

int main() {
    vector<string> words = randomGen();
    Timer timer;
    timer.start();

#ifdef _MAP
    map<string, int> wc;
    int cnt = 0;
    string best;

    for (int i = 0; i < (int)words.size(); i++) {
        if (wc.count(words[i])) {
            int n = ++wc[words[i]];

            if (n > cnt) {
                cnt = n;
                best = words[i];
            }
        }
        else
            wc[words[i]] = 1;
    }
    cout << cnt << ": " << best << endl;
#endif

#ifdef _HASH
    hashMap<string, int> wc(stringHash, WORDS_NUM);
    int cnt = 0;
    string best;

    for (int i = 0; i < (int)words.size(); i++) {
        if (wc.hashElement(words[i])) {
            int n = ++wc.getElement(words[i])->second;

            if (n > cnt) {
                cnt = n;
                best = words[i];
            }
        }
        else
            wc.insert(words[i], 1);
    }
    cout << cnt << ":" << best << endl;
    wc.printHashBalance();
#endif

    timer.end();

    return 0;
}

以下比較結果。HashMapの_sizeは入力文字数と同じサイズにした。一つのハッシュ値に格納される入力値の種類は平均1.45くらいだった。

WORDS_NUM	100,000	1,000,000	10,000,000
Map [s]	0.346	4.999	66.71
HashMap [s]	0.2	2.052	20.274

もっと速くなると思ったが、どこがネックなんだろう。オーダー的にはハッシュ値の計算が重いけど、実際はvectorのメモリ関連処理が遅そうな気がする。プロファイルとってみるかー。

Haskellで有名アルゴリズム実装（5）

最近英語力が落ちてる感があるので、英語で書く。

The fifth trial of implementing a classic algorithm with Haskell!

Today I'll try Pascal' Triangle.

pascalTriangle x = x : pascalTriangle y
             where y = zipWith (+) (0:x) (x++[0])

Nothing special. Thanks to lazy evaluation, you can define an infinite list of Pascal's Triangle.

Here's a similar problem:

Let's consider how many times it lands heads if you toss a coin. More specifically, what the probability of the coin landing heads k times out of n tosses?

It seems it's nothing to do with Pascal's Triangle. But it actually is.

prbTree x = x : prbTree y
       where y = zipWith (+) (map (*0.5) (0:x)) (map (*0.5) (x++[0]))

The tree structure is just like Pascal's Triangle.

In the above source, I used floating point numbers. So they don't sum up to 1.0 due to the limitation of accuracy of the type.

You can also use rational numbers, which is in Data.Ratio module, instead of floating point numbers. In that case, they sum up to 1.0 when added together.

import Data.Ratio

half = (*) (1%2)
prbTree x = x : prbTree y
          where y = zipWith (+) (map half $ 0:x) (map half $ x++[0])

Go to Haskell's page.

クイックソート VS マージソート

今日は、軽い実験。Haskellの練習がてら、quick sortとmerge sortを比較する。

まず、quick sort。

-- quick sort
qsort :: (Ord a) => [a] -> [a]
qsort []     = []
qsort (x:xs) = let left  = filter (<x) xs
                   right = filter (>=x) xs
               in qsort left ++ [x] ++ qsort right

Haskellは美しい。次に、merge sort。

msort :: (Ord a) => [a] -> [a]
msort [x] = [x]
msort [x, y] = [min x y, max x y]
msort xs = let len = length xs `div` 2
               pre = msort $ take len xs
               suf = msort $ drop len xs
           in merge pre suf

merge [] [] = []
merge x []  = x
merge [] y  = y
merge (x:xs) (y:ys)
     | x <= y    = x : merge xs (y:ys)
     | otherwise = y : merge (x:xs) ys

もうちょっと綺麗に書けそうだが、こんなもんだろう。

適当に乱数をソートしてみる。

getRandoms :: Int -> [Int]
getRandoms n = take n $ randomRs (1, 100000) $ mkStdGen 0

n=100,000で比較。

quick sort 0.653[s]

merge sort 0.891[s]

quick sortの弱点を突いてみる。ソート済みのリストを入力として与えてみる。[1..100000]で比較。

quick sort 815.877[s]

marge sort 0.465[s]

quick sortはpivot選択すれば改良はできるけど、それでも理論上の最悪計算量はO(n^2)。merge sortは常に平衡二分木を形成するように分割するので、計算量は常にO(n log n)。ランダムなデータに対してはquick sortが速い場合が多いようだが、実用上はmerge sortを使用した方が安心な気がする。

モナドを実装する

　今日はモナド実装にチャレンジ。

モナドの機能は、contextualな値に関数を適用すること。

functorやapplicativeも似たようなことを行うが、モナドのいいところは演算子のつなぎ目をきれいに書けるところ。

　とりあえず、書いてみる。まずモナドクラスのインスタンスにする型を定義する。

Maybe型を応用して、即成功、即失敗、判断を後回しにする、という3つのcontextを持つ型を定義する。

data Check a = YES | NO | Unknown a deriving (Show, Eq)

　この型を利用するために以下の関数を定義する。

passTest f x  = if f x == True then YES else Unknown x
failTest f x  = if f x == True then NO else Unknown x
determine f x = if f x == True then YES else NO

　さらに、この型をモナドインスタンスとして実装する。

instance Monad Check where
 return x = Unknown x
 YES >>= _ = YES
 NO >>= _ = NO
 Unknown x >>= f = f x

　これを使うと、ややこしい制御フローがきれいに書ける。例えば、うるう年判定をするプログラム。やたらとif-elseが多い。

isLeapYear x = if x `mod` 400 == 0
              then True
              else if x `mod` 100 == 0
                      then False
                      else if x `mod` 4 == 0
                              then True
                              else False

　モナドを使うときれいに書ける。

divisable y x = x `mod` y == 0
isLeapYear x = return x >>= passTest (divisable 400)
                      >>= failTest (divisable 100)
                      >>= determine (divisable 4)

　モナドは、制御フローを表に出さずに演算を繋げることができるデザインパターンと考えることもできる。見てのとおり制御フローをcomposableな要素として部品化することができ、if-else-break文のボイラープレート除去を行うことができる。条件が増えたら、新しいバインドで関数をつないであげれば、ifとかelseとかbreakとか書かなくてもよい。（composableうんぬんは、tanakhさんがつぶやいていた話です。多分もっと高尚なお話だったと思いますが、自分が理解した限りだとこんな話。。のはず。）

functorやapplicativeと比べて、モナドがやたらと注目されている理由は、この3つの中で最も抽象度が高いからという理由もあるが、上のようなデザインパターンとして使用できるからというのが大きいのだと思う。

最後に、今回実装したモナドがモナド則を満たしていることを確認する。

適当に値を与えて確認する。１つ１つやってもいいけど、せっかくなので最近覚えたlist applicativeを使って複数個のパターンを一気に確認する。

nums  = [-5, 0, 5]
ck    = [YES, NO, Unknown (-5), Unknown 0, Unknown 5]
funcs = [passTest (divisable 400),
        failTest (divisable 100),
        determine (divisable 4)
       ]

assertLeftIdentity  = (\f x -> (return x >>= f) == f x) <$> funcs <*> nums
assertRightIdentity = (\m -> (m >>= return) == m) <$> ck
assertAssociativity = (\f g m -> ((m >>= f) >>= g)
                         == (m >>= (\x -> f x >>= g))) <$> funcs <*> funcs <*> ck

イケてた。モナドイケイケ～～。

Haskellで有名アルゴリズム実装（4）

第4弾。今日は、素数列挙のアルゴリズム。

ある数が素数かどうかの判定は、擬多項式時間でできます。ある区間[1, n]内の素数を列挙するには、普通にやるとO(n*sqrt(n))かかります。

エラトステネスの篩というアルゴリズムを使用すると、O(n)で計算できます。

では、実装。

Haskellのエラトステネスは美しいです。

eratos []     = []
eratos (x:xs) = x : eratos [y | y <- xs, mod y x /= 0]

これで十分な場合がほとんどですが、100万までの素数列挙とかすると、上の実装ではたえられません。エラトステネスが速いのは、ある素数pの倍数を2*p, 3*p,..とpおきに消していけるからで、上の方法だとすべての要素をなめているので遅いです。加えてリストは参照にO(n)かかるので遅いです。

ある素数に対して必要な要素だけ参照する+高速に参照することができるデータ構造に変換する。というのが必要です。100万くらいの素数列挙なら、Setを使えば十分でしょう。あとは、thunkがメモリを圧迫するのを防ぐ+遅延評価によるスタックオーバーを防ぐため正格評価を用います。

primes :: [Integer] -> [Integer]
primes = toList . (seive 2) . fromList

seive :: Integer -> Set Integer -> Set Integer
seive p xs
     | p*p > findMax xs = xs
     | otherwise        = seive (p+1) $ foldl' (\acc x -> delete x acc) xs [2*p,3*p..q]
                          where q = findMax xs

ただ、上の実装でも本来のO(n)の速度には届かないので、より大きな範囲（[1, 10000000]など）の素数列挙をしたい場合は、工夫が必要です。Setは参照がO(log n)なので、ランダムアクセスが可能なデータ構造を利用する + なるべく遅延評価の機構をいれないようにするで達成できるかと思います。

Haskellのページに戻る。

Haskellでメモ化する（中級）

　前回のpostで、Haskellでメモ化する方法について書きました。そのとき、疑問点として、「なんで、メモ化と不動点が繋がってるねん！なんでやねん！」というのがありました。その疑問が解決したので中級編を書きます。

　どうやら不動点コンビネータを利用してやりたいことは、「任意の関数が与えられたときに、その計算をメモ化する関数を作成する」ということのようです。この発想にはびっくりです。どれだけ抽象度の高いことをやるんだよ。という感じです。imparative、あんど、object-orientedな思想に育てられた（毒された？）普通のプログラマーである私には思いもよらない発想でした。。

　まずは、メモ化と、不動点コンビネータについて少しおさらいをします。フィボナッチ数列をメモ化再帰で解くプログラムを書いてみます。

memoFib =
   let fib 0 = 1
       fib 1 = 1
       fib n = memoFib (n-1) + memoFib (n-2)
   in (map fib [0..] !!)

これは、list-comprehensionを用いて以下のように書くこともできます。（こっちの方が分かりやすいかも？？）

memoFib = [fib x | x <- [0..]]
         where fib 0 = 1
               fib 1 = 1
               fib n = memoFib !! (n-1) + memoFib !! (n-2)

と、とりあえず、関数がgivenであれば、その関数をメモ化してあげることは出来ます。

　次に、不動点コンビネータについて。Haskellでは、Control.Monad.Fixに不動点コンビネータが定義されています。これを使ってもいいですが、勉強のために自分で実装してみます。再帰関数を不動点コンビネータに書きかえる方法についてはここのページがとても分かりやすいです。

フィボナッチ関数を不動点コンビネータで書きなおします。フィボナッチ数列は、ラムダ式を用いて

fib = (\x -> if x <= 1 then 1 else fib (x-1) + fib (x-2))

と書けます。さらに、再帰で呼び出される関数を引数で与えることにすると、

fib = (\f x -> if x <= 1 then 1 else f (x-1) + f (x-2)) fib

と書けます。上の式をみると、フィボナッチ関数は

(\f x -> if x <= 1 then 1 else f (x-1) + f (x-2))

の不動点となっています。つまりフィボナッチ関数は、

fix (\f x -> if x <= 1 then 1 else f (x-1) + f (x-2))

と書けるというわけです。（詳しくは上のお勧めページ参照。）

　と、準備はここまで。いよいよ任意の関数をメモ化関数に書きかえる関数を定義するわけですが、任意の関数というのは、Haskellビギナーの私にはちと難しいので、(Int -> Int)型の任意の関数をメモ化することにします。

　まず、引数にそのままメモ化したい再帰関数を渡してしまうと意味がない（メモ化機構を入れられない）ので、メモ化したい関数を不動点にとるような関数を引数として渡してあげます。

fix f = f (fix f)

とすると、fを渡すという意味です。

ここまで来ると、なんとなくできそうな感じがします。なんとなく実装したのが、これ。

memorizeCalc f = ([f (memorizeCalc f) x | x <- [0..]] !!)

計算自体は正しいですが、メモ化が機能していないようです。前にも書きましたが、Haskellの場合、グラフ簡約のメモ化（ここでいうメモ化は、グラフ簡約されて同じ計算がされないという意味のメモ化）が働く場合と働かない場合があります。この形は前に見たことがあります。おそらく、引数をやめて実体を無限リストとして定義してあげればグラフ簡約が働くはず。で、書いたのがこちら。

memorizeCalc f = memo
       where memo = ([f memo x | x <- [0..]] !!)

これで、「与えられた関数をメモ化関数に変換する関数」が出来ました。実際に確かめてみます。

fib' = (\f x -> if x <= 1 then 1 else f (x-1) + f (x-2))
trippleFib' = (\f x -> if x <= 2 then 1 else f (x-1) + f (x-2) + f (x-3))

main = do print $ memorizeCalc fib' 300
         print $ memorizeCalc trippleFib' 300

値も正しく、高速化もされていることが確認できました。

　次のテーマは、

• (Int -> Int)型以外の関数のメモ化
• 高速化

です。高速化については、無限平衡二分木を用いた方法が紹介されています。また、Stateモナドを用いた方法があるようです。他にもData.ArrayのArrayが使えるんじゃないかなあ、とにらんでいます。。この辺が出来たら、上級編を書きたいです。