Haskell模写修行（2）

問題: Codeforces Beta Round #84 (Div. 2 Only) B. Lucky String

自分が書いたコード:

import Control.Monad

main = do n <- liftM read getLine
          putStrLn . take n . cycle $ "abcd"

お手本コード（bjin氏参考）:

main = interact$(`take`cycle"abcd").read

interact使って書くとごちゃごちゃしそうだったので、doブロックで書いた。Control.MonadのliftMは便利。適用する関数の引数の数に応じて、liftM2, liftM3, liftM4, liftM5 などもある。

bjinさんはゴルフ好きな感じである。takeをinfix operationにしてセクション（二項演算子の部分適用をするための手法）を実現してるのは面白いテクニックだなーと。

Haskell模写修行（1）

"New Zealand Autumn" by Abaconda

　

　モナドを実装したり、メモ化したり、コンビネータを書いたりとかしたけど、Haskellらしいコードが全然書けない。一応は動くものは書けるけど、どうも命令型で関数型っぽくない。ということで、Haskellらしいコードが書けるように修行する。

　"らしい"コードを書けるようになるためには、素晴らしいプログラマが書いた"らしい"コードをたくさん読んで吸収していくのが一番じゃないかなと思う。ということで、しばらく以下をやっていこうと思います。

1. codeforceの問題をHaskellで解く
2. ratingの高いHaskellerが書いたコードと自分のコードを比べる
3. いいところを吸収する

問題: Codeforces Beta Round #89 (Div. 2) A. String Task

自分が書いたコード:

import Char

isVowel 'a' = True
isVowel 'i' = True
isVowel 'u' = True
isVowel 'e' = True
isVowel 'o' = True
isVowel 'y' = True
isVowel 'A' = True
isVowel 'I' = True
isVowel 'U' = True
isVowel 'E' = True
isVowel 'O' = True
isVowel 'Y' = True
isVowel x = False

addDot   [] = []
addDot (x:xs) = '.' : x : addDot xs

main = do 
         input <- getLine
         putStrLn . map toLower . addDot . filter (not . isVowel ) $ input
         

お手本コード（bjin氏参考）:

import Char
g x | x < '0' || elem x "aoyeui" = ""
    | 1 > 0           = '.':[x]
main = interact $ concatMap (g.toLower)

ポイント:

1. elem
2. concatMap
3. interact

　”Haskellらしい”うんぬんより、まず便利な関数を覚えていくところからかな。

UNIX環境でサービスのポート番号を調べる方法

教えてもらったので忘れないうちにメモ。

例えば、MySQLのポートを調べたいときは、以下のようにします。

$ sudo lsof -P | grep mysql | grep TCP

　lsofはプロセスによって開かれているファイルの情報を出力するらしいです。ここでいう"ファイル"とは、通常のファイルやディレクトリ、さらにはライブラリ、ストリーム、ネットワークファイルまで含まれます。

-Pオプションは、ポート番号をポートの名前（サービス名）に変換するのを抑制するという意味です。-Pオプションを付けない場合は

TCP *:mysql (LISTEN)

と表示されるのが、-Pをつけると、

TCP *:1234 (LISTEN)

のようにポートの番号で表示されます。

ちょっと調べてみたところ、-iオプションというものもあって、これを使うとネットワーク関連の情報だけを取り出せます。上のコマンドは以下のようにしてもいいです。

$ sudo lsof -i tcp | grep mysql

　lsofの他にも、socklistというコマンドがあるみたい。このコマンドはオプションなしで

$ sudo socklist

とすると、lsofと同様に通信プロトコル、ポート番号、プロセスID、サービス名などが表示されます。こっちの方が楽ですねー。

MCMCで二次元正規分布に従う乱数を生成する

"The Borneo office" by angela7dreams

　マルコフ連鎖モンテカルロシミュレーション（MCMC）を用いて、事象(x, y)の発生確率が二次元正規分布

に従う事象をサンプリングしてみます。この分布はグラフに書くと下のようになります。

さて、この確率分布に従う乱数を発生させたいと思います。今回は一番簡単なメトロポリス法を用いて推移確率を決定します。

Haskellで書きました。自分の実力足らずでいろいろとやばいです。IOモナドのところが何かぎこちないです。実行速度も遅いです。Haskellもっと勉強しないとっ。。

import System.Random
import Text.Printf

select :: (Double, Double) -> IO (Double, Double)
select (x, y) = do
                    dx <- getStdRandom $ randomR (-5,5)
                    dy <- getStdRandom $ randomR (-5,5)
                    return (x+dx, y+dy)

prob :: (Double, Double) -> Double
prob (x, y) = 1/(2*pi) * exp(-(x*x+y*y)/2)

nextStatus :: (Double, Double) -> IO (Double, Double)
nextStatus (x, y) = do
                 (xt, yt) <- select(x, y)
                 let p = min 1  (prob(xt, yt) / prob(x, y))
                 q <- getStdRandom $ randomR (0, 1)
                 if p >= q then return (xt, yt) else return (x, y)

markovChain :: (Double, Double) -> [IO (Double, Double)]
markovChain (x, y) = iterate (>>= nextStatus) (return (x, y))

printPoints :: (Double, Double) -> IO()
printPoints (x, y) = printf "%f %f\n" x y

main = mapM_ (>>= printPoints) . take 3000 . drop 1000 . markovChain $ (0, 0)

それっぽい形になっています。E(x) = -0.00368、E(y) = -0.0121、V(x) = 1.03、V(y) = 1.02でした。試行回数を増やせば、それぞれ0, 0, 1, 1に近づいていくでしょう。

MCMC入門まとめ

マルコフ連鎖モンテカルロシミュレーション（MCMC）の入門資料をいくつか読んだので、分かったことと考えたことをメモしておく。

［分かったこと］

■アルゴリズムの流れ

1. 不定分布が分かっているものから、マルコフ性を仮定して遷移確率を推定する。

2. 適当な初期条件を選んで、1.の推移確率を用いて推移させる。

