平方数について

平方数について簡単なまとめ。

1.平方数とは？

整数の二乗で表すことのできる数のこと。0, 1, 4, 9, 16, ...と続く。

2. ２つの平方数の和

２つの平方数の和に関しては、フェルマーのクリスマス定理が有名。この定理では奇素数が２つの平方数の和で表すことができるか否かについて述べられているが、

(x₁² + y₁²)(x₂² + y₂²) = (x₁x₂ - y₁y₂)² + (x₁y₂ + y₁x₂)²

の公式より、一般の自然数nについて以下のことが言える。

n = n₁² * n_2, ただし、n₂は法4の下で3と合同な奇素数を因数にもたない。

と分解できる場合、nは２つの平方数の和で表すことができる。上記のように分解できない場合は２つの平方和で表すことはできない。

3. ３つの平方数の和

３つの平方和ですべての自然数を表すことができるか？答えはできない。法8の下で7と合同な自然数は３つの平方数の和で表すことはできない。なぜなら、平方数は8を法として{0, 1, 4}のいずれかと合同になるが、これらの任意の３つの組み合わせの和は7と合同になりえないからである。

4. ４つの平方数の和

4つの平方数の和を用いれば、2. および3.で表すことができなかった自然数を表すことができる。

7 ≡ 3 (mod 4)および7 ≡ 7 (mod 8)より、7は２つの平方数の和でも３つの平方数の和でも表すことができない。しかし、

7 = 1+1+1+4

と4つの平方数の和で表すことができる。

すべての自然数は、４つ以下の平方数の和で表すことができる。定理の名前は知らなかったが、「ラグランジュの四平方定理」というらしい。

SRM 527 DIV1 450 P8XMatrixRecovery

2部グラフのマッチングまでは、気付いてた。あと一歩だった。
"ある条件を満たす辞書順最小のものを求めよ。"という問題はgreedyで解ける場合が多いような気がする。

int r;
int c;
bool edge[60][60];
bool used[60];
int pr[60];

class P8XMatrixRecovery {
    bool dfs(int s) {
        used[s] = 1;

        REP(v, 2*c) {
            if (!edge[s][v])
                continue;
            int w = pr[v];
            if (w == -1 || (!used[w] && dfs(w))) {
                pr[s] = v;
                pr[v] = s;
                return true;
            }
        }
        return false;
    }

    int match(vector<string> &rows, vector<string> &cols) {
        memset(edge, 0, sizeof(edge));
        REP(i, c) REP(j, c) {
            edge[i][j+c] = edge[j+c][i] = 1;
            REP(k, r) {
                if (cols[j][k] != '?' && rows[k][i] != '?' && rows[k][i] != cols[j][k]) {
                    edge[i][j+c] = edge[j+c][i] = 0;
                    break;
                }
            }
        }

        int ret = 0;
        memset(pr, -1, sizeof(pr));
        REP(i, 2*c) {
            if (pr[i] >= 0)
                continue;
            memset(used, 0, sizeof(used));
            if (dfs(i))
                ++ret;
        }

        return ret;
    }

public:
    vector<string> solve(vector<string> rows, vector<string> columns) {
        r = rows.size();
        c = rows[0].length();

        REP(i, r) REP(j, c) {
            if (rows[i][j] == '?') {
                rows[i][j] = '0';

                if (match(rows, columns) < c)
                    rows[i][j] = '1';
            }
        }

        return rows;
    }
};

メルセンヌ素数と完全数の関係

メルセンヌ素数と完全数には面白い関係がある。

1. メルセンヌ素数とは？

　メルセンヌ素数とは、メルセンヌ数かつ素数であるような数のこと。メルセンヌ数とは、2^k-1という形で表される数のことで、2進数で表したときに、1のみで構成される数のことである。

2. 完全数とは？

完全数とは、その数のproper divisor（自分以外の約数）の和が自分自身と等しいような数のこと。

6 = 1 + 2+ 3

28 = 1 + 2 + 4 + 7 + 14

が最初の2つの完全数である。

3. 面白い関係って？

　メルセンヌ素数と完全数の間には以下の定理が成り立つ。

メルセンヌ数2^k-1が素数であるとき、2^k-1(2^k-1)は完全数である。逆に、すべての偶数の完全数は、2^k-1(2^k-1)という形（2の冪乗とメルセンヌ素数の積）で表現できる。（※注意）

　証明のために関数σを導入します。σ(n)は、nの約数（properな約数ではなく、自分自身も含むことに注意）の和を表します。例えば、

σ(6) = 1 + 2 + 3 + 6 = 12

です。nが完全数の場合は、σ(n) = 2nです。

　

　この関数σについて以下の補題を考えます。

N = ab（a, bは互いに素）と表せるとき、σ(N) = σ(a) * σ(b)である。

　証明)　σ(a) * σ(b) = (a1 + a2 + a3 + ... ) * (b1 + b2 + b3 + ... ) = Σ(aの約数 * bの約数)である。ここでa1,a2,..はaの約数、b1,b2,..はbの約数である。よって、（aの約数 * bの約数）がNの約数となることと、Nの任意の約数dがある一組の(i, j)についてd = ai * bjとなることを示せばよい。

　まず前半だが、a, bが互いに素であるため、a x + b y = 1となるような整数(x, y)が存在する。両辺にNをかけると、a x N + b y N = Nとなる。aの約数 | a かつ bの約数 | Nなので、（aの約数 * bの約数）はa x Nを割りきる。同様に（aの約数 * bの約数）はb y Nを割りきる。よって（aの約数 * bの約数）はNを割りきる。

　次に後半を示す。a, bの素因数分解を

a = p₁^k1p₂^k2p₃^k3...p_n^kn

b = q₁^t1q₂^t2q₃^t3...q_m^tm

と表記する。ki > 0, i = 1,2,...n、およびtj > 0, j = 1,2,..mである。また、aとbが互いに素なので、p_i≒q_j, i=1,2,..n, j =1,2,...mである。

N= a * bなので、Nの素因数分解は

N = p₁^k1p₂^k2p₃^k3...p_n^kn q₁^t1q₂^t2q₃^t3...q_m^tm

となる。ここで、Nの任意の約数は、

x = p₁^k'1p₂^k'2p₃^k'3...p_n^k'n q₁^t'1q₂^t'2q₃^t'3...q_m^t'm

と表すことができる。ただし、0 ≦k'i ≦ki, i = 1,2,...n、0 ≦ t'j ≦tj, j=1,2,...,m.である。

