約数が1000個ある最小の数は？

どこかの酔っぱらいに問題を出されたので、解いてみた。

「約数が1000個ある最小の自然数はいくつか？但し答えはintの範囲におさまる。」

答えは、810810000。あってますかね？

ある自然数nの約数の数は、nの素因数p1, p2, ..., pmを用いて

n = p1^(k1+1) * p2^(k2+1) * p3^(k3+1) *....., pm^(km+1)

と書けるので、これを使ってメモ化再帰した。

view plainprint?

1. const double inf = (1e+300);  
2. double dp[1024][1024];  
3. bool isPrime[1<<20];  
4. vector<int>primes;  
5. vector<int>facs;  
6.   
7. void init() {  
8.     REP(i, 1<<20) isPrime[i] = true;  
9.     isPrime[0] = isPrime[1] = false;  
10.   
11.     for (int i = 2; i * i <= (1<<20); i++) {  
12.         if (!isPrime[i])  
13.             continue;  
14.         for (int j = i*i; j < (1<<20); j += i)  
15.             isPrime[j] = false;  
16.     }  
17.     primes.clear();  
18.     REP(i, 1<<20) if (isPrime[i]) {  
19.         if (primes.size() > 1024)  
20.             break;  
21.         primes.push_back(i);  
22.     }  
23. }  
24.   
25. double rec(int pw, int pr) {  
26.     if (pw == 1)  
27.         return 1;  
28.   
29.     if (dp[pw][pr] > 0)  
30.         return dp[pw][pr];  
31.   
32.     double ret = inf;  
33.     REP(i, facs.size()) {  
34.         if (pw % facs[i] == 0)  
35.             ret = min(ret, pow(primes[pr], facs[i]-1) * rec(pw/facs[i], pr+1));  
36.     }  
37.   
38.     return dp[pw][pr] = ret;  
39. }  
40.   
41. int main() {  
42.     int n;  
43.     cin >> n;  
44.   
45.     init();  
46.     FOR (i, 2, 1001) {  
47.         if (n % i == 0)  
48.             facs.push_back(i);  
49.     }  
50.   
51.     REP(i, 1024)REP(j, 1024) dp[i][j] = -1;  
52.     cout << (int)rec(n, 0) << endl;  
53. }  

const double inf = (1e+300);
double dp[1024][1024];
bool isPrime[1<<20];
vector<int>primes;
vector<int>facs;

void init() {
    REP(i, 1<<20) isPrime[i] = true;
    isPrime[0] = isPrime[1] = false;

    for (int i = 2; i * i <= (1<<20); i++) {
        if (!isPrime[i])
            continue;
        for (int j = i*i; j < (1<<20); j += i)
            isPrime[j] = false;
    }
    primes.clear();
    REP(i, 1<<20) if (isPrime[i]) {
        if (primes.size() > 1024)
            break;
        primes.push_back(i);
    }
}

double rec(int pw, int pr) {
    if (pw == 1)
        return 1;

    if (dp[pw][pr] > 0)
        return dp[pw][pr];

    double ret = inf;
    REP(i, facs.size()) {
        if (pw % facs[i] == 0)
            ret = min(ret, pow(primes[pr], facs[i]-1) * rec(pw/facs[i], pr+1));
    }

    return dp[pw][pr] = ret;
}

int main() {
    int n;
    cin >> n;

    init();
    FOR (i, 2, 1001) {
        if (n % i == 0)
            facs.push_back(i);
    }

    REP(i, 1024)REP(j, 1024) dp[i][j] = -1;
    cout << (int)rec(n, 0) << endl;
}

この結果は面白くて、約数の個数が素数の場合は、その数は合成数個の約数を持つ数と比べて大きな数になることを意味している。

例えば、11個の約数を持つ数。

11はこれ以上因数分解できないため、11個の約数を持つ最小の自然数は上の約数の個数の式より、

2^(10) = 1024となる。

これに対して、12 = 3 * 2 * 2と因数分解できるので、12個の約数を持つ最小の自然数は、

2^2 * 3^1 * 5^1 = 60

となる。

約数の数が1000個である最小の自然数と、約数の数が101個である最小の自然数はどちらが大きいでしょうか？

という問いを考えると、約数の個数が多い1000個の方が大きく思える。しかし実際は、

1000個の場合は、810810000。（答えがあってれば。）

101個の場合は、2^100 ≓　10^30

と101個の方がケタ違いに大きな数となる。

文字列に関するcoding tips

ちょっとしたtipsを紹介します。

1. accumulate()で文字結合
2. mapは使わない
3. 文字をint型に代入する

まず、1.番目。vectorで与えられた入力をconcatenate して新しい文字列を作成せよ。よくありますね。これをどう書くか。

view plainprint?

1. vector<string> x;  
2.   
3. void init() {  
4.     x.push_back("acchan");  
5.     x.push_back("yuko");  
6.     x.push_back("yukirin");  
7.     x.push_back("mayuyu");  
8. }  
9.   
10. int main() {  
11.     init();  
12.   
13.     string y = "";  
14.     REP(i, x.size())  
15.         y += x[i];  
16.   
17.     cout << y << endl;  
18.     return 0;  
19. }  

vector<string> x;

