コンビネーション,二項係数の求め方

n≒1000くらいまでコンビネーションnCkを計算できるプログラムを作ろうと思っています。
階乗を使った公式では直ぐに破綻してしまうので
7C3=(7・6・5)/(3・2・1)とやるようなプログラムを組んだのですが
希望より小さなnで破綻してしまいます。
とりあえず、今は7C3＝(7/3)・(6/2)・(5/1)とやるような計算法で凌いでいますが
途中で実数計算(整数計算でないという意味)をせざるを得ず、ちょっと気持ち悪いです。
究極のプログラムを組もうという気は無いのですが
もう少し現状を改善したいと思っています。

良きアドバイスをいただけたら幸いです。

No.8 ベストアンサー 回答者： Ishiwara 回答日時： 2007/11/27 12:26

確かに実数（小数点以下）が入ってくるのは気分が悪いですよね。

一番いいのは、分子・分母別に「出てくるたびに素因数分解」してしまい、配列に格納することではないでしょうか。

分子をa(素因数)、分母をb(素因数)とし、
7/3 が出たら a(7)=a(7)+1, b(3)=b(3)+1 とし
次に 6/2 がでたら、6=2x3 と分解して、
a(2)=a(2)+1, a(3)=a(3)+1, b(2)=b(2)+1 とします、そして直ちに a(2),b(2)が「両者ゼロ」でないことをとがめて相殺（約分）してしまいます。

このようにして行けば、Nが1000ぐらい楽勝と思います。

No.11 回答者： Ishiwara 回答日時： 2007/11/30 12:49

＃８,９,１０です。

「分子なら＋１、分母なら－１」
すばらしいアイディアじゃないですか。もうできたも同然でしょう。がんばってください。

No.10 回答者： Ishiwara 回答日時： 2007/11/28 13:13

＃８,９です。

再帰型のプログラムを作って実験しましたが、せいぜいｍ＜２０ぐらいまでしか使いものになりませんでした。理由は、１つのFunctionがそれぞれ２つのFunctionを呼び出すので「下請け」の回数が膨大になってしまうことです。
ご自分でプログラムを新規に構築されるなら、やはり＃８の方法がお勧めだと思います。

この回答への補足

何度も回答ありがとうございました。
No.10で書かれた通り、再帰タイプは組むのは簡単ですが、実行は大変のようです。

n=1000にもなると何百桁にもなることから
通常の整数型変数では対処しきれないので
No.8で教えていただいたように素因数分解を用いる方法で挑戦したいと思います。
ただ思ったのですが、分子も分母も同じ配列で良いのではないでしょうか?
分母に現れたら-1して、分子に現れたら+1すればいいわけなので…。
(組む人の好みに依ると思いますが…)

また、No.7で教えていただいたことを組み合わせると符号無し変数も使えそうです。
(EXCEL VBAでは実装していませんが)

問題は関数から抜ける時に、どうやって値を渡すかが次の問題になりそうですが、
とりあえず行けるところまで行ってみたいと思います。
ありがとうございました。

補足日時：2007/11/29 00:55

No.9 回答者： Ishiwara 回答日時： 2007/11/28 01:18

＃８です。

システムが再帰呼び出しをサポートしていれば、それが一番簡単です。
5C2=4C2+4C1 ですから、5C2を自分で計算せず、仕様書を渡すだけで、２人の下請けに4C2と4C1を計算させます。どの下請けも、なるべくさらに下請けに出そうとします。ただし「発注者」と「受注者」は、実は全部同一人物（同一プログラム）であって、自分で仕様書を書いて自分で受け取り、結果が出たら自分に納品する、という仕組みです。つまり、プログラマーは、mCn=(m-1)Cn+(m-1)C(n-1)というプログラムを１個書くだけでおしまいなのです。ただし、そのプログラムの中で、もしkC0 やkCkを受注したら、もう下に出せないので、値1を上に返します。
ちょっとややこしい感じがしますが、実体は、ほんの数行です。一度作って味を占めると、次回からはすごく簡単です。ぜひお試しください。

No.7 回答者： kmvblxc 回答日時： 2007/11/25 00:16

>とりあえず、今は7C3＝(7/3)・(6/2)・(5/1)とやるような計算法で凌いでいますが

>途中で実数計算(整数計算でないという意味)をせざるを得ず、ちょっと気持ち悪いです。

次のアルゴリズムを使えば、途中で整数計算でない計算をしないですむ。
以下は、フリーソフト Risa/Asir での計算結果。
comb(10000,3000)の計算を一瞬でやってくれる。


[0] def comb(A,B){for(I=1, C=1; I<=B; I++) C *= (A-I+1)/I; return C;}$
[1] [2] comb(10000,3000);
77578708414622781919184089365853843098672731962523453409968751787633293257691437

77964834277963683734825766127802024646897842450448606376443199326481209988685564

28053665674682245456870011002895215023109478980481892649071494308290750871267577