Nの任意の約数を選ぶことと、(k'1,k'2, ..., k'n, t'1,t'2,..., t'm)の値を決めることは一対一に対応している。(k'1,k'2, ..., k'n, t'1,t'2,..., t'm)が決まれば、(k'1,k'2, ..., k'n)、(t'1,t'2,..., t'm)が決まり、これは、対応するaの約数、bの約数が一意に決まることを意味する。□

　それでは、メルセンヌ素数と完全数に関する定理の証明に入ります。まず、最初の部分を証明します。

2^k-1は素数なので、σ(2^k-1) = 1 + 2^k-1 = 2^kです。2^k-1は、2のべき乗なので約数は1, 2, 2², .... 2^k-2, 2^k-1である。よって、σ(2^k-1) = 1 + 2 + 2² + ... + 2^k-1 = 2^k-1。2^k-1（2の冪乗）と2^k-1（奇数）は互いに素なので補題より、σ(2^k-1(2^k-1)) = σ(2^k-1) * σ(2^k-1) = (2^k-1) * 2^k = 2 * 2^k-1(2^k-1)となる。

σ(n) = 2nなので、2^k-1(2^k-1)は完全数である。

　次に後半を示す。偶数の完全数を2^k-1mと表記する。ただし、k≧2、mは奇数。

補題よりσ(2^k-1m) = σ(2^k-1) * σ(m) = (2^k-1) σ(m)。

また、2^k-1mは完全数であることから、σ(2^k-1m) = 2 * 2^k-1m = 2^km。

よって、

2^km = (2^k-1) σ(m)。

ここで、(2^k-1) | 2^km、かつ2^kと(2^k-1)は互いに素であることから、(2^k-1) | m。

よって、m = (2^k-1) Mとして先ほどの式に代入すると

2^k (2^k-1) M = (2^k-1) σ(m)となり、両辺を (2^k-1) で割ると、

2^kM = σ(m)。

m + M = (2^k-1) M + M = 2^kMより、σ(m) = m + Mとなるが、これはmが素数で、M=1であることを意味する。よって、m = 2^k-1（素数）となり、もともとの完全数は、2^k-1(2^k-1)と表される。□

（※注意）奇数の完全数が存在するか否かは不明です。現時点で発見されている完全数はすべて偶数です。

参考：

[1] http://primes.utm.edu/notes/proofs/EvenPerfect.html
[2] http://ja.wikipedia.org/wiki/%E3%83%A1%E3%83%AB%E3%82%BB%E3%83%B3%E3%83%8C%E6%95%B0

n進数表記の数を(n-1)で割れるかどうか

面白い問題があったので、紹介。

http://poj.org/problem?id=1152

　n進数表記の数値を(n-1)で割り切れるかどうかは、すべての桁の数値を足した数が(n-1)で割れるかどうかを調べればいいです。例えば、10進数表記の111111111、45、1233、123453は桁の合計値が9で割れるので9の倍数です。

　n進数表記のm桁の数 A = A_(m-1).. A_2 A_1 A_0について考えます。つまり、

A = A_(m-1) * n^(m-1) + ... A_2 * n^(2) + A_1 * n^(1) + A_0

です。 n = 1 (mod (n-1))なので、

A = A_(m-1) + ... + A_2 + A_1 + A_0 (mod (n-1))

です。

ということで、Aの各桁の合計値を出してそれが(n-1)で割れるかどうかをみればOKです。

　C++のスパゲッティなソースを張っておきます。（dividable()の中で毎回合計値計算してるのは無駄ですね。。最初まじめにHormerの公式を使って余りを計算して解こうとしていたのの名残です。）

int num[128];

bool dividable(int k, const string &str) {
    int len = str.length();
    int digsum = 0;
    REP(i, len) digsum += num[(int)str[i]];

    return digsum % (k-1) == 0;
}

int main() {
    string str;

    for (int i = '0'; i <= '9'; i++) num[i] = i - '0';
    for (int i = 'A'; i <= 'Z'; i++) num[i] = i - 'A' + 10;
    for (int i = 'a'; i <= 'z'; i++) num[i] = i - 'a' + 36;

    while (cin >> str) {
        int len = str.length();
        int base = 0;
        REP(i, len) base = max(base, num[(int)str[i]]);
        base = max(2, base+1);

        int k;
        for (k = base; k <= 62; k++)
            if (dividable(k, str))
                break;

        if (k <= 62) cout << k << endl;
        else cout << "such number is impossible!" << endl;
    }

    return 0;
}

最大クリーク問題

http://poj.org/problem?id=3692

を解いていて、「最大クリーク問題」という問題について知った。

　与えられたグラフの部分グラフのうち完全グラフとなるものをクリーク（clique）と呼ぶ。節点数が最大となるクリークを求める問題を「最大クリーク問題」という。この問題はNP完全らしい。

　完全グラフの補グラフは、null graph（枝集合が空集合であるグラフ）であるので、あるグラフの最大クリークを求めることと、その補グラフの最大独立集合を求めることは同値である。

　二部グラフにおいては最小点カバーと最大マッチングが一致すること、および、一般のグラフにおいて最小点カバーと最大独立集合の和が節点数と一致することから、対象の補グラフが二部グラフである場合には最大クリークを簡単に求めることができる。

　”クリーク”という言葉は、社会学から生まれたらしい。ソーシャルグラフの中からお互いが全員知り合いであるような最大集合の大きさを求めたいというのが研究のはじまり。まさに、PKUの問題もこれと同様の問題。実際のソーシャルグラフの補グラフは必ずしも二部グラフになるとは限らないため、さまざまな近似アルゴリズムが研究されているらしい。

　PKUの問題の解答例を載せておきます。

int n;
bool edge[400][400];
bool used[400];
int pr[400];

bool match(int s) {
    used[s] = true;
    REP(i, n) {
        if (!edge[s][i]) continue;
        if (pr[i] == -1 || (!used[pr[i]] && match(pr[i]))) {
            pr[s] = i;
            pr[i] = s;
            return true;
        }
    }

    return false;
}