void init() {
    x.push_back("acchan");
    x.push_back("yuko");
    x.push_back("yukirin");
    x.push_back("mayuyu");
}

int main() {
    init();

    string y = "";
    REP(i, x.size())
        y += x[i];

    cout << y << endl;
    return 0;
}

でもいいけど、accumulate()を使うとちょっとだけかっこよく書けます。accumulate()はdefaultでは演算子+を用いる。stringに対する+は文字の結合なのでそのまま使えます。

view plainprint?

1. int main() {  
2.     init();  
3.   
4.     string y = accumulate(ALL(x), string());  
5.   
6.     cout << y << endl;  
7.     return 0;  
8. }  

int main() {
    init();

    string y = accumulate(ALL(x), string());

    cout << y << endl;
    return 0;
}

次に、2.番目。文字の出現回数を数えています。char -> intへの写像なので、mapを使うのが自然に見えます。

view plainprint?

1. int main() {  
2.     string x = "riverponiiteerutoshushuheavyrotationeverydaykachusha";  
3.   
4.     map<char, int> cnt;  
5.     REP(i, x.size())  
6.         cnt[x[i]]++;  
7. }  

int main() {
    string x = "riverponiiteerutoshushuheavyrotationeverydaykachusha";

    map<char, int> cnt;
    REP(i, x.size())
        cnt[x[i]]++;
}

でも、charなんて高々[0, 256)までなので、配列に入れるのもありです。C++の配列はランダムアクセスなので、こっちの方が計算量は安いです。[0, 256)と書きましたが、ASCIIコード表を見ると、実際に使われるコードは[0, 128)で十分のようです。AKBのメンバーの名前を覚える暇があったら、ASCIIコードを覚えましょう。（→自分）

view plainprint?

1. int main() {  
2.     string x = "riverponiiteerutoshushuheavyrotationeverydaykachusha";  
3.   
4.     int cnt[128];  
5.     memset(cnt, 0, sizeof(cnt));  
6.     REP(i, x.size())  
7.         cnt[x[i]]++;  
8. }  

int main() {
    string x = "riverponiiteerutoshushuheavyrotationeverydaykachusha";

    int cnt[128];
    memset(cnt, 0, sizeof(cnt));
    REP(i, x.size())
        cnt[x[i]]++;
}

最後に、3.番目。これは、文字をインクリメントして走査したい場合です。それぞれの数字がdigit中に何回出てくるかを標準出力に表示します。

view plainprint?

1. int main() {  
2.     string digit = "1234254723075027659279743057423953";  
3.   
4.     for (char i = '0'; i <= '9'; i++)  
5.         cout << count(ALL(digit), i) << endl;  
6. }  

int main() {
    string digit = "1234254723075027659279743057423953";

    for (char i = '0'; i <= '9'; i++)
        cout << count(ALL(digit), i) << endl;
}

とか、

view plainprint?

1. int main() {  
2.     string digit = "1234254723075027659279743057423953";  
3.   
4.     REP(i, 10)  
5.         cout << count(ALL(digit), i+'0') << endl;  
6. }  

int main() {
    string digit = "1234254723075027659279743057423953";

    REP(i, 10)
        cout << count(ALL(digit), i+'0') << endl;
}

とかやっていましたが、以下のようにできますね。プロコン界伝統の悪しきマクロで書いてますが、iはintで持っています。charではなく、intで持っていた方が何かと便利です。

view plainprint?

1. int main() {  
2.     string digit = "1234254723075027659279743057423953";  
3.   
4.     FOR (i, '0', '9'+1)  
5.         cout << count(ALL(digit), i) << endl;  
6. }  

int main() {
    string digit = "1234254723075027659279743057423953";

    FOR (i, '0', '9'+1)
        cout << count(ALL(digit), i) << endl;
}

ベイズ推定と最尤推定の違い

ベイズ推定と、最尤推定の違いについて。

まず、両者の違いを説明する前に、確率で使われる用語の説明を少しします。

観測データ: 実際に観測された事象のことです。

例）サイコロを振ったら1が出たとか、5が出たとか。

仮説: 観測されたデータが満たしていたであろうと考えられる条件。

例）投げたサイコロは、すべての目がでる確率が均等だった。

実は偶数の目が出やすいように細工されたサイコロだった。

以下では、観測データをD, 仮説をHi, i = 1,2,..., m.とおきます。

尤度: ある仮定Hiのもとで、事象Dが起こる確率です。条件付き確率で、P(D | Hi)と書きます。

　　　例）すべての目が均等に出るサイコロを振って6が出る確率。

偶数の方がでやすいように細工されたサイコロを振って6が出る確率。など。

事後確率：　ある事象Dが発生した場合、仮説Hiが正しい確率。条件付き確率で、P(Hi | D)と書きます。

例）サイコロを振って6がでたとき、そのサイコロは均等に目がでるサイコロであった確率。

サイコロを振って6がでたとき、そのサイコロは偶数が出やすいようにいかさまされたサイコロであった確率。

この尤度と事後確率の違いは大事です。

最尤推定は、尤度を最大化する仮説を選択する推定法です。これに対して、ベイズ推定は、事後確率を最大化する仮説を選択する推定法です。

ベイズの定理より、

P(Hi | D) = P(D | Hi) * P(Hi) / P(D)

