De Andreas' answer:
Aquí está un antiguo algoritmo que es exacta y no se desborde a menos que el resultado es demasiado grande para un long long
unsigned long long
choose(unsigned long long n, unsigned long long k) {
if (k > n) {
return 0;
}
unsigned long long r = 1;
for (unsigned long long d = 1; d <= k; ++d) {
r *= n--;
r /= d;
}
return r;
}
Este algoritmo está también en Knuth "El arte of Computer Programming, 3rd Edition, Volume 2: Algoritmos Seminumerical "Creo.
ACTUALIZACIÓN: Hay una pequeña posibilidad de que el algoritmo se desbordará en la línea:
r *= n--;
para muy n grande. Un límite superior ingenuo es sqrt(std::numeric_limits<long long>::max()) lo que significa un n menos de 4,000,000,000.
Considere n == 67 yk == 33. El algoritmo anterior se desborda con un 64 bit sin signo de larga duración. Y, sin embargo, la respuesta correcta es representable en 64 bits: 14,226,520,737,620,288,370. Y el algoritmo anterior no dice nada acerca de su desbordamiento, seleccione (67, 33) devuelve:
8.829.174.638.479.413
Una respuesta creíble pero incorrecta.
Sin embargo, el algoritmo anterior se puede modificar ligeramente para que no se desborde nunca, siempre que la respuesta final sea representable.
El truco está en reconocer que en cada iteración, la división r/d es exacta. Temporalmente la reescritura:
r = r * n/d;
--n;
Para que esto sea exacta, que significa que si se expandió R, N y D en sus factores primos, entonces uno podría fácilmente anular d, y se quedará con un valor modificado para n, llamada que t, y luego el cálculo de r es simplemente:
// compute t from r, n and d
r = r * t;
--n;
una manera rápida y fácil de hacer esto es encontrar el máximo común divisor de I + D, lo llaman g:
unsigned long long g = gcd(r, d);
// now one can divide both r and d by g without truncation
r /= g;
unsigned long long d_temp = d/g;
--n;
Ahora podemos hacer la s ame cosa con d_temp yn (encuentre el divisor común más grande). Sin embargo, como sabemos a priori que r * n/d es exacto, también sabemos que gcd (d_temp, n) == d_temp, y por lo tanto no es necesario que lo calculemos. Así podemos dividir n por d_temp:
unsigned long long g = gcd(r, d);
// now one can divide both r and d by g without truncation
r /= g;
unsigned long long d_temp = d/g;
// now one can divide n by d/g without truncation
unsigned long long t = n/d_temp;
r = r * t;
--n;
Limpieza:
unsigned long long
gcd(unsigned long long x, unsigned long long y)
{
while (y != 0)
{
unsigned long long t = x % y;
x = y;
y = t;
}
return x;
}
unsigned long long
choose(unsigned long long n, unsigned long long k)
{
if (k > n)
throw std::invalid_argument("invalid argument in choose");
unsigned long long r = 1;
for (unsigned long long d = 1; d <= k; ++d, --n)
{
unsigned long long g = gcd(r, d);
r /= g;
unsigned long long t = n/(d/g);
if (r > std::numeric_limits<unsigned long long>::max()/t)
throw std::overflow_error("overflow in choose");
r *= t;
}
return r;
}
Ahora se puede calcular elegir (67, 33) sin rebosadero. Y si intentas elegir (68, 33), obtendrás una excepción en lugar de una respuesta incorrecta.
para el uso grande int véase:
[http://stackoverflow.com/questions/1188939/representing-128-bit-numbers-in-c ] (http://stackoverflow.com/questions/1188939/representing-128-bit-numeros-in-c)
[http://stackoverflow.com/questions/1055661/bigint-bigbit-library](http: //stackoverflow.com/questions/1055661/bigint-bigbit-library)
[http://stackoverflow.com/questions/238343/big-number-in-c](http://stackoverflow. com/questions/238343/big-number-in-c)
– SjB