int main() {
    int b, g, m;

    int t = 0;
    while (cin >> b >> g >> m) {
        if (!b && !g && !m) break;

        int x, y;
        n = b + g;
        REP(i, n) REP(j, n) edge[i][j] = 1;
        REP(i, b) REP(j, b) edge[i][j] = 0;
        REP(i, g) REP(j, g) edge[i+b][j+b] = 0;

        REP(i, m) {
            cin >> x >> y;
            --x, --y;
            edge[x][y+b] = edge[y+b][x] = 0;
        }

        int ret = n;
        memset(pr, -1, sizeof(pr));
        REP(i, n) {
            if (pr[i] == -1) {
                memset(used, 0, sizeof(used));
                if (match(i)) --ret;
            }
        }

        printf("Case %d: %d\n", ++t, ret);
    }

    return 0;
}

参考：
1. simezi_tanの日記
2. Clique problem - wikipedia

サーバー奮闘記（18） gitを使う

"Cowboy Cafe - Dubois - Wyoming - July 2005" By GeordieMac Pics

　結構前にgitはローカルとvpsともに入れていたのですが、GitHubにソースを上げるときくらいしか使ってませんでした。

　vpsにマスターのレポジトリを作って、ローカルにいくつかレポジトリを作ってpushしたりpullしたりして遊んでみました。簡単なメモを残します。

1. sshのconfigファイルの設定

　sshを使ってリモートログインしたいサーバーが複数台ある場合は、configファイルを書いておくと便利です。サーバーのアドレスや、sshのポート番号、どのrsa秘密鍵を使うのか、などを記述しておくと、sshコマンドを使う際にいろいろとオプションを書かなくてよくて楽です。

~/.ssh/config

HOST vps

　HostName 123.12.1.123

　Port 4321

　IdentityFile ~/.ssh/id_rsa_vps

　User taro

　それぞれのサーバーに付き、上のような設定を書いておくと、

ssh vps

のようにHOST名を書くだけでログインできるようになります。

2. サーバー側にbareリポジトリを作る

　調べてみた感じだと、2つやり方があるようです。1つ目は、サーバー上で普通のレポジトリを作成して、それをbareリポジトリにcloneして、さらにそれをローカルにcloneするという方法。2つ目は、サーバー側にbareレポジトリを作成して、クライアント側でリモートサーバーを登録する方法です。

1つ目の方法だと、cloneした時点で、push/pull先が自動でremoteに登録されるので明示的にリモートサーバー名を追加する必要はありません。ただし、サーバー側でcloneする前に、適当なファイルを作ってcommitしておかないといけません。

参考: http://d.hatena.ne.jp/Naruhodius/20110415/1302832214

2つ目の方法は、サーバー側では

$ git --bare init

とするだけです。ただし、1つ目の方法と違って、remoteサーバー名が自動で登録されないので、サーバーに名前を付けたい場合（※）は、

$ git remote add origin ssh://vps/~taro/git/repos/hello

のようにremote サーバーを明示的に追加する必要があります。

参考: http://blog.champierre.com/670

※remoteサーバー名は必ず登録する必要はないですが、毎回push/pullする際に

$ git push ssh://vps/~taro/git/repos/hello

とするのは面倒なので、登録しておいた方が便利。

おそらく2つ目のやり方の方が一般的かと思います。

SRM 525 DIV2 950 MagicalSquare

DPで解く問題。どのような状態を保持するかが悩ましかった。最初、dp[r][i][j][k]としてr行目まで見たときに、各カラムが先頭からi, j, k文字まで埋まっている場合のパターン数を保持するようにした。しかしよく考えると、kはi, jが決まれば一意にきまるので、dp[r][i][j]とした。

　計算量は、(行数 4) * (カラムの埋まっている文字 50^2) * (現在注目している行の文字の分割パターン 50^2) * (文字列比較 50)であやしかったけど、投げてみる。TLE。ループの中にbreakを入れて枝狩りをしてみる。Accept。

　んー、通ったものの、計算量はもう一段階落ちると思う。文字列の比較をハッシュとか使ってやればいいのかなー。Editorialに期待します。

long long int dp[4][51][51];

class MagicalSquare {
public:
    long long int getCount(vector<string> rowStrings, vector<string> columnStrings) {
        memset(dp, 0, sizeof(dp));

        vector<string> rows(4);
        rows[0] = "";
        REP(i, rowStrings.size()) rows[i+1] = rowStrings[i];

        dp[0][0][0] = 1;
        int rLen[4];
        int cLen[3];

        REP(i, 4) rLen[i] = rows[i].length();
        REP(i, 3) cLen[i] = columnStrings[i].length();

        if (rLen[1] + rLen[2] + rLen[3] != cLen[0] + cLen[1] + cLen[2])
            return 0;

        int s = 0;
        REP (r, 3) {
            s += rLen[r];
            REP (i, cLen[0]+1) REP(j, cLen[1]+1) {
                if (!dp[r][i][j]) continue;

                int k = s - i - j;

                REP(p, rLen[r+1]+1) {
                    string s1 = rows[r+1].substr(0, p);
                    if (columnStrings[0].substr(i, s1.length()) != s1) break;

                    FOR (q, p, rLen[r+1]+1) {
                        string s2 = rows[r+1].substr(p, q-p);
                        string s3 = rows[r+1].substr(q);

                        if (columnStrings[1].substr(j, s2.length()) != s2) break;
                        if (columnStrings[2].substr(k, s3.length()) != s3) continue;

                        dp[r+1][i+s1.length()][j+s2.length()] += dp[r][i][j];
                    }
                }
            }
        }

        return dp[3][cLen[0]][cLen[1]];
    }
};

SRM 525 DIV1 250 dropcoins

完敗(^O^)／

どうせ探索問題だろうとおもったのが、すべての敗因でした。このやり方知ってるのに出て来なかったのは悔しい。ばか、もうほんとばか。

int cnt[31][31];
class DropCoins {
public:
    int getMinimum(vector<string> board, int K) {
        int r = board.size();
        int c = board[0].length();

        int ret = INF;
        memset(cnt, 0, sizeof(cnt));
        FOR (i, 1, r+1) FOR (j, 1, c+1) {
            cnt[i][j] = cnt[i-1][j] + cnt[i][j-1] - cnt[i-1][j-1] + (board[i-1][j-1] == 'o');
        }

        FOR (i, 1, r+1) FOR (j, 1, c+1) FOR (it, i, r+1) FOR (jt, j, c+1) {
            int x = cnt[it][jt] + cnt[i-1][j-1] - cnt[it][j-1] - cnt[i-1][jt];

            if (x == K) {
                int ans = (i-1)+(j-1)
                        + (r-it)+(c-jt)
                        + min(i-1, r-it)
                        + min(j-1, c-jt);
                ret = min(ret, ans);
            }
        }

        if (ret == INF) return -1;
        return ret;
    }
};