が成り立ちます。この式から、事後確率とは、尤度に事象Hiが起こる事前確率（※）をかけた値に比例することが分かります。P(H1) = P(H2) = ... = P(Hn)が成り立っている場合、 最尤推定とベイズ推定は同じ結果を返します（後半の例題で具体的に説明します。）。このP(Hi), i=1,2,..n.をどうとるかというのは難しい問題で、主観的な判断に任せるをえない場合があります。

※事前確率とは、他になんの情報もない状況で、ある事象が観測される確率のことです。これに対して、事後確率は、ある事象が観測されたことをヒントにして事前確率をより精度の高いものに修正したものであると考えることができると思います。例えば、モンティーホールとかはこのことを実証する良い例だと思います。

実際に例を見てみましょう。参考サイト2に簡単な例題があったのでそれを使います。

部屋Aと部屋Bがあります。

部屋Aには、男1人、女1人がいます。部屋Bには女が2人います。

部屋の外にいるあなたが、「おーい」と呼びかけたら、女の人の声がしました。あなたの呼びかけに答えたのはAの部屋の人？Bの部屋の人？（どちらの方が確率が高いでしょう？）

まず、最尤推定。

女の人が返事をするという事象をDとします。返事をしたのがAの部屋の人であるという仮説をHa, Bの部屋の人であるという仮説をHbとします。

それぞれの尤度は、

P(D|Ha) = 0.5

P(D|Hb) = 1

です。尤度が高いのは、仮説Hbなので、あなたの呼びかけに答えたのは、Bの部屋の人であると判断します。

すごい簡単な例ですが、おもしろいです。

何も情報が与えられなかったら、どちらの部屋の人が答えたのか分かりません。強いて答えろと言われたら、ランダムにAかBを選ぶでしょう。（この場合は、P(Ha) = P(Hb) = 0.5）

そこに、"女の人が答えた"という事象が観測され、それをヒントにすることで、この事前確率P(Ha), P(Hb)をより正しい（尤もらしい）方向に修正できたと言えます。

次に、ベイズ推定で解いてみます。

P(Ha), P(Hb), P(D|Ha), P(D|Hb)は最尤推定で出したので、P(D)を求めます。（実はP(D)はiによらない値なので、本来求めなくてもいいです。が、ここでは計算しておきます。）

P(D) = P(D|Ha) * P(Ha) + P(D|Hb) * P(Hb) = 0.5*0.5 + 1*0.5 = 0.75

この場合は4人の人それぞれが返事をする確率が一様なので、部屋のことは忘れて単純に4人いてそのうち3人が女性なので、3/4と出してもいいです。

Ha, Hbの事後確率を求めます。

P(Ha) = P(D|Ha) * P(Ha) / P(D) = 0.5 * 0.5 / 0.75 = 1/3

P(Hb) = P(D|Hb) * P(Hb) / P(D) = 1 * 0.5 / 0.75 = 2/3

で、ベイズ推定の場合も仮説Hbを解に選択します。

上の式を見て分かるように、P(D)はどちらにもかかってるfactorなので、無視できます。また、P(Ha) = P(Hb) = 0.5なので、このfactorも無視できます。

すると、事後確率と尤度が同じ（スケールが違うだけ）になっています。これ（事前確率が一様であること）が、最尤推定とベイズ推定の結果が同じになる条件です。

じゃあ、事前確率が一様じゃないときは？やってみましょう。

Aの部屋の人が答える確率と、Bの部屋の人が答える確率を無意識のうちに0.5ずつとしていましたが、以下の条件が与えられたらどうでしょう？

Aの部屋の2人はとても社交的である。Bの部屋の2人は内向的で部屋に引きこもりがちで居留守をよく使うらしい。

この条件をもとに、以下のような事前確率を設定します。（このように事前確率の設定は主観的なものになりがちです。このことが最尤派とベイズ派の論争の中心となってるようです。[参考サイト1]）

P(Ha) = 0.7

P(Hb) = 0.3

最尤推定は、事前確率を判断材料に用いないので前回の結果と変わりません。

ベイズ推定の場合は、結果が変わります。

P(D) = P(D|Ha) * P(Ha) + P(D|Hb) * P(Hb) = 0.5*0.7 + 1*0.3 = 0.65

P(Ha) = P(D|Ha) * P(Ha) / P(D) = 0.5 * 0.7 / 0.65 = 0.54

P(Hb) = P(D|Hb) * P(Hb) / P(D) = 1 * 0.3 / 0.65 = 0.46

で、P(Ha)の方が大きいので、Aの部屋の人が返事をしたと判断します。このように、事前確率が一様でない場合は、最尤推定とベイズ推定の結果は異なります。ベイズ推定の方がより多くの情報を活かした推定となっています。しかし、今回のように主観的な確率を用いた場合は、その結果が本当に客観的に確かなのか疑問が残ります。また、事前確率が統計的に求められる場合にも、このfactorが問題を引き起こす場合もあるようです[参考サイト3]。

ちょっと長いpostになってしまいましたが、読んでいただきありがとうございます。間違ってる場所がありましたら、ご指摘いただけると助かります。

参考サイト：

1. 「ベイズ統計入門＆モンテカルロ法と逆問題」

