Самый быстрый способ генерации биномиальных коэффициентов

Ответ 1

Если вам нужны полные разложения при больших значениях n, свертка БПФ может быть самым быстрым способом. В случае биномиального разложения с равными коэффициентами (например, серии штрафных монетных монет) и четного порядка (например, числа бросков) вы можете использовать симметрии, таким образом:

Теория

Представляют результаты двух монетных томов (например, половину разницы между общим числом голов и хвостов) с выражением A + A * cos (Pi * n/N). N - количество выборок в вашем буфере - биномиальное расширение четного порядка O будет иметь коэффициенты O + 1 и требует буфера из N >= O/2 + 1 выборок - n - это номер образца, который генерируется, а A - масштабный коэффициент, который обычно будет либо 2 (для генерации биномиальных коэффициентов), либо 0,5 (для генерации биномиального распределения вероятности).

Обратите внимание, что по частоте это выражение напоминает биномиальное распределение этих двух монетных тонов - в положениях, соответствующих числу (head-tails)/2, есть три симметричных пика. Поскольку моделирование общего распределения вероятностей независимых событий требует свертывания их распределений, мы хотим свернуть наше выражение в частотной области, что эквивалентно умножению во временной области.

Алгоритм

• выберите размер буфера, N, то есть, по крайней мере, O/2 + 1 и может быть удобно DCTed
• инициализируйте его выражением pow (A + A * cos (Pi * n/N), O/2)
• применить прямой DCT-I
• считывает коэффициенты из буфера - первое число является центральным пиком, где head = tails, а последующие записи соответствуют симметричным парам последовательно дальше от центра.

Точность

Там предел того, насколько высокий O может быть до накопления ошибок округления с плавающей запятой, лишает вас точных целочисленных значений для коэффициентов, но я бы предположил, что число довольно велико. Плавная точка с двойной точностью может представлять собой 53-битные целые числа с полной точностью, и я буду игнорировать потери округления, связанные с использованием pow(), потому что генерирующее выражение будет иметь место в регистрах FP, что даст нам дополнительный 11 бит мантиссы, чтобы поглотить ошибку округления на платформах Intel. Поэтому, предполагая, что мы используем DCT-I 1024-х точек, реализованные через FFT, это означает, что точность 10 бит равна погрешности округления во время преобразования и не намного больше, оставляя нам ~ 43 бита чистого представления. Я не знаю, какой порядок биномиального расширения генерирует коэффициенты такого размера, но я смею сказать, что он достаточно велик для ваших нужд.

Асимметричные разложения

Если вам нужны асимметричные разложения для неравных коэффициентов a и b, вам нужно будет использовать двухсторонний (сложный) ДПФ и сложную функцию pow(). Выразите выражение A * A * e ^ (- Pi * я * n/N) + A * B + B * B * e ^ (+ Pi * я * n/N) [используя сложную функцию pow() для повышения это до половины порядка расширения] и DFT. То, что у вас в буфере, опять-таки, центральная точка (но не максимальная, если A и B очень разные) при смещении нуля, за ней следует верхняя половина распределения. Верхняя половина буфера будет содержать нижнюю половину распределения, соответствующую отрицательным значениям heads-минус-хвостов.

Обратите внимание, что исходные данные являются эрмитовыми симметричными (вторая половина входного буфера представляет собой комплексное сопряжение первого), поэтому этот алгоритм не является оптимальным и может быть выполнен с использованием комплексно-сложного БПФ, равного половине требуемого для оптимальной эффективности.

Излишне говорить, что все сложное возведение в степень будет переваривать больше времени процессора и повысить точность по сравнению с чисто реальным алгоритмом для симметричных распределений выше.

Ответ 2

Вы используете динамическое программирование для создания биномиальных коэффициентов

Вы можете создать массив и использовать цикл O (N ^ 2), чтобы заполнить его

C[n, k] = C[n-1, k-1] + C[n-1, k];

где

C[1, 1] = C[n, n] = 1

После этого в вашей программе вы можете получить значение C (n, k), просто просматривая ваш 2D-массив по индексам [n, k]

UPDATE 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];

где N, K - максимальные значения ваших n, k

Ответ 3

Если вам нужно вычислить их для всех n, ответ Ribtoks, вероятно, лучший. Для одного n вам лучше делать следующее:

C[0] = 1
for (int k = 0; k < n; ++ k)
    C[k+1] = (C[k] * (n-k)) / (k+1)

Деление точнее, если сделано после умножения.

И будьте осторожны с переполнением с помощью C [k] * (n-k): используйте достаточно большие целые числа.

Ответ 4

Это моя версия:

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

Ответ 5

Недавно я написал фрагмент кода, который должен был вызвать бинарный коэффициент около 10 миллионов раз. Таким образом, я использовал комбинацию lookup-table/calculate, которая все еще не слишком расточительна для памяти. Вы можете найти это полезным (и мой код находится в общественном достоянии). Код находится в

http://www.etceterology.com/fast-binomial-coefficients

Было предложено включить код здесь. Большая контрольная таблица для отслеживания кажется пустой тратой, поэтому здесь конечная функция и Python script, которая генерирует таблицу:

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))

Ответ 6

Если вам действительно нужен только тот случай, когда n намного больше p, можно было бы использовать формулу Стирлинга для факториалов. (если n → 1 и p - порядок один, Stirling приближенно n! и (n-p)!, держите p! как и т.д.)

Ответ 7

Самое быстрое разумное приближение в моем собственном бенчмаркинге - это аппроксимация, используемая библиотекой Apache Commons Maths: http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/special/Gamma.html#logGamma(double)

Мы с коллегами пытались понять, можем ли мы победить, используя точные вычисления, а не приближенные. Все подходы оказались неудачными (многие заказы медленнее), кроме одного, что было в 2-3 раза медленнее. В наиболее эффективном подходе используется https://math.stackexchange.com/a/202559/123948, вот код (в Scala):

var i: Int = 0
var binCoeff: Double = 1
while (i < k) {
  binCoeff *= (n - i) / (k - i).toDouble
  i += 1
}
binCoeff

Очень плохо подходит, когда различные попытки реализации Pascal Triangle используют хвостовую рекурсию.

Ответ 8

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 )

Если мы обрезаем одни и те же члены числителя и знаменателя, мы остаемся с минимальным умножением. Мы можем написать функцию в C, чтобы выполнить 2p умножения и 1 деление, чтобы получить 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;
}

Ответ 9

Я искал то же самое и не мог его найти, поэтому написал сам, который кажется оптимальным для любого биномиального Coeffcient, для которого endresult вписывается в Long.

// 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;
}

Ответ 10

Если я понимаю обозначение в вопросе, вы не просто хотите nCp, вы действительно хотите все nC1, nC2,... nC (n-1). Если это правильно, мы можем использовать следующее соотношение, чтобы сделать это довольно тривиальным:

• для всех k > 0: nCk = prod_ {из я = 1..k} ((n-i + 1)/i)
• то есть. для всех k > 0: nCk = nC (k-1) * (n-k + 1)/k

Вот фрагмент питона, реализующий этот подход:

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

Если вам нужны все коэффициенты до некоторого k > потолка (n/2), вы можете использовать симметрию, чтобы уменьшить количество операций, которые вам нужно выполнить, остановив при коэффициенте для потолка (n/2), а затем просто засыпая насколько вам нужно.

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