3. 十分な回数ループしたら終了。そうでなければ、2.へ戻る。

　1. では、ある特定の不定分布におちつく遷移確率行列は複数あることに注意。推移確率の選定には｢詳細釣り合い｣の条件（状態遷移図を書いたときに、確率のフローのようなものが釣り合っていると考える）が一般的に使われる。詳細釣り合いだけでは、自由度が残り値がfixしないため以下いずれかの方法で値を決める。

• メトロポリス法
• ギブスサンプラー
• メトロポリスヘイスティング（M-Hアルゴリズム）

　2. では、いろいろな初期解に対するある時刻における各事象の発生確率（パターン数を無限にする）と、ある初期解について無限回推移させた場合の各事象の発生確率が等しくなる（エルゴード性）という特性が活かされている。

■用途

1. 任意の分布に従うデータをサンプルできる
2. モンテカルロ積分

1. は、最初よく分からなかったが、良く考えると任意の確率分布にしたがう乱数をつくるのは結構むずかしい（[参考文献3] 参照）。普段よく利用するのは一様分布にしたがう乱数だが、正規分布N(0, 1)にしたがう乱数とか、もっとわけの分からない適当なp(x)にしたがう乱数発生器を作れて言われると難しい。MCMCを用いると、簡単に任意の乱数発生器を作成できる。

2.は、1.の応用のように思える。例えば、事象xが起こる確率をp(x)としてf(x)の期待値を求めよ。という問題を考える。∫ f(x) p(x) d(x)を計算することになるが、xの次元が大きいと計算が困難なためMCMCを利用したモンテカルロ積分がしばしば利用される。

［考えたこと］

　MCMCはマルコフ性のある事象に対するシミュレーション手法だと勘違いしていたが、どちらかというと任意の確率分布にしたがうサンプリングデータを生成できるということに利用価値がありそうだ。

　MCMCでマルコフ性のある事象をシミュレーションする場合は、確率遷移行列の選び方がとても重要だと思う。パラメータの取り方に自由度があるメトロポリスヘイスティングが一番使えそうだが、提案密度の選び方にはどうしても主観が入ってしまいそう。MCMCはベイズ識別にも応用されているらしいので、そのあたりを詳しく調べる必要がありそうだ。（というかそこがそもそも知りたかったところなので。）

[参考文献]

1. マルコフ連鎖モンテカルロ法入門-1
2. マルコフ連鎖モンテカルロ法入門-2
3. Markov chain Monte Carlo

確率遷移行列の形と定常状態

マルコフ性（単純マルコフ過程）をもつ事象がある確率推移のもと、どのような定常状態に落ち着くのか実験してみた。

　Haskellで書きました。こういうのはHaskellが綺麗に書けますねー。

genNextPrb a b c d = \(x, y) -> (a*x+b*y, c*x+d*y)
main = mapM_ print . take 200 . iterate nextPrb $ (1, 0)
       where nextPrb = genNextPrb 0.2 0.2 0.8 0.8

分かったこと。確率遷移行列をP, 事象の分布の初期値をx0, 収束時の事象の分布をxnとおきます。

1. Pが単位行列の場合は、xn = x0となる。
2. すべての行iにおいてPi0 = Pi1 = Pi2 =...の場合（列の値が同じ）、xn = [Pi0, Pi1, Pi2, ...]となる。
3. Pが対象行列の場合、xnはすべての要素の値は等しい。（一様分布になる）
4. 1.以外の場合は、xnはx0の値によらない。（エルゴード性）

　少しだけ補足。

　1.はベクトルに単位行列かけるだけなので確率分布がかわらないのは当たり前。

　2.は、次に特定の事象がおこる確率は、現在の状況に依存しないことを表している。例えば、サイコロをふる場合など。これはマルコフ性を持つ事象とは呼ばないのかもしれない。

　3.はおもしろい。x ∈ R^2のときは感覚的に分かりやすい。

例えば、天気の推移にマルコフ性があると仮定して

• 晴れの次が晴れの確率を0.9
• 晴れの次が雨の確率を0.1
• 雨の次が晴れの確率0.1
• 雨の次が雨の確率0.9

とする。このとき確率推移行列は対称行列である。一旦晴れになるとなるべくその状態を保とうとするし、一旦雨になるとなるべくその状態を保とうとする。よって仮にx0 = [1, 0]、つまり晴れの確率100%などとしても、次の日には雨が10%になって、その10%は晴れ側にほとんど戻ってこないのでそれが蓄積されていって最終的にxn = [0.5, 0.5]で落ち着く。これは感覚的に分かりやすい。

　では一般にx ∈ R^nの場合はどうか？これは対称行列を写像として見たときの意味を考えるか、対角化して極限を取るみたいなことをして数学的に証明できるのかもしれないが、大変そうなので直感的に理解する方法を考えてみる。オートマトンみたいな状態遷移図を書けば分かりやすそう。。3次元の場合も、以下の図をみると直感的には定常状態では均等に落ち着きそうです。

たくさん持っているところからはたくさん出て行き、少ししかないところからは少ししかでていかない。同じ色の線は同じ遷移確率を持るのがポイントです。（電子回路のオペアンプをバッファとして使う場合に定常状態がどうなるかを考えるのに似ているかも。）

4. は、確率遷移行列が単位行列以外の場合は、定常状態が初期状態によらないことを意味している。初期状態が何でもいいというのは、シミュレーションをする場合にとてもありがたい特性である。また、初期状態によらないということは、マルコフ性を持つ確率モデルは遷移行列によってのみ特徴づけられることを意味している。固有値とか固有ベクトルとかがここでも使えそうである。

sedの最短一致

最近知ったけど、sedはstream editorの略らしい。

それはどうでもよくて、sedの正規表現のマッチングの小技的なものをメモとして書いておく。

HTMLをパースする場合を考える。

　

　例えば、