60970600128136307191477472836543587696792303287128651789677102700955487241433699

25219040594723598959713315458952496604157370703092408567053250339982272723332954

95090803099389306787186346799897659251373046947323832760862334897746475792126428

85578743907836354019701783883048023872206622806636352045276069616724975579837711

77536628952968902911778940157712556142108685234016571884302818227186854838921026

56238550371710479704719161489083941430788365654578114839371386886584288433093076

48201538869279089975997921408491177067667611792575679011657080668803280848189309

51557800873241241913727923481130817780397935214520329518215582779431369109297672

97917309857981118602214902115337340123842561040687972572992706363241895569394746

29411781084725425913515068172351801315118208125402060325840725900500604082184865

35106084042110995262883556887775255353024188124937610410154541587073021928598892

58219479036015569951701536201695291500782833561272918912889102668085475897203941

05911344257379277304549575049354285576843730220694295275435485251033346293408317

92791337167735991209574134284900979037388725415353002644086675253706300756958689

01491155314398185889794446701668340271033348025086189433897878546599910935246670

93991139230457632005685109254296862862826621766815964756364190444633813355612593

36735740136912220940851507605550235053600960655904111979726170852945765906231731

71523179442047734527549751878920299318799059678333466774145141112529093977474288

26679839888112128175195018049063902123230517964865415017075715173083126695073092

12926852059136998254713532203720878594474642442600068850530397063561306733997993

99531720578192075968428145233725903775783028769106314532758299198525616678156266

90325720395145676223982087510693989808055283808938171380053608991466724348163261

76042231138572997365314617417029504686844779359477617159502969359631317307660994

04108941955509732724735406216005215402013744496820941747633265342650204486695473

04892857748211091972028737309834508503333952623674119987540629556228390883273717

15611479953835617054258066986050637436743293286345363059606749473161357960484850

45627870210064752650297291710729496150779459221910583125069854687964764841900837

17582220249017921906106011522284298198930211499275805765495010980922961271399465

49285526172579286638223966962211318582738278666968096448819053381901740491909522

27635744996678531849775951330445119253429762987839242945447998153030655886592465

95943283200

この回答への補足

回答ありがとうございます。
＞[0] def comb(A,B){for(I=1, C=1; I<=B; I++) C *= (A-I+1)/I; return C;}$

この順に行っていけば、分母の因数が、
それ以前に必ず分子に含まれるわけですね。
こんなに簡潔な方法で“気持ち悪さ”が除去できるとは…。
感動しました。ありがとうございました。

補足日時：2007/11/27 23:08

No.6 回答者： kabaokaba 回答日時： 2007/11/24 23:34

n=1000の二項係数ですか・・・・

工夫をしないと，まともな時間で計算は不可能です．
それと，定義にそのまま当てはめると
まず間違いなく誤差が生じます．

>組み込み関数として導入したいのです。
エクセルが持ってるなら，もう組み込み関数なのでは？
VBAから呼び出せばよいのでしょう．
ただし，桁が大きいので，マイクロソフトがどのように
実装してるか，どう最適化してるのかによっては使い物にならないかも．

で，比較的単純かつ簡単な実装方法．
この手の組合せの数を計算させるのには
（1）再帰
（2）キャッシュ
を使うのが，一番簡単な定石でしょう．

nCk を求める関数を Comb(n,k)とすると
Comb(n,k) = Comb(n-1,k)+Comb(n-1,k-1)
任意のkに対して，Comb(0,k)=Comb(k,k)=1
であることを利用して，計算させるのが再帰．
これは No.2 さんのいってる通り

けど，これだけじゃ実際は計算が遅くてやってられません．
同じ計算を何度も繰り返すことになるのです．
そこで，一度計算した値は二度と計算しないですむように
記録させておきます．これがキャッシュ

昔つくったPerlのスクリプトを晒しておきます．
VBAとはかなり違う文法ですが，厄介なことはしてないので
アルゴリズムそのものを追うことは比較的容易だと思います．
要は，再帰で計算させるときに
その計算結果がすでに分かってるならばそれを流用し，
分かってないときだけ計算させ，さらにその結果を
記録させておくということをしているだけです．

use strict;
use warnings;
no warnings "recursion";
use bigint;

{my %cache;
sub comb_cache{
my ($n, $r)=@_;
$cache{"$n,$r"} = 0, return 0 if $r<0 || $n<$r ;
$cache{"$n,$r"} = 1, return 1 if $r==0 && $n==0;
$cache{"$n,$r"} = comb_cache($n-1,$r)+comb_cache($n-1,$r-1)
unless $cache{"$n,$r"};
return $cache{"$n,$r"}
}
}

#Comb(1000,200)の計算例
print comb_cache(1000,200);

_END_

6617155560659303656271633461324588318973
2170301763866936478813470889179595672641
1057801285583913163781806953211915554723
3739314518470598302521758877124573075476
4935413546061929638388295789716188963628
0577155889117185

