ポワソン分布

今日、ポアソン分布という分布を知った。

高校の物理でやったような記憶もあるが・・・(ポアソン方程式だった。。しかもやったのは大学。)

まず、以下の問題を考えてみてください。

ある女の人は一日平均2回ナンパされるそうです。（とっても美人です。）

この女の人が、今日3回ナンパされる確率を求めなさい。

この問題、簡単そうで普通には解けません。（少なくとも私には解けませんでした。）

一日にx回ナンパされる回数をP(x)とおくと、

とりあえず、ここまでは分かる。

で、P(3)は？？・・・・・・

そこで登場するのがポアソン分布。ポアソン分布が使えるのは以下の条件（ポアソン過程）を満たす場合のみらしい。

1. 事象はいかなる時点でもランダムに発生しうる。
2. 与えられた時間区間での事象の発生は，それと重複しない他の区間に対して独立である。
3. 微小時間⊿tにおける事象の発生確率は⊿tに比例して小さくなっている。
4. 微小時間⊿tの間に事象が２回以上発生する確率は無視できる。
5. 時間tの間に当該事象が発生する平均発生回数λがおおむね5以下である。

でこれがポアソン分布の式。

数学的な裏付けは分からないが、この式を使用すれば単位時間あたり平均λ回発生する事象が単位時間内にx回発生する確率を求めることができる。

この式考えた人天才すぎます。。。

先ほどのナンパの問題だと、

となります。すごい、すごすぎる、この女。じゃなくてMr. ポアソン。

いろいろなλの値に対するポアソン分布のグラフがwikipedhiaにありましたので、下記の参考サイトからご覧ください。

引用サイト：

http://next1.cc.it-hiroshima.ac.jp/MULTIMEDIA/prob/node14.html

参考サイト：

http://ja.wikipedia.org/wiki/%E3%83%9D%E3%82%A2%E3%82%BD%E3%83%B3%E5%88%86%E5%B8%83

TopCoder惨敗とPKU祝杯！

昨日開催されたTopCoder SRM 483（division 2）.

Easyに手間取り8分くらい使ってしまった。良く見ると、途中で変数の名前が入れ替わってた（笑）焦りは禁物ですね。

Middleは、20分かからないくらいで解けた。これは調子がいい。

初めての本格的なHardの挑戦。なんか行ける気がする。。結局行けなかった・・・・。

そして、チャレンジフェーズ。まさかのMiddle撃沈。ショックを受けている間に、次々と撃沈されるソース達。いったい何が・・

最終的には、部屋の中の上2人以外はMiddleが撃沈されている始末！これは前代未聞の出来事では。。

撃沈されていない人のソースをみると、、、、なるほど1人のときは行と列のうち片方見るだけでいいのか。。確かに。。次回からは、例外っぽいケースはきちんと確認するようにしよう。

そして、PKU。

今日、ついに念願の100問、突破～～～。。

意外と早く達成できて驚いている。次の目標は150問。しかし、そろそろ簡単な問題が無くなってきている。ここからの問題はなかなかの猛者達ばかりだろう。

ワンピースで言うと『グランドライン』到達くらい。。

がんばろう！！targetは東北大学のmatsuさん。。すぐに抜かす。

エラトステネスの脅威

今日は、素数を列挙するアルゴリズムを考えます。

最初は基本的な実装を行い、徐々に高速化していきます。最終的には『エラトステネスの篩』を実装します。

また、それぞれのアルゴリズムのオーダーを見積り、実際に実行速度を計測します。

まずは、基本的な実装から。これ。

void getPrime1(int n) {
    VI prime;

    FOR (i, 2, n) {
        bool isPrime = true;

        FOR (j, 2, i) {
            if (i % j == 0) {
                isPrime = false;
                break;
            }
        }
        if (isPrime)
            prime.PB(i);
    }

    printf("%d\n", prime.size());
}

まー、だれでも思いつくような方法。

少し考えてみると、整数iが素数か否かを判定するのにi-1までの割り算をしなくても、sqrt(i)までの割り算をすればよいことに気付く。これをふまえて、改良するとこうなる。

void getPrime2(int n) {
    VI prime;

    FOR (i, 2, n) {
        bool isPrime = true;

        FOR (j, 2, (int)sqrt(i) + 1) {
            if (i % j == 0) {
                isPrime = false;
                break;
            }
        }
        if (isPrime)
            prime.PB(i);
    }

    printf("%d\n", prime.size());
}

ちょっと変わっただけですが、かなり速くなりそうです。

