組み合わせ順列

質問

nＣrを求める関数combination(n,r)をC言語で作りたいのですが、どうすればよいか教えてください。また、参考となるようなサイトを教えてください。僕の作った関数だと、すぐに桁あふれになってしまいます。そのことを考慮して、桁あふれにになりにくい関数もつくりました。これは パスカルの三角形の関係を使ったnＣr=n-1Ｃr +n-1Ｃr-1の関係を使っての再帰関数です。しかし、これだと結果を出すのに時間がかかってしまいます。僕がつくった関数をいくつか出しておきますのでいい考えがあれば教えてください。64Ｃ32が高速に正確に出れれば最高です。
long combination(int n,int r)
{
int i, a, b;
long c;
if(r>n/2) r=n-r;
a= n;
b= 1;
c= 1L;
for(i=0 ;i< r ;i++){
c= c* a/ b;
a--;
b++;
}
return c;
}
これはすぐに桁あふれになってしまう。
long combination(int n,int r)
{
int i, a, b;
double c;
if(r>n/2) r=n-r;
a= n;
b= 1;
c= 1.0;
for(i=0 ;i< r ;i++){
c= c/b*a;
a--;
b++;
}
return (long)c;
}
これはｃをdoubleにして計算する分、丸めこみが生じ誤差がでる。
long combination(int n,int r)
{
long c;
if(r> n-r) r=n-r;

if(r==0) return 1;
else if(r==1) return n;
else{
c= combination(n-1,r-1)+ combination(n-1,r);
return c;
}
}
これは誤差は出ず正確であるがいかんせん遅い！

回答No.9

επιστημη（@episteme） · ベストアンサー

お詫びに:

#include <cstdio>

__int64 combination(int n, int r) {
__int64* work = new __int64[n+1]; /* n は 64 まで */
for (int i = 0; i <= n; i++) {
for (int j = i - 1; j > 0; j--) {
work[j] += work[j - 1];
}
work[i] = 1;
}
__int64 result = work[r];
delete[] work;
return result;
}

namespace std {}
using namespace std;

int main() {
for ( int i = 0; i < 33; ++i )
printf("64C%d = %I64d\n", i, combination(64,i));
return 0;
}

----- 結果 -------
.....
64C25 = 401038568751465792
64C26 = 601557853127198688
64C27 = 846636978475316672
64C28 = 1118770292985239888
64C29 = 1388818294740297792
64C30 = 1620288010530347424
64C31 = 1777090076065542336
64C32 = 1832624140942590534

だそうです。# VC6でコンパイル/実行

お礼

なんどもお返事ありがとうございます。大変感謝しております。64C32 = 1832624140942590534 と分かっただけでも大変参考になりました。これから自分がつくらないといけないものは1024C512なのでガンバってみます。これからも何かあればよろしくお願いします。

回答No.8

επιστημη（@episteme）

あああ、失礼しましたっ。

ダイナミックレンジ上げるためにunsigned long
にしちまったのが敗因でした。

# __int64にしとけばよかった...

回答No.7

noname#30727

