2

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

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

#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;
}

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

6
  • 64C32を手元で計算すると"1832624140942590534"になるのですが、その解答は本当にあってますか?
    – raccy
    Commented 2019年11月20日 12:05
  • 得られた計算結果の下9桁を出力するプログラムなので間違いないです。わかり辛くて申し訳ないです。 Commented 2019年11月20日 12:11
  • 下9桁のみ表示すればいいと言うことですね?10桁以上の数値は求めなくても良いと言うことであってますね?
    – raccy
    Commented 2019年11月20日 12:16
  • 1
    計算の順序を変えれば long long int に収まりますが… combi(60, 28)/(5*29)*(2*21*61) ま、参考までにどうぞ。
    – user39889
    Commented 2019年11月20日 13:16
  • 3
    この問題はどこから持ってきたものですか?プログラミングコンテストの問題っぽくもあるのですが、そういうところでは(特にコンテスト中の)他者との相談を禁止しているところが多いです。できれば出典を明記していただけませんか?
    – letrec
    Commented 2019年11月20日 13:16

4 件の回答 4

8

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

(n - r + 1) * combi(n, r - 1) / r

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

33 * combi(64, 31) / 32

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

したがって、このオーバーフローが起こらないように計算式を修正する必要があります。剰余がある元で二項係数を求める方法は色々ありますが、今回の場合はパスカルの三角形で求めるのが簡単そうです。つまり、

nCr = (n-1)Cr + (n-1)C(r-1)

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

3

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

//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;
}

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


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

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;
}

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

comb(1234567, 5000) #=> 772988800
1
  • stackoverflow 英語版に 似たようなコードが ありました。 stackoverflow.com/a/58558434 約分を使う方法が一番 最速で効率的ですよね。 私のは今一歩でした。
    – user18851
    Commented 2019年11月21日 23:38
2

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

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

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;
}

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

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;
}

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

//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;
}

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

//デバッグ用、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];
}

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

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

      1
     1 1
    1 2 1
   1 3 3 1
  1 4 6 4 1
 1 5 10 10 5 1

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


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

1

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;
と表します。

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

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

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) {
                // 例:3の5倍 (15) に 2 が 入っている場合 x[3] と x[5] に 2 を足す;
                x[i] += x[ij];
                x[j] += x[ij];
                x[ij] = 0;
            }
        }
    }
}

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


/** 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;
}

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

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);
}

main 関数と、定数宣言等

#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;
}

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

C# で BigInteger を使えるなら

// 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;

たったこれだけです。

482219923991114978843459072919892677776312893440000000 / 263130836933693530167218012160000000 = 1832624140942590534

となります。

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

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。