¿Cómo calcular multinomios de forma eficaz?

Question

Estoy intentando reproducir MULTINOMIAL de Excel en C# de modo que

$${MULTINOMIAL(a,b,..,n)} = \frac{(a+b +\cdots +n)!}{a!b! \cdots n!}$$

¿Cómo puedo hacerlo sin provocar un desbordamiento debido a los factoriales?

Answer 1

Otra forma más eficiente de calcular exactamente los coeficientes, que generalmente no debería desbordarse a menos que el resultado lo haga, es utilizar la caracterización del coeficiente multinomial como un producto de coeficientes binomiales: $${a+b+c+\cdots+n\choose a\;b\;c\cdots\;n} = {a+b\choose b}{a+b+c\choose c}\cdots{a+b+c+\cdots +n\choose n}$$ Esto es fácil de demostrar multiplicando el lado derecho y observando que los factores de $(a+b)!,\;(a+b+c)!,\;\ldots$ cancelar. Esto simplifica el problema de calcular un coeficiente multinomial a sólo calcular uno binario - y hay varias maneras directas de hacerlo sin desbordamiento. La más sencilla puede ser utilizar la fórmula recursiva ${n+1\choose k+1} = \frac{n+1}{k+1}{n\choose k}$ junto con una manipulación cuidadosa de los resultados (exactos) para asegurarse de que siempre está dividiendo los DGC; para obtener más información sobre este enfoque, consulte http://blog.plover.com/math/choose.html . Alternativamente, puedes usar el teorema de Kummer de que la potencia de un primo $p$ dividiendo ${a+b\choose b}$ es el número de transportes cuando $a$ y $b$ se escriben como base- $p$ números y añadidos. Haciendo esto para un intervalo adecuado de números primos, se puede hallar la factorización en primos del resultado y, por tanto, el propio resultado (de hecho, se podría simplemente acumular el resultado hallando el exponente divisor de todos los primos menores que el máximo de su coeficiente binomial por turno).

Answer 2

Si tu ordenador puede manejar \logaritmos y exponenciales de números, podrías calcular lo siguiente: $$\log(\frac{(a+b +\cdots +n)!}{a!b! \cdots n!})=\log(a+b +\cdots +n)!)-\log(a!b! \cdots n!)=$$ $$\log(2)+\log(3)+\dots+\log(a+b +\cdots +n)-\log(a!)-\log(b!)-\ldots -\log(x!)=$$ $$\log(2)+\log(3)+\dots+\log(a+b +\cdots +n)-\log(2)-\ldots-\log(a)-\log(2)\ldots-\log(b)\ldots \log(2)-\ldots \log(x)$$

Cuando haces todas estas sumas, exponencias el resultado usando la base con la que empezaste (como está escrito la base es diez, pero podrías elegir otra base).

Como nota histórica, esto es esencialmente lo que la gente hacía antes de los ordenadores. Utilizaban $\log$ mesas. No soy programador informático, así que no sé hasta qué punto esto es factible en la práctica o si existen métodos mejores. También se pueden multiplicar números enormes más rápidamente con una transformada rápida de Fourier, pero eso podría ser excesivo.

Answer 3

Se trata más bien de programación, pero no de matemáticas :)

Existen cuatro soluciones (véase también código en github ) :

1. Utilice BigInteger para números grandes. Véase también esta pregunta sobre SO .
2. Usa mi código de abajo. No es el más eficiente, pero bastante bueno y parece ser mejor que la implementación mazo.
3. Calcular valores para ciertos argumentos en Excel y almacenarlos en una tabla para su uso posterior. Pero tal vez no todos los casos se pueden calcular debido a su gran cantidad.
4. Utiliza los logaritmos como se describe en la respuesta de @Baby Dragon (mi código de logaritmos también está abajo).

Mi aplicación:

public static ulong Mutinomonal(params uint[] numbers)
{
    uint numbersSum = 0;
    foreach (var number in numbers)
        numbersSum += number;

    uint maxNumber = numbers.Max();
    var denomFactorPowers = new uint[numbers.Max() + 1];
    foreach (var number in numbers)
        for (int i = 2; i <= number; i++)
            denomFactorPowers[i]++;
    for (int i = 2; i < denomFactorPowers.Length; i++)
        denomFactorPowers[i]--; // reduce with nominator;

    uint currentFactor = 2;
    uint currentPower = 1;
    double result = 1;
    for (uint i = maxNumber + 1; i <= numbersSum; i++)
    {
        uint tempDenom = 1;
        while (tempDenom < result && currentFactor < denomFactorPowers.Length)
        {
            if (currentPower > denomFactorPowers[currentFactor])
            {
                currentFactor++;
                currentPower = 1;
            }
            else
            {
                tempDenom *= currentFactor;
                currentPower++;
            }
        }
        result = result / tempDenom * i;
    }

    return (ulong)result;
}

Aplicación ingenua:

public static ulong Mutinomonal(params uint[] numbers)
{
    uint numbersSum = 0;
    foreach (var number in numbers)
        numbersSum += number;

    ulong nominator = Factorial(numbersSum);

    ulong denominator = 1;
    foreach (var number in numbers)
        denominator *= Factorial(number);

    return nominator / denominator;
}

public static ulong Factorial(ulong n)
{
    ulong result = 1;
    for (ulong i = 1; i <= n; i++)
        result = result * i;
    return result;
}

He comprobado ambas implementaciones en $1, 2, 3, 4, 5, 6$ números y en la última implementación se ha producido desbordamiento, a diferencia de la primera.

ACTUALIZACIÓN

Estimé la respuesta de @Baby Dragon y escribí esta implementación del logaritmo en C#:

public static List<double> Logarithms = new List<double>();

