フェルマーの小定理

LayCurseさんの日記を見て面白そうなのでチャレンジした問題。

http://www.codeforces.com/contest/57/problem/C

①重複組み合わせ（Homogeneous Combination）と

②包除原理（Inclusion-exclusion principle）

の練習にと解いてみたが、見事にhogeりました。。

LayCurseさんのソースと睨めっこすること約2日、どうやら素数pと自然数aに対して

a^p = a (mod p)

が成り立つのではないか？？と気付きました。

これは「フェルマーの小定理」と呼ばれているようです。

別の書き方をすれば、

a^(p-1) = 1 (mod p)や

a^(p-2) = 1/a (mod p)

です。

上の問題は、この定理を用いて重複組み合わせを計算します。

普通、組み合わせはパスカルの三角形を利用してDPで出せますが、この問題の場合は抽出対象の個数が大きいためメモリーオーバーとなり、NGです。

フェルマーの小定理を利用すると、法の世界の割り算ができます。

普通に

nCk = n! / (n-k)!k!

を計算しましょう。

pが大きいので、累乗の計算は繰り返し二乗法で行いましょう。

以下ソースです。

const int MOD = 1000000007;
long long int f[200001];
long long int rf[200001];

long long int pw(long long int x, long long int y) {
    long long ret = 1LL;

    REP(mask, 50) {
        if (y >> mask & 1)
            ret = ret * x % MOD;
        x = x * x % MOD;
    }
    return ret % MOD;
}

long long int homoComb(int x, int y) {
    int n = x + y - 1;
    int k = y;

    rf[0] = f[0] = 1;
    FOR (i, 1, n + 1) {
        f[i] = f[i-1] * i % MOD;
        rf[i] = pw(f[i], MOD-2);
    }

    long long int ret = f[n];
    ret = ret * rf[k] % MOD;
    ret = ret * rf[n-k] % MOD;

    return ret;
}

int main() {
    int n;

    cin >> n;
    memset(f, 0x00, sizeof(f));
    memset(rf, 0x00, sizeof(rf));

    long long int ret = homoComb(n, n) % MOD;
    ret = 2 * ret % MOD;
    ret -= n;
    if (ret < 0) ret += MOD;

    cout << ret << endl;

    return 0;
}

Google+Facebook=??

Facebookアプリを作成しよう！と思いついたが、どうやら自分でサーバーを立てないとダメみたい。

ホスティングサービスしてくれ。。と思いつつも、

よく考えると、、、

Google App Engineがあるじゃん！！

これに気付いた自分はかなりの天才では！？

と思いましたが、すでにやってる人結構いるみたいですね。。

実は、GAEは昔ちょっとやってました。

そのときはJavaでやってたのですが、FacebookのAPIが簡単に使えるPythonに乗り換えてチャレンジしました。

PythonのGAEは、簡単すぎます。。。

必要な設定ファイルは、一つだけ。

使用するコマンドは、

• テスト用にローカルサーバーをあげて画面確認する
• GAEにソースをアップロードする

の2つだけ。

これは、楽すぎるのでは。

あとは、FacebookでアプリケーションIDを取得し、自作のGAEサイトをコールバックするように設定するだけ。

今日作ったのはここまで（少っ。。笑）

http://apps.facebook.com/kenjih-test-eleven/

VMwareでUbuntu

最近VMwareをインストールして、Ubuntuを試している。

VMwareとはOSのエミュレータ。windows上でunixやMac OSを動かすことができます。

cygwinでもいいですが、いろいろなOSを試せるというのがVMwareの利点。

下から、VMware playerをインストールできます。

もちろん無料。

http://www.vmware.com/jp/products/player/

そして、VMware上で、OSをエミュレートするにはゲストOSと呼ばれる"動かしたいOS"をダウンロードしなければいけません。

今回は、ubuntuを入れてみました。

下のサイトに、VMware用のubuntuがありましたのでそれを落としました。