<tr><td>id</td><td>data</td></tr><tr><td>1</td><td>AKB48</td></tr><tr><td>2</td><td>perfume</td></tr><tr><td>3</td><td>bump.y</td></tr>
…

のようなデータからHTMLのタグだけを消して、意味のある値だけを取り出したいとする。つまり、

id data

1 AKB48

2 perfume

3 bump.y

…

のような出力をえたい。

　このときに、

sed s/\<.*\>//g test.html

みたいな書き方をすると悲しいことになる。（行全体にマッチしてしまい、すべて空行になってしまう。）

これは、複数のマッチングが考えられる場合なるべく長いものにマッチする（最長一致）からである。この性質はsedに限ったことではなくて正規表現全般に言えることで、デフォルトでは最長一致である。

　最短一致を利用するためのオプションを持つ言語もあるが、sedには無いので何らかの工夫が必要となる。今回の場合は、閉括弧が現れた時点で一旦マッチングをやめるようにすればいいので、以下のように書くと思惑のデータがえられる。

sed s/\<[^\>]*\>//g test.html

　データの間にスペースを入れたい場合は、

sed 's/<[^>]*>/ /g' test.dat

とする。（スペースを用いるので''で第一引数を囲む。''を使うとエスケープ文字の扱いが変わるので<>はエスケープしない。※注）

※注)

　シングルクオーテーションで囲んだ場合、シングルクオーテーション以外の文字は普通の文字と同じ扱い。ダブルクオーテーションで囲んだ場合は、"$`\以外の文字は普通の文字と同じ扱いになるそうです。

http://genki01.cc.hokudai.ac.jp/reo/pukiwiki/index.php?shellscripts

エピメニデスのパラドックス

昔、クレタ人のエピメニデスはこう言った。

「すべてのクレタ人は嘘つきである。」

これはパラドックス（矛盾）である。

　もし、彼の言ったことが本当だとする。本当のことを言った彼は嘘つきではない。しかし、彼はクレタ人は嘘つきだと言っているので、彼自身は嘘つきということになる。これは矛盾。

　もし、彼の言ったことが嘘だったとする。彼は嘘を言ったので嘘つきである。彼は嘘をいっているので真実は「すべてのクレタ人は正直者である」となる。ということは彼は正直者であるはずである。これは矛盾。

　しかし、上の論法には引っかかるところがある。「すべてのクレタ人は嘘つきである。」の反対は、「すべてのクレタ人が嘘つきだとは限らない。」である。これをふまえて彼が嘘をついた場合をもう一度考えてみる。

　もし、彼の言ったことが嘘だったとする。彼は嘘を言ったので嘘つき。彼は嘘を言っているので真実は「すべてのクレタ人が嘘つきだとは限らない。」となる。これは彼が嘘つきだという仮定に矛盾しない。ということは、上記の命題は「偽」ということで実は矛盾ではない。

　ちょっと元の命題を変えてみる。エピメニデスは言った。

「私は嘘つきだ。」

この場合は、矛盾である。エピメニデスが正直者なら「私は正直者だ。」というはずだし、嘘つきなら嘘を付いて「私は正直者だ。」というはずである。さてさて、ここでまた面白い問題が浮かんでくる。

嘘つき（嘘しかつかない人）と正直者（嘘は絶対つかない人）しかいない村があります。通りがかりの人に質問をして、その人が嘘つきか正直者かを確かめたいです。あなたは何と質問すればいいでしょう？

（「あなたは正直ものですか？」と聞いてもどちらも「はい。」と答えます。）

答えは、

・・・・

（自分で考えたい人は、ここから先は後で読んでください。）

・・・・

答えは、2とおりあります。

1. 『お前は人間か？』と質問をする。
2. 『「あなたは正直者か？」という問いに「はい」と答えますか？』と質問する。

Visualize Apache Log on Google Map

I'm lately playing by mapping Apache log on Google Map. This is it, named Access Visualizer. (I recommend you see it with Google Chrome. Some browsers cannot handle some features in the app.)

If you put your mouse cursor on a marker, you can see the next information:

• ISP
• Organization
• Country
• Region
• Coordinate (latitude, longitude)

I used a witty service that can convert IP addresses into location information. And this service allows you to use their functionality by API (XML-based REST).

What I did to create the app is:

1. consider the DB structures and create DB tables.
2. write a shell script and some PHP scripts to parse Apache log, run queries, and store the data into MySQL database.
3. write a PHP script that automatically generate Javascript code which utilize Google Map API.

Not so difficult tasks. Let's go deeper into each section.

1. is an interesting part. Since the query number is limited to 100 times a day, I think it's good to store IP-location data into my own database so that I can re-use the same queries previously used. In the current IP address scheme, there's 2^32 IP addresses, which is about 4*10^9 in decimal. But an organization has a several, or hundreds of global IP addresses. So considering an efficient database table structure is interesting!

2. is kind of tedious. Just wrote a shell script and PHP scripts.

3. is the first time for me to use Google Map API. I actually tried to use the API before, but I gave up it since the use was restricted and some registration was needed at the time. But now its version is updated, and now, you can utilize their APIs without any registration and restriction!! Yikes, thanks Google!!

The future trials:

With lots of markers, there are some parts you cannot see. Japan is hidden and unseen due to the markers now lol So I'm considering to create some layers that control whether show a marker or not. Some markers are to be unseen when zoom-out, but be seen when zoom-in.

Plus, public location information about IP is not so accurate unfortunately . For example, KDDI, University of Tokyo and other organizations are located in the exactly same point according to utrace. Due to this, you cannot see the marker for University of Tokyo and KDDI. I have to scatter the points at the exactly same point in some manner.

オイラーの定理についてちょっとだけ

Google Code Jam Japan決勝の問題B「バクテリア」が非常におもしろかったので、オイラーの定理について以下の本で復習したことをまとめておきます。（数論の入門書としてかなりお勧めです。オイラーの定理、フェルマーの小定理、拡張ユークリッド互除法、中国の剰余定理などプロコンに頻出の定理はほぼ網羅しています。）