次に、「エラトステネスの篩」を用いた実装。（アルゴリズムの考え方については、参考ページ参照のこと。）ソースを見てもらえば何をしているか分かると思います。

void getPrime3(int n) {
    bool *isPrime = new bool[n];
    VI prime;

    memset(isPrime, 1, sizeof(bool) * n);
    for (int i = 2; i < sqrt(n) + 1; i++) {
        if (!isPrime[i])
            continue;

        for (int j = i + i; j < n; j += i)
            isPrime[j] = 0;
    }
    FOR(i, 2, n)
        if (isPrime[i])
            prime.PB(i);

    delete [] isPrime;
    printf("%d\n", prime.size());
}

ここまで、3つの方法を記述しました。それぞれの計算量を見積ってみましょう。そして実際にn=1,000,000のときの実行速度を計測してみましょう。

#include <sys/time.h>

inline double gettimeofday_sec() {
    struct timeval t;
    gettimeofday(&t, NULL);
    return (double) t.tv_sec + (double) t.tv_usec * 1e-6;
}

int main() {
    double s = gettimeofday_sec();
    getPrime1(1000000);
    printf("%lf\n", gettimeofday_sec() - s);

    s = gettimeofday_sec();
    getPrime2(1000000);
    printf("%lf\n", gettimeofday_sec() - s);

    s = gettimeofday_sec();
    getPrime3(1000000);
    printf("%lf\n", gettimeofday_sec() - s);

    return 0;
}

結果は以下のとおりになりました。

アルゴリズム	オーダー	実行時間[s]
getPrime1	n^2	157.9
getPrime2	n^1.5	6.943
getPrime3	n	0.016

やっぱりエラトステネス速いですね。。getPrime1, getPrime2のオーダーは見てのとおりですが、getPrime3（エラトステネス）のオーダーは少し分かり難いので少し補足します。

まず、一番外側のループは見てのとおりO(n^0.5)です。次に内側のループですが、n/i回まわります。

iによって値が違うので平均値を考えます。i は 2 から n^0.5 までなので、内側のループの回数は、それぞれのiに対して、n/2, n/3, n/4, ...... n/n^0.5となります。この平均は、およそn/{(2+n^0.5)/2}なので内側のループの回数はO(n^0.5)です。よって2重ループ部分のオーダーは、O(n^0.5) * O (n^0.5) = O(n)です。最後にisPrimeを線形探索する部分もO(n)なので、最終的にエラトステネスのオーダーは、O(n)となります。

参考ページ:

エラトステネスの篩について　http://ja.wikipedia.org/wiki/%E3%82%A8%E3%83%A9%E3%83%88%E3%82%B9%E3%83%86%E3%83%8D%E3%82%B9%E3%81%AE%E7%AF%A9

多倍長演算 --doubleの精度!?

今日は、ついに多倍長演算デビュー。

PKU4問解いてやった。C++で自力で実装しようと思ってたが結局Javaでやることにした。信頼性の高い既成品があれば、当然そちらを使うべき。（あー、いつの間にかビジネスチックな思考になってしまった。。）

Javaは2種類の多倍長演算に対応する型を装備している。

まずは、BigInteger。その名の通り大きな整数。これはC++でいうlong long intの最大値より大きな値を扱いたい場合に使用する。long long intの最大値は、(2^64) / 2 - 1でだいたい8.0 * 10^18。

話はちょっとそれるが、『2^10 ≒ 10^3』と覚えておくと累乗の計算はすぐできるので便利。

2つ目は、BigDecimal。これはdoubleより精度の高い小数を扱うときに用いる。

BigIntegerの概念は分かりやすいが、BigDecimalは・・・。はて？？doubleの精度ってそもそもどのくらいだったっけ？

ということで、今日は、

1. doubleの最大値、最小値
2. doubleの精度
3. 丸め誤差

について、実際にPC上でプログラムを走らせながら検証します。（ちょっと長くなりそうだな～）

その前に、doubleの内部構造についておさらいします。doubleは1ビットの符号部、11ビットの指数部、52ビットの仮数部から成り立っています。

（詳しくは、参考ページをご参照ください。）

簡単に言うと10進数の指数表記（2.999 * 10e+10みたいなやつ）と同様のことを2進数でやってるだけです。

では、準備もできたので、まずは、doubleの最大値、最小値を求めてみましょう。指数部が11ビットなので概算では大体2^1024くらいがdoubleで表すことのできる最大値です。だいたい10^300くらいですねー。最小値は、2^(-1024)なので、10^(-300)程度だと思います。