http://www.ubuntulinux.jp/products/JA-Localized/vmware

ダウンロードしたファイルを解凍すると、「Ubuntu.vmx」というファイルが出来ます。このファイルをVMwareにドラッグすると、Ubuntu OSが立ちあがります。

ちょー簡単。。

ちなみに外観はこんな感じになります。

いやー便利ですねー。

僕のようなunix初心者にはありがたい限りです。

How simple it could be!?

Lately I strongly think that things should be as simple as possible.

The simpler, the better.

It's true for computer programming.

If the source code is simple, there's fewer bug contained.

And to help you make things simple, there are some great stuff called "STL."

See the problem below:

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

I solved similar problems several times. Guess it's got to be a quite famous problem, and I solved it correctly. Then again, I was not sure what approach was appropriate to this problem.

I used to solve it with some queues.

But it hit me today.

It can be solved with priority_queue in much easier way. Once I noticed it, things went simple, neat and tidy!

Pretty simple, huh??

int main() {
    int p, q, r, k;

    cin >> p >> q >> r >> k;
    priority_queue<LL>que;

    que.push(-1);

    int cnt = 0;
    while (!que.empty()) {
        LL n = que.top();
        while (!que.empty() && que.top() == n)
            que.pop();

        if (cnt++ == k) {
            cout << -n << endl;
            break;
        }

        que.push(n*p);
        que.push(n*q);
        que.push(n*r);
    }

    return 0;
}

Karp-Rabin Algorithm

今日はKarp-Rabinの文字列アルゴリズム（通称KR法）について。

Karp-Rabin法は、大量のテキストデータから対象の文字列を検索するアルゴリズムである。

1,000,000文字のテキストから、1,000文字のキーワードを抽出したいシチュエーションを考えよう。

以下のbrute-forceな解法では、1,000,000 * 1,000 = 10^9の計算時間がかかってしまう。

実際に下のコードを実行すると、10秒くらいかかった。

string text;
string target;

bool doit1() {
  int L = text.length();
  int l = target.length();

  REP(i, L-l+1) {
      bool ck = true;
      REP(j, l) {
          if (text[i+j] != target[j]) {
              ck = false;
              break;
          }
      }
      if (ck) return true;
  }
  return false;
}

int main() {
  // case 1
  text.assign(1000000, 'a');
  target.assign(999, 'a');
  target.append("b");

  cout << doit1() << endl;
}

Karp-Rabin法では、比較文字列をハッシュ値に変更することで、文字列の比較をO(1)で実施するようにしている。

ここでポイントとなるのは、ハッシュ値の算出方法である。通常、長さlの文字列のハッシュ値を計算するためにはO(l)必要だが、「ローリングハッシュ」とよばれる窓関数のような特性をもったハッシュ関数を用いることでハッシュ値の計算をO(1)で実施することができる。

まずは、簡単なハッシュ関数として文字列内の文字のアスキーコードをすべて足す関数を用いてみる。

bool doit2() {
  int L = text.length();
  int l = target.length();

  int textH = 0, tarH = 0;
  REP(i, l) {
      textH += text[i];
      tarH += target[i];
  }

  REP(i, L-l+1) {
      if (textH == tarH) {
          bool ck = true;
          REP(j, l) {
              if (text[i+j] != target[j]) {
                  ck = false;
                  break;
              }
          }
          if (ck) return true;
      }
      textH = textH - text[i] + text[i+l];
  }
  return false;
}

これで、先ほどの例題は、1秒以下で解くことができる。

Karp-Rabin法の弱点は、ハッシュ値の衝突が起きた際に文字列の比較をしなければいけないという点にある。つまり、最悪の場合（ほとんどの場合衝突が起こる場合）の計算量は、brute-forceの場合と変わらない。

例えば、このようなインプットに対しては、先ほどのハッシュ値ではほぼ毎回衝突が起きてしまう。計算時間も約10秒かかってしまった。

