組合せ $_n C _r$ を求める関数がC++には無いようなので、自分ように作ってみた。
Qiitaへの投稿、Markdownの記述が初めてなので、誤りが含まれているかもしれません。
目標
$n=50$ で、最大になる $r = 25$ のときの $_{50} C _{25} = 126,410,606,437,752 $ をオーバーフローなしで返す関数の作成。
作成例1
組み合わせの下の式を使う。
_n C _r \ (n \geq r \geq 0) = \frac{_n P _r}{r!}
= \frac{m \times (m - 1) \times \cdots \times (n - r + 1)}{n \times (n - 1) \times \cdots \times 2 \times 1}
- 分子を計算。
- 分母を計算、分母と分子が 2 か 3 の倍数ならxの倍数のxで分母と分子を割る。
- 分子 ÷ 分母 を返す。
ソースコード
comb1.cpp
using ull = unsigned long long;
ull comb(ull Left, ull Right) {
ull top = 1, bottom = 1, cnt = 0;
while(cnt < Right) {
top *= Left - cnt;
cnt++;
}
while(Right > 0) {
bottom *= Right;
if(Left % 2 == 0 && Right % 2 == 0) {
top /= 2;
bottom /= 2;
}
if(Left % 3 == 0 && Right % 3 == 0) {
top /= 3;
bottom /= 3;
}
Right--;
}
return top / bottom;
}
void output_comb(ull a, ull b) {
cout << "comb(" << a << ", " << b << ")"
<< " = " << comb(a, b) << "\n\n";
}
int main() {
ull n, r;
n = 4;
r = 2;
output_comb(n, r);
n = 50;
r = 10;
output_comb(n, r);
}
実行結果
result
comb(4, 2) = 6
comb(50, 10) = 10272278170
問題点
$n = 50, \ r = 12$ 以降で $_n P _r$ がオーバーフローする。
comb1.cpp
int main() {
ull n, r;
n = 50;
r = 12;
output_comb(n, r);
}
result
comb(50, 12) = 5867193126
本来は $121,399,651,100$ となる必要がある。
作成例 2
改善点
- 分母・分子を同時に計算する。
- 配列に素数を 100 個登録しておき、一回掛けるるごとに分母・分子の公約数で割る。
- $_n C _r \ = \ _n C _{n-r}$ を用いて
n / 2 < rのときr = n - rにする。
ソースコード
comb2.cpp
using ull = unsigned long long;
ull comb(ull Left, ull Right) {
if(Left / 2 < Right)Right = Left - Right;
ull top = 1, bottom = 1, cnt = 0;
ull prime[100] = { 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,
73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,
179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,
283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,
419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541 };
while(cnt < Right) {
// 分母に n-cnt 分子に r-cnt を掛ける。
top *= Left - cnt;
bottom *= Right - cnt;
// 公約数で割る
for(int i = 0; i < 100; i++) {
if(top % prime[i] == 0 && bottom % prime[i] == 0) {
top /= prime[i];
bottom /= prime[i];
}
}
cnt++;
}
return top / bottom;
}
void output_comb(ull a, ull b) {
cout << "comb(" << a << ", " << b << ")"
<< " = " << comb(a, b) << "\n";
}
int main() {
ull n, r;
n = 50;
r = 11;
while(r <= 39) {
output_comb(n, r);
r++;
}
}
実行結果
result
comb(50, 11) = 37353738800
comb(50, 12) = 121399651100
comb(50, 13) = 354860518600
comb(50, 14) = 937845656300
comb(50, 15) = 2250829575120
comb(50, 16) = 4923689695575
comb(50, 17) = 9847379391150
comb(50, 18) = 18053528883775
comb(50, 19) = 30405943383200
comb(50, 20) = 47129212243960
comb(50, 21) = 67327446062800
comb(50, 22) = 88749815264600
comb(50, 23) = 108043253365600
comb(50, 24) = 121548660036300
comb(50, 25) = 126410606437752
comb(50, 26) = 121548660036300
comb(50, 27) = 108043253365600
comb(50, 28) = 88749815264600
comb(50, 29) = 67327446062800
comb(50, 30) = 47129212243960
comb(50, 31) = 30405943383200
comb(50, 32) = 18053528883775
comb(50, 33) = 9847379391150
comb(50, 34) = 4923689695575
comb(50, 35) = 2250829575120
comb(50, 36) = 937845656300
comb(50, 37) = 354860518600
comb(50, 38) = 121399651100
comb(50, 39) = 37353738800
目標の $_{50} C _{25} = 126410606437752$ まで返すことができた。
$n=51$ 以降の数では prime[99] == 541 以上の素数での除算漏れが原因のオーバーフローが発生する可能性がある。
おわりに
最後まで読んでいただきありがとうございました。
これよりも高速かつ効率の良い方法は大量にあるとおもうので、改善案や記事のミスがあればコメントで教えていただけると幸いです。
参考にしたサイト・記事
順列・組合せ - ke!san
コンビネーション（mCn）の計算方法 - 具体例で学ぶ数学
1000までの階乗表
素数生成機