2. 「事後確率と尤度――系統推定における最尤法とベイズ法の最前線」

3. 「ベイズ推定　臨床検査における偽陽性」

C言語の4つのメモリー領域

ちょっと気になることがあったので勉強しなおした。プロコンやってる人はヒープか、スタックなのかはよく気にすると思います。ヒープはメモリ領域が広く、スタックはメモリ領域が狭いからです。スタックのメモリオーバーしてクラッシュするとかよくあります。

　例えば、1000*1000のint[]配列をmain()関数内で宣言すると、プログラム実行時にクラッシュします。そこで、この1000*1000の配列を外部変数にしてあげます。すると、プログラムはクラッシュしません。

これは、メモリ領域が広いヒープに格納されているからだと思っていましたが、グローバル変数が格納される場所は”静的領域”と呼ばれるそうです。

ということで、勉強したことを書きます。

まず、Cには4つの記憶領域があります。

• プログラム領域
• 静的領域
• ヒープ
• スタック

　プログラム領域は、プログラム（マシン語）が格納される場所です。「なんじゃそりゃ？」と思われる方もいるかと思いますが、ノイマン型コンピュータでは、プログラム内蔵方式でプログラムが実行されます。つまり、命令自体もメモリに格納されていて、プログラムカウンタという現在実行している命令が格納されているアドレスを保持するものが参照している場所の命令が実行されます。

メモリに命令を書くと言う感覚は分かりにくいですが、アセンブリをやってみると分かります。

　次に、静的領域です。静的領域は、その名のとおりプログラムの実行中に変化しない変数を格納する領域です。具体的には、グローバル変数やstatic変数が置かれます。この静的領域の大きさはプログラム実行時に確保され固定サイズです。（プログラム実行中にサイズが大きくなったり、小さくなったりしません。）

　そして、ヒープ。ここは、new とか malloc() とかで動的確保するところです。ヒープのサイズはプログラム実行中も変化します。ただし、一般に１つのアプリケーションが使用できるメモリは2GBまでですので、これ以上のサイズのメモリを確保することは普通はできません。

　最後に、スタック。これは局所変数が置かれるところです。スタックのサイズは固定です。このポストの冒頭で書いたように大きなメモリをスタック上に確保すると、すぐにクラッシュします。また、一般にスタックのメモリは大きい方から小さいほうに伸びます。逆にヒープの場合は、小さい方から大きいほうに伸びます。具体的にいうと、スタックの場合は、関数内でint x, int yを宣言した場合、後に宣言したyの方が若いアドレスを持っています。

　んー、ただこれはC言語の構造上のメモリ分類であって、それが物理的にハードウェア上でどのように実現されているかというところまで抑えないと、なんか本質を理解してないような気がする。。

でも、少なくとも上の言葉の分類を知っていれば、技術者同士のコミュニケーションが円滑に進みます。そういう意味では、言葉を覚えるというのは大事ですね。

P2P通信NAT超え調査メモ

前回、Javaでソケット通信をして遊んだ。

簡単なP2Pチャット（片方のノードがグローバルIPを持っている OR LAN内での使用限定）を作ってみたり、接続したリモートPCを遠隔操作して遊んだりしたけど、やはりNATを超えたい衝動が抑えられないので、ちょっと勉強してみた。

調べた感じだと、たぶん自分で実装できそう。

今後のためにメモしておく。

[リンク集]

NATの仕組み、パケットのヘッダデータ書き換え

http://www.n-study.com/network/nat.htm

IPマスカレード(NAPT)

http://www.asi.co.jp/techinfo/unix/nat.html

NATトラバース(STUN)

http://www.itmedia.co.jp/enterprise/articles/0505/30/news070_3.html

ルータのタイプいろいろ

http://www.skame.nu/P2Pmeeting/0409sunouchi.pdf

Full Cone 開いたポートに送られるパケットは送信元がどこでも受け取る

Restricted Cone 送ったことのあるIPアドレスから送られるパケットのみ受け取る

Port Restricted Cone 送ったことのあるIPアドレス/ポートから送られてくるパケットのみ受け取る

Symmetric 最初に送った宛先以外からは受け取らない （同一ソケットから別の宛先に送った場合は、ルータに別のポートがマップされる）

[考えたこと]

P2PでNATの奥の端末同士が通信している間も、STUNとの通信はオープンになっている？

一度、STUNサーバーにつなげば、通信したい相手側のルーターのglobal IP、portが分かるので、

そこに捨てパケットをなげてやる。これで、次回からは通信可能となる。

-> Restricted Cone、Port Restricted ConeでもNAT超えられる。

P2P間で接続ができたら、その時点でSTUNサーバーとのコネクションはクローズしてもいいように思える。

開けたままだと、通信の内容がSTUNサーバーにも漏れてしまうことになりそうなので、

プライバシーの問題からおそらくクローズするの妥当だろう。この辺は実験して確認する。symmetricでは、NAT traverseは無理。

蛙逢引 ～一次不定式で愛を探す～

蛙が逢引するロマンティックな問題。この問題は、問題文的にも、内容的にも好きだ。

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

ここに訳してくれている人がいる。http://d.hatena.ne.jp/kabus/20060728

以下問題の解き方。