216桁，当然，No.3さんと一致します．
私の機械では上のわずかなプログラムで
約78.7秒でこの答えが出てきます．
ちなみにこの関数は322001回呼び出されてました．

No.5 回答者： kumipapa 回答日時： 2007/11/24 21:32

＞「recursive call」というのがよくわからないのですが

＞　自分自身の関数を呼ぶ、という意味でしょうか?
＞　それは、どのような言語で可能ですか?
＞　EXCELのVBAで可能でしょうか?
はい、自分自身の関数を呼ぶという再帰呼び出しと呼ばれるもので、EXCELのVBAでも可能なようです。数値計算ではあまり利用されることは無いのかもしれませんが、ループで記述するとやっかいな論理問題などをキレイに記述できるので、覚えておかれると便利でしょう。
・・・何せ大昔にやってた話なので、懐かしい限りです。

この回答への補足

(フォローも含め)回答ありがとうございます。
EXCELのVBAでも動作しました。
50C10位で、いかんともしがたくなってしまいました。

とは言え
「recursive call」という言葉も教えていただき勉強になりました。
ありがとうございました。

補足日時：2007/11/27 22:59

No.4 回答者： kumipapa 回答日時： 2007/11/24 21:17

#2です。

ごめんなさい、細かな話ですが、r = 0 のときが抜けてました。
nＣ0 = 1
でしたね。一応念のため。

No.3 回答者： info22 回答日時： 2007/11/24 20:20

nCr=n!/{r!*(n-r)!}はrについて対称ですからrについて半分だけ計算すれば良いですね。

nの偶数、奇数でnCrの項数がそれぞれ奇数、偶数となることを考慮にいれてプログラムを組めばいいですね。
計算機で計算できる有効桁数が限られている場合は
＞今は7C3＝(7/3)・(6/2)・(5/1)とやるような計算法で凌いでいます
の計算法は有効ですね。

でも有効桁数が限られている場合は限界があります。
１）出来るだけ有効桁数の多い計算機を使うか、
２）有効桁数に限度のない数式処理ソフトで計算する
ことが考えられます。

２）はMathmaticaやMapleがありますがいずれも有料ソフトで、多くの大学や高専で導入されていて、学生は無料で使えます。また学生用価格でも低価格でも提供されています。このようなソフトを使えば数式通りの整数型の計算でオーバーフローすることもなく計算できます。1000!/{200!*800!}の計算も瞬時に出てきます。↓
6617155560659303656271633461324588318973217030176386693647881347088
9179595672641105780128558391316378180695321191555472337393145184705
9830252175887712457307547649354135460619296383882957897161889636280
577155889117185
プログラムも組めますが、桁数が多くなりますので大変です。

１）の場合でもフリーのソフトのMaxima(Win PC版もある)などを使えば　
69!までの階乗に相当する整数まで桁落ちしないで扱えます。
例えば
300!/{200!*100!}
の計算が正確に計算できます。
逆に
{200!*100!}/300!
も正確に上記の逆数の分数で計算結果が出てきます。

以下のどこかのサイトからWindows版をダウンロードしてインストールすればMaximaが使えます。
http://www.google.co.jp/search?num=100&hl=ja&q=M …

この回答への補足

回答ありがとうございます。

とりあえず今はEXCELのVBAで考えています。
EXCELの関数にコンビネーションがあることは知っていますが
コンビネーションを求めること自体が目的ではなく
組み込み関数として導入したいのです。

(2)につきましては関数電卓やEXCELがありますし
maximaも使っております。
＞69!までの階乗
今のところ、目安はn=1000ですので…。

補足日時：2007/11/24 21:09

No.2 回答者： kumipapa 回答日時： 2007/11/24 20:00

専門家ではないのでおかしな答えかも知れませんが、recursive call が許されるなら、それを利用する手もあるのではないでしょうか。

nＣr = (n-1)Ｃr + (n-1)Ｃ(r-1)
nＣ1 = n
nＣn = 1
を利用します。

プログラミング自体は良く知らないので細かなところは見逃して頂くとして、

combination(n,r)
{
　　if ( n < r )
　　　　return(0);
　　else if ( n == r )
　　　 return(1);
　　else　if ( r == 1 )
　　　　return(n);
　　else {
　　　　return(combination(n-1,r)+combination(n-1,r-1));
　　}
}
のようにしますと、割り算が入らず整数計算だけでＯＫですし、掛算もなく、途中の計算値が求める値を越えることはないでしょう。rが適当に小さくなったところで計算してしまっても良いのかも知れません。
究極のプログラムを求めているわけではないとのことですので、気軽に回答してみました。

この回答への補足

回答ありがとうございます。
補足質問させてください。

「recursive call」というのがよくわからないのですが
自分自身の関数を呼ぶ、という意味でしょうか?
それは、どのような言語で可能ですか?
EXCELのVBAで可能でしょうか?

補足日時：2007/11/24 21:13