C++の関数オブジェクト

前から気になっていたのですが、C++にgreater<int>というものがあります。例えば、vector<int> xsを降順にソートしたい場合に

sort(xs.begin(), xs.end(), greater<int>());

のようにして使います。また、priority_queue<int>のデータ構造でtopに最小のものがくるようにしたい場合、

priority_queue<int, vector<int>, greater<int> >que;

のようにします。

　ここで気になることがあります。sort()の場合は、greater<int>()と括弧つきで使用して、priority_queueの場合は括弧なしで使用しています。この違いは一体・・？そもそもこれは関数なのか、関数ポインタなのか？

ということで調べてみました。

http://www.geocities.jp/ky_webid/cpp/library/015.html

　これは、関数オブジェクト（ファンクタ）と呼ばれるらしいです。operator()がオーバーロードされていてオブジェクトを通常の関数と同様の方法で呼び出し可能にする仕組みらしいです。ということでgreater<int>()は以下のように使うことができます。

cout << greater<int>()(3,1) << endl;        // true
cout << greater<int>()(1,2) << endl;        // false
cout << greater<int>()(-1,-1) << endl;      // false

または、こんなかんじ。こっちの書き方がオブジェクトを関数のように使ってるというのが分かりやすいかも。

greater<int> g;

cout << g(3,1) << endl;        // true
cout << g(1,2) << endl;        // false
cout << g(-1,-1) << endl;      // false

　でも何だかもやもやしています。これって何が嬉しいんだろう。関数ポインタと一緒じゃ。。ということでまた調べてみました。

http://d.hatena.ne.jp/pknight/20090818/1250584515#c

　関数ポインタと比べたときの関数オブジェクトの利点は以下の2つだそうです。

• コンパイラがインライン展開しやすい
• メンバ変数を持てるので状態を持った関数を作れる

　なるほど。。2つ目の利点が面白そうで、これってクロージャとかジェネレータとかがつくれるってことなんじゃ？と思って書いてみました。

　まず、クロージャ的なやつ（実行時に入力された数値によって挙動の異なる関数を生成したように見せかける）。

template <typename T>
class Closure {
    T pw;

public:
    T operator() (T x) {
        T ret = 1;
        for (int i = 0; i < pw; i++)
            ret *= x;

        return ret;
    }

    Closure(T pw) {
        this->pw = pw;
    }
};

int main() {
    int n, m;

    cin >> n >> m;
    Closure<int> f(n), g(m);

    cout << f(10) << endl;           // 10^n
    cout << g(10) << endl;           // 10^m

    return 0;
}

　次に、ジェネレータ的なやつ（フィボナッチ数列生成器）。

class Gen {
    int a, b;

public:
    int operator() () {
        int tmp = a;

        a = b;
        b += tmp;

        return tmp;
    }

    Gen() {
        a = 1;
        b = 1;
    }
};

int main() {
    Gen gen;

    for (int i = 0; i < 30; i++)
        cout << gen() << endl;

    return 0;
}

なんかここまでくると、関数型言語みたいなことができそうじゃなイカ！？と思ってfunction composition的なものを書いてみます。最初templateをプレースホルダーのように使って書いてみました。

template <typename _T, typename _OuterFunction, typename _InnerFunction>
class Composite {
public:
    _T operator() (_T __x) {
        _OuterFunction __g = _OuterFunction();
        _InnerFunction __f = _InnerFunction();
        return __g(__f(__x));
    }
};

template <typename _T>
class dbl {
public:
    _T operator() (_T __x) {
        return 2*__x;
    }
};

template <typename _T>
class succ {
public:
    _T operator() (_T __x) {
        return __x + 1;
    }
};

int main() {
    Composite<int, succ<int>, dbl<int> > comp1;
    Composite<int, dbl<int>, succ<int> > comp2;

    cout << comp1(10) << endl;     // ((+1).(*2)) 10
    cout << comp2(10) << endl;     // ((*2).(+1)) 10

    return 0;
}

そのあと、動的に渡すような形で書けないかなーと試行錯誤して以下のようなものを書きました。

template <typename _T, typename _OuterFunction, typename _InnerFunction>
_T composite(_OuterFunction __f, _InnerFunction __g, _T __x) {
    return __f(__g(__x));
}

template <typename _T>
class dbl {
public:
    _T operator() (_T __x) {
        return 2*__x;
    }
};

template <typename _T>
class succ {
public:
    _T operator() (_T __x) {
        return __x + 1;
    }
};

int main() {
    cout << composite(succ<int>(), dbl<int>(), 10) << endl;    // ((+1).(*2)) 10
    cout << composite(dbl<int>(), succ<int>(), 10) << endl;    // ((*2).(+1)) 10

    return 0;
}

　書いた後、気付いたのですが、compositionの上のバージョンはpriority_queueと同様、下のバージョンはsort()と同様の書き方になっていました。こうやって見ると、priority_queueの場合は関数オブジェクトのクラスを、sort()の場合は関数オブジェクトのインスタンスを渡しているということに気付きます。たしかに、priority_queue<>の<>内で指定するものは"型"、sort()の引数で渡すのは実体を持った"値、または、インスタンス"であることを考えると、どっちに括弧付けないどっちに付けないのかってのは納得できます。実際にSTLのライブラリのソースを読んでみると、priority_queueの場合、第三引数で渡す_Compare（greater<int>に対応するもの）は"型"として使われていて、こういう使い方は関数ポインタだと書きづらいのかなと思いました。

木構造のメモリ確保

トライ木を使って以下の問題にトライしたところ、TLEを食らいました。

http://poj.org/problem?id=3630

　最初に書いたコードはこれです。トライ木を構築して、語尾にあたる節点が子を持っていたらエラーとします。このエラー検出と同時に動的確保したメモリを解放しています。

struct Trie {
    bool tail;
    Trie *next[10];
    Trie() {
        tail = false;
        fill(next, next+10, (Trie*)0);
    }
};

Trie *find(char *t, Trie *r) {
    for (int i = 0; t[i]; ++i) {
        int c = t[i]-'0';
        if (!r->next[c]) r->next[c] = new Trie;
        r = r->next[c];
    }
    r->tail = true;
    return r;
}