x + mt = y + nt (mod L)より、

x + mt = y + nt + kL

(n-m) t + Lk = x - y となるような整数(t, k)の組を見つければいい。（問題では、最小の正の整数t*を求めることになるが、これは求まったtに対する代表元を考えればいい。）

このような解が存在する必要十分条件は、(x-y)がgcd(n-m, L)で割り切れることである。

まず、必要条件の証明。

もし、(x-y)がgcd(n-m, L)で割り切れない場合に、(t, k)が存在したとする。

このとき、

左辺 ≡ 0 (mod gcd(n-m, L))に対して、右辺 > 0 (mod gcd(n-m, L))となり矛盾。よって(x-y)がgcd(n-m, L)で割り切れない場合は、(t, k)の組は存在しない。

次に、十分条件の証明。

両辺をgcd(n-m, L)で割る。

(n-m)/gcd(n-m, L) をA、L/gcd(n-m, L)をB、(x-y)/gcd(n-m, L)をCとおくと、

At + Bk = Cと書ける。このとき、AとBは互いに素なので、Aを法にしたときのBの逆元が存在する。

つまり、y ≡ 1/B (mod A)となるようなyが存在する。

よって、1 - By ≡ 0 (mod A)。これは、(1 - By)がAで割り切れることを意味しているので、

x = (1 - By) / Aとなる整数xが存在する。これを変形すると、Ax + By = 1となり、この式を満たす整数(x, y)が存在する。これを両辺C倍してやると、A(Cx) + B(Cy) = Cとなり、(Cx, Cy)が(t, k)に他ならない。

ということで、解が存在するか否かは、gcdを出すだけで良いです。加えて、(n-m) t + Lk = x - y を満たす(t, k)の組は、(n-m)t' + Lk' = gcd(n-m, L)を満たす(t', k')を拡張ユークリッド互除法で出し、gcd(n-m, L)と(x - y)の比に合わせてスケーリングしなおせばOK。

最後に、問題となるのは、求まったtに対して、それと合同な最小の正の整数t*をどう出すか？そもそも何を法にして考えるか？これは、目盛りがL個付いている時計を(n-m)目盛りずつジャンプしていくと、何回ジャンプすると、同じ場所に戻ってくるかを考えればいいです。L/gcd(L, n-m)が答えです。ということで、tの法(L/gcd(L, n-m))における代表元を出せばOK。

以下回答ソース。

view plainprint?

1. LL extgcd(LL a, LL b, LL &x, LL &y) {  
2.     if (b == 0) {  
3.         x = 1;  
4.         y = 0;  
5.         return a;  
6.     }  
7.   
8.     int d = extgcd(b, a%b, y, x);  
9.     y -= (a/b) * x;  
10.     return d;  
11. }  
12.   
13. LL solve(LL x, LL y, LL m, LL n, LL L) {  
14.     if (m == n)  
15.         return -1;  
16.   
17.     LL k, t;  
18.     LL g = extgcd(n-m, L, k, t);  
19.     if ((x-y) % g)  
20.         return -1;  
21.   
22.     LL ans = (x - y) / g * k;  
23.     L /= g;  
24.   
25.     return (ans % L + L) % L;  
26. }  
27.   
28. int main() {  
29.     LL x, y, m, n, L;  
30.   
31.     cin >> x >> y >> m >> n >> L;  
32.   
33.     LL ret = solve(x, y, m, n, L);  
34.     if (ret < 0)  
35.         puts("Impossible");  
36.     else  
37.         printf("%I64d\n", ret);  
38.   
39.     return 0;  
40. }  

LL extgcd(LL a, LL b, LL &x, LL &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }

    int d = extgcd(b, a%b, y, x);
    y -= (a/b) * x;
    return d;
}

LL solve(LL x, LL y, LL m, LL n, LL L) {
    if (m == n)
        return -1;

    LL k, t;
    LL g = extgcd(n-m, L, k, t);
    if ((x-y) % g)
        return -1;

    LL ans = (x - y) / g * k;
    L /= g;

    return (ans % L + L) % L;
}

int main() {
    LL x, y, m, n, L;

    cin >> x >> y >> m >> n >> L;

    LL ret = solve(x, y, m, n, L);
    if (ret < 0)
        puts("Impossible");
    else
        printf("%I64d\n", ret);

    return 0;
}

PKU Candy Distribution

This is about number theory and I like these kinds of problems:

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

Since in the problem, N is quite large, you need some analysis to solve it in a given time.

The first observation is how the number would increment if you thought it in the linear system, not in the circular system.

0, 1, 3, 6, 10, 15, 21, 28, ....

It's obvious that the sequence above is sums of an arithmetic sequence, more specifically you'll notice that p(i) = i * (i+1) / 2.

The second observation is its cycle. In modular arithmetic, there's almost always some cycle. Notice that the differences of the sequence above is 1, 2, 3,...., so the minimum possible cycle, --what we are discussing here is the cycle of jumping--, is n, since 0, 1, 2,...n-1 is equivalent to n, n+1, n+2, ... , 2n-1 in the modulo n.

Now the hardest part -- at least in my case--. The cycle of jumping is n as I explain above. So what's the cycle of position. The cycle of this system is the gcd of these two cycles.