int main() {
  // case 2
  text.assign(1000000, 'b');
  target.assign(998, 'b');
  target.append("ac");

  cout << doit2() << endl;
}

Karp-Rabin法では、如何にローリングハッシュ関数を決めるかが重要である。例えば、アスキーコードを掛け合わせてハッシュ値を作るというのはいいかもしれない。但し、個の場合、値が非常に大きくなるため適当な素数のMODをハッシュ値として用いなければならない。（また、次の部分文字列のハッシュ値を計算する際、割り算が必要だが割り算は合同式の分配法則は成り立たないので工夫が必要である。）以下にサンプルソースを記載する。この場合は、上記のサンプルでも1秒程度で計算できた。

bool doit3() {
   const int MOD = 1000007;
   int L = text.length();
   int l = target.length();

   int textH = 1, tarH = 1;
   REP (i, l) {
       textH = textH*text[i]%MOD;
       tarH = tarH*target[i]%MOD;
   }

   REP(i, L-l+1) {
       if (textH == tarH) {
           bool ck = true;
           REP(j, l) {
               if (text[i+j] != target[j]) {
                   ck = false;
                   break;
               }
           }
           if (ck) return true;
       }
       while (textH % text[i])
           textH += MOD;
       textH = (textH/text[i]*text[i+l]) % MOD;
   }
   return false;
}

※合同式の割り算の部分のwhile()が何回くらい実行されるのか気になるが、最大でもtext[i]回である。ほぼ定数オーダーと考えていいと思う。

CDTでデバッグ

最近仕事でデバッガの重要性を認識し、家でプログラミングするときもデバッガが必要なのではと思い始めた。

仕事のプログラミング（特に保守プログラミング）では、デバッグがないと仕事にならない。

なぜなら、人が書いたコードは分からないから。

家でやるプログラミングは、デバッガ不要だと思っていた。

しかし、よく考えてみると、難しい問題を解く場合はcoutやcerrなどで途中結果を表示させて確認することが多い。（特に文字列処理の問題では、よくやってる気が・・・）

これは、デバッガでやるべきだろう。。

で、試してみました。

開発環境は、eclipse+CDTなので、gdbを入れてeclipse上で使えるようにした。

これは簡単。

F11　デバッガ開始

F8 次のブレークポイントまで実行

F5　ステップイン（次の行に進む）

F6 ステップオーバー

F7 ステップリターン

という感じで使えます。

気になった点は、インスペクタがちょっと独特です。

「変数」ウィンドウでスタック変数の内容を確認できますが、なぜがグローバル変数は参照できません。

「式」ウィンドウというのがあるので、そこに監視する式を入力すれば、その値を参照できます。

グローバル変数を参照したいときは、この「式」ウィンドウを使うと良いみたい。

あーー、やっぱeclipse便利やな～～。

SRM BrushUp: ReversalChain

We -- A friend of mine and I try the past SRMs every Sunday-- have started challenging division 1. It seems 500-pointer problems in div.1 is almost as same level as 1000-pointer problems in div.2. And I don't think they're far beyond our abilities.

Today I solved a 500 problem we were not able to solve the last weekend.

Once I read the editorial, there's nothing difficult. All you have to do is recursive and memorize return values not to exceed the time limit.

The problem is like this:

You are given two strings, init and goal.

You can choose a pair of indices{i, j}, and reverse the string init within the range of {i, j}.

Starting from init, how many times you have to do the operation above to get the string goal. Return the minimum possible value. If it's impossible return -1.

The approach to the problem is quite simple. Just check the head and the tail of each string.

If none of them have their equivalent, you cannot make the objective string.

If some have their equivalent, then move to the next recursion.