bool validate(Trie *r) {
    bool hasCld = false;
    bool ret = true;
    for (int i = 0; i < 10; i++) {
        if (r->next[i]) {
            hasCld = true;
            ret &= validate(r->next[i]);
        }
    }

    if (r->tail && hasCld)
        ret = false;

    delete r;
    return ret;
}

int main() {
    int t, n;
    char num[12];

    scanf("%d", &t);
    while (t--) {
        scanf("%d", &n);
        Trie *rt = new Trie;
        REP(i, n) {
            scanf("%s", num);
            find(num, rt);
        }

        if (validate(rt))
            puts("YES");
        else
            puts("NO");
    }

    return 0;
}

　なんとか高速化する方法はないかと調べていたら、komiyamさんの日記におもしろいアイディアがあったので、参考にさせていただきました。毎回メモリを動的確保するのではなく、予め必要な分だけ静的領域に確保しておき、節点を生成するときは予め確保しておいたメモリのアドレスを渡すという方法です。これにより、メモリの確保/解放に係るオーバーヘッドがなくなります。また、メモリの解放をせずにすむために枝狩りができます。私の最初のソースでは、エラーを検出してもトライ木全体をトラバースしていました。これはメモリを解放するためです。メモリの解放がいらなくなると、エラーを検出した時点で探索をやめることができます。

struct Trie {
    bool tail;
    Trie *next[10];
    Trie() {
        tail = false;
        fill(next, next+10, (Trie*)0);
    }
};

Trie nodes[100000];
int ptr;

Trie *find(char *t, Trie *r) {
    for (int i = 0; t[i]; ++i) {
        int c = t[i]-'0';
        if (!r->next[c]) {
            nodes[++ptr] = Trie();
            r->next[c] = nodes + ptr;
        }
        r = r->next[c];
    }
    r->tail = true;
    return r;
}

bool validate(Trie *r) {
    for (int i = 0; i < 10; i++) {
        if (r->next[i]) {
            if (r->tail)
                return false;
            if (!validate(r->next[i]))
                return false;
        }
    }
    return true;
}

int main() {
    int t, n;
    char num[12];

    scanf("%d", &t);
    while (t--) {
        scanf("%d", &n);
        ptr = 0;

        nodes[ptr] = Trie();
        Trie *rt = &nodes[ptr];
        REP(i, n) {
            scanf("%s", num);
            find(num, rt);
        }

        if (validate(rt))
            puts("YES");
        else
            puts("NO");
    }

    return 0;
}

　以上がメイントピックですが、3点ほどちょっとしたTipsを書いておきます。

[C++のstructについて]

　C++のstructはCの構造体とは違います。C++の場合はstruct内にメンバ関数を持つことができます。classの場合アクセス修飾子を省略するとprivateになりますが、structの場合はpublicになるそうです。

http://www2s.biglobe.ne.jp/~ragnarok/program/cpp/struct.htm

[文字から数値への変換]

'0'　- '9'のcharを0 - 9のintに変換するとき、komiyamさんはおもしろいやり方をしています。

c & 15

でchar -> intに変換されます。確かに、'0'は16進数だと0x30なので、0x0fでマスクしてあげるとintに変換できます。

[ローカル変数のdelete]

最初、rootの接点だけローカル変数で宣言していて、validateでメモリ解放するとクラッシュしました。当然っちゃ当然ですがスタックに積まれた変数をdeleteしようとすると予期せぬ結果が発生します。再帰的処理で一括してメモリ解放するときなどに、ある特別なデータだけスタック変数だったとか結構ありそうな落とし穴なので注意。

http://stackoverflow.com/questions/441831/c-calling-delete-on-variable-allocated-on-the-stack

SRM 524 DIV2 1000 multipleswithlimit

modular arithmeticの絡んだ数論チックな問題（好きなタイプの問題）。余り毎にメモ化しておけばいいのはすぐ気付いたけど、最初のsubmissionでは、dp[2][10000]みたいにして、rolling table/flying tableを使って解いて、TLE。前回新しく更新されたところだけをpick upすればいいので、priority_queueを利用したダイクストラ法のようなやり方をすればいいのに気付く。これでAC。

（追記）editorial読んだら、queueでいいっぽい。最初に埋まったやつがbestになるから、確かにqueueでOK。

typedef pair<int, pair<string, int> >  ISI;
string dp[10000];

class MultiplesWithLimit {
    string format(string x) {
        if (x.length() < 9)
            return x;

        char ret[512];
        sprintf(ret, "%s...%s(%d digits)", x.substr(0, 3).c_str(), x.substr(x.length()-3).c_str(), x.length());

        return ret;
    }

    bool contain(vector<int> vec, int n) {
        REP(i, vec.size())
            if (vec[i] == n)
                return true;

        return false;
    }

public:
    string minMultiples(int N, vector<int> forbiddenDigits) {
        vector<int> nums;

        REP(i, 10) {
            if (!contain(forbiddenDigits, i))
                nums.push_back(i);
        }

        priority_queue<ISI, vector<ISI>, greater<ISI> > Q;
        Q.push(MP(0, MP("", 0)));
        string ret = "";

        REP(i, 10000) dp[i] = "";
        while (!Q.empty()) {
            string str = Q.top().second.first;
            int rm = Q.top().second.second;

            Q.pop();
            if (rm == 0 && str != "") {
                ret = str;
                break;
            }
            if (dp[rm].length() < str.length() ||
                (dp[rm].length() == str.length() && dp[rm] < str))
                continue;

            REP(i, nums.size()) {
                int rm_t = (10 * rm + nums[i]) % N;
                string str_t = str + (char)(nums[i] + '0');

                if (str_t == "0") continue;
                if (dp[rm_t] == "" ||
                    str_t.length() < dp[rm_t].length() ||
                    (str_t.length() == dp[rm_t].length() && str_t < dp[rm_t])
                ) {
                    dp[rm_t] = str_t;
                    Q.push(MP(str_t.length(), MP(str_t, rm_t)));
                }
            }
        }

        if (ret == "")
            return "IMPOSSIBLE";

        return format(ret);
    }
};

SRM 523 DIV2 1000 smallbricks31

　On top of the basement of size 1*1*w, we like to construct a building object using three types of bricks, whose sizes are 1*1*1, 1*1*2, 1*1*3, respectively. The height of the building is at most h. How many different type of buildings could you construct under some constraints? -- A concise summery is like that.

At first glance, I was like, "Umm, this is a simple bit DP problem." But it took me longer to implement the idea than I thought. Guess I have to solve this level of problem within 30 minutes to be a yellow coder. Just need more practices!!