Let's consider it by separating two situations: when n is odd, and n is even.

When n is odd, p(n) = n * (n+1) / 2. Since (n+1) is dividable by two, p(n) = 0 (mod n).

While when n is even, p(n) = n * (n+1)/2 = n/2 * n + n/2 = n/2 (mod n).

So the cycle of position is n and 2*n, for an odd number n and an odd number n, respectively.

So it's enough to think just 2*n moves, but it's still a big number -- the biggest possible number is 10^9--, and more analysis is required.

Let's consider when n is odd. p(0) = 0, and p(n-1) = (n-1) * n/2 = 0 (mod n). So for i in {0, 1, 2, ..., n-1}, two number, 0 and n-1, have the same value when mapped onto p. It's impossible to cover all numbers on p with {0, 1, ... n-1}. And since the cycle for odd numbers is n, we can say that the answer for test cases where n is an odd number is "No".

How about when n is even? As discussed in the previous paragraph, p(n) = n/2 (mod n). So for i in {n, n+1, n+2, ....., 2*n-1}, p(i) covers the opposite side of number set against {p(i) | i in {0, 1, ...n-1} }. Therefore it's possible to cover all the numbers if p(n/2) covers all the numbers.

We come to a conclusion that the answer ret(n) is:

"Yes" when n is 1,

"No" whe n is an odd number,

ret(n/2) when n is an even number.

Which comes down to the fact the answer is "Yes" if and only if n is two to the power of i, i = 0, 1, 2,.... . The implementation is pretty easy. Just count the bit set to 1. You can use __builtin_popcount() or n & (n-1) (*).

(*) n & (n-1) sets the LSB of n to zero. In this manner, n & (n-1) is zero iff n is two to the power of i.

Javaで簡易チャットを作る

Javaのソケット通信で遊んでみた。

以前Cでソケット通信をやってみたけど、手続きがややこしくて挫折してしまった。

Javaの場合は、Cと比べ簡単にソケット通信ができる。

ということで、簡易チャットを作って遊んでみた。

ソケットの読み/書きを、SocketReader、SocketWriterというクラスで実装した。チャットの読み/書きは非同期に独立して行われるべきなので、この2つのクラスはスレッドにした。

ソースはこちら。

実行ファイルはこちら。

ネットワークプログラミングはあまりやったことがなかったけど、とても楽しかった。上のプログラムは、いろいろな応用テーマが考えられる。

1. ソケットの暗号化　シーザー符号とかでもいいけど、フェルマーの小定理を用いた暗号化がおもしろそう。　具体的には、両者間で素数pとべき数pwを決めておく。送信者側で、f : x -> x^pw (mod p)の変換をしてデータを送り、受信者側では、f : x -> x^(p-pw) (mod p)でデータを復元する。あとは、RSAとかを自分で実装してみるとか。(※)
2. リモートマシンの制御　ソケットでコマンドを送信してあげて、受信側で別プロセスを立ち上げてコマンドを実行する。コマンドの結果をソケットに流すようにすると、別のマシンを制御することができる。これは、自分のサーバーとローカルマシンで試したみたんだけど、いろいろと危ない。
3. NAT超え　通信を行うコンピュータが両方ともグローバルIPを持っていれば問題なく通信ができるが、そうでない場合は、ルータでポートフォワーディングの設定が必要となる。ポートフォワーディングの設定をしたとしても、ルータ側のIPを知っている必要がある。仲介サーバとかをつかえば出来そうだけど、勉強中。

（※ちょっと勉強したので追記）

RSAは、フェルマーの小定理というよりそれを一般化したオイラーの定理を利用している。上の自分が考えたやり方だと、公開鍵方式にはならない。。

合成数(x = p1 * p2)に対してオイラー関数（totient functionのこと）が、(p1-1) * (p2-1)となるのが嬉しい。p1とp2が分からないと、a^pw = a (mod x)となるpwが分からないってところがポイント。

参考：サルにもわかるRSA暗号

C++のConcept Checkについて

STLのソースを読んでいて、Concept Checkという概念について知った。

C++でジェネリックプログラミングをするときは、templateを使って引数の型に依存しない関数を定義する。この際、templateとして使用できる型が満たすべき条件を"concept"と呼ぶ。C++には、明示的なconcept checkの機能がないため、template関数に予期しない型を使用した場合、コンパイルが吐き出すエラーが意味をなさないものになりがちである。

例えば、以下の例を考える。

view plainprint?

1. int main()  
2. {  
3.     vector<std::complex<float> > v;  
4.     stable_sort(v.begin(), v.end());  
5. }  

int main()
{
    vector<std::complex<float> > v;
    stable_sort(v.begin(), v.end());
}

上のソースをコンパイルすると、長くて意味不明なエラーが出力される。（参考ページ1.参照）

これには、ちょっとびっくりさせられる。

エラーが発生した理由は、complexが演算子<をサポートしていないからである。しかし、templateとは単なるplaceholderにすぎないので、stable_sort()の引数であるiteratorが満たすべきconcept "LessThanComparable"をチェックすることができない。結果、深いところまで潜っていって矛盾がでたところでよく分からないあまり意味のないエラーが大量に吐き出される。

これを解決するには、コンパイル時にconcept checkを有効にしてあげるといい。