では実際に計算しましょう。

#include <float.h>

#define EXP_MAX 2046     // 2047 is booked for exception procedure
#define EXP_MIN 1        // 0 is booked for exception procedure
#define EXP_BIAS 1023    // The biased exponent range: [-1022, 1023]
#define SIG_BIT 52       // significand bits

int main() {
   // max
   double dblmax = pow(2.0, EXP_MAX - EXP_BIAS);
   double significand = 1.0;
   REP(i, SIG_BIT)
       significand += pow(2.0, -(i+1));
   dblmax *= significand;
   printf("dblmax = %e, DBL_MAX = %e\n", dblmax, DBL_MAX);

   // min
   double dblmin = pow(2.0, EXP_MIN - EXP_BIAS);
   printf("dblmin = %e, DBL_MIN = %e\n", dblmin, DBL_MIN);

   return 0;
}

んー、実際は、最大値が1.797693e+308、最小値が2.225074e-308でした。ほぼほぼ予想通りです。指数部の値は、0と1024が例外処理のため予約語となっているようです。（参考ページ参照のこと。）やられた。これに気付かず1時間ぐらい悩んでしまった。。。

次に、doubleの精度を考えてみましょう。指数部が数値のスケール、仮数部が数値の精度を表しているので仮数部に注目すればOKです。52ビットなので、およそ10^15程度（有効数字15桁）の精度は持っているのではないでしょうか。（経験的にはそんなに無い気が・・・）実際にPC上で確認してみましょう。

int main() {
   double min_significand = pow(2.0, -SIG_BIT);
   printf("%.16f\n", min_significand);

   double x = 1e15 + 1.0;
   double y = 1e16 + 1.0;
   printf("%.1lf\n", x);
   printf("%.1lf\n", y);
}

んーと、仮数部の最小値は2.0 * 1e-16程度でした。スケールに対して2.0 * e-16程度の大きさの値までしか表現できないという意味なので、有効数字は15桁と考えていいでしょう。

実際に上のソースで検証しましたが、1.0 e+16に1を加算しても無視されてしまいます。

あれ、そうか！これが世間で言われている『情報落ち』ってやつか！！

最後に、「丸め誤差」について。言わずと知れた「10進数を2進数に変換するときにでる誤差」のことです。代表的な例に、0.1を10回足しても1にならないというのがあります。

今回は、0.1をdoubleで表したときにどのような表現になるのか、またその値は10進数ではどの程度なのか調べたいと思います。標準関数frexp()を使えば簡単にできるようですが、なかなか面白そうなテーマなので、今回は自力で実装します。

ポイントは、『仮数部は52ビットあるので全件探索ではNG。何らかの”工夫”が必要』ってことです。まーでもこの場合は工夫と言ってもアレしかないですね（笑）

double getSigVal(LL sigbit) {
   double ret = 1.0;

   REP(i, SIG_BIT) {
       if ((sigbit >> SIG_BIT - 1 - i) & 1)
           ret += pow(2.0, -(i+1));
   }
   return ret;
}

double getDoubleVal(int exp, LL sigbit) {
   return getSigVal(sigbit) * pow(2.0, exp - EXP_BIAS);
}

int main() {
   double x = 0.1;
   LL sig, sigl, sigr;
   int expo = 1;

   while (pow(2.0, expo + 1 - EXP_BIAS) <= x)
       expo++;

   sigl = 0LL;
   sigr = (1LL << 52) - 1LL;
   sig = (sigl+ sigr) / 2;

   REP(i, 200) {
       double val = getDoubleVal(expo, sig);
       if (ABS(val - x) < 1e-18)
           break;
       else if (val > x)
           sigr = sig;
       else
           sigl = sig;
       sig = (sigl + sigr)/2;
   }


   printf("significand = %lld, exponent = %d\n", sig, expo);
   printf("0.1 is expressed something like %.18lf on the PC\n", getDoubleVal(expo, sig));
   return 0;
}

52ビット（10進数で10^15程度）を線形探索していたのでは、時間がありません。Binary Searchで解きましょう。精度は1e-18程度でいいでしょう。ループは200回くらい回せば探索範囲が1e-60倍くらいに縮小するので200回くらい回せばこの場合は十分です。（Binary Searchはバグを生みやすいので、無限ループはやめて必ず有限回の探索で打ち切るようにしましょう。）

参考ページ:

doubleの内部構造

http://www.cc.kyoto-su.ac.jp/~yamada/pB/float.html

指数部の例外処理予約値

http://www.asahi-net.or.jp/~uc3k-ymd/Lesson/Section01/section01_06.html

