コンビネーションnCrの計算

[keito]

今回が初投稿です。（C言語歴5ヶ月ほどの初心者、for、do~while文、文字列の簡単な処理は理解している位）

コンビネーション（組合せ）nCrを計算するプログラムを作ったのですが、nとrが大きな値になると答えが合わなくなってしまいます。

↓ソースコード

コード:

//Microsoft Visual Studio Express 2015 Windows for Desktop
//Visual C++

#include<stdio.h>
#include<stdlib.h>
double combination(int, int);

int main(void) {
	int n = 0, r = 0;
	char str[8];

	do {								//n,rの入力、但し0<=r<=n
		printf("n=");
		fgets(str, sizeof(str), stdin);
		n = atoi(str);
	} while (n < 0);
	do {
		printf("r=");
		fgets(str, sizeof(str), stdin);
		r = atoi(str);
	} while (n < r);

	printf("%dC%d = %.0f\n", n, r, combination(n, r));
	getchar();
	return 0;
}

double combination(int n, int r) {		//nCr の計算
	int i;
	double C = 1;
	for (i = 0; i < r; i++) {
		C *= n - i;
		C /= i + 1;
	}
	return C;
}

100C50をこれでやってみると、100891344545564184479528910848　とでました。
（実際は 100891344545564193334812497256 なので、約9兆ほどズレてる。
他の場合も自分のPCで試した結果、nが60以下の時は問題無いが、nが61を超えると、rの値によっては答えがずれてしまう）

自分が思いつく限りのソースコードの修正をしてはみましたが、これ以上改善しませんでした・・・。

最終的にはnが100以下ならば絶対答えが合うようにしたいのですが、分かる方はアドバイスをお願いします。

[みけCAT]

double型の有効桁数を超えていますね。
ライセンス的に問題が無ければ、The GNU MP Bignum Libraryなどの多倍長計算ライブラリを使うのが簡単でしょう。

[みけCAT]

この場合は、nが100までなら「1桁」を10000くらいに設定すれば比較的めんどくさい複数「桁」同士の掛け算や割り算は発生しないので、自前で多倍長計算を実装してもいいでしょう。

[keito]

ありがとうございます！
有効桁数の問題でしたか・・・なんでもかんでもdouble型に持っていけば何とかなると勝手に思っていました・・・。

では多倍長計算ライブラリを利用する方法を採用して解決してみます。

[かずま]

こんなのはどうでしょうか？

コード:

#include <stdio.h>

#define N 100

unsigned long long comb[2*N+1][2*N+1][2];

int main(void)
{
    int i, j, n, r, carry;

    for (i = 0; i <= N; i++)
        comb[i][0][0] = comb[0][i][0] = 1;

    for (n = 2; n <= 2*N; n++)
        for (j = 1; j <= n; j++) {
            i = n - j;
            if (i > 0 && i <= 2*N && j > 0 && j <= 2*N) {
                comb[i][j][0] = comb[i-1][j][0] + comb[i][j-1][0];
                carry = 0;
                if (comb[i][j][0] >= 1000000000000000000ULL) {
                    comb[i][j][0] -= 1000000000000000000ULL;
                    carry = 1;
                }
                comb[i][j][1] = comb[i-1][j][1] + comb[i][j-1][1] + carry;
            }
        }

    while (printf(">> "), scanf("%d%d", &n, &r) == 2) {
        i = n - r;
        if (comb[i][r][1])
            printf("%llu%018llu\n", comb[i][r][1], comb[i][r][0]);
        else
            printf("%llu\n", comb[i][r][0]);
    }
    return 0;
}

[かずま]

ちょっと無駄があったので修正。

コード:

#include <stdio.h>

#define N 100

long long comb[N+1][N+1][2];