g++ -o hogehoge -D_GLIBCXX_CONCEPT_CHECKS hogehoge.cpp

-Dxxxxxとは、マクロxxxxを定義してプリプロセッサを実行するという意味。concept checkを有効にすると、LessThanComparableに関する内容がエラーに出力される。

参考ページ:

1. The Boost Concept Check Library (BCCL)

2. gccのオプション

便利なマクロの使い方

STLライブラリのソースを読んでたら、見慣れないマクロがいくつかあったので、勉強しなおした。

便利そうなのをまとめておく。

マクロ名	展開される値
#	変数名
__FILE__	ソースファイル名
__func__	関数名
__LINE__	ファイル内での行数

#はSRMでよく使用する。他の3つはより大規模な開発で重宝する。さらに、マクロは引数を(...)とすると、可変長引数を使用することができる。その際、可変長引数の参照は、__VA_ARGS__で行う。

以下、サンプル。

まずは、#を使った例。

view plainprint?

1. #define SHOW(x) cout << #x << ": " << x << endl;  
2.   
3. int main() {  
4.     int x = 10;  
5.   
6.     ++x;  
7.     SHOW(x);      // x: 11  
8. }  

#define SHOW(x) cout << #x << ": " << x << endl;

int main() {
    int x = 10;

    ++x;
    SHOW(x);      // x: 11
}

次に、ログ出力。

view plainprint?

1. #define DEBUG_PRINT(fmt, ...) \  
2.     printf("In the function %s(), at line %d in %s\n", __func__, __LINE__, __FILE__); \  
3.     printf(fmt, __VA_ARGS__)  
4.   
5. int doit(int x, int y) {  
6.     int z = x + y;  
7.   
8.     DEBUG_PRINT("  x=%d, y=%d\n", x, y);    // In the function doit(), at line 12107 in ..\pku.cpp  
9.     return z*z;                             //   x=10, y=10  
10. }  
11.   
12. int main() {  
13.     int x = 10;  
14.     int y = 10;  
15.   
16.     int z = doit(x, y);  
17.     cout << z << endl;  
18. }  

#define DEBUG_PRINT(fmt, ...) \
    printf("In the function %s(), at line %d in %s\n", __func__, __LINE__, __FILE__); \
    printf(fmt, __VA_ARGS__)

int doit(int x, int y) {
    int z = x + y;

    DEBUG_PRINT("  x=%d, y=%d\n", x, y);    // In the function doit(), at line 12107 in ..\pku.cpp
    return z*z;                             //   x=10, y=10
}

int main() {
    int x = 10;
    int y = 10;

    int z = doit(x, y);
    cout << z << endl;
}

3つ目の例では、ログ出力のオン、オフを切り替えれるようにします。マクロの内容そのものの定義を、DEBUG_MODEが定義されているか否かで切り替えています。この書き方はかっこいいですねー。しびれます。

view plainprint?

1. #define DEBUG_MODE  
2.   
3. #ifdef DEBUG_MODE      
4. #define DEBUG_PRINT(fmt, ...) \  
5.     printf("In the function %s(), at line %d in %s\n", __func__, __LINE__, __FILE__); \  
6.     printf(fmt, __VA_ARGS__)  
7. #else  
8. #define DEBUG_PRINT(fmt, ...)  
9. #endif  
10.   
11. int doit(int x, int y) {  
12.     int z = x + y;  
13.   
14.     DEBUG_PRINT("  x=%d, y=%d\n", x, y);    // In the function doit(), at line 12107 in ..\pku.cpp  
15.     return z*z;                             //   x=10, y=10  
16. }  
17.   
18. int main() {  
19.     int x = 10;  
20.     int y = 10;  
21.   
22.     int z = doit(x, y);  
23.     cout << z << endl;  
24. }  

#define DEBUG_MODE

#ifdef DEBUG_MODE    
#define DEBUG_PRINT(fmt, ...) \
    printf("In the function %s(), at line %d in %s\n", __func__, __LINE__, __FILE__); \
    printf(fmt, __VA_ARGS__)
#else
#define DEBUG_PRINT(fmt, ...)
#endif

int doit(int x, int y) {
    int z = x + y;

    DEBUG_PRINT("  x=%d, y=%d\n", x, y);    // In the function doit(), at line 12107 in ..\pku.cpp
    return z*z;                             //   x=10, y=10
}

int main() {
    int x = 10;
    int y = 10;

    int z = doit(x, y);
    cout << z << endl;
}

本当にユーザーが欲しいものは何か？

youtubeで興味深いビデオを見つけた。

簡単に要約すると、こんかんじ。

「アイディアはいろんなところからやってくる。ユーザー、技術者、役員、分析から。

そして、アイディアが出たら、次はプロトタイプだ。きちんとアイディアを実証できるものがつくれるか？どうやってユーザーのイメージを正確にとらえるか。

イノベーションをもたらすのに一番大切なのは、イタレーションである。何かを試す、フィードバックをもらう、また新しい何かを試す。これの繰り返し。そして、ユーザーが本当に求めているものに近づいていく。

googleが成功している理由は、いいアイディアを持っているから、いいアイディアをうまく実行しているから、そしてとても小さなチームで働いているからだ。

