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

Мне нужно рассчитать комбинации для числа.

Какой самый быстрый способ вычислить nCp, где n & gt; & gt; p?

Мне нужен быстрый способ для генерации биномиальных коэффициентов для полиномиального уравнения, и мне нужно получить коэффициент всех слагаемых и сохранить его в массиве.

(a+b)^n = a^n + nC1 a^(n-1) * b + nC2 a^(n-2) * ............ +nC(n-1) a * b^(n-1) + b^n

Каков наиболее эффективный способ расчета nCp?

 Pieter Geerkens22 мар. 2013 г., 23:33
@SteveJessop: у вас должно быть физическое или инженерное образование.
 Burhan Khalid14 июн. 2012 г., 14:12
whathaveyoutried.com
 Steve Jessop14 июн. 2012 г., 16:46
@ rao_555: хорошо, еслиn 10 ^ 10 иp являетсяn/2тогда nCp имеет приблизительно 45 миллиардов десятичных цифр. Лично я бы попытался изменить требования на этом этапе.
 Steve Jessop14 июн. 2012 г., 14:31
@ rao_555: какие размерыn а такжеp? А что за язык? Например, в языках, которые имеют существенно более быструю арифметику для небольших чисел, это имеет значение для ответа, даже стоит ли этоtrying вычислить результат вunsigned long long (пример взят из C) и избегайте переполнения с помощью хитрых уловок, или же результат намного больше, чем этот, так что нет смысла пытаться.
 Steve Jessop14 июн. 2012 г., 16:51
Извините, 10 миллиардов цифр. Еще.

Ответы на вопрос(11)

омиальных коэффициентов.

Вы можете создать массив и чем использовать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), просто взглянув на двумерный массив с индексами [n, k]

UPDATE что-то подобное

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

 06 февр. 2014 г., 17:07
Код @ScottNguyen показан в наиболее общей форме. Вы можете изменить N и K в соответствии с вашими пожеланиями
 06 февр. 2014 г., 15:18
не должны ли N и K иметь одинаковые значения? Я не думаю, что есть смысл запоминать непропорциональный треугольник Паскаля.
 15 июл. 2014 г., 13:07
+1. Хорошее решение Просто чтобы быть конкретным, вместо использования O (n * 2) памяти, его можно оценить с использованием O (N) памяти, так как в соответствии с определением проблемы в вопросе, он просто заинтересован в коэффициентах N.
 02 авг. 2018 г., 14:17
@ Мечтатель, да O (n) памяти на вопрос, заданный ОП, как вы и Скотт указываете, но все же O (n ^ 2) времени.

когда n намного больше, чем p, одним из способов будет использованиеФормула Стирлинга для факториалов. (если n & gt; & gt; 1 и p - порядок один, Стирлинг приближается к n! и (n-p) !, оставьте p! как есть и т. д.)

Космическая сложность: O (1)

public class binomialCoeff {
    static double binomialcoeff(int numerator, int denominator) 
    { 
        double res = 1; 
        //invalid numbers
        if (denominator>numerator || denominator<0 || numerator<0) {
            res = -1;
            return res;}
        //default values
        if(denominator==numerator || denominator==0 || numerator==0)
            return res;


        // Since C(n, k) = C(n, n-k) 
        if ( denominator > (numerator - denominator) ) 
            denominator = numerator - denominator;


        // Calculate value of [n * (n-1) *---* (n-k+1)] / [k * (k-1) *----* 1] 
        while (denominator>=1)
        { 

        res *= numerator;
        res = res / denominator; 

        denominator--;
        numerator--;
        } 

        return res; 
    } 

    /* Driver program to test above function*/
    public static void main(String[] args) 
    { 
        int numerator = 120; 
        int denominator = 20; 
        System.out.println("Value of C("+ numerator + ", " + denominator+ ") "
                                + "is" + " "+ binomialcoeff(numerator, denominator)); 
    } 

}

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

for all k>0: nCk = prod_{from i=1..k}( (n-i+1)/i ) i.e. for all k>0: nCk = nC(k-1) * (n-k+1) / k

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

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 & gt; потолок (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
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;
}
 08 сент. 2017 г., 09:18
num может превышать 32/64 бит, даже если результат не может быть. Функция практически безопасна только приp <= 12.

ляется приближение, используемое библиотекой 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

Действительно плохие подходы, когда различные попытки реализации треугольника Паскаля с использованием хвостовой рекурсии.

 19 дек. 2016 г., 00:32
Реализовал это в общих чертах - спасибо
Решение Вопроса

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

Theory

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

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

Другими словами, повышая выражение косинуса для результата двух бросков до степени (например, чтобы смоделировать 500 бросков, увеличьте его до степени 250, поскольку оно уже представляет пару), мы можем организовать биномиальное распределение для большого номер для отображения в частотной области. Поскольку все это реально и даже, мы можем заменить DCT-I на DFT для повышения эффективности.

Algorithm

decide on a buffer size, N, that is at least O/2 + 1 and can be conveniently DCTed initialise it with the expression pow(A + A*cos(Pi*n/N),O/2) apply the forward DCT-I read out the coefficients from the buffer - the first number is the central peak where heads=tails, and subsequent entries correspond to symmetrical pairs successively further from the centre

Accuracy

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

Asymmetrical expansions

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

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

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

ответ Ribtoks, вероятно, будет лучшим. Для одного n вам лучше сделать так:

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

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

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

который должен вызывать двоичный коэффициент около 10 миллионов раз. Поэтому я применил комбинированный подход «таблица поиска / расчет», который все еще не слишком расточителен для памяти. Вы можете найти это полезным (и мой код находится в свободном доступе). Код находится на

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

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

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

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

// 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;
}
 17 апр. 2016 г., 21:18
@ Атом, можете ли вы объяснить, почему этот цикл недопустим, указав, для какого входа ответ неверен? Я только что проверил это снова, он находится в maven io.github.htools в lib.MathTools и, кажется, работает для меня ...
 16 сент. 2014 г., 12:16
Хорошая точка зрения. Однако пытались ли вы использовать метод с плавающей запятой для очень больших чисел? Интересно, не теряет ли он точность в какой-то момент? Этот метод все еще может вычислять точные результаты, используя BigDecimal вместо longs.
 11 авг. 2014 г., 22:34
Довольно неплохо, но вложенный в цикл, предназначенный для избежания арифметики с плавающей точкой, кажется более дорогим, чем просто использование арифметики с плавающей точкой.
 17 апр. 2016 г., 21:29
@JeroenVuurens Для некоторых комбинацийn а такжеk break утверждение никогда не достигается и, следовательно, разделитель не устраняется.
 17 апр. 2016 г., 01:24
Этот алгоритм неверен. Петляfor (int divider: dividers) ... является недействительным.

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

Ваш ответ на вопрос