※一部自分の考えも入っているため、間違えている箇所がありましたらご指摘ねがいます。

　a^N を nで割った余りについて考えます。n = p1^k1 * p2^k2* .... pm^kmと素因数分解できるとします。このとき、a = A * p1^a1 * p2^a2 * .... pm^amとおきます。(Aはnと互いに素な数)。考えられるパターンは、以下の3つです。

1. a1 = a2 = .... am = 0のとき（aとnが互いに素なとき）
2. a1 >0, a2 > 0, ... am >0のとき(aがnの素因数をすべて含むとき)
3. その他のとき(aがいくつかのnの素因数を含むとき)

　1. はオイラーの定理より、a^φ(n) = 1なので、φ(n)はa^Nをnで割った余りの列の周期になります。

　2. のとき。a1 * N >= k1, a2 * N >= k2, ..., am * N >= kmとなるようなNに対してa^Nをnで割ったあまりは0となります。十分に大きなNに対しては、φ(n)はa^Nをnで割った余りの列の周期になります。"十分に大きな"とはどれくらいか？式を見た感じだと、かなり悪く見積もっても最悪でN = nとすればよさそうです。

　3. のとき。これは長くなるので書けませんが、

φ(n) = φ(p1^k1) * φ(p2^k2) * .... φ (pm^km)

となるのを利用して、mod(pi^ki), i = 1,2, ..., mにおいてa^Nを評価して結果を中国の剰余定理で統合するようなことをしています。この場合もφ(n)はa^Nをnで割った余りの列の周期になります。

断言はできませんが、こちらも式を追う限りだと、最悪でもN=nとすれば周期的になりそうです。

　プログラムでテストしてみます。

const int MAX_N = 1000;

int gcd(int a, int b) {
    if (b == 0)
        return a;
    return gcd(b, a%b);
}

int euler(int x) {
    int ret = x;

    for (int i = 2; i * i <= x; i++) {
        if (x % i == 0)
            ret = ret * (i-1) / i;
        while (x % i == 0)
            x /= i;
    }
    if (x > 1)
        ret = ret * (x-1) / x;

    return ret;
}

int modPow(int x, int p, int mod) {
    int ret = 1;
    while (p > 0) {
        if (p & 1)
            ret = ret * x % mod;
        x = x * x % mod;
        p >>= 1;
    }
    return ret;
}

bool cover(int x, int n) {
    for (int i = 2; i * i <= n; i++) {
        if (n % i == 0 && x % i > 0)
            return false;
        while (x % i == 0)
            x /= i;
        while (n % i == 0)
            n /= i;
    }
    if (n > 1 && n != x)
        return false;
    return true;
}

bool assert1(int x, int mod) {
    int v = euler(mod);

    if (modPow(x, v, mod) != 1) {
        fprintf(stderr, "Error1 %d, %d", x, mod);
        return false;
    }
    return true;
}

int assert2(int x, int mod) {
    int y = 1;
    int i = 0;
    for (; y > 0; i++)
        y = y * x % mod;

    if (i > mod) {
        fprintf(stderr, "Error2 %d, %d", x, mod);
            return -1;
    }
    return i;
}

int assert3(int x, int mod) {
    set<int> used;
    int y = x % mod;

    int i = 0;
    for (; !used.count(y);i++) {
        used.insert(y);
        y = y * x % mod;
    }
    if (i > mod) {
        fprintf(stderr, "Error3 %d, %d", x, mod);
            return -1;
    }

    return i;
}

int main() {
    // consider sequences i^k mod j

    int passed1 = 0;
    int passed2 = 0;
    int passed3 = 0;

    int loop2 = 0;
    int loop3 = 0;
    for (int i = 1; i < MAX_N; i++) {
        for (int j = 2; j < MAX_N; j++) {
            int g = gcd(i, j);

            if (g == 1)
                passed1 += assert1(i, j);    // eurler's theorem
            else if (cover(g, j)) {          // j has all prime factor of i
                int loop = assert2(i, j);
                passed2 +=  loop > 0;
                loop2 += loop;
            }
            else {
                int loop = assert3(i, j);
                passed3 += loop > 0;
                loop3 += loop;
            }
        }
    }

    int total = passed1 + passed2 + passed3;
    printf("case1 event prob: %.2lf\n", 1.*passed1/total);
    printf("case2 event prob: %.2lf average loop: %.2lf\n", 1.*passed2/total, 1.*loop2/passed2);
    printf("case3 event prob: %.2lf average loop: %.2lf\n", 1.*passed3/total, 1.*loop3/passed3);

    return 0;
}

以下、実験結果。ランダムに2つの正の整数を取ったとき、半分以上の確率でお互いに素になるようです。iがjの素因数をすべてカバーするケースはまれでMAX_Nが大きくなるに従って減っていくように見えます。また、i^nがmod jにおいて周期的になるのはn < MAX_Nと見積もっておけば十分だといえます。（case 2については、式と結果から考えてlogオーダーっぽい。）

MAX_N = 100

case1 event prob: 0.61

case2 event prob: 0.09 average loop: 2.07

case3 event prob: 0.30 average loop: 7.57

MAX_N = 1000

case1 event prob: 0.61

case2 event prob: 0.02 average loop: 2.70

case3 event prob: 0.37 average loop: 46.79

k-smooth numberの分布を考える

そろそろ100桁くらいの数の素因数分解でもチャレンジしたいなと思ってるけど、general number field sieveは難しすぎてまだ無理。ちょうど手ごろな感じのquadratic sieveというのがあったので、これなら実装できそうである。110桁程度までだと、GNFSよりQSの方が効率がいいらしい。

　さて、このquadratic sieveだが、大抵の数は小さな素因数に分解できるというところが大きなポイントとなる。一般に、最大素因数の数がk以下の正の整数をk-smooth numberと呼びます。では実際に、smooth numberってどれくらいあるのでしょうか？

　以下の問題設定をします。

