Open9

C++で組み合わせ数を計算するコードを考える

log5log5

https://mathtrain.jp/permucombi

AtCoderの問題を解くときに必要になったので...

案1:パスカルの三角形を使う

#include <bits/stdc++.h>

using namespace std;
using ll = int64_t;

void comb(vector<vector<ll> > &v) {
    for (int i = 0; i < v.size(); i++) {
        v[i][0] = 1;
        v[i][i] = 1;
    }
    for (int k = 1; k < v.size(); k++) {
        for (int j = 1; j < k; j++) {
            v[k][j] = (v[k - 1][j - 1] + v[k - 1][j]);
        }
    }
}
ll ncr1(ll n, ll r) {
    vector<vector<ll>> v(n + 1, vector<ll>(n + 1, 0));
    comb(v);
    return v[n][r];
}

案2: 公式を使う

#include <bits/stdc++.h>

using namespace std;
using ll = int64_t;

ll ncr2(ll n, ll r) {
    if (r > n) return 0;
    if (r * 2 > n) r = n - r;
    if (r == 0) return 1;

    ll result = n;
    for (ll i = 2; i <= r; ++i) {
        result *= (n - i + 1);
        result /= i;
    }
    return result;
}
log5log5

結果確認用に以下を実行してみた
https://paiza.io/projects/267jeWJF0uKRLGsiw42vYw

int main() {

    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);
    ll MOD = 1000000007;
    for (int i = 100; i < 1000; ++i) {
        for (int j = 0; j < i + 1; ++j) {
            ll c1 = ncr1(i, j);
            ll cc = ncr2(i, j);
            if (c1 != cc) {
                printf("Result not match! on C(%d, %d) -- expected:%lld, actual:%lld \n", i, j, c1, cc);
                exit(0);
            }
        }
    }
    cout << "check finished!" << endl;
}

結果

Result not match! on C(100, 16) -- expected:1345860629046814650, actual:192939124439967674 

Process finished with exit code 0
log5log5

なお _{100}C_{16}=1~345~860~629~046~814~650 なので ncr1 が正しく計算できている。

ただし、
1~345~860~629~046~814~650-192~939~124~439~967~674\\=1~152~921~504~606~846~976=2^{60}
なので、 ncr2 の計算方法が誤っているとは言えなさそう...

なんかコンピュータの都合かなにかで切り取られてしまった印象...

log5log5

for 文の範囲が広すぎるかなと思って縮めてみた
https://paiza.io/projects/eMiC1dYgksg_w8Vou8BFWw

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);
    ll MOD = 1000000007;
    for (int i = 1; i < 100; ++i) {
        for (int j = 0; j < i + 1; ++j) {
            ll c1 = ncr1(i, j);
            ll cc = ncr2(i, j);
            if (c1 != cc) {
                printf("Result not match! on C(%d, %d) -- expected:%lld, actual:%lld \n", i, j, c1, cc);
                exit(0);
            }
        }
    }
    cout << "check finished!" << endl;
}

結果は

Result not match! on C(62, 28) -- expected:349615716557887465, actual:-309196571788882235 
log5log5

_{62}C_{28}=349~615~716~557~887~465 なのでまたも ncr1 が正しいが...

349~615~716~557~887~465+309~196~571~788~882~235\\=658~812~288~346~769~700\\=100100100100100100100100100100100100100100100100100100100100_2\\=44444444444444444444_8

なので、 何かありそう な気はする...
(ラップアラウンドは言うまでもないとして)

log5log5

なんか ncr2 の実装おかしいんですかね...
(forループが2重になっていないので ncr2 に期待したいのに...)

log5log5

ncr2、おかしくはないんだろうけど、掛け算が壁にぶち当たってる感じがあるなぁ...

ためしに int64_t のかわりに uint64_t を使ってみたら、 _{63}C_{28} 自体は正しく計算できているので、やはり「壁」の突破が問題かな...
https://paiza.io/projects/N12z9nFr1qsIwhfLwDmC-A

log5log5

ncr2、もしかして逆側から掛けていったほうがよさそう?

いや、大きい数で割り切れなくなるから無理ぽ..

log5log5

こういう仕様として受け入れたほうが良い気がしてきた。

派生タスクとして1000000007等で割った値を上手く計算できるように改造するのが良さそう