const LL MOD = 1000000007LL;

LL dp[10+1][1<<10];
LL ptr[1<<10][1<<10];

void init(int x, int y, int b) {
    if ((1<<b) > x) return;

    init(x, y, b+1);
    if (x >> b & 1) {
        init(x, y | 1<<b, b+1);
        ++ptr[x][y | 1<<b] %= MOD;
    }
    if ((x >> b & 1) && (x >> (b+1) & 1)) {
        init(x, y | 3<<b, b+2);
        ++ptr[x][y | 3<<b] %= MOD;
    }
    if ((x >> b & 1) && (x >> (b+2) & 1)) {
        init(x, y | 7<<b, b+3);
        ++ptr[x][y | 7<<b] %= MOD;
    }
}

class SmallBricks31 {
public:
    int countWays(int w, int h) {
        memset(dp, 0, sizeof(dp));
        memset(ptr, 0, sizeof(ptr));

        REP(i, 1<<10) {
            if (i > 0) {
                int b = 0;
                while (!(i & 1<<b))
                    ++b;
                init(i, 0, b);
            }
        }

        dp[0][(1<<w)-1] = 1;
        REP(i, h) REP(j, 1<<w) {
            if (!dp[i][j]) continue;

            REP(k, 1<<w) {
                if (ptr[j][k]) {
                    LL tmp = dp[i][j] * ptr[j][k] % MOD;
                    dp[i+1][k] = (dp[i+1][k] + tmp) % MOD;
                }
            }
        }

        LL ret = 0;
        FOR (i, 1, h+1) REP(j, 1<<w) ret = (ret + dp[i][j]) % MOD;

        return ++ret % MOD;
    }
};

サーバー奮闘記（17） Hadoopインストール

"stone setting" by angela7dreams

　Hadoop（最新安定ver 0.20.203）試してみました。Pseudo-Distributedモードでexampleをいくつか動かしてみました。ほとんど公式ホームページにあるとおりにやればOKですが、いくつかひっかかったところがあったので、メモしておきます。

[SSHのオプション]

SSHのポートをデフォルトの22以外にしている場合は設定ファイルにsshのオプションを設定。

（例：ポート1234の場合）

conf/hadoop-env.sh に以下追記。

export HADOOP_SSH_OPTS="-p 1234"

[Webインターフェース]

NameNode、JobTrackerのステータス確認用のwebインターフェース

http://localhost:50070/

http://localhost:50030/

があって、プラスα的な機能かなと勝手に思い込んでいたら、このページを開かないとhdfsがきちんと機能しないらしい。

[データ保存領域の設定]

Ubuntuだと再起動後に、dfsのテンポラリーディレクトリに保存したファイルが消えてしまいます。保存先をデフォルトのtmpから変更してあげればOKです。

conf/core-site.xmlに以下の設定を追記。

<property>
    <name>hadoop.tmp.dir</name>
    <value>/home/kenji/hadoop-data</value>
    <description>hadoopのデータ保存領域</description>
</property>

[example実行]

bin/hadoop jar hadoop-examples-0.20.203.0.jar xxxx

とする。xxxxは実行したいサンプル。

wordcountを実行したい場合は、

bin/hadoop jar hadoop-examples-0.20.203.0.jar wordcount input output

のようにする。

サーバー奮闘記（16） SSHのRSA鍵設定、rm、emacsの設定など

vps初めて半年くらい経ちました。当初借りたは良いけど、何をしたらいいのか分からない状態でしたが、適当にwebサービスを作ったりしてまあまあ活用できている気がします。最近気に成っていたことがいくつかあったので、まとめて片づけました。

1. SSHのRSA鍵設定
2. aliasでrmをrm -i にする
3. emacsのバックアップファイルを作成しないようにする

　まず1。プレインテキストでユーザー情報とかパスワード流すのやばいだろう。。（※訂正：teratermのオプションにプレインテキストって書いてあったので平文で流れているのかと思いましたが、mojavy先生曰くsshではすべての通信が暗号化されているため、プレインテキストでパスワードが流れるなんてことは無いそうです。）vpsで誰かがtcpdumpとかしてたら漏れ漏れじゃないですか@_@　と前々から思っていて、いつの間にか忘れていて、セキュリティの読み物見ていたら思い出して、即実行しました。

　http://d.hatena.ne.jp/elwoodblues/20070620/1182330984

を参考にしました。

　上のサイトにあるようにローカルで鍵ペアを作成して、サーバーに公開鍵を送るようにしましょう。サーバーで作って、秘密鍵をローカルに転送とかやってたら本末転倒なので。

　あと、SSH認証を使いますっていうのをtera termの初期設定ファイルに保存しておくとログインがらくです。「設定」 -> 「SSH認証」 と進み、

1. 「ユーザ名」にユーザ名を入力
2. 「RSA/DSA鍵を使う」にチェック
3. 「秘密鍵」項目で秘密鍵のファイルを選択
4. 「OK」ボタン押下

とし、「設定」 -> 「設定の保存」を選んで、INIファイルに設定を保存します。

　2. はtypoでrm *とかやったら泣けるなーと思ってその予防策。

~/.bashrcに以下追記。

alias rm='rm -i'

　3. は、毎回チルダから始まるバックアップファイルを作っているのがemacsだと分かったので作らないように設定。

~/.emacsに以下の設定追記。

(setq make-backup-files nil)

ヨセフス数（Josephus Number）の求め方

ヨセフスの問題（Josephus problem）という有名問題があります。これの一般解を求める式があるらしいので調べてみました。

　以下では、円を組む人の数をn、スキップ数をkとします。本来は、n番目（つまり最後に）処刑される人の番号を求めればよいのですが、より汎用的にs番目に処刑される人の位置を求めることにします。

　ここのページにあるように、以下のソースコードでs番目に処刑される人の位置がわかります。

int josephus(int n, int k, int s) {
    int ret = k * s;

    while (ret > n)
        ret = ((ret - n) * k - 1) / (k - 1);

    return ret;
}

　えーーっと、よくわかりません。ということで、何故上のコードで答えがでるのか考えてみました。（結構がんばりました（汗）。）

まずは、正の整数xに対して、f_k: x -> (x * k - 1) / (k - 1)という写像が何をするのか考えます。

k = 2のとき

f_2(1) = 1

f_2(2) = 3

f_2(3) = 5

f_2(4) = 7

