60706


**行いたいこと**  
組み合わせの式 64C32 の結果の下9桁を出力するプログラムをC言語で書きたい

**試したこと**  
以下のようなコードで実行を試みたが、オーバーフローが生じてしまう

[code]

    #include<stdio.h>
    #define ll long long
    
    ll combi(ll n, ll r){
      if (r == 0) {
          return 1;
      }
      return (n - r + 1) * combi(n, r - 1) / r;
    }
    
    int main(){
      printf("%09lld\n",combi(64 ,32 ) );
      return 0;
    }
    
[/code]

解答は 942590534 なのですが、-275878929539273115 のような結果になってしまいます。  
どのようにコードを書けばいいでしょうか。


CC BY-SA 4.0
2019-11-20T11:04:09.067
60698
2019-11-22T07:04:11.800
2019-11-22T00:56:34.773
3060
34550
question
c アルゴリズム
組み合わせ計算で 64C32 の結果がオーバーフロー


BigInteger が使えると簡単なんですけどね・・。  
独自で BigInteger の割り算を実装しようかと思ったのですが、  
分数の約分で実装してみました。

int の配列を用意して 2 * 3 * 4 * 5 を  
int[] a = {0,0,1,1,1,1}  
と表します。  
このa を 3 倍する時は a[3] = a[3] + 1;  
5倍する時は a[5] = a[5] + 1;  
と表します。

[code]

    /** 掛け算 */
    void mul(int ar[], int n) {
        if (n > 0 && n < MAX_AR) {
            ar[n] ++;
        }
    }
    
[/code]

また 約分をするための関数を作ります。  
64 は  
a[64] = 1;  
で表し  
これを 約分すると 2 * 2 * 2 * 2 * 2 * 2 なので 2 を 6 回 掛ける事なので  
a[2] = 6 と表します。

[code]

    void yakubun(int x[]) {
        int prime[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47 };
        for (int i = 0; i < 15 ; i++) {
            int pi = prime[i];
            for (int j = (MAX_AR-1)/i ; j > 1; j --) {
                int ij = i * j;
                if (x[ij] > 0) {
                    // 例：３の５倍 （１５） に 2 が 入っている場合 x[3] と x[5] に 2 を足す;
                    x[i] += x[ij];
                    x[j] += x[ij];
                    x[ij] = 0;
                }
            }
        }
    }
    
[/code]

次に 割り算を作ります。  
分数の約分を使って、  
必ず割り切れるという前提での実装です。  
汎用的には使えません。

[code]

    /** 2つの数字の割り算 */
    long long div(int b[], int a[]) {
        int c[MAX_AR];
        for (int i = 0; i < MAX_AR; i++) {
            c[i] = 0;
        }
        // b / a を実施する 分母と分子で同じ物を約分する 
    
        // 配列に入っている数字を約分する
        yakubun(a);
    
        yakubun(b);
    
        for (int i = 0; i < MAX_AR; i++) {
            if (b[i] >= a[i]) {
                c[i] = b[i] - a[i];
            }
        }
    
        // 最終的な計算結果
        long long s = 1;
        for (int i = 1; i < MAX_AR; i++) {
            if (c[i] > 0) {
                for (int j = 0; j < c[i]; j++) {
                    s = s * i;
                    // 掛け算の余りを使うとオーバーフローせず 正しい答えが出せる
                    if (s >= 1000000000) {
                        s = s % 1000000000;
                    }
                }
            }
        }
        return s;
    }
    
[/code]

最後に combi は 再帰を使わず 掛け算と割り算に展開します。  
`nCr = ( n から 数を下げながら r 個の掛け算) / (1 から 数を上げながら r 個の掛け算)`

[code]

    long long combi(long l, long r) {
    
        int a[MAX_AR];
        int b[MAX_AR];
    
            // 変数初期化
        for (int j = 0; j < MAX_AR; j++) {
            a[j] = 0;
            b[j] = 0;
        }
    
        for (int i = 0; i < r; i++) {
                    // 1 から r までの掛け算
            mul(a, i + 1);  // a = a * (i + 1);
                    // l から r + 1 までの掛け算
            mul(b, l - i); // b = b * (64 - i);
        }
    
        return div(b, a);
    }
    
    
[/code]

main 関数と、定数宣言等

[code]

    #define MAX_AR 65
    
    int main()
    {
        char ans[32];
        char last9[32];
    
        sprintf_s(ans, "%lld", combi(64, 32));
    
        strcpy_s(last9 , &ans[strlen(ans) - 9]);
    
        printf("%s\n", last9);
    
        return 0;
    }
    
[/code]

割り算を 手本なしに作る自信がなかったので 約分という方法を使ってみました。  
かける数も 64 までの 小さな数という事が前提なので 汎用的な BigInteger の掛け算  
を作っていません。

C# で BigInteger を使えるなら