My source is like this: (It takes lots from ACRush's source code.)

class ReversalChain {
public:
    int solve(string init, string goal) {
        if (init.length() == 1) {
            if (init == goal) return 0;
            else return INF;
        }
        if (memo.count(init+"-"+goal))
            return memo[init+"-"+goal];

        int ret = INF;
        char i0 = init[0], in = init[init.length()-1];
        char g0 = goal[0], gn = goal[goal.length()-1];

        if (i0 == g0)
            ret = min(ret, solve(init.substr(1), goal.substr(1)));
        if (in == gn)
            ret = min(ret, solve(init.substr(0, init.length()-1),
                    goal.substr(0, goal.length()-1)));

        string rev = init;
        reverse(ALL(rev));
        if (i0 == gn)
            ret = min(ret, 1+solve(rev.substr(0, rev.length()-1),
                    goal.substr(0, goal.length()-1)));
        if (in == g0)
            ret = min(ret, 1+solve(rev.substr(1), goal.substr(1)));

        return memo[init+"-"+goal] = ret;
    }

    int minReversal(string init, string goal) {
        memo.clear();
        int ret = solve(init, goal);
        return ret==INF ? -1 : ret;
    }
private:
    map<string, int>memo;
};

ポインタの参照渡し

今日は、ポインタとアドレスの違いについて。

まず、両者の違いについて。

ポインタ（以下では、アドレスとの違いを明確にするため”ポインタ変数”と呼ぶ）はあくまでも、アドレスを指し示す入れ物である。

ポインタ変数自体はアドレスではなく、intやcharと同じ変数である。

これに対してアドレスは、固定番地である。メモリの物理的な位置である。

ちょっと突っ込んで考えると、ポインタ変数にもアドレスが存在するし（アドレスを指し示す変数のアドレスというイメージ）。

よってポインタ変数を関数に参照渡ししたい場合は、ポインタ変数のアドレスを渡す必要がある。

具体的な例を見てみよう。

main()でchar型のポインタ変数を宣言し、setStr()という関数にポインタ変数を渡して、文字列を取得する。

まずは、ダメな例から。

void setStr(char *a) {
   a = new char [8];
   strcpy(a, "Hello");
}

int main() {
   char *a = NULL;

   setStr(a);
   cout << a << endl;
   return 0;
}

上の書き方では、setStr()を抜けた後もaはNULLのままです。

setStr()に渡したのは、ポインタ変数です。

setStr()では、仮引数aのコピーを作り、関数内だけで有効なaに対して処理をします。よって、setStr()をreturnするときには、aに対する変更はmainに反映されません。

次に、良い例。

下の書き方では、ポインタ変数のアドレスを渡しています。

void setStr(char **a) {
    *a = new char [8];
    strcpy(*a, "Hello");
}

int main() {
    char *a = NULL;

    setStr(&a);
    cout << a << endl;
    return 0;
}

これだと、期待した結果が得られます。

くどいようですが、動的メモリの確保をする場合は「ポインタ変数とアドレスは異なるものだ」としっかり理解しておくことが大切です。

アルゴリズマー V.S. システムエンジニア

アルゴリズムコンテストに参加するようなアルゴリズマーがIT関連の職場に出たとして、

果たして彼らは活躍できるのか？

一般の人たちからすれば同じプログラミングだろうが、実際に要求されるスキルはまったく異なる。正直に言うと、現場の職場ではあまり高いレベルのプログラミングを書くことはない（SIer、IT consulting firmなどの場合）。

やることといえば、

①画面からの入力を受けて、

②その入力に応じてDBにアクセスしたり、

③計算したりして、

④画面に出力を出す。

簡単に言うとこれだけ。③はせいぜい掛け算、足し算くらい。素数を求めたり、因数分解したり、バックトラッキングしたりとかは、普通の現場ではまず無いと言っていいだろう。

しかし、現場のプログラミングが簡単かと言われると決してそうではない。技術力以外に必要な要素が非常に多い。

今日は、自分なりに違いをまとめてみました。興味があればどうぞ。

	アルゴリズムコンテスト	現場プログラミング
ひらめき	必要	あまりいらない
数学的思考能力	必要	あまりいらない
論理的思考能力	必要	少しだけ必要
体力	ちょっとだけ必要	必要不可欠
精神力	たまに必要	絶対必要
デバッガ	デバッグなどしませんが。。	無いと無理
妥協	自分の裁量で	多いに必要
コミュニケーション能力	あまりいらない	必要
チームワーク	場合によりけり	絶対必要
英語力	必須	いらない
ショートコーディング	大事	やったら怒られる
環境構築能力	ほぼいらない	大事

と見てみると、アルゴリズマーには知力が、システムエンジニアには人間力が必要なようだ。

どっちも備わった人間は神ですね。。

One of the Best Hackers in the World!

Facebook HackerCupのsemifinalが迫ります。

finalには行けなくてもT-shirtが欲しいです。

ロゴはこんな感じ。

"You are one of the best hackers in the world.

certified by Facebook"

やばい。。

今日は1Aの問題を解いてみた。「Wine Tasting」という問題。

g本ワインを試飲して、c本以上正しく正解する場合のパターン数を求めよ。という問題。

Ai = {1,2,3,4, ... g}と定義すれば、Ai == iとなる要素数がc以上になる場合の数列のパターン数を求めよ。

という問題に帰着できる。

SjをAi==iとなる要素数がj個である場合のパターン数と定義すると、求める値が

Sum_{c <= j <= g} Sj

となるのは自明。

あとは、Sjの求め方を考えればよい。

まず、Ai==iとなる場所を選ぶパターンは、gCjである。

次に、残ったg-j個をAi == iとならないように並べる。

このパターン数をptn(j)とおくと、

ptn(j) = (j-1) * (ptn(j-1) + ptn(j-2))という漸化式が成り立つ。

これは、難しいので解説。（これに気付くのに2日かかった。。。この問題解いた人は閃いたのか？それとも知っていたのか？？）

j-1個のAi==iとならないパターン{A1, A2, ..., Aj-1}の後ろにAjをおく。

Ajと任意のAk, 1<=k<=j-1を交換してもこの数列は、Ai == i, 1<=i<=jとならない。

よって、(j-1) * ptn(j-1)のパターンを数え上げた。

次に、Ajとスワップする要素に関しては、最初の段階でAi==iとなっていても問題ない（スワップするため）。Ai==iとする場所をj-1個の要素から選び、残りのj-2個を並び変えるパターンは、

(j-1) * ptn(j-2)

よって、前述を2つを足して、

ptn(j) = (j-1) * (ptn(j-1) + ptn(j-2))

となるわけである。

むずっ。。

ここまで出来れば、実装は簡単。

コンビネーションは「パスカルの三角形」で計算すればOK。

残ったものの並び変えも、DPで計算すればOK。

あとは、MODULOの比較的大きいため、64bitに入れないと掛け算でオーバーフローするので注意。

ソースはこちら。（今回は何回もチャレンジして書きなおしたため、最終的なソースは結構きれい。。）

#define MOD 1051962371
uint64_t cmb[126][126];
uint64_t ptn[126];

int main() {
   cmb[0][0] = 1;
   FOR (i, 1, 126) cmb[i][0] = 1;
   FOR (i, 1, 126)FOR (j, 1, i+1)
       cmb[i][j] = (cmb[i-1][j-1] + cmb[i-1][j]) % MOD;

   ptn[0] = 1, ptn[1] = 0;
   FOR (i, 2, 126) {
       ptn[i] = (i-1) * (ptn[i-1] + ptn[i-2]);
       ptn[i] %= MOD;
   }

   int n, g, c;
   cin >> n;
   while (n--) {
       cin >> g >> c;
       uint64_t ret = 0;
       FOR(i, c, g+1) {
           int t = cmb[g][i] * ptn[g-i] % MOD;
           ret = (ret + t) % MOD;
       }
       cout << ret << endl;
   }

   return 0;
}