タイトルが紛らわしいですが、車輪の再発明をしてみたという話です。
この記事は
多倍長ができないC++で普通に組み合わせの計算をしようとすると、階乗でオーバーフローするのはよく知られていると思いますが、その上でパスカルの三角形的に計算しようと実装したら、自力では全然速度が出なかったので、次回以降のため、自分用のメモを含めて記録しようという記事です。
自力コード１
まずは組み合わせの式に基本的に成り立つ漸化式？と、基本の値
{}_nC_r = {}_{n-1}C_{r-1} + {}_{n-1}C_{r}\\
{}_nC_n = {}_nC_0 = 1
(漸化式にタイポがあったので修正しました)
を見て、再帰関数にしたらめちゃくちゃ簡単に書けるんじゃないかと思って書いたのがこちら
long long combi(long long n, long long k) {
if (n == k || k == 0)
return 1;
else {
return combi(n - 1, k - 1) + combi(n - 1, k);
}
}
実際これはnが小さいうちは当然問題なく動作したが、n=30あたりからかなり辛くなってきた。
繰り返し計算している部分が多いはずなので、メモ化したらいいのか？と思って（雑に）メモ化したのが次の自力コード２になる。
自力コード２
long long combi(long long n, long long k) {
std::map<std::string, long long> memo;
if (n == k || k == 0)
return 1;
else {
std::string key = std::to_string(n) + std::to_string(k);
if (memo.count(key) != 0) {
return memo[key];
} else {
ll out = combi(n - 1, k - 1) + combi(n - 1, k);
memo[key] = out;
return out;
}
}
}
雑なメモ化で早くなった試しがないのでだめだろうなと思いつつも回してみると、案の定大幅な改悪に成功した。(n=25ですでにn=30の自力コード１よりも遅い)実際はkについて計算したらn-kについても計算したことになるので、そちらもメモしたら多少マシになるのかなとも思ったが、すでに改悪されているので手間の割に効果は無いだろうと実装はしなかった。
2018/3/23/7:36追記
キーの競合が発生するので実際は
to_string(n)+ ","+ to_string(k);
とするべきでした。
2018/3/23/8:58 追記
まさかのmemoがスコープ的に共有されていないというありえないミスをしていたことに気づいたので、修正したところメモ化すれば再帰してもなんとか実用レベルの速さになりました（60C30で4000 microseconds）
自力ではもう思いつかなかったので、google先生に頼って最終的に行き着いたのが次のコードになる。
最終コード
このコードをかくにあたって以下のサイトを参考にさせて頂いた（というかほとんどそのままお借りする形になってしまった。）
caprestのブログ C++ vectorクラスで二次元配列 & Combination の計算
http://caprest.hatenablog.com/entry/2016/05/29/181102
std::vector<std::vector<long long>> comb(int n, int r) {
std::vector<std::vector<long long>> v(n + 1,std::vector<long long>(n + 1, 0));
for (int i = 0; i < v.size(); i++) {
v[i][0] = 1;
v[i][i] = 1;
}
for (int j = 1; j < v.size(); j++) {
for (int k = 1; k < j; k++) {
v[j][k] = (v[j - 1][k - 1] + v[j - 1][k]);
}
}
return v;
}
返り値に対して、v[n][r]が${}_nC_r$になる。元プログラムでは余計な代入などを行っていないので、私が書いたものよりも更に高速に動作する（すごい）が、これでも自力のよりものすごく早い。
時間測定
自力２は測定が不可能なくらい遅いので省略(自力１も死ぬほど遅いけど)
| 式 | 自力１(microsecond) | 最終(microsecond) |
|---|---|---|
| 30_C_15 | 701192 | 40 |
| 35_C_17 | 20369929 | 24 |
| 60_C_30 | - | 86 |
2018/3/23/7:36 追記
元々70C35が表にありましたが、この時点でオーバーフローしているので消去しました。
最初に書いた再帰のプログラムも最終的なプログラムも計算しているのはたかだかパスカルの三角形の範囲に収まる節点の数なはずなのだが、こうも差が出ると、最初のプログラムの計算量が考えているよりも多いのか、はたまた関数の呼び出しとif文の判定にかなり時間が取られてしまっているのか。
あとでゆっくり考えてみたい。