・・・・・

k = 3のとき

f_3(1) = 1

f_3(2) = 2

f_3(3) = 4

f_3(4) = 5

f_3(5) = 7

・・・・・

ふーむ、どうやらf_kは、自然数からなる点列(1,2,3,...)からkの倍数を取り除くような効果があるようです。本当か？ちょっと数学的に考えてみます。

f_k(x) = (k * x - 1) / (k - 1)なので、f_k(1) = (k*1 - 1) / (k - 1) = 1です。

f_k(2) = (k*2 - 1) / (k - 1) はコンピュータの整数演算では2ですが、数学上は、2余り1です。

同様に、

f_k(3) = 3余り2

f_k(4) = 4余り3

・・・・

f_k(k) = k余り(k-1)

です。f(k)は(k-1)余ってますが、(k-1)で割るのでこの余りは消えて商が1増えます。ということで、f(k) = k+1となります。これで、(1,2,3,.....k)という点列が(1,2,3,....,k-1, k+1)となりました。さらに続きを調べるために、g_k(x) = (k * x - 1) を (k-1)で割った余りが0になる場所を調べます。

g_k(1) = 0 (mod (k-1))

g_k(k) = 0 (mod (k-1))

g_k(2*k-1) = 0 (mod (k-1))

・・・・

なので、y = 1 + i * (k-1), i = 0,1,2,3,...なる正の整数yに対してg_k(y) = 0 (mod (k-1))となります。先の例で正の整数の点列からkが除去されることを確認しました。2*k-1という場所は本来2*kがある場所です。（kが無くなってるから1引かれる）。ここで同様に商の繰り上がりがおこることから、2*kは点列から除去されます。同様に3*k-2に対応する3*kも除去されます。以降同様にして、kの倍数がすべて除去されるというわけです。

　ここまでで、写像f_kがkの倍数だけを除いた整数点列を生成できることが分かりました。ソースコードを見ると、まず

ret = k * s;

としているので、retは、(1,2,3...n)の点列を拡張して、(1,2,3,....,n,n+1,n+2,....... k*s)にしたものを考えていることが分かります。ここから逆算して、k*sが区間[1,n]の整数のどれに対応しているかをみていこうというわけです。(ret - n)でnからはみ出ている分が分かります。このはみ出ている部分に対して、写像f_kを適用することで、retは消えたものを除いて数えた場合にそれがあるべき場所へとマッピングされます。これをret <= nになるまで繰り返しているというわけです。ちょっと説明足らずな感じなので、表で示します。n = 5, k = 3, s = 5とした場合、

k*s	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
k*s - n	-	-	-	-	-	1	2	3	4	5	6	7	8	9	10
f_k(k*s - n)	-	-	-	-	-	1	2	4	5	7	8	10	11	13	14

となります。例えば、処刑される順番は、k*s = 3, 6, 9, 12, 15の順です。それぞれ見て行きます。

1番目に処刑される人は、3番目の人です。

2番目に処刑されるのは、6ですが、これはnより大きいので、もともといる位置に変換します。表のf_k(k*s - n)を見ると1なので、1番目の人です。

3番目に処刑されるのは、9ですが、これも同様に対応位置に変換して5になります。

4番目に処刑されるのは、12 -> 10 -> 7 -> 2と変換されるので、2番目の人です。

同様に、最後に処刑されるのは、15 -> 14 -> 13 -> 11 -> 8 -> 4で4番目の人です。

ということで、(3, 1, 5, 2, 4)の順で処刑されるということです。4の位置に居れば助かりますね！

アセンブラを出力して遊ぶ

簡単なプログラムをアセンブラで出力して、スタックの動きを追ってみます。簡単なプログラムですが、結構面白いです。環境は、Ubuntu、gcc 4.4.3です。

　まず、Cのソースコード（test.c）。

#include <stdio.h>

main()
{
  int i, j, k;

  i = 1;
  j = 2;
  k = i + j;
  printf("%d %d %d\n", i, j, k);
}

　そして、アセンブラ（test.s）。

$ gcc -S test.c

でアセンブラのコードを出力できます。

    .file    "test.c"
    .section    .rodata
.LC0:
    .string    "%d %d %d\n"
    .text
.globl main
    .type    main, @function
main:
    pushl    %ebp
    movl    %esp, %ebp
    andl    $-16, %esp
    subl    $32, %esp
    movl    $1, 28(%esp)
    movl    $2, 24(%esp)
    movl    24(%esp), %eax
    movl    28(%esp), %edx
    leal    (%edx,%eax), %eax
    movl    %eax, 20(%esp)
    movl    $.LC0, %eax
    movl    20(%esp), %edx
    movl    %edx, 12(%esp)
    movl    24(%esp), %edx
    movl    %edx, 8(%esp)
    movl    28(%esp), %edx
    movl    %edx, 4(%esp)
    movl    %eax, (%esp)
    call    printf
    leave
    ret
    .size    main, .-main
    .ident    "GCC: (Ubuntu 4.4.3-4ubuntu5) 4.4.3"
    .section    .note.GNU-stack,"",@progbits

　いろいろ分かることがあります。

　まず、スタックのアドレスが番地の大きい方から小さいほうへと広がっているところです。これはx86マシンでは共通だと昔聞いたことがあります。

　それから、andl $-16, %esp があるので、16byteのaddress alignmentを使っているようです。FFFFFFF0でマスキングして16の倍数になる次の数をスタックポインタに代入しています。これだとベースポインタとスタックポインタの整合性があわないようにみえるんですけど、結局ベースポインタ使ってないし、リターンする時ベースポインタをスタックポインタに戻す処理も省略されてるので別にいいような気もします。

　意外だったのは、スタック内の変数にスタックポインタからの差分でアクセスしているところです。スタックポインタは実行時に移動するから、ベースポインタからの差分で変数にアクセスするのが普通かと思いましたが、そうでもないみたいです。

　あと、最後にprintfの可変長引数のところの仕組みがおもしろいです。第一引数の%dが3つあるんで、第一引数をpopしたあと、3回popすればいいんですねー。当然といえば当然ですが、実際コード上でみるとちょっとした感動がありました。printfがcallされるときのstackの状態は下のようになっています。

素数とは宇宙の創造主が人間に与えた暗号である

"On the trail" by Mr eNil

最近、素数に関するNHKのドキュメンタリーを見た。それそれは、もう、感動でした。宇宙の神秘、数学の美に魅せられました。元ネタは、『魔性の難問 ～リーマン予想・天才たちの闘い～』です。要約を書きます。