[1, N]までの正の整数に含まれるk-smooth numberの数がN/2を超えるのはkがいくらのときでしょうか？

このkが小さければ、小さいほどQSを実装する側としては嬉しいわけです。他にも面白い解法がありそうですが、取りあえずパソコンのパワーを信じてがりがりエラトステネスで篩にかけます。普通のエラトステネスの場合は、素数か合成数かをマークしますが、今回は最大の素因数をキャッシュしていきます。

const int MAX_N = 1000;
int fac[MAX_N+1];
map<int, int> smooth_pop;

int main() {
    memset(fac, 0, sizeof(fac));
    fac[1] = 1;

    for (int i = 2; i <= MAX_N; i++) {
        if (fac[i] > 0)
            continue;
        for (int j = i; j <= MAX_N; j += i)
            fac[j] = i;
    }

    /* for debug
    for (int i = 1; i <= MAX_N; i++)
        cout << i << ": " << fac[i] << endl;
    */

    smooth_pop.clear();
    for (int i = 1; i <= MAX_N; i++)
        ++smooth_pop[fac[i]];

    int cnt = 0;
    for (map<int, int>::iterator itr = smooth_pop.begin(); itr != smooth_pop.end(); itr++) {
        cnt += itr->second;
        printf("%d %.3lf\n", itr->first, 1.*cnt/MAX_N);
    }

    return 0;
}

以下、Nの値と、k-smooth numberの数がN/2以上になるkの値です。

N	k
10^3	43
10^6	3257
10^8	53773

　10^8程度の数をquadratic sieveで分解したい場合は、53773-smooth numberまで調べれば、2個に1個は使える数が入っているということです。篩にかける際に、実際に必要となるメモリのサイズは、p(53773)の二乗程度です。p(x)はxが何番目の素数かを表します。p(x)はここでは素数定理を用いてx / log(x)で概算します。p(53773) = 4937です。4937 * 4937程度の行列なら特別なこと（非ゼロの要素だけもつようなデータ構造）をしなくても、普通にメモリに格納できます。この感覚だと100桁程度の数をQSで分解しようとすると、スパースな行列を扱うためのアルゴリズムが必要になりそうです・・・。

以下、k-smooth numberが{i | i = 1,2 ... N}の集合全体の何%を占めるかを示します。横軸がk、縦軸が占める割合です。予想どおりのグラフの形になりました。上から順にN=1,000、 1,000,000、 100,000,000のときです。

apacheのログを見て遊ぶ

最近、apacheのログの見方が分かるようになりました。いろいろ整形して遊んでみると面白いです。デフォルトでは、/var/log/apache2にログが吐かれます。まだ基本的なことしかできないですが、一応まとめ。

1. 直近のアクセスログを見る

tail access.log

2. 自分のローカルマシン以外のアクセスログを見る（自分のマシンのIPがxx.xx.xx.xxの場合）

grep -v xx.xx.xx.xx access.log

3. ユニークユーザー数を見る

cat access.log | awk '{print $1}' | sort | uniq | wc -l

4. アクセス元を訪問数の多い順にソートして表示する

cat access.log | awk '{print $1}' | sort | uniq -c | sort -r

5. リファラー（どこからこのページにやってきたか）を多い順にソートして表示する

cat access.log | awk '{print $11}' | sort | uniq -c | sort -r

4., 5. で$1、$11などと書きましたが、これはapacheの設定によって違うので、自分の環境にあわせて適当に変えてください。

あと、IPから場所を特定できるサイトがあるので、これと組み合わせるとおもしろいです。Google Analyticsだと都市名までしか出てきませんが、このサービスを使うと、大学名や企業名まで分かります。有名な海外大学からアクセスがあったりするとテンションがあがります。このサイトからは、APIが公開されていてHTTPでxmlデータを取り込むことができます。非商用のクエリは100件まで無料ということです。アクセス元のIDをためていって、適当なタイミングでAPIを流すようにして、アクセス元IP - 場所の対応データを溜めていくのもおもしろいかもしれません。

Google Code Jam Japan 2011決勝

113位でGoogle T-shirtsゲットです！！Aが早く解けてしまって、暫く10位台だったので、かなり興奮してましたが、その後何もできませんでした。BとCトライしましたが、WAでした。やったことを書きます。

A. アンテナ修復

　角度が等しいので、三角形の面積は辺の長さによってのみ決まります。長い辺同士を掛けた方がお得なので、greedyに入れて行きます。アンテナが円の形をしているので、双方向キュー（STLのdeque）を使って実装しました。本番は直観のまま突き進んで、今回は運よくACだったけど、証明せずにsubmitっていうのはあまりよくない習慣です。ちょっと証明を考えてみました。

　一番長いエレメントをx1, 二番目に長いエレメントをx2とおく。このとき、私の解法では最適解においてx1 * x2の項が存在しています。これを背理法で証明してみます。

もし、x1*x2の項がないと仮定する。このとき、ある項xp * xqのxqとx2を交換する。すると、x1*xq、xp*x2という項ができる。もとの最適値に対して、この解は、-(x2 - xq) * x1 + (x2- xq) * xp = (-x1 + xp) (x2 - xq)だけ変化したことになります。この変化分は負なので、x1*x2という組み合わせは最適値に含まれます。

以下同様にしてgreedyにとっていけばOKです。

B. バクテリアの増殖

　ただの繰り返し二乗法じゃない？と思ってやったけど最後のサンプルがあわない。何度やってもあわない。Aの採点方式が間違っているというアナウンスがあったため、Bのサンプルおかしいんじゃない？と疑いだす始末。暫くしてレッドコーダーの人達がACし出すので、おそらく何か勘違いしてるんだろうなと思いだす。コンテスト終わって、modのところで大きな勘違いをしているのに気付きました。

　A ^p (mod N)　は、 A^(p mod N) (mod N)とは違いますね。。何やってんだか。。この問題、オイラーの定理使って周期出すだけじゃんと思いましたが、AとCがco-primeじゃない場合もあるので、難しいですねー。。十分に大きな値に対してはφ関数の値が周期になるはずなので、これを使えばいいっぽいです。十分に大きなって微妙な表現だけど、どれくらい大きければいいかは数論の本を見て復習しなければ・・。