配列を関数に渡す

C/C++を使用する人なら、値渡しとアドレス渡しがあるのは知っていると思います。

しかし、配列は値渡しをしたつもりでもアドレス渡しになってしまうようです。

また、配列を関数に渡す場合には、どうすればいいか迷うことも多いでしょう。特に多次元配列の場合はそうです。

まず、配列が他の型とは違う挙動をすることを実際にアドレスを表示させることで確認してみましょう。

/*
* int address
*/
void func1(int n) {
   printf("The address of input is %d.\n\n", &n);
}

/*
* buffer address
*/
void func2(int n[]) {
   printf("The address of input is %d.\n\n", n);
}

/*
* STL address
*/
void func3(vector<int> n) {
   printf("The address of input is %d.\n\n", &n);
}

/*
* class address
*/
class Human{
public:
   int age;
   string name;
};
void func4(Human n) {
   printf("The address of input is %d.\n\n", &n);
}

int main() {
   int a;
   printf("The address of input is %d.\n", &a);
   func1(a);

   int b[10];
   printf("The address of input is %d.\n", b);
   func2(b);

   vector<int> c;
   printf("The address of input is %d.\n", &c);
   func3(c);

   Human d;
   printf("The address of input is %d.\n", &d);
   func4(d);


   return 0;
}

配列の場合のみ、呼び出し元と呼び出し先で同じアドレスが出力されます。つまり、配列は値渡ししたつもりでも実際はアドレス渡しになるということです。もちろん関数のプロトタイプ宣言部のint n[]をint *nと書くことも出来ます。こちらの方がアドレスを渡していると明示的に分かるので良いかもしれません。

では、次の問題。多次元配列を関数に渡す時はどうするか。

多次元配列は実際は一次元方向にメモリーが展開されるということを知っておくことが大切です。実際に見てみましょう。

int main() {
  int x[3][3][3];

  REP(i, 3)
      REP(j, 3)
          REP(k, 3)
              printf("%d\n", &x[i][j][k]);

  return 0;
}

どうでしょうか？

メモリが32bit intの環境では4byteずつ増えていると思います。

つまりメモリ上は多次元配列は一次元配列と同じ形をしているということです。

多次元配列を宣言した関数内では、それぞれが何次元で構成されているか分かりますが、関数に渡してしまうと呼び出しさきでは、何次元目がどれくらいの単位で区切られているか分かりません。

よって、関数に渡す際には、呼び出し先の関数のプロトタイプに配列のサイズを明記しなければなりません。

こんな感じ。

void func5(int n[][10]) {
  REP(i, 10)
      REP(j, 10)
          printf("%d%c", n[i][j], (j==9) ? '\n' : ' ');

  return;
}

int main() {
  int buf[10][10];
  memset(buf, 0, sizeof(buf));
  func5(buf);

  return 0;
}

今日の話は基本中の基本って感じですが、やっぱり基本はしっかり押さえておきたいものです。

IFRS対応でITバブル！？

今日は巷を賑す（！？）「IFRS」について書こうと思います。

まず「IFRS」って何ぞや？ということで、簡単に説明します。

International Financial Reporting Standards（ 国際会計基準）の略で、「アイファース」と読まれることが多いです。

読んで字の如く、「企業の財政状態や経営成績を表すための会計基準の世界標準」のことです。

現在、日本の企業は日本独自の会計基準に基づいてレポーティングを実施していますが、近い将来このIFRSが日本の上場企業に強制適用される可能性があります。

2012年を目途に、金融庁によって上場企業へのIFRS強制適用の是非が判断され、強制適用が決定した場合は早ければ2015年よりIFRSに準拠したレポーティングが必須となります。

お試し期間（？）として、2010年3月末より、一定の要件を満たす上場企業の連結財務諸表についてIFRSを任意に適用できるようになりました。既にIFRS取り入れている企業もあるみたいです。。

IFRSの適用されることになれば、決算・財務報告のプロセス及び関連する情報システムの大規模な改革が必要となります。

しかし、会計系システム以外にも影響を受ける部分は多くあります（下表参考）。

影響を受ける業務	変更が必要なシステム
売上計上基準の変更	販売管理システム
棚卸資産の範囲/原価算定方法の変更	購買管理システム/在庫管理システム
固定資産の減価償却の扱いの変更	固定資産管理システム
リース取引の根本的な変更	リース管理システム
未消化有給休暇の費用計上	人事給与システム