[code]

    // 64C32 を計算する
    
                BigInteger a = 1, b = 1;
    
                for(int i=0;i<32;i++)
                {
                    a = a * ( i + 1);
                    b = b * (64 - i);
                }
                BigInteger c = b / a;
    
[/code]

たったこれだけです。

`482219923991114978843459072919892677776312893440000000 /
263130836933693530167218012160000000 = 1832624140942590534`

となります。

上記プログラムで MAX_AR を大きくして それに合わせた素数を計算する事で  
かなり大きな数でも 高速に計算できることを確認しました。


CC BY-SA 4.0
2019-11-20T14:36:10.873
60705
2019-11-22T05:20:21.877
2019-11-22T05:20:21.877
18851
18851
60698
answer


オーバーフローが起こっているのはこの計算をするときです。

[code]

    (n - r + 1) * combi(n, r - 1) / r
    
[/code]

たとえば `combi(64, 32)` の呼び出しでは以下の計算をすることになりますが、

[code]

    33 * combi(64, 31) / 32
    
[/code]

₆₄C₃₁ = 1,777,090,076,065,542,336 なので、これに 33 を掛けると 2⁶⁴ - 1 =
9,223,372,036,854,775,807 を超えます。

したがって、このオーバーフローが起こらないように計算式を修正する必要があります。剰余がある元で二項係数を求める方法は[色々あります](https://w.atwiki.jp/uwicoder/pages/2118.html)が、今回の場合は[パスカルの三角形](https://ja.wikipedia.org/wiki/%E3%83%91%E3%82%B9%E3%82%AB%E3%83%AB%E3%81%AE%E4%B8%89%E8%A7%92%E5%BD%A2)で求めるのが簡単そうです。つまり、

[code]

    nCr = (n-1)Cr + (n-1)C(r-1)
    
[/code]

という式を元に順番に足し合わせることで求めることができます。これなら途中式が足し算だけなので計算途中に剰余をとることができます。


CC BY-SA 4.0
2019-11-20T14:50:37.427
60706
2019-11-22T06:56:06.220
2019-11-22T06:56:06.220
19110
19110
60698
answer


64ビット符号付整数の最大値は2⁶³-1 = 9223372036854775807 なので、nekketsuuu さんの書かれたように `33 *
combi(64, 31)`を計算する際にこの値を超えてオーバーフローしてしまいます。

ですが、64C32の値は 1832624140942590534 なので、計算の順番を工夫してやれば、ギリギリ正しい結果は出せそうです。

[code]

    ll gcd(ll a, ll b) {
        ll temp;
        while (b != 0) {
            temp = a % b;
    
            a = b;
            b = temp;
        }
        return a;
    }
    
    ll combi1(ll n, ll r){
        if (r == 0) {
            return 1;
        }
        //掛け算を先にするとオーバーフローする
        //return (n - r + 1) * combi(n, r - 1) / r;
        //と言うわけで掛ける前に割っておく
        ll n1 = n - r + 1;
        ll n2 = combi1(n, r - 1);
        ll d = r;
        //n2/d だけでは割り切れない可能性があるので、最大公約数を求める
        ll div = gcd(n2, d);
        n2 /= div;
        d /= div;
        //組み合わせは必ず整数になるからもう割り切れる
        n1 /= d;
        //割り算が終わったら掛ける
        return n1 * n2;
    }
    
[/code]

下9桁が欲しい場合には 1000000000 (0が9個)で割った余りを求めます。

[code]

    int main(int argc, const char * argv[]) {
        //参考
        printf("INT64_MAX=%lld\n", INT64_MAX);
    
        ll c = combi1(64, 32);
        printf("%dC%d=%lld\n", 64, 32, c);
    
        //下9桁のみの出力
        printf("%09lld\n", c % 1000000000);
    
        return 0;
    }
    
[/code]

「結果の下9桁」だけを求めれば良いと言う出題の条件を利用するには、nekketsuuu
さんが書かれたように途中計算を足し算だけにすると言うのは有効なのですが、`nCr = (n-1)Cr +
(n-1)C(r-1)`を素朴に再起的に表現すると、計算量が爆発して64C32程度でも実用的な時間では解けません。

[code]

    //nCr = (n-1)Cr + (n-1)C(r-1) を素朴に実装する
    ll combi2(ll n, ll r) {
        if( r > n - r ) {
            return combi2(n, n - r);
        }
        if( r == 0 ) {
            return 1;
        }
        return (combi2(n - 1, r) + combi2(n - 1, r - 1)) % 1000000000;
    }
    
[/code]

「順番に足し合わせる」と言う趣旨の部分を活かす、つまりモロにパスカルの三角形をメモリ上で計算しながら作る、と言う方針で行くと、こんなコードになります。

[code]

    //デバッグ用、nが大きい値の時は綺麗に出力されない
    void dumpPascalTriangle(ll ** pt, ll n) {
        for( int i = 0; i <= n; ++i ) {
            for( int j = 0; j < n - i; ++j ) {
                putchar(' ');
            }
            for( int j = 0; j <= i; ++j ) {
                printf(" %lld", pt[i][j]);
            }
            putchar('\n');
        }
    }
    //パスカルの三角形をもろに計算する
    ll combi3(ll n, ll r) {
        ll **pascalTriangle = (ll **)malloc((n + 1) * sizeof(ll *));
        for( int i = 0; i <= n; ++i ) {
            pascalTriangle[i] = (ll *)malloc((i + 1) * sizeof(ll));
        }
        pascalTriangle[0][0] = 1;
        for( int i = 1; i <= n; ++i ) {
            for( int j = 0; j <= i; ++j ) {
                ll n1 = j - 1 < 0 ? 0: pascalTriangle[i - 1][j - 1];
                ll n2 = j > i - 1 ? 0: pascalTriangle[i - 1][j];
                pascalTriangle[i][j] = (n1 + n2) % 1000000000;
            }
        }
        //dumpPascalTriangle(pascalTriangle, n);
        return pascalTriangle[n][r];
    }
    
[/code]

これだと計算量はたかだかnの2乗程度、特に工夫しなくても64C32程度ならすぐに計算が終わります。

5C2の時に`dumpPascalTriangle(pascalTriangle, n);`をコードインした場合の出力例:

[code]

          1
         1 1
        1 2 1
       1 3 3 1
      1 4 6 4 1
     1 5 10 10 5 1
    
[/code]

(左から3つ目の 10 が5C2の値を表している。)

* * *

「計算の順番を工夫」の方は、頑張ればもう少し簡単に出来そうですが、とりあえずすぐ思いつく範囲だとこんなもんですね。


CC BY-SA 4.0
2019-11-20T22:45:49.647
60709
2019-11-20T22:45:49.647
13972
60698
answer


C++のコードですがこのようにすることもできます。

[code]

    //aとbの最大公約数
    int gcd(int a, int b) {
        switch (b) {
        case 0: return a;
        default: return gcd(b, a % b);
        }
    }
    long long int comb(int n, int r) {
        //n!
        std::vector<int> numerator;
        for (int i = n; i > 0; --i)
            numerator.push_back(i);
    
        //numeratorをr!で割る
        for (auto i = 1; i <= r; ++i) {
            int d = i;
            for (int j = 0; j < numerator.size(); ++j) {
                int g = gcd(d, numerator[j]);
                d /= g;
                numerator[j] /= g;
            }
        }
        //numeratorを(n-r)!で割る
        for (int i = 1; i <= n - r; ++i) {
            int d = i;
            for (int j = 0; j < numerator.size(); ++j) {
                int g = gcd(d, numerator[j]);
                d /= g;
                numerator[j] /= g;
            }
        }
    
        //numeratorを掛け合わせる。
        long long int result = 1;
        for (int i = 0; i < numerator.size(); ++i) {
            result = (result * numerator[i]) % MOD;
        }
        return result;
    }
    
[/code]

基本的な考えは、分子の階乗部分と分母の階乗部分を掛け合わせる前に約分していくというものです。

* * *

今回の問題に関しては明らかに不要な効率化ですが、2と5の素因数を含まない数についてオイラーの定理等を使って計算して、後で2と5を掛けるという方法もあります。

[code]

    long long int power(long long int base, long long int exp) {
        switch (exp) {
        case 0: return 1;
        case 1: return base;
        default: return power(base * base % MOD, exp / 2) * power(base, exp % 2) % MOD;
        }
    }
    //n!をfactorで割ることができる回数
    int count_factor(int n, int factor) {
        switch (n) {
        case 0: return 0;
        default: return n / factor + count_factor(n / factor, factor);
        }
    }
    int remove2_5(int n) {
        if (n % 2 == 0) return remove2_5(n / 2);
        if (n % 5 == 0) return remove2_5(n / 5);
        return n;
    }
    long long int comb(int n, int r) {
        long long int numerator = 1, denominator = 1;
        for (auto i = 1; i <= n; ++i) {
            numerator = numerator * remove2_5(i) % MOD;
        }
        for (auto i = 1; i <= r; ++i) {
            denominator = denominator * remove2_5(i) % MOD;
        }
        for (auto i = 1; i <= n - r; ++i) {
            denominator = denominator * remove2_5(i) % MOD;
        }
        return numerator * power(denominator, 400000000 - 1) % MOD *
            power(2, count_factor(n, 2) - count_factor(r, 2) - count_factor(n - r, 2)) % MOD *
            power(5, count_factor(n, 5) - count_factor(r, 5) - count_factor(n - r, 5)) % MOD;
    }
    
[/code]

これだとn = 10^6ぐらいでも高速に求められます。

[code]

    comb(1234567, 5000) #=> 772988800
    
[/code]


CC BY-SA 4.0
2019-11-21T05:53:19.810
60723
2019-11-21T08:55:09.340
2019-11-21T08:55:09.340
33033
33033
60698
answer