algorithm - teorema - La forma más rápida de generar coeficientes binomiales
teorema del binomio (10)
Necesito calcular combinaciones para un número.
¿Cuál es la forma más rápida de calcular nCp donde n >> p?
Necesito una forma rápida de generar coeficientes binomiales para una ecuación polinómica y necesito obtener el coeficiente de todos los términos y almacenarlo en una matriz.
(a + b) ^ n = a ^ n + nC1 a ^ (n-1) * b + nC2 a ^ (n-2) * ............ + nC (n-1) ) a * b ^ (n-1) + b ^ n
¿Cuál es la forma más eficiente de calcular nCp?
Esta es mi versión:
def binomial(n, k):
if k == 0:
return 1
elif 2*k > n:
return binomial(n,n-k)
else:
e = n-k+1
for i in range(2,k+1):
e *= (n-k+i)
e /= i
return e
Estaba buscando lo mismo y no pude encontrarlo, por lo que escribí uno que parece óptimo para cualquier Coeficiente Biólogo para el cual el resultado final se ajusta a un Largo.
// Calculate Binomial Coefficient
// Jeroen B.P. Vuurens
public static long binomialCoefficient(int n, int k) {
// take the lowest possible k to reduce computing using: n over k = n over (n-k)
k = java.lang.Math.min( k, n - k );
// holds the high number: fi. (1000 over 990) holds 991..1000
long highnumber[] = new long[k];
for (int i = 0; i < k; i++)
highnumber[i] = n - i; // the high number first order is important
// holds the dividers: fi. (1000 over 990) holds 2..10
int dividers[] = new int[k - 1];
for (int i = 0; i < k - 1; i++)
dividers[i] = k - i;
// for every dividers there is always exists a highnumber that can be divided by
// this, the number of highnumbers being a sequence that equals the number of
// dividers. Thus, the only trick needed is to divide in reverse order, so
// divide the highest divider first trying it on the highest highnumber first.
// That way you do not need to do any tricks with primes.
for (int divider: dividers) {
boolean eliminated = false;
for (int i = 0; i < k; i++) {
if (highnumber[i] % divider == 0) {
highnumber[i] /= divider;
eliminated = true;
break;
}
}
if(!eliminated) throw new Error(n+","+k+" divider="+divider);
}
// multiply remainder of highnumbers
long result = 1;
for (long high : highnumber)
result *= high;
return result;
}
Hace poco escribí un código que requería un coeficiente binario de aproximadamente 10 millones de veces. Así que hice una combinación de tabla de búsqueda / enfoque de cálculo que aún no es demasiado desperdicio de memoria. Puede que le resulte útil (y mi código es de dominio público). El código está en
http://www.etceterology.com/fast-binomial-coefficients
Se sugirió que incluya el código aquí. Una gran tabla de búsqueda de bocinazos parece una pérdida, así que aquí está la función final, y una secuencia de comandos de Python que genera la tabla:
extern long long bctable[]; /* See below */
long long binomial(int n, int k) {
int i;
long long b;
assert(n >= 0 && k >= 0);
if (0 == k || n == k) return 1LL;
if (k > n) return 0LL;
if (k > (n - k)) k = n - k;
if (1 == k) return (long long)n;
if (n <= 54 && k <= 54) {
return bctable[(((n - 3) * (n - 3)) >> 2) + (k - 2)];
}
/* Last resort: actually calculate */
b = 1LL;
for (i = 1; i <= k; ++i) {
b *= (n - (k - i));
if (b < 0) return -1LL; /* Overflow */
b /= i;
}
return b;
}
#!/usr/bin/env python3
import sys
class App(object):
def __init__(self, max):
self.table = [[0 for k in range(max + 1)] for n in range(max + 1)]
self.max = max
def build(self):
for n in range(self.max + 1):
for k in range(self.max + 1):
if k == 0: b = 1
elif k > n: b = 0
elif k == n: b = 1
elif k == 1: b = n
elif k > n-k: b = self.table[n][n-k]
else:
b = self.table[n-1][k] + self.table[n-1][k-1]
self.table[n][k] = b
def output(self, val):
if val > 2**63: val = -1
text = " {0}LL,".format(val)
if self.column + len(text) > 76:
print("/n ", end = "")
self.column = 3
print(text, end = "")
self.column += len(text)
def dump(self):
count = 0
print("long long bctable[] = {", end="");
self.column = 999
for n in range(self.max + 1):
for k in range(self.max + 1):
if n < 4 or k < 2 or k > n-k:
continue
self.output(self.table[n][k])
count += 1
print("/n}}; /* {0} Entries */".format(count));
def run(self):
self.build()
self.dump()
return 0
def main(args):
return App(54).run()
if __name__ == "__main__":
sys.exit(main(sys.argv))
La aproximación razonable más rápida en mi propia evaluación comparativa es la aproximación utilizada por la biblioteca Apache Commons Maths: http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/special/Gamma.html#logGamma(double)
Mis colegas y yo tratamos de ver si podíamos superarlo, usando cálculos exactos en lugar de aproximaciones. Todos los enfoques fallaron miserablemente (muchas órdenes más lentas) excepto uno, que fue 2-3 veces más lento. El mejor enfoque de rendimiento utiliza https://math.stackexchange.com/a/202559/123948 , aquí está el código (en Scala):
var i: Int = 0
var binCoeff: Double = 1
while (i < k) {
binCoeff *= (n - i) / (k - i).toDouble
i += 1
}
binCoeff
Los enfoques realmente malos fueron varios intentos de implementar el Triángulo de Pascal usando recursividad de cola.
Si desea expansiones completas para valores grandes de n, la convolución de FFT podría ser la forma más rápida. En el caso de una expansión binomial con coeficientes iguales (por ejemplo, una serie de lanzamientos de monedas) y un orden par (por ejemplo, número de lanzamientos), puede explotar las simetrías de esta manera:
Teoría
Representa los resultados de dos lanzamientos de monedas (por ejemplo, la mitad de la diferencia entre el número total de caras y cruces) con la expresión A + A * cos (Pi * n / N). N es el número de muestras en su buffer - una expansión binomial de orden par O tendrá coeficientes O + 1 y requerirá un buffer de N> = O / 2 + 1 muestras - n es el número de muestra generado, y A es un factor de escala que generalmente será 2 (para generar coeficientes binomiales) o 0.5 (para generar una distribución de probabilidad binomial).
Observe que, en frecuencia, esta expresión se asemeja a la distribución binomial de esos dos lanzamientos de monedas: hay tres picos simétricos en las posiciones correspondientes al número (caras-colas) / 2. Dado que el modelado de la distribución de probabilidad global de eventos independientes requiere la convolución de sus distribuciones, queremos convolucionar nuestra expresión en el dominio de la frecuencia, que es equivalente a la multiplicación en el dominio del tiempo.
En otras palabras, elevando nuestra expresión de coseno para el resultado de dos lanzamientos a una potencia (por ejemplo, para simular 500 lanzamientos, elevándola a la potencia de 250 dado que ya representa un par), podemos organizar la distribución binomial para una gran número para aparecer en el dominio de la frecuencia. Como todo esto es real e incluso, podemos sustituir el DCT-I por el DFT para mejorar la eficiencia.
Algoritmo
- decidir sobre un tamaño de búfer, N, que es al menos O / 2 + 1 y se puede DCTed convenientemente
- Inicialízalo con la expresión pow (A + A * cos (Pi * n / N), O / 2)
- aplicar el DCT-I hacia adelante
- lea los coeficientes del buffer - el primer número es el pico central donde heads = tails, y las entradas subsiguientes corresponden a pares simétricos sucesivamente más lejos del centro
Exactitud
Hay un límite de cuán alto O puede ser antes de que los errores acumulados de redondeo de punto flotante te roben los valores enteros precisos para los coeficientes, pero supongo que el número es bastante alto. El punto flotante de precisión doble puede representar enteros de 53 bits con total precisión, y voy a ignorar la pérdida de redondeo involucrada en el uso de pow () porque la expresión generadora tendrá lugar en los registros FP, lo que nos da un 11 adicional. bits de mantisa para absorber el error de redondeo en las plataformas Intel. Entonces, suponiendo que usemos un DCT-I de 1024 puntos implementado a través de la FFT, eso significa perder precisión de 10 bits en el error de redondeo durante la transformación y no mucho más, dejándonos con ~ 43 bits de representación limpia. No sé qué orden de expansión binomial genera coeficientes de ese tamaño, pero me atrevo a decir que es lo suficientemente grande para sus necesidades.
Expansiones asimétricas
Si quiere las expansiones asimétricas para coeficientes desiguales de ayb, necesitará usar un DFT de dos lados (complejo) y una función pow () compleja. Genere la expresión A * A * e ^ (- Pi * i * n / N) + A * B + B * B * e ^ (+ Pi * i * n / N) [utilizando la función compleja pow () para aumentar a la potencia de la mitad del orden de expansión] y DFT lo. Lo que tienes en el búfer es, de nuevo, el punto central (pero no el máximo si A y B son muy diferentes) en el desplazamiento cero, y va seguido de la mitad superior de la distribución. La mitad superior del búfer contendrá la mitad inferior de la distribución, que corresponde a los valores de cabezas-menos-colas que son negativas.
Observe que los datos de origen son herméticos simétricos (la segunda mitad del búfer de entrada es el conjugado complejo del primero), por lo que este algoritmo no es óptimo y puede realizarse utilizando una FFT compleja a compleja de la mitad del tamaño requerido para obtener el óptimo eficiencia.
Huelga decir que toda la exponenciación compleja masticará más tiempo de CPU y dañará la precisión en comparación con el algoritmo puramente real para las distribuciones simétricas anteriores.
Si entiendo la notación en la pregunta, no solo quieres nCp, realmente quieres todo nC1, nC2, ... nC (n-1). Si esto es correcto, podemos aprovechar la siguiente relación para hacer esto bastante trivial:
- para todo k> 0: nCk = prod_ {de i = 1..k} ((n-i + 1) / i)
- es decir, para todo k> 0: nCk = nC (k-1) * (n-k + 1) / k
Aquí hay un fragmento de python que implementa este enfoque:
def binomial_coef_seq(n, k):
"""Returns a list of all binomial terms from choose(n,0) up to choose(n,k)"""
b = [1]
for i in range(1,k+1):
b.append(b[-1] * (n-i+1)/i)
return b
Si necesita todos los coeficientes hasta un cierto límite de k> (n / 2), puede usar simetría para reducir el número de operaciones que necesita realizar deteniéndose en el coeficiente de techo (n / 2) y luego simplemente rellenando hasta donde necesitas.
import numpy as np
def binomial_coef_seq2(n, k):
"""Returns a list of all binomial terms from choose(n,0) up to choose(n,k)"""
k2 = int(np.ceiling(n/2))
use_symmetry = k > k2
if use_symmetry:
k = k2
b = [1]
for i in range(1, k+1):
b.append(b[-1] * (n-i+1)/i)
if use_symmetry:
v = k2 - (n-k)
b2 = b[-v:]
b.extend(b2)
return b
Si necesita calcularlos para todo n, la respuesta de Ribtoks es probablemente la mejor. Para una sola n, es mejor que hagas esto:
C[0] = 1
for (int k = 0; k < n; ++ k)
C[k+1] = (C[k] * (n-k)) / (k+1)
La división es exacta, si se hace después de la multiplicación.
Y tenga cuidado con el desbordamiento de C [k] * (nk): use números enteros suficientemente grandes.
Si realmente solo necesitas el caso donde n es mucho más grande que p, un camino a seguir sería usar la fórmula de Stirling para los factoriales. (si n >> 1 yp es orden uno, Stirling se aproxima a n! y (np)!, mantenga p! como es, etc.)
Ustedes usan programación dinámica para generar coeficientes binomiales
Puede crear una matriz y luego usar el loop O (N ^ 2) para llenarlo
C[n, k] = C[n-1, k-1] + C[n-1, k];
dónde
C[1, 1] = C[n, n] = 1
Después de eso, en su programa puede obtener el valor C (n, k) simplemente mirando su matriz 2D en índices [n, k]
ACTUALIZAR así smth
for (int k = 1; k <= K; k++) C[0][k] = 0;
for (int n = 0; n <= N; n++) C[n][0] = 1;
for (int n = 1; n <= N; n++)
for (int k = 1; k <= K; k++)
C[n][k] = C[n-1][k-1] + C[n-1][k];
donde N, K - valores máximos de su n, k
nCp = n! / ( p! (n-p)! ) =
( n * (n-1) * (n-2) * ... * (n - p) * (n - p - 1) * ... * 1 ) /
( p * (p-1) * ... * 1 * (n - p) * (n - p - 1) * ... * 1 )
Si podamos los mismos términos del numerador y el denominador, nos quedamos con una multiplicación mínima requerida. Podemos escribir una función en C para realizar multiplicaciones 2p y 1 división para obtener nCp:
int binom ( int p, int n ) {
if ( p == 0 ) return 1;
int num = n;
int den = p;
while ( p > 1 ) {
p--;
num *= n - p;
den *= p;
}
return num / den;
}