いや、これは大変ですね。。(ちょっと難しい経済用語が並んでいるので、その説明については参考サイトを参照下さい。)

上場企業にとっては、お金も時間もかかる一大事ですが、System IntegraterやIT Cosulting Firmにとっては大きなビジネスチャンスとなります。

実際に私の会社でもIFRSは今後もっとも重要視すべきビジネスチャンスとして捉えています。

2010年-2015年にかけてITバブルの再来を期待します。給料あがらないかなー（笑）

引用サイト：

http://www.fsa.go.jp/news/21/sonota/20100423-2.html

http://otona.yomiuri.co.jp/life/kabunushi/kabu100319.htm

http://www.ntts.co.jp/publish/column/ifrs/02/

参考サイト：

棚卸資産とは？http://tax.xrea.jp/auction/tanaoroshi.html

原価とは？http://www.sk-genka.com/kanri/kan02_genka.html

減価償却とは？http://itpro.nikkeibp.co.jp/word/page/10011680/

知ってると便利なSTL(2) pair

intとintのペアを扱いたいとき、ついつい自分でクラスを作ってしまっていた。

何が面倒かというとソート。いちいちoperatorを定義してあげないとsort()が使えない。

最近、pairというSTLが便利だということを知った。

2つの型のペアを格納でき、さらに比較演算子が定義されているため、sort()が使える。かなり便利。

比較演算では、pairの第一要素の比較結果が返される。第一要素が同じ値だった場合は、第二要素の比較結果が返される。

唯一不便そうな、make_pairも下のようなマクロを定義してやれば楽々！！

#define MP make_pair

まー、いつもの如くコードで説明します。誕生日を格納するために、pairを使用しています。そして月日の順でソートをして表示します。

int main() {
vector<pair<int, int> >birthday(20);

REP(i, 20)
 birthday[i] = MP(rand() % 12 + 1, rand() % 30 + 1);

SORT(birthday);

REP(i, 20)
 printf("%02d/%02d is My Birthday.\n", birthday[i].first, birthday[i].second);
return 0;
}

いやー、これはかなり便利だねー。

更に、性と名をpairに格納することも可能。所謂”名前順”ってやつでソートしてみましょう！

int main() {
vector<pair<string, string> >name(10);
name[0] = MP("Tanaka", "Ai");
name[1] = MP("Yamada", "Taro");
name[2] = MP("Yamada", "Hanako");
name[3] = MP("Okada", "Yusuke");
name[4] = MP("Yamaguchi", "Satoshi");
name[5] = MP("Suzuki", "Ichiro");
name[6] = MP("Suzuki", "Jiro");
name[7] = MP("Suzuki", "Saburo");
name[8] = MP("Suzuki", "Shiro");
name[9] = MP("Oda", "Nobuo");

SORT(name);

EACH(itr, name)
  cout << "I'm " << itr->first << " " << itr->second << endl;
return 0;
}

何か気付きましたか？

そうです、鈴木が多いです。

じゃなくて、stringを使っています。const char*ではダメです。

なぜなら、比較演算子が定義されていないからです。（const char*の場合はstrcmp()で比較しますよね！）

TOEIC受験

今日TOEICを受けた。

1年半ぶりだが、そのときよりも良い結果がでそう。。でも正直出来はいまいちだった。。

まず、リスニング。

問題が上から下に設問が書いてあるのに、左から右に設問を読んでしまい、3問ミス笑。。

いやーヘタレすぎる。。

部分的には対策をしたが、きちんと過去問解いておくべきだった。

そして、リーディング。

かなり読むのが速くなってた。でも、最後まで解けなかった。

多分原因は、選択肢をすべて確認していたこと。

選択肢Aが正解だと分かっても、確認のためB, C, Dも吟味していた（笑）

これが無かったら全部解けたかも。。

同じく、過去問を解いておくべきだった。

今回学んだことは、

『何事もリハーサルが必要である。そして、そのリハーサルは可能な限り本番と同じ環境を作り上げて行うこと。』

それから、TOEICのTipsとして以下2点。

1. ちゃんと問題の番号を確認する！！（超基本（笑））
2. これが正解と思ったらすぐにマークする！！

まー、900点取れてなかったら、来年もう一回受けてみようかと思います。

コンマ演算子

世の中には「ショートコーダー」と言われる人達がいる。

コードを短く書くことすべての情熱を注ぐ変態天才たちである。

彼が重宝するテクニックが、「コンマ演算子」。

これを使うと、

1. ループの終了判定が標準入力から入力したものに依存するとき、breakを使わずに簡単に書ける。
2. ループの括弧を省略することができる。（Pythonみたい。）