C. ワイルドカード

　largeは諦めて、small狙いでいきました。単純な全探索でいいのですが、C++で正規表現関数を自作して、sampleは何とか通せたけど、submitしたらWA食らいました。暫くバグ取りと格闘して、Pythonでやろうと迷いましたが、暫くPython書いてないので頑張ってC++でごりごりやってたら時間切れ。言語によって差が出てしまうような問題はあまり好きじゃないですが、問題に応じて適切な言語を使うことができるのも技術者として重要な能力の一つと考えれば文句は言えません。

　コンテスト後、Pythonで復習しました。

import re

def smaller(x, y):
    if len(x) < len(y):
        return True
    if len(x) > len(y):
        return False
    if x.count("*") < y.count("*"):
        return True
    if x.count("*") > y.count("*"):
        return False
    return x < y

def solve():
    A = raw_input()
    B = raw_input()
    
    ret = "*" * 1024;
    L = len(A)
    for mask in range(1<<L):
        At = ""
        for i in range(L):
            if (mask >> i & 1):
                At += A[i]
            else:
                if len(At) == 0 or At[-1] != '*':
                    At += "*"
            
        if not re.match("^" + At.replace("*", ".*") + "$", B):
            if smaller(At, ret):
                ret = At
           
    return ret        

T = int (raw_input())
for t in range(T):
    print "Case #%s: %s" % (t+1, solve())
    

プロシージャやサブクエリでLIMIT句を使う場合の注意点

　MySQLでストアードプロシージャ作っていて詰まったのでメモ。limitで指定する値をプロシージャの引数で与えようとしましたが、エラーが出ました。

[ダメな例]

delimiter //

DROP PROCEDURE IF EXISTS `LRU_WORDS`//

CREATE PROCEDURE LRU_WORDS(IN p1 INT)

BEGIN

SELECT * FROM MyDictionary ORDER BY last_delivered LIMIT p1;

END

//

delimiter ;

LIMITは定数じゃないとダメらしいです。

http://dev.mysql.com/doc/refman/4.1/en/select.html

　prepared statementを使う場合は、以下のように変数を利用できます。

一旦ユーザー変数に代入しないとダメなようです。なぜ引数から直で行けないのか・・・

[良い例]

delimiter //

DROP PROCEDURE IF EXISTS `LRU_WORDS`//

CREATE PROCEDURE LRU_WORDS(IN num INT)

BEGIN

SET @limit_num = num;

PREPARE SET_STMT FROM 'SELECT * FROM MyDictionary LIMIT ?';

EXECUTE SET_STMT USING @limit_num;

END

//

delimiter ;

これで解決。

　上のような単純な構文だとOKだけど、LIMITを用いたクエリをサブクエリとして使いたいときにはさらにひと癖あるみたいです。LIMIT句を含むクエリはサブクエリとして利用できないため、以下のように仮想テーブルに一旦落とす必要があります。

http://okwave.jp/qa/q4856233.html

What numbers can Ax + By represent?

It's well known that, when two integers A and B are given, N can be represented as Ax + By when N can be dividable by (A, B). So when A and B are relatively prime, that is (A, B) = 1, All the integers can be represented as Ax + By for some integer (x, y).

But what if when a restriction is given that x and y are positive integers?

For instance, let's consider (A, B) = (4, 3). Since two numbers are co-prime, you can represent 1 as 4 * (-2) + 3 * 3 = 1. But here y is a negative integer. Other than (x, y) = (-2, 3), many other candidate solutions will hit your mind-- (-11,15), (-8,11), (-5,7), (1,-1), (4,-5), (7,-9), and so forth--. But all of them have a negative integer.

In general, the following holds:

When A and B are relatively prime, A x + B y can represent all the numbers above and equal to 2AB for some positive integers (x, y).

The following graphs are self-explainary.

The first one is plotting 4x+3y = 0. The second 4x+3y = 1, the third 4x+3y=11, a case where a feasible solution exists. And the last one shows why Ax + By can represent all the numbers over 2AB.

サーバー奮闘記（15） php CLI版インストール

phpをコマンドラインから使おうと思ったけど、使えなかった。

CLI（Command Line Interface）版のphpがインストールされていなかったことが原因だった。

ということで、CLI版のphpをインストールした。

$ sudo apt-get install php5-cli

これで行けるはず。以下のコマンドでバージョンが表示されればOKです。

$ which php

CLI版のphpは何かと重宝しそうです。perlでもいいんですけど、個人的にperlは象形文字みたいで気持ち悪くてあまり好きじゃないんです。あと今回やりたかったテーマは、MySQLからデータを引っ張ってきてバッチ処理をするだったので、すでにCGIではMySQLの設定が出来てるphpでやりたかったのです。

とりあえず、データはCGIでやるような感じで引っ張ってこれました。

#! /usr/bin/php -q

<?php

// connect to MySQL Server
$link = mysqli_connect("host", "usr", "passwd", "db");
if (mysqli_connect_errno()) {
  printf("Connect failed: %s\n", mysqli_connect_error());
  exit();
}

// send a query
$query = 
"SELECT 
   dic.word, 
   sample.sentence 
 FROM 
   SampleSentences sample right outer join MyDictionary dic 
 ON 
   sample.word_id = dic.id";

$result = mysqli_query($link, $query);
if (!$result) {
  print("An error occurred while processing SQL queries.<br>\n");
  exit();
}

// show word lists
$word = "";
while ($row = mysqli_fetch_array($result))
  print $row['word'] . " " . $row['sentence'] . "\n";
?>

CLI版のphpの説明はここが詳しそうですかね。

http://www.spencernetwork.org/memo/tips-6.php

群、環、体について