>しくしく、吹っ飛びます。
>
>> for (i = 0; i <= n; i++) {
>> for (j = i - 1; j > 0; j--) {
>
>ココで j が -1 になるので、
>work[-1] += work[-2]; // ぶっとびー

for の条件式 j > 0 によって、負のインデックスでのアクセスはありえません。

回答No.6

επιστημη（@episteme）

> それを以下のようにコーディングします。

しくしく、吹っ飛びます。

> for (i = 0; i <= n; i++) {
> for (j = i - 1; j > 0; j--) {

ココで j が -1 になるので、
work[-1] += work[-2]; // ぶっとびー

回答No.5

noname#30727

パスカルの三角形の再帰を展開してみたらどうでしょうか？
配列内の値がどのように変化するかを 7Cr まで書きます。
0Cr: 1
1Cr: 1 1
2Cr: 1 2 1
3Cr: 1 3 3 1
4Cr: 1 4 6 4 1
5Cr: 1 5 10 10 5 1
6Cr: 1 6 15 20 15 6 1
7Cr: 1 7 21 35 35 21 7 1

それを以下のようにコーディングします。
ただし、64C32 だと 32bit の long ではオーバーフローするので、n の範囲を制限するか、より大きな型を使用してください。

long combination(int n, int r)
{
long work[65]; /* n は 64 まで */
int i, j;

for (i = 0; i <= n; i++) {
for (j = i - 1; j > 0; j--) {
work[j] += work[j - 1];
}
work[i] = 1;
}
return work[r];
}

お礼

大変ありがとうございます。やはり、誤差がでないのはパスカルの三角形ですね。ほかにいろいろかんがえてみます。大変参考になります。また、何かあればよろしくお願いします。

回答No.4

nitscape

約分だけでなく#3で言われているように素因数分解もつけてみました。しかしやはり桁あふれをおこします。

除算を複数回行うとどうしても誤差が生じてしまうので、誤差なしで求めるには無限桁計算のルーチンが必要になりますね。ということで乗算部分までは作ってみました。あとは無限桁の除算ルーチンを加えれば完成（？）です。
ただ無限桁の除算ルーチンは今まで作れたことがないので私にはできません。

void Soinsuu(int nData,CDWordArray& adwElement)
{
int i;
int j;
int n;
bool b;

if(nData == 1)
return;

b = false;
n = int(sqrt(nData) + 1);

for(i = 2; i <= n; i += 2)
{
if(i == 4) //2で１回、次からは2n+1で処理を
i--; //行なうように4のときに1を減算

j = 0;
while(1)
{
if(nData % i == 0)
{
j++;
nData /= i;
adwElement.Add(i);
b = true;
}
else
{
break;
}
}
}
if(b == false)
adwElement.Add(nData);

return;
}

void exee(CDWordArray& adwData1,CDWordArray& adwData2)
{
int i;
int j;

for(i = 0; i < adwData1.GetSize(); i++)
{
for(j = 0; j < adwData2.GetSize(); j++)
{
if((adwData2[j] % adwData1[i]) == 0)
{
adwData2[j] /= adwData1[i];
adwData1[i] = 1;
break;
}
}
}
}

//adwData は0～9999まで！
void Mul(CDWordArray& adwResult,CDWordArray& adwData)
{
int i;
int j;
DWORD dwUp;

adwResult.RemoveAll();
adwResult.Add(1);

dwUp = 0;
for(i = 0; i < adwData.GetSize(); i++)
{
if(adwData[i] != 1)
{
for(j = 0; j < adwResult.GetSize(); j++)
{
adwResult[j] *= adwData[i];
adwResult[j] += dwUp;
dwUp = 0;
if(adwResult[j] >= 10000)
{
dwUp = adwResult[j] / 10000;
adwResult[j] %= 10000;
}
}
if(dwUp > 0)
{
adwResult.Add(dwUp);
dwUp = 0;
}
}
}
for(i = adwResult.GetSize() - 1; i >= 0; i--)
{
TRACE("%d ",adwResult[i]);
}
TRACE("\n");
}

void main()
{
int i;
int n;
int r;
DWORD dwResult;
DWORD dwBuff;
CDWordArray adwUp;
CDWordArray adwDown;
CDWordArray adwUp2;
CDWordArray adwDown2;

n = 64;
r = 32;

for(i = n - r + 1; i <= n; i++)
adwUp.Add(i); //　n(n-1)(n-2)...(n-r+1)

for(i = 1; i <= r; i++)
adwDown.Add(i); //　r！

exee(adwUp,adwDown); //素因数分解する前に簡単な約分

for(i = 0; i < adwUp.GetSize(); i++)
Soinsuu(adwUp[i],adwUp2); //分子の素因数分解
for(i = 0; i < adwDown.GetSize(); i++)
Soinsuu(adwDown[i],adwDown2); //分母の素因数分解

exee(adwUp2,adwDown2); //素因数分解してから約分

dwResult = 1;
for(i = 0; i < adwUp2.GetSize(); i++)
{
if(adwUp2[i] != 1)
{
dwResult *= adwUp2[i];
TRACE("%d * ",adwUp2[i]);
}
}
TRACE("\n");

dwBuff = 1;
for(i = 0; i < adwDown2.GetSize(); i++)
{
if(adwDown2[i] != 1)
{
dwBuff *= adwDown2[i];
TRACE("%d * ",adwDown2[i]);
}
}
TRACE("\n");
dwResult /= dwBuff; //求め終わり．．．but、桁あふれあり
TRACE("%dC%d = %d\n",n,r,dwResult);

adwUp.RemoveAll();
Mul(adwUp,adwUp2); //分子の無限桁計算（桁あふれなし）

adwDown.RemoveAll();
Mul(adwDown,adwDown2); //分母の無限桁計算（桁あふれなし）

//多桁でadwUp/adwDownの計算ができれば誤差なしで求められる
}

お礼

僕が最後に考えたのは多倍長演算です。最終的には実験に1024C512が要りそうなので、多倍長演算で計算を実行してみます。

回答No.3

επιστημη（@episteme）

nCr = n! / ((n-r)! * r!)
だから、分子 n! と 分母 (n-r)! : r! それぞれを
素因数分解し、分子,分母から共通項を消去すれば、
かなりイケそうな気がします。

# 気がするだけ^^; すんません。

お礼

ちょっと、やってみます。ありがとうございます。

回答No.2

noname#11476

一つの方法は、tera242さんのdoubleを使うルーチンで、最後に丸め処理( c + 0.5 してから整数化する＝四捨五入）方法です。
それだと丸め誤差は出てこないと予想されます。
ただ、計算の方法として、

a = n;
b = r;
c = 1.0;
for(i=0 ;i< r ;i++){
c= c/b*a;
a--;
b--;
}
return (long)(c + 0.5);

としたほうが賢明かもしれません。（小数部を有効に使うことにより、cの数が大きくなりすぎるのを防ぐ）
（ただ返す数字はlongでdoubleを考えると多分こんなことはしなくても大丈夫と思いますが）

他の方法は、ln(n!) を計算する、factln(n) (精度double)を用意して、

c = factln(n) - factln(n-r) - factln(r)
return (long)(c + 0.5);

とするやり方も考えられます。

参考になりそうな数学的な話は、
http://mathworld.wolfram.com/BinomialCoefficient.html
http://mathworld.wolfram.com/Factorial.html
で、あとは上記の参考文献などを見ればもっと効率の良い方法はありそうです。
が、、、詳細は見たことが無いので分かりません＿ｏ＿

どれも確かめたわけではありませんので、自信なしとしておきます。

ご参考になれば。

お礼

四捨五入は大変参考になりました。それを使うことによって他のプログラミングに生かせそうです。

回答No.1

nitscape

CDWordArrayはDWORD型の自由に配列サイズを変えられる配列です。そこにどんどん乗算すべき値を代入しています。dwUpは分子、dwDownは分母としています。数値が大きくなりすぎないように約分してから最後に一回だけ除算を行っています。あとTRACEはテスト用に入れたものでメッセージを表示するのに使いました。意味はありません。
適当に作ったので非常に効率は悪いですが、こういう人間的なアルゴリズムもあるという参考程度にどうでしょうか？＆もしかしたらnCrの定義を間違えて計算しているかもしれません。

int i;
int j;
int n;
int r;
DWORD dwResult;
DWORD dwBuff;
CDWordArray dwUp;
CDWordArray dwDown;

n = 64;
r = 32;

for(i = n - r + 1; i <= n; i++)
{
dwUp.Add(i); //　n(n-1)(n-2)...(n-r+1)
}

for(i = 1; i <= r; i++)
{
dwDown.Add(i); //　r！
}

for(i = 0; i < dwDown.GetSize(); i++)
{
for(j = 0; j < dwUp.GetSize(); j++)
{
if((dwUp[j] % dwDown[i]) == 0)
{
dwUp[j] /= dwDown[i];
dwDown[i] = 1;
break;
}
}
}

dwResult = 1;
for(i = 0; i < dwUp.GetSize(); i++)
{
dwResult *= dwUp[i];
TRACE("%d * ",dwUp[i]);
}
TRACE("\n");

dwBuff = 1;
for(i = 0; i < dwDown.GetSize(); i++)
{
dwBuff *= dwDown[i];
TRACE("%d * ",dwDown[i]);
}
TRACE("\n");

dwResult /= dwBuff;

TRACE("%dC%d = %d\n",n,r,dwResult);

お礼

回答大変ありがとう。ございます。参考にさせていただきます。