という利益を享受できる。

すばらしい。。。

百聞は一見に如かずということでソースコードを。

例によって、PKUの問題から。

私の回答はこんな感じ。

int main() {
int L, M, s, e;
while (scanf("%d %d", &L, &M), L++ && M) {
REP(i, M)
    scanf("%d %d", &s, &e),
    L -= e-s+1;
printf("%d\n", L);
}
return 0;
}

注目すべきは、「コンマ演算子」。

まず、whileの中について。ここでコンマを使用することにより、複数の処理を判定条件の中に埋め込むことができます。判定条件として判断されるのは、一番最後の部分（この場合は、L++ && M部分）です。

次に、REP()の中ループについて。ループの内容は2行ですが、括弧を使用していません。

これも「コンマ演算子」によってなせる技。。コンマで繋がれた処理は1行の処理とみなされます。

これは、是非ともマスターしたいテクニック！！！

でもやり過ぎは厳禁。業務でこんなの使いまくってたら殺されますね（笑）

最小公倍数と最大公約数

自数数a, b (a > b)の最小公倍数と最大公約数を求めるプログラムを作成せよ。

んーー。。

最大公約数（greatest common devisor）は言わずと知れたユークリッドの互除法で解けます。

まー一応ソースを載せるなら、

こんな感じ。

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

じゃー、最大公約数（least common multiple）は？

・・・・・

ちょっと迷った（笑）

でもちょっと考えてみると、a * bを出してa と bの共通の約数で割っていくと、lcmが出ることに気付く。

ということは、下のコードでイケますね！

int lcm(int a, int b) {
    return a*b/gcd(a,b);
}

今日は、若干手抜きな感がありますが、悪しからず。

PKU Top10,000入り！！

ついにPKU上位10,000に入りました。

思ったより早く達成できました。

やっぱり、地道に解くことが大切だと痛感しています。

日々の小さな積み重ねが山となるのです。

まあ簡単な問題しか解いてないような気がしますが・・・・

これから更に上に行こうと思うなら、簡単な問題を解くだけでは厳しいでしょう。

次の目標は、上位5000位！！

頑張ります！

最近思うのは、問題を解くだけではダメだということです。

やはり、解くだけだと理解が自分の中で閉じてしまいます。

自分の解いた問題の解説や問題から学びとったTips、テクニックを記録して公開することで、より理解が深まると思います。

例えば、この人。すごいです。

私も解いた問題の難易度・概要・感想などをまとめているので、そのうちデータベース化してGoogle App Engine上にでも公開しようかと思います。

シフト演算って何に使う？(2)

前回に引き続き、シフト演算について。

今回は、その活用法についてです。

シフト演算を使うことで簡単に実装できる例題を2つ取り上げます。

まず、シフト演算は、数値を2進数と考えて処理をしているということから、2進数に関連する演算に使えそうです。ということで、『シフト演算を用いた10進・2進数変換』を紹介します。

いろいろなサイトに10進数を2進数にしたり、その逆をやったりする場合のコードがありますが、どれも・・・・な感じです。(イケてないです。。)

まあ普通にやるなら、

10進数 -> 2進数変換は、再帰を使ってやります。

2進数 -> 10進数変換は、Hormerの公式を使ってやるでしょう。

今回は、シフト演算を使って実装します。

こんなかんじ。

const int B_SIZE = 32;
char bin[B_SIZE + 1];

void ToBase2(int n, char *bin) {
REP(i, B_SIZE)
bin[B_SIZE - 1 - i] = (n >> i & 1) + '0';
bin[B_SIZE] = '\0';
}
int ToBase10(char *bin) {
int ret = 0;
REP(i, B_SIZE)
ret += bin[i] - '0' << B_SIZE - 1 - i;
return ret;
}

すごいシンプルです。気を付けないといけない点は、シフト演算の優先順位です。+, -よりも優先度が低いということに注意して下さい。

そして、2つ目はべき乗の計算。

シフトが何故べき乗に関係あるのか・・・。ちょっと説明します。

簡単な例として、3^128 を計算してみましょう。

3 * 3 * 3* .... (128回掛算して) ... * 3

なんてやってたらダメです。

3^128 = (3^64)^2 = ((3^32)^2)^2 = ....と分解して行けば、n乗の計算量はO(log(n))まで縮みます。

今回は、128（2の累乗）という数字だったので簡単でしたが、例えば、3^11なんてのはどうしましょう？