数体ふるい法のアルゴリズムを読んでいると分からない言葉が出てきてすぐ挫折します。。1行読むのにも一苦労。とりあえず、数学力が足りません。ということでちょっとずつ勉強していこうと思う。

今日は、群、環、体について。英語では、group、ring、fieldとよばれます。1つずつ概念を抑えていこうと思います。ここでは具体例を考えながら、その性質を考えていきます。以下では、ここに記載されている定義にしたがって議論を進めます。

まず、全体的な話。これらの概念は、集合+その集合に対する演算をひとくくりにした概念です。集合だけではなく、それにたいする演算も統一的に扱って議論を進めたい場合にこのgroupとかring、filedといった概念が必要になります。

group

　群では、集合S上の要素x, yについて二項演算opを施した結果もまた集合Sに含まれます。たとえば、整数の集合Zを考えます。Zの要素に対して二項演算+を実行します。すると結果は整数になります。

"整数 + 整数 = 整数"なんて当たり前じゃないか、と思うかもしれないですが、これはよく考えるととても面白いことです。

群の定義を見ると、いくつかの例は簡単に思いつきます。例えば、整数の集合と演算子+は群を成します。しかし、正の整数と+は群を成しません。単位元が存在しないからです。正の整数 + {0}の集合にすると、単位元は存在します。しかし逆元が存在しないため、群にはなりません。

　次に、整数の集合と*を考えてみます。単位元は1です。しかし、逆元が存在しません。一般に整数nの逆元は1/nとなりますが、これは必ずしも整数にはなりません。それでは、整数じゃなくて実数集合を考えてみます。これは、群になるような感じがしますが、0に対する逆元は存在しません。よってこれも群ではありません。実数集合から0を除いたものと演算*は群をなします。

　上記では、無限集合を考えましたが、有限集合の場合を考えてみます。{0, 1}という集合を考えます。これに対して演算ORを考えます。これは群をなすでしょうか？

　まず、演算ORは集合に対してclosedです。単位元は0です。演算はassociativeです。但し、1に対する逆元が存在しません。よって、{0, 1}と演算ORは群をなしません。では演算ANDに対してはどうでしょうか？演算はclosed、単位元は1、演算はassociativeですが、0に対する逆元がありません。有限集合で群をなすものはそう簡単には見つけられないのでしょうか。。

　法3の世界を考えます。剰余類の完全代表系は、{0, 1, 2}です。この法の世界において演算+を考えます。演算はclosed、単位元は1、演算はassociativeです。すべての要素について逆元は存在するでしょうか？(0, 1, 2)に対して、(0, 2, 1)が逆元になります。よって、この系は群をなします。

ring

　環は集合Sと二つの演算add、mulを統一的に扱うための概念です。群では考える演算子は1つでしたが環では2つです。それぞれ、足算と掛け算に似たような性質を持つためここではadd、mulと表記します。定義を見ると、(S, add, mul)が環であるとき、(S, add)は群をなすようです。つまり、(S, add)が群であるということは、(S, add, mul)が環であるための必要条件らしいです。

　まず、対象の集合が無限集合である環を考えてみます。整数の集合と+はアーベル群をなすので、これに*を加えると環になるか考えてみます。アーベル群とは、演算がcommutativeな群のことです。

　ちょっと話を脱線します。普通演算子はcommutativeな気がします。commutativeじゃない群（アーベル群じゃない群）ってどういうのがあるんでしょうか？

　"-"はcommutativeじゃありません。整数集合と"-"は群を成すでしょうか？演算は集合に対してclosed、単位元は存在しません。n - 0 = nですが、0 - n = -nとなってしまうので演算-に対する単位元は存在しません。んー、これは困った。除算についても同様の議論で、群になりません。

　行列の演算はどうでしょうか？行列の掛け算はcommutativeじゃないです。正方行列と演算*は群をなすでしょうか？演算はclosedです。単位元は単位行列です。行列の積はassociativeです。逆元の存在はややこしい。行列式が0だと逆行列が存在しないのでダメでは・・・？ぐぬぬ。。ということで集合の設定を変えます。

「正則行列集合と演算子*は群を成すか？」正則行列同士の積は正則行列になります。ということで、演算は集合に対してclosed。単位元は単位行列。正則行列なので、逆行列が存在し、これが逆元になります。行列の積はassociative。よってこれは群です。しかし、行列の積はcommutativeじゃないので、これはアーベル群ではありません。出ました！これが群だけどアーベル群じゃない例ですね。

　話を元に戻します。環の話です。整数集合と+, *は環をなすか？まず、整数集合と+はアーベル群を成します。*はcommutativeです。*は+に対してdistributiveです。ということで、これは環ですね。こうやってみると、環は演算子addについては厳しい制約を課しますが、mulについては緩いです。

　

　では、有限集合の環を考えてみます。群のところでみたように、法3の世界でその剰余類の完全代表系と加算+は群をなしました。これに乗算*を加えると、これらは環をなすでしょうか？まず、{0, 1, 2}と+はアーベル群をなします。*はcommutativeで、+に対してdistributiveです。よってこれは環です。

field

　最後、体です。定義によると、体は環です。さらに、体では、演算addに対する単位元以外の要素集合とmulが群をなすことを要求します。環でゆるゆるだったmulに対する条件を厳しくしたものが体であるといえます。

　無限集合の体を考えてみます。整数集合と+, *は環でした。では、これは体でしょうか？つまり、0以外の整数と*は群を成すか？答えはNoです。逆元1 / nが整数になるとは限らないからです。よって、整数集合と+, *は環ですが、体ではありません。

　では、実数の集合と+, *を考えてみましょう。実数集合と+はアーベル群です。また、*はcommutative、+に対してdistributiveなので、これは環です。では、これは体でしょうか？+の単位元0以外の実数集合と*は群を成すでしょうか？これは群のところで議論したように、答えはYesです。よって、実数集合と+, *は体です。