public static ulong Multinominal(params uint[] numbers)
{
    uint numbersSum = 0;
    foreach (var number in numbers)
        numbersSum += number;

    double sum = 0;
    for (int i = 2; i <= numbersSum; i++)
    {
        if (i - 2 < Logarithms.Count)
            sum += Logarithms[i - 2];
        else
        {
            var log = Math.Log(i);
            Logarithms.Add(log);
            sum += log;
        }
    }

    foreach (var number in numbers)
        for (int i = 2; i <= number; i++)
            sum -= Logarithms[i - 2];

    return (ulong)Math.Exp(sum);
}

Answer 4

Otra implementación (en pari/GP) evita los logaritmos y funciona completamente en modo entero. Es una extensión del método para calcular el binomio secuencialmente : [update] $\binom {a+b}{a} = {a+1\over 1} \cdot{a+2\over 2} \cdot{a+3\over 3}\cdot{a+4\over 4}\cdots{a+b\over b} = (((((a+1)/1 \cdot (a+2) )/2 \cdot (a+3) )/ 3 \cdots(a+b) / b$ . (Esto también está en la respuesta de Steven) [/actualización]
Aquí está el código:

{mulbin(v)=local(pr,g,maxf);
 pr = vector(#v,r,1); 
 g  = vector(#v); g[1]=1;for(k=2,#v,g[k]=g[k-1]+v[k-1]);
 maxf=1; for(k=2,#v,maxf=max(maxf,v[k]));
for(f=1,maxf,
 for(k=2,#v , if(f>v[k],next());
    pr[k] *= g[k]; pr[k] /=f;
    g[k] ++ ;
  );
 );
 return(prod(k = 2,#v,pr[k])); }

 gp> mulbin([8,2,3])
 %293 = 12870

Es una optimización proporcionar primero los elementos del vector con el número más alto; porque el bucle exterior (maxf) es el máximo de las entradas penúltimas en el vector solamente.

El método tiene además la ventaja de que cuando uno decide utilizar los logaritmos de todos modos: entonces debido a la organización del bucle externo y el interno no necesita calcular el logaritmo de un número dos veces (que es una operación costosa). El código cambiaría simplemente a

{logmulbin(v)=local(logpr,g,maxf,logf);
 logpr = vector(#v,r,0); 
 g  = vector(#v); g[1]=1;for(k=2,#v,g[k]=g[k-1]+v[k-1]);
 maxf=1; for(k=2,#v,maxf=max(maxf,v[k]));
for(f=1,maxf, logf=log(f);
 for(k=2,#v , if(f>v[k],next());
    logpr[k] += log(g[k]); logpr[k] -= logf;
    g[k] ++ ;
  );
 );
 return(sum(k = 2,#v,logpr[k])); }

gp> exp(logmulbin([8,3,2]))
%294 = 12870.0000000

Answer 5

La respuesta que dio Steven Stadnicki es sencilla de aplicar utilizando el lenguaje de programación Python; véase este solución (y se acaba de actualizar en enero / 2020).

Cada coeficiente multinomial puede descomponerse en un producto de coeficientes binomiales. Si es fundamental optimizar el tiempo de ejecución, se puede utilizar un fino/esparcido/diagonal tabla de coeficientes binomiales podría almacenarse en una tabla de consulta almacenada en memoria con la especificación de alto nivel

$\tag 1 {\text{BC-Table}_{\;n,k} =\displaystyle {n \choose k}} \text{ where } k \approx \frac{n}{2}$

El programa emplearía la lógica para desglosar $\text{MULTINOMIAL}(a,b,..,n)$ en el producto de al menos un coeficiente binomial y otros dos coeficientes multinomiales, intentando que el factor binomial esté en la tabla; si el algoritmo controlador "renuncia a una coincidencia en la tabla", el trabajo sobre la pieza binomial se pasaría a un algoritmo que utilice el método de Steven.

El programa también emplearía métodos de programación paralela/concurrente en las piezas generadas.

---

Nota: Uno (de los muchos) desgloses completos de $\text{MULTINOMIAL}(a_1,a_2,\dots,a_k)$ en coeficientes binomiales viene dada por

$\tag 2 \displaystyle{\binom{n}{a_1,a_2,\dots,a_k} = \binom{n}{a_1}\binom{n-a_1}{a_2}\binom{n-a_1-a_2}{a_3}\cdots \binom{a_k}{a_k}}$

Utilizando $\text{(2)}$ podemos escribir

$\quad \displaystyle{\text{MULTINOMIAL}(1,2,3,4,5,6,7) = \binom{28}{1}\binom{27}{2}\binom{25}{3}\binom{22}{4}\binom{18}{5}\binom{13}{6}\binom{7}{7}}$

Pero nuestro (hipotético) programa comenzaría con

$\tag 3 \displaystyle{\text{MULTINOMIAL}(1,2,3,4,5,6,7) =}$

$\quad \quad \quad \displaystyle{\text{MULTINOMIAL}(14,14) \times \text{MULTINOMIAL}(1,2,5,6) \times \text{MULTINOMIAL}(3,4,7)}$

Si necesita comprobar $\text{(3)}$ véase wolfram .

Completar la descomposición en coeficientes binomiales,

$\tag 4 \displaystyle{\text{MULTINOMIAL}(1,2,5,6) =}$

$\quad \quad \quad \displaystyle{\text{MULTINOMIAL}(7,7) \times \text{MULTINOMIAL}(1,6) \times \text{MULTINOMIAL}(2,5)}$

y

$\tag 5 \displaystyle{\text{MULTINOMIAL}(3,4,7) =}$

$\quad \quad \quad \displaystyle{\text{MULTINOMIAL}(7,7) \times \text{MULTINOMIAL}(3,4) \times \text{MULTINOMIAL}(7)}$