3^11 = 3^(8+2+1)と考えましょう。指数部分を2進数に分解しているようなイメージです。そして、逐次2乗を計算していき、2進数の桁数が1の場合は、答えに積算といった感じです。

『百聞は一見に如かず』ということでコードをどうぞ。

double powerSht(double x, int n) {
double ret = 1.0;
double tmp = x;

REP(i, 32) {
if (n >> i & 1)
ret *= tmp;
tmp *= tmp;
}
return ret;
}

実際に普通にn回掛け算する場合と速度を比較しましたが、1.000001^10000000を計算した場合、

約2000倍ほど差がでました。

シフト演算って何に使う？

シフト演算って何に使ってますか？

とりあえず、簡単に分かるのは、『シフト演算は、加減乗除の演算より速い』ということ。それから、2進数に対してシフトを行うので、『被演算子は、左へ1シフトすれば2倍に、右に1シフトすれば2倍になる』ということ。

じゃー、実際にどれくらい速いのだろうということで、実験をしてみます。

ソースは以下のとおり。

int main() {
int x;
double t;

t = gettimeofday_sec();
REP(i, 1e8)
x = i + i;
printf("%lf\n", gettimeofday_sec() - t);

t = gettimeofday_sec();
REP(i, 1e8)
x = 2* i;
printf("%lf\n", gettimeofday_sec() - t);

t = gettimeofday_sec();
REP(i, 1e8)
x = i << 1;
printf("%lf\n", gettimeofday_sec() - t);

}

同様の処理を行う3つのループ。それぞれ、足算、掛算、シフト演算で計算を行っている。処理速度を見てみると・・・・。。

以下実行結果。

0.354000

0.274000

0.252000

うーーん、コンパイル時に掛算はシフト演算に置き換えられているのかもしれない。

（ちなみに、コンパイラはgcc3.4.5）

足算より掛算の方が遅くなるはずだが・・・

何はともあれ、シフト演算が速いのは間違いなさそう。そして、さり気にシフト演算を使っていると何となく”出来そう”に見える。

ということで、次回は、シフト演算を使ってどんな事ができるのか見て行きましょう。

2次元座標処理Tips

2次元平面上で何かをするような問題はよくある。

コンピュータ上で、実装する場合は、だいたい2次元配列を用意して処理をする。x[i][j]を座標(i, j)とみなして2次元平面を表現するわけだ。

例えば、”8個の接している点の中で自身と同じ値を持っているものの個数を数える”とか”上下左右の接している4点のみ移動できる時○○せよ。”なんて問題はよく見かける。

私も数問解いた記憶がある。。。

しかし、どうも上手く実装できない。

何故かコードが長くなってしまう。

関数化できない。。

こんなヘタレな悩みが昨日解決された。やっぱり人のソースコードを見るのは大事ですね。。

例によって、PKUの問題で説明します。

今回の問題は、Red and Black。まあよくある問題です。

で、これが私が書いたヘタレコード。。笑

class Plot {
public:
int x;
int y;
Plot(int x, int y) {
this->x = x;
this->y = y;
}

};

int main() {
int w, h;
for (;;) {
cin >> w >> h;
if (!w && !h)
break;

VS tiles;
REP(i, h) {
string tileLine;
char c;
REP(j, w) {
 cin >> c;
 tileLine += c;
}
tiles.PB(tileLine);
}

queue<Plot> prc;
REP(i, h) {
REP(j, w) {
 if (tiles[i][j] == '@') {
  prc.push(Plot(i, j));
  break;
 }
}
}

int ret = 0;
bool done[h][w];
memset(done, 0, sizeof(done));

while (prc.size()) {
int x = prc.front().x;
int y = prc.front().y;
prc.pop();
++ret;

if (x > 0) {
 if (tiles[x - 1][y] == '.' && !done[x - 1][y]) {
  prc.push(Plot(x - 1, y));
  done[x - 1][y] = true;
 }
}
if (x < h - 1) {
 if (tiles[x + 1][y] == '.' && !done[x + 1][y]) {
  prc.push(Plot(x + 1, y));
  done[x + 1][y] = true;
 }
}
if (y > 0) {
 if (tiles[x][y - 1] == '.' && !done[x][y - 1]) {
  prc.push(Plot(x, y - 1));
  done[x][y - 1] = true;
 }
}
if (y < w - 1) {
 if (tiles[x][y + 1] == '.' && !done[x][y + 1]) {
  prc.push(Plot(x, y + 1));
  done[x][y + 1] = true;
 }
}
}
cout << ret << endl;
}

return 0;
}