int main(void)
{
    int i, n, r, carry;

    for (i = 0; i <= N; i++)
        comb[i][0][0] = comb[0][i][0] = 1;

    for (n = 2; n <= N; n++)
        for (r = 1; r < n; r++) {
            i = n - r;
            comb[i][r][0] = comb[i-1][r][0] + comb[i][r-1][0];
            carry = 0;
            if (comb[i][r][0] >= 1000000000000000000LL) {
                comb[i][r][0] -= 1000000000000000000LL;
                carry = 1;
            }
            comb[i][r][1] = comb[i-1][r][1] + comb[i][r-1][1] + carry;
        }

    while (printf(">> "), scanf("%d%d", &n, &r) == 2) {
        i = n - r;
        if (comb[i][r][1])
            printf("%lld%018lld\n", comb[i][r][1], comb[i][r][0]);
        else
            printf("%lld\n", comb[i][r][0]);
    }
    return 0;
}

[かずま]

結局、パスカルの三角形を作るだけのことでした。

コード:

#include <stdio.h>

#define N 100

long long comb[N+1][N+1][2];

int main(void)
{
    unsigned n, r, carry;

    for (n = 0; n <= N; n++)
        comb[n][0][0] = comb[n][n][0] = 1;
    for (n = 2; n <= N; n++)
        for (r = 1; r < n; r++) {
            comb[n][r][0] = comb[n-1][r-1][0] + comb[n-1][r][0];
            carry = comb[n][r][0] >= 1000000000000000000LL;
            if (carry) comb[n][r][0] -= 1000000000000000000LL;
            comb[n][r][1] = comb[n-1][r-1][1] + comb[n-1][r][1] + carry;
        }
    while (printf(">> "), scanf("%d%d", &n, &r) == 2 && n <= N && r <= N)
        if (comb[n][r][1])
            printf("%lld%018lld\n", comb[n][r][1], comb[n][r][0]);
        else
            printf("%lld\n", comb[n][r][0]);
    return 0;
}

long long の変数を 2つ使うだけの、多倍長(2倍長)です。

[keito]

ライブラリのビルド方法がよく分からず・・・。という訳で今回ライブラリを使用する方法は諦めます・・・。

かずまさんの方法のが一番簡潔でソースも短いですね、こちらを参考にしてみます。

[かずま]

元のアルゴリズムを多倍長で実装してみました。

コード:

#include <stdio.h>
#include <stdlib.h>

#define N 100 // max n of nCr
#define M 8   // multiple length of int

void print(const int *c)
{
    int i = M;
    while (c[--i] == 0) ;
    printf("%d", c[i]);
    while (--i >= 0) printf("%04d", c[i]);
    printf("\n");
}

void mul(int *c, int n)
{
    int i, carry = 0;
    for (i = 0; i < M; i++) {
        c[i] = c[i] * n + carry;
        carry = c[i] / 10000;
        c[i] %= 10000;
    }
}

void dvd(int *c, int n)
{
    int i = M, borrow = 0;
    while (--i >= 0) {
        c[i] += borrow * 10000;
        borrow = c[i] % n;
        c[i] /= n;
    }
}

void combination(int n, int r, int *c)
{
    int i;
    c[0] = 1;
    for (i = 1; i < M; i++) c[i] = 0;
    for (i = 0; i < r; i++) mul(c, n - i), dvd(c, i + 1);
}

int main(void)
{
    int n, r, c[M];
    while (printf(">> "), scanf("%d%d", &n, &r) == 2 && n >= 0 && n <= N && r <= n) {
        combination(n, r, c);
        printf("C(%d,%d) = ", n, r);
        print(c);
    }
    return 0;
}

int の変数を 8つ使う多倍長です。

[かずま]

combination が 0 になることはないのですが、
print が 0 を表示できないバグがあったので修正します。
また、mul がオーバーフローを返し、dvd が余りを返すようにしました。

コード:

void print(const int *c)
{
    int i = M;
    while (i > 0 && c[--i] == 0) ;
    printf("%d", c[i]);
    while (--i >= 0) printf("%04d", c[i]);
    printf("\n");
}
 
int mul(int *c, int n)
{
    int i, carry = 0;
    for (i = 0; i < M; i++) {
        c[i] = c[i] * n + carry;
        carry = c[i] / 10000;
        c[i] %= 10000;
    }
    return carry;
}
 
int dvd(int *c, int n)
{
    int i = M, borrow = 0;
    while (--i >= 0) {
        c[i] += borrow * 10000;
        borrow = c[i] % n;
        c[i] /= n;
    }
    return borrow;
}