　まず、素数とは何か？それ自身と1以外に約数を持たない正の整数のことです。2,3,5,7,11,13....と無限に続きます。この素数列の並びはランダムに見えます。しかし、2人の有名な数学者が素数と自然界においてとても重要な数を結びつけました。レオンハルト・オイラーは、素数と円周率πの間に以下の関係を見出しました。

まったくランダムに見える素数に関する式を無限個掛けあわせると、なんと世界で最も美しい数πがでてくるのです。これによって、オイラーは素数と円や球体が何らかの関係を持っていることを発見しました。

　それから、カール・フリードリヒ・ガウスは、ある正数xまでに現れる素数の数は以下の式で近似できることを発見しました。いわゆる素数定理ってやつですね。

ln()は自然対数で、その底はeです。eは、カタツムリの渦巻き、台風、銀河なのどの自然界で多くみられる螺旋形状と関連のある数字です。

　オイラーとガウスの手によって、ランダムな素数が数のキングおよびクイーンと関連付けられました。そして、素数と自然界には何らかの関係があるのではと考えられるようになったのです。その考えはベルンハント・リーマンによって『何らかの関係がある』から以下のような数学的な問題へと変換されました。

ゼータ関数

の非自明な零点は一直線上に並ぶ。

ゼータ関数が零になるようなx（複素数）を複素平面上にプロットすると、なんとそれが一直線上に並ぶというのです。衝撃的ですね！

　上記のリーマン予想は、現在証明されておらず、N ≓　NP予想などと同様にミレニアム懸賞問題の一つです。今まで、ジョン・ナッシュやアラン・チューリングといった偉大な数学者、暗号科学者たちがその証明を試みましたが、失敗に終わっています。

　時は流れて、違ったアプローチでリーマン予想にチャレンジする数学者が現れました。ヒュー・モンゴメリーは、直線状にならぶ零点の間隔に興味をもちました。そして、たまたま出会った物理学者フリーマン・ダイソンに零点の間隔の式を見せた時に、ダイソンは言いました。「それはすごい！その式はウランなどの原子核のエネルギーレベルの間隔を表す式と同じような形をしている！」

素数と量子力学の世界、まったく異なる分野が一つになった瞬間でした。

　この発見を機に、多くの物理学者がゼータ関数の零点と同じ間隔をもった原子を探すようになりました。すると今度は、非可換幾何学というまったく異なる学問を研究しているアラン・コンヌが閃いたのです。「モンゴメリーの発見を発端として、物理学者たちは素数と関連をもつ空間をさがしている。そして、その空間とは非可換幾何学と完璧に一致する空間である。」当時、非可換幾何学は、量子力学を書き変え森羅万象の原理を説明する究極の法則『万物の理論』の基礎になると期待されていた学問でした。宇宙の創生と素数が繋がっていることを示す非常に興味深い話です！

　ある数学者はいいます。

「非可換幾何学を使って素数の謎がとかれるとき、森羅万象を説明する万物の理論--いわば創造主による宇宙の設計図--も明らかになるだろう。」と。

Google Developer Day Japan 2011

　『Webの世界はすさまじい速度で進化してきた。次世代のWebがどのようなものになるのか私は分からない。しかし確かなことがある。次世代の新しいWebはこの部屋にいる技術者たちによって作られだろう。そして未来のビジネスはあなた方がつくった新しい仕組みのうえで展開される。』

　―先日、Google Developer Day 2011が横浜で開催されました。Googleの製品へ高い貢献をしているDeveloper、およびDevQuizで高得点を取り招待状をゲットしたエンジニア、約2000人が一同に会しました。

　Google Developer Dayは午前の部と午後の部に大きく分かれていて、午前はKeynote Speech、午後は複数のSessionが同時進行するという形で構成されていました。主に、

• Android
• Chrome
• App Engine
• HTML5
• Google+

などの技術の最新動向が発表されました。また、技術者同士が会話をしたり、Googleのエンジニアとコミュニケーションを取ったりできるような工夫が施されており、とても楽しいイベントでした。私自身も、缶バッジ交換のおかげでゲーム業界のエンジニアの人とお話することができ非常にいい経験となりました。

　特に私の印象に残ったのは以下の3点です。

1. ScalabilityとAgility
2. Googleマインド
3. エンジニアとしての社会貢献

　1. 『これからのWebのサービスにおいて重要なものは、どうスケールするか、そして如何に柔軟な開発スタイルを取り入れていくかだ。』特に、アジャイル開発の重要性に関しては、Googleの東日本大震災に対する取り組みの中でも語られました。

　アジャイル開発を行うためには、プロトタイピング力を持ったエンジニアの存在が絶対条件だと思います。しかし、浮かんだアイディアをすぐに実装できるエンジニアは日本には決定的に少ないように感じます。これはそもそも、日本のシステム業界（とりわけSIer業界?）が要件定義/設計/開発/テスト/保守のように縦割りになっていることが原因じゃないかと。。要件定義はできるけど、fizzbuzzも書けないエンジニア。開発はできるけど要件がまったく分からないエンジニア。アジャイル開発を取り入れるためには、まずこのあたりの体質を変えていく必要があるのかなと思いました。

　2. これは、単純にいいな～と（羨ましいなと）。。会社のスローガンって見栄えだけ綺麗で中身が伴っていないことがしばしばありますが、Googleが大切にしているという以下の三箇条は本物だなと感じました。

• なにごともエンジニアありき
• 百聞は一デモに如かず
• 日本で「イケる！」と思ったら、世界のみんなも同感するかも

　3. 「Googleのクライシスレスポンス」というSessionで、東日本大震災に対するGoogleの取り組みが語られました。『すべてが自動化された完璧なシステムを届ける必要はない。途中に人の手が入っていてもかまわない。』『当初、被災地の人達が欲していたのは、スマートフォン上で起動するアプリでもクラウド上で動作するアプリでもなかった。エクセルやアクセスといったごく単純な技術だった。』

　技術者として最新技術を追うことや完璧を求めることも大事だけど、今自分が持っている技術がどう世の中に役に立つのかと考えて実践するというのは本質的により重要なことだなと改めて考えさせられました。自分の技術なんて大したことないやと卑屈にならずに、少しでも世の中が変わる、少しでも世の中に役に立つはずと信じてサービスを提供していきたいなと思いました。

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())