なんとか、4つの方向に同様の処理をしている部分を関数化したかったが、関数化するにも、2次元配列使っているので面倒臭いことになるなと考えて挫折。。。

で、他の人のソース見てたら、シンプルな方法があったので、頂きました。。（ごちそうさま。）

こんな感じ。

const int dx[] = { -1, 1, 0, 0 };
const int dy[] = { 0, 0, -1, 1 };

int main() {
int w, h;
for (;;) {
cin >> w >> h;
if (!w && !h)
break;

VS tiles;
REP(i, h) {
string tileLine;
char c;
REP(j, w) {
 cin >> c;
 tileLine += c;
}
tiles.PB(tileLine);
}

queue<Plot> prc;
REP(i, h) {
REP(j, w) {
 if (tiles[i][j] == '@') {
  prc.push(Plot(i, j));
  break;
 }
}
}

int ret = 0;
bool done[h][w];
memset(done, 0, sizeof(done));

while (prc.size()) {
int x = prc.front().x;
int y = prc.front().y;
prc.pop();
++ret;
REP(i, 4) {
 int xx = x + dx[i];
 int yy = y + dy[i];

 if (xx < 0 || xx > h - 1 || yy < 0 || yy > w - 1)
  continue;
 if (tiles[xx][yy] == '.' && !done[xx][yy]) {
  prc.push(Plot(xx, yy));
  done[xx][yy] = true;
 }
}
}
cout << ret << endl;
}

return 0;
}

いや、だいぶシンプル。

始めに、変化量を予め定数の配列で確保しておけばいいですね！！ここでポイントは、2次元方向の変化量を次元毎に個別に持たせて2重ループで回そうって考えないことです。変化パターンをまとめて書きだして、それぞれのx変化量、y変化量を配列dx, dyに格納しましょう。

フィボナッチあれこれ

世の中には、フィボナッチ数列がたくさん。。

例えば、うさぎのつがいの問題。これはあまりにも有名。

うさぎのつがいがいる。このつがいは1か月経つと、大人になる。そして2か月経つと子供を産み始める。さて、0か月目に1組のつがいが居たとして、nか月後にはつがいは何組になるでしょうか？

最初、見たとき意味分かりませんでした。そもそも”つがい”って何？？笑

雄と雌の１組のことを”つがい”と言うらしいです。（英語ではbreeding pairです。）

では、さっそく考えてみましょう。

0か月目: 1組。

1か月目: 1組が大人になる。でもまだ1組のまま。

2か月目: 大人になったつがいが子供を産むので2組。

3か月目：2か月目に生まれたつがいが大人になる。最初からいるつがいは子供を産むので、3組。

4か月目: 2か月目に生まれたつがいと最初からいるつがいが子供を産むので2組増えて5組。

……

はい、nか月目にいるつがいの数をn=0から列挙すると、

1,1,2,3,5,...

うむむ。なんかフィボナッチっぽい香りがしますね。。

それでは、より一般化して考えてみましょう。nか月目にいるつがいの数をf(n)とします。

f(n) = nか月目にいるつがいの数

= (n-1)か月目の時点で存在しているつがいの数 + nか月目に生まれるつがいの数

　　 = (n-1)か月目の時点で存在しているつがいの数 + (n-1)か月目の時点で大人であるつがいの数

=(n-1)か月目の時点で存在しているつがいの数 + (n-2)か月目の時点で存在しているつがいの数 　　= f(n-1) + f(n-2)

上の考え方は、かなり大事です。

トップダウン的な思考で、より一般的な現象を把握します。

上の式より、

f(n) = f(n-1) + f(n-2)

なので、フィボナッチですね。プログラムで解くときは、メモ化再帰か動的計画です。個人的にはフィボナッチは、メモ化再帰で解くのが好きです。フィボナッチに限らず漸化式系の問題一般に適用できるからです。

他にも、フィボナッチ数列になる現象はいろいろあります。

例えば、

あなたは、1段、または、2段ずつ階段を上れます。n段目まで上がる場合何通りの上り方があるでしょうか？

2*nの長方形の形をした部屋があります。この部屋に1*2の畳を敷き詰めるとき、敷き詰め方は何通りあるでしょうか？

どちらの問題も有名どころです。。アルゴリズマーの中では、常識の範疇です。

上記２つの問題も最初問題を把握するために簡単な例を解いてみて（ボトムアップ的思考）、そのあと、問題の意味を把握したら、n番目の状態とn-1番目、n-2番目の関係を考えてみる（トップダウン的思考）というやり方をするとすぐに分かると思います。