Воспользуйтесь соотношениями:
C(k, k) = 1
C(n, k) = n * C(n - 1, k) / (n - k)
Получается что какое количество дробей ни возьми произведение всегда целое:
k + 1 k + 2 k + 3 k + i
----- * ----- * ----- * ... * ----- * ...
1 2 3 i
Что приводит к такой процедуре:
unsigned long long cnk(unsigned long long n, unsigned long long k) {
unsigned long long p = 1;
for (unsigned long long i = 1; i <= n - k; ++i) {
p = p * (k + i) / i;
}
return p;
}
Переполнения не будет пока n <= 62
. Это решение можно улучшить до оптимального: если результат помещается в разрядную сетку мы его вычислим.
Проблема кода выше - это переполнение произведения p * (k + i)
. Избавимся от него. Пусть известно что a * b / c
- целое, тогда
a * b / c = (a / НОД(a, c)) * (b / (c / НОД(a, c)))
В правой части все частные целые и нет переполнения (мы предполагаем что результат левой части в тип помещается).
Получается такое решение:
unsigned long long gcd(unsigned long long a, unsigned long long b) {
while (b > 0) {
unsigned long long c = a % b;
a = b;
b = c;
}
return a;
}
// a * b / c
unsigned long long mul_div(unsigned long long a, unsigned long long b, unsigned long long c) {
unsigned long long ac = gcd(a, c);
return (a / ac) * (b / (c / ac));
}
unsigned long long cnk(unsigned long long n, unsigned long long k) {
unsigned long long p = 1;
for (unsigned long long i = 1; i <= n - k; ++i) {
p = mul_div(p, i + k, i);
}
return p;
}
Оно может вычислить Cnk
для любых комбинаций n
и k
для которых результат помещается в разрядную сетку. Например любые n <= 67
:
cnk(67, 33) = 14226520737620288370
Последняя вещь которую нужно сделать - уменьшить число итераций если n
очень велико, а k
мало. В таком случае будем заменять k
на n - k
:
unsigned long long cnk(unsigned long long n, unsigned long long k) {
if (k < n - k) {
k = n - k;
}
unsigned long long p = 1;
for (unsigned long long i = 1; i <= n - k; ++i) {
p = mul_div(p, i + k, i);
}
return p;
}
Программа максимум достигнута - если результат помещается в 64 бита он будет вычислен.
Максимальные n
для небольших k
:
cnk(18446744073709551615, 0) = 1
cnk(18446744073709551615, 1) = 18446744073709551615
cnk( 6074001000, 2) = 18446744070963499500
cnk( 4801280, 3) = 18446738006366306560
cnk( 145056, 4) = 18446483332847246040
cnk( 18580, 5) = 18442234518422931216
cnk( 4868, 6) = 18426048264474788224
cnk( 1913, 7) = 18399302838933135756
cnk( 967, 8) = 18419736117819661560
cnk( 577, 9) = 18348006354228436600
cnk( 386, 10) = 17991165343481265936