　いよいよ終盤戦です。有限集合の体を考えてみましょう。法3の下で、{0, 1, 2}, *, +は環でした。これは体でしょうか？{1, 2}は法3の下で*に対して群をなすでしょうか？まず、演算が集合に対してclosedか否か。{1, 2}から任意の2つの要素をとって*を適用しても(mod 3)で0にはなりません。よってclosedです。単位元は1で、すべての要素に対する逆元も存在します。また、*はassociativeです。よって、これは体であるといえます。

　どうやら、法の世界での体には面白い法則があるようです。今度は法4の下で同様の議論をしてみます。まず、{0, 1, 2, 3}と+について考えます。以下の表に演算+に対する結果をまとめます。

演算は集合に対してclosedで、単位元0を集合内にもっています。また、表から明らかなようにすべての要素が逆元をもっています。+はassociativeかつ、commutativeです。よって、これはアーベル群です。

+	0	1	2	3
0	0	1	2	3
1	1	2	3	0
2	2	3	0	1
3	3	0	1	2

　

　では、{0,1,2,3}, +, *は環でしょうか？*はassociativeかつ+に対してdistributiveです。よってこれは環です。ここまでは法3のときと同じです。

　それでは、{0, 1, 2, 3}, +, *は体でしょうか？つまり、{1,2,3}, *は群でしょうか？法4のもとでの乗算の結果を表に書いてみます。

*	0	1	2	3
0	0	0	0	0
1	0	1	2	3
2	0	2	0	2
3	0	3	2	1

　

　対象となる集合は{1,2,3}なので、0のところは一応書きましたが無視してください。ここで、1の行と3の行は、2の行と違っています。よく見てください。1の行は、(0,1,2,3)となっていて、0-3までの数字がすべて出て来ます。3の行も(0,3,2,1)となっていて同様です。しかし、2の行は、(0,2,0,2)となっていて、{1,3}が出てきません。これは法の4と2が共通因数を持っているからです。これが体とどう関係があるのか見てみます。

　{1, 2, 3}と*が群をなすか？まず、演算がclosedになっていません。2*2 = 0 (mod 4)なので、集合{1,2,3}に対して閉じていません。単位元はあります。1です。associativityもOKです。逆元はどうでしょう？2に対する逆元はありません。これは、2と4が互いに素ではないためです。

　

　なぜaとbが互いに素でない場合、ax = 1 (mod b)となるようなxは存在しないのでしょうか？これは、目盛りがb個ついている時計を考えれば分かります。この時計の針を毎回a目盛りずつ進めると、どんなにまわしてもgcd(a, b)単位でしか針はまわらないのは明らかです。よってgcd(a, b) = 1出ない限り、ax = 1 (mod b)となるようなxは存在しません。

　ということで、{0, 1, 2, 3}と+, *は環ですが、体ではありません。上記の逆元の存在の理論から、一般に、法pにおける剰余類の完全代表系と+, *からなる系は、pが素数なら体になる。pが素数でなければ、環だが、体ではない。となります。おもしろいですねー。

二項関係とグラフ

"transitive disclosure"（推移閉包）について調べていたら、グラフは二項関係の宝庫だなーと思ったのでメモしておく。

例えば、グラフの頂点集合Vについて二項関係R1を頂点間に辺があると定義する。

例えば上のグラフだと、R1 = {(1,4),(1,6),(2,3), (2,4), (3,2), (3,5), (3,6), (4,1), (4,2), (4,5), (5,3), (5,4), (6,1), (6,3)}である。ここで、集合と順序対の表記の違いについて書いておく。集合は、{}で書かれる。順序対は()で書かれる。集合の場合は、要素の順序や重複は関係ない。つまり、{1,2,3}は、{3,2,1}と同じであり、{1,3,5}は、{1,1,3,3,3,5,5}と同じである。これに対しては、順序対は要素の順序、重複ともに意味がある。

さて、二項関係には特徴的な法則がいくつかある。（詳細はこちら）

• Reflexivity
• Symmetricity
• Anti-Symmetricity
• Transitivity

これらの法則は満たされる場合もあるが、満たされない場合もある。例えば、上記のR1の場合。Symmetricityが満たされる。頂点(v, w)間に枝がある場合は、頂点(w, v)間にも枝がある。無効グラフなので当然である。

次に、以下の有向グラフを考える。

頂点vからwに有向枝が存在するという二項関係をR2と定義する。R2 = {(1,2), (2,3),(3,4)}である。このR2はtransitivity（x R y ⋀ y R z →　x R z）を満たしていない。しかし、(1, 3), (2, 4), (1, 4)の3つの順序対をR2に追加すると、これはtransitivityを満たす。

このように、二項関係Rに要素を追加してある特性pが満たされるようになった場合、この拡張された集合R'をRの特性pに対する閉包と呼ぶ。但し、R'を閉包と呼ぶためには、拡張される要素集合R' - Rは最小である必要がある。

上の例だと、R2がtransitivityを満たすようにするには、最低でもV = {(1, 3), (2, 4), (1, 4)}を追加しないといけないので、R3 := R2 ∨ VはR2の推移閉包であると言える。さてこのR3は(v, w)間に有効路が存在する（v -> wに到達可能）であるという関係になっている。

プログラミング上では、R2はグラフの隣接行列で表される。R3はR2からWarshall-Floydで出すことができる。

const int N = 128;   // number of  verticies
bool conn[N][N];
void transitive_disclosure() {
    for (int k = 0; k < N; k++)
        for (int i = 0; i < N; i++)
            for (int j = 0; j < N; j++)
                conn[i][j] |= conn[i][k] && conn[k][j];
}

この推移閉包には、おもしろい性質がある。ある二項関係Rの推移閉包は以下のように表さる。

transitive disclosure(R) = R^1 ∨ R^2 ∨ R^3 …

ここで、R^1 = R、R^{i+1} = R^i ◦ Rです。この◦は、compositionの記号で、Haskellとかでやるfunction compositionと同様のものです。二項関係におけるcompositionはここに定義されているとおり。Warshall-Floydのアルゴリズムもよく考えるとこの性質を使っている。