その小さなチームの中で、技術者たちは、決断をする。何がこのサービスの一番の売りなのか？

本当にユーザーが必要としているのは何か？最高の製品をどのようにつくりあげるのか？と。」

• ユーザーが本当に欲しいものをきちんと考える。
• 小さいチームが、決断権を持って仕事を進めていく

というのは非常に考えさせられるものがあった。

以前いたプロジェクトで以下のような話を聞いた。

「お客さんがどうしても作ってほしいっていうシステムがあった。難しい案件だったけど、なんとかがんばって作った。そのシステムの構築には、数千万円かかったらしい。でも、そのシステムができて、1年くらい経つけど、実際に使われたのって2回しかないらしいんだよね」

先輩の話によると、お客さん側のシステム部とユーザーの間で認識が違っていたらしく、システム部や経営陣が作って欲しいと言っていたシステムは、実はユーザーからすると本当は必要なシステムではなかったらしい。

おそらくシステム屋さんであれば、よく聞く話だと思う。

• 実際にものを作る人
• 実際にものを使う人
• 決定を下す人

これらの間で、綿密なやりとりが無ければ、作ってみて実は要りませんでした～ってことに成りかねない。そして、これらの3種類の人間が一堂に会する機会というのはほとんどないというのが現実だと思う。

例えば、実際にものをつくる私たちシステム屋が、お客さんのユーザーと話す機会や、お客さんの意思決定権を持つひとと話すことは、ほぼない。お客さん側のシステム部という部署の人たちから要件をいただく。

細かいやりとりを抜かせば、コミュニケーションの大きな流れは以下の２つだと思う。

• ユーザー -> (開発依頼) -> システム部 -> (意思決定依頼) -> 経営陣 -> (意思決定) -> システム部 -> (開発依頼)->開発者
• 開発者 -> (要件、仕様確認) -> システム部 -> (確認) -> ユーザー -> (返答) -> システム部 -> (返答) -> 開発者

上記では、システム構築者を開発者としてまとめたが、実はこのシステム構築者も一次受け、二次受け、三次受け、・・・、オフショアなどと細かく階層化されている。

問題点は、２つある。１つは、情報の伝達が遅くなってしまうこと。２つ目は、ユーザーの真意が開発者に直接伝わらないことだ。大きい企業だと、このように階層化された情報伝達は仕方ないのだろうか？

私は、学生時代にアルバイトで社内システム（住所録管理、在庫管理、部品発注システムなど）を作ったことがある。

そのときは、「窓がぱっと出てきて、ボタンぽんって押したら、シュッてこの項目が別画面に飛んでボタン押したら、こうなって、・・・みたいなの作って欲しんだけど？」と言われ、「はい、分かりました。」というかんじで要件が確定し、早ければその日のうちに試作バージョンを作って、「こんな感じだけどどうでしょう？」と聞いて、「うーん、ここはもうちょっとこうなってた方がいいなー」と言われ、また作りなおして、・・

という感じで進めていた。

このやり方だとスピード感が全然違う。開発者である私と、ユーザーと、そして意思決定権のある社長。３人が画面を見て、どうするのかをその場で決める。そして、要件定義書や設計書というイメージしにくい文書で話をするのではなくて、実際にものを作って、それを改良していくので、ユーザー、開発者間の認識の差異が生まれにくい。

”アジャイル開発”という言葉が先進的な開発スタイルみたいな感じで語られることが多いが、私はこれが本来あるべき普通のシステム開発の姿のように思う。

グラフ理論の紛らわしい言葉

歩道とか経路とか小道とか紛らわしい用語の整理。こういう学術的な専門用語は日本語で覚えるより英語で覚えた方がいいと思っているので（英語で覚えると世界中のエンジニアと会話ができるから）、英語で書きます。

まず、walkというのが大きな集合としてあります。walkとは、接点と枝の列です。そしてwalkには、特に重要な以下の4つを形があります。

• trail
• path
• circuit
• cycle

まず、始点と終点が同一か否かで、trail, pathとcircuit、cycleに分かれます。問題は、trailとpathの違い、circuitとcycleの違いです。

まず、trailとpathの違いについて。trailは、含まれる枝がdistinctなものです。pathは含まれる接点がdistinctなものです。

上の図だと、1->2->3->4->2はtrailです。しかし、2の節点を2度通っているため、pathではありません。1->2->3はpathです。

　ここで思い浮かぶのは、shortest pathという言葉です。あれは、pathですよね？ダイクストラをやるときに、一度通った点はもう通らないというような処理をしているかと思います。また、Euler trailという言葉があります。あれは、同じ接点を複数回通るのでpathじゃなくて、trailなんです。同様に、Hamiltonian pathとう言葉があります。あれは、trailじゃなくて、pathですよね。

次に、始点と終点が同一になる、circuitとcycleについて。この2つも同様の基準でカテゴライズされています。circuitは含まれている枝がdistinctであるもの、cycleは含まれている節点がdistinctであるものを指します。ループが一つだけのもの（始点と終点だけ）がcycle、他にもループを持っているものをcircuitと呼ぶと考えればいいです。

以上から、集合論上は、以下の関係が成り立つかと思います。

path ⊂ trail ⊂ walk

cycle ⊂ circuit ⊂ walk