数値の組み合わせを計算する必要があります。
n >> pの場合にnCpを計算する最速の方法は何ですか?
多項式の二項係数を生成する高速な方法が必要であり、すべての項の係数を取得して配列に格納する必要があります。
(a + b)^ n = a ^ n + nC1 a ^(n-1)* b + nC2 a ^(n-2)* ............ + nC(n-1 )a * b ^(n-1)+ b ^ n
nCpを計算する最も効率的な方法は何ですか?
数値の組み合わせを計算する必要があります。
n >> pの場合にnCpを計算する最速の方法は何ですか?
多項式の二項係数を生成する高速な方法が必要であり、すべての項の係数を取得して配列に格納する必要があります。
(a + b)^ n = a ^ n + nC1 a ^(n-1)* b + nC2 a ^(n-2)* ............ + nC(n-1 )a * b ^(n-1)+ b ^ n
nCpを計算する最も効率的な方法は何ですか?
二項係数を生成するために動的計画法を使用します
配列を作成し、O(N ^ 2)ループを使用して配列を埋めることができます
C[n, k] = C[n-1, k-1] + C[n-1, k];
どこ
C[1, 1] = C[n, n] = 1
その後、プログラムで[n、k]インデックスの2D配列を見るだけでC(n、k)値を取得できます。
そのような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の最大値
すべてのnについてそれらを計算する必要がある場合は、Ribtoksの答えがおそらく最良です。単一のnの場合、次のように実行することをお勧めします。
C[0] = 1
for (int k = 0; k < n; ++ k)
C[k+1] = (C[k] * (n-k)) / (k+1)
乗算後に行われる場合、除算は正確です。
また、C [k] *(nk)でオーバーフローすることに注意してください。十分な大きさの整数を使用してください。
nの値が大きい場合に完全な展開が必要な場合は、FFT畳み込みが最速の方法である可能性があります。係数が等しく(たとえば、一連の公正なコイントス)、順序が偶数(たとえば、トスの数)の二項展開の場合、対称性を利用できます。
仮説
式A+A * cos(Pi * n / N)を使用して、2回のコイン投げの結果(たとえば、頭と尾の総数の差の半分)を表します。Nはバッファー内のサンプル数です-偶数次Oの二項展開はO+1係数を持ち、N> = O / 2 +1サンプルのバッファーを必要とします-nは生成されるサンプル数であり、Aは通常は2(二項係数を生成する場合)または0.5(二項確率分布を生成する場合)のいずれかのスケール係数。
頻度では、この式はこれら2つのコイントスの二項分布に似ていることに注意してください。数(頭-尾)/2に対応する位置に3つの対称的なスパイクがあります。独立したイベントの全体的な確率分布をモデル化するには、それらの分布を畳み込む必要があるため、時間領域での乗算に相当する周波数領域で式を畳み込みます。
言い換えると、2回のトスの結果の正弦式を累乗することで(たとえば、500回のトスをシミュレートし、すでにペアを表しているため、250の累乗に上げる)、大規模な二項分布を調整できます。周波数領域に表示される番号。これはすべて現実的で均一であるため、効率を向上させるためにDFTの代わりにDCT-Iを使用できます。
アルゴリズム
正確さ
累積された浮動小数点の丸め誤差が係数の正確な整数値を奪う前に、Oがどれだけ高くなるかには限界がありますが、その数はかなり高いと思います。倍精度浮動小数点は53ビット整数を完全な精度で表すことができます。生成式はFPレジスタで行われるため、pow()の使用に伴う丸め誤差は無視し、11を追加します。 Intelプラットフォームでの丸め誤差を吸収するための仮数のビット。したがって、FFTを介して実装された1024ポイントのDCT-Iを使用すると仮定すると、変換中の丸め誤差で10ビットの精度が失われ、それ以外の場合はほとんど失われず、約43ビットのクリーンな表現が残ります。二項式展開の順序がそのサイズの係数を生成するかどうかはわかりませんが、それはあなたのニーズに十分な大きさであるとあえて言います。
非対称展開
aとbの係数が等しくない場合に非対称展開が必要な場合は、両側(複素数)DFTと複素数pow()関数を使用する必要があります。式A*A * e ^(-Pi * i * n / N)+ A * B + B * B * e ^(+ Pi * i * n / N)を生成します[複雑なpow()関数を使用して発生させます膨張次数の半分の累乗]とDFTします。バッファにあるのは、ここでも、オフセットゼロの中心点(ただし、AとBが大きく異なる場合は最大値ではありません)であり、その後に分布の上半分が続きます。バッファーの上半分には、負のheads-minus-tails値に対応する、分布の下半分が含まれます。
ソースデータがエルミート対称であることに注意してください(入力バッファーの後半は最初の複素共役です)。したがって、このアルゴリズムは最適ではなく、最適化に必要なサイズの半分の複素数から複素数へのFFTを使用して実行できます。効率。
言うまでもなく、すべての複雑なべき乗は、上記の対称分布の純粋に実際のアルゴリズムと比較して、より多くのCPU時間を消費し、精度を損ないます。
これは私のバージョンです:
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
私は最近、二項係数を約1,000万回呼び出す必要のあるコードを作成しました。そこで、メモリをあまり浪費しないルックアップテーブル/計算の組み合わせアプローチを実行しました。あなたはそれが役に立つと思うかもしれません(そして私のコードはパブリックドメインにあります)。コードはにあります
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))
nがpよりもはるかに大きい場合のみが本当に必要な場合、1つの方法は、階乗にスターリングの公式を使用することです。(n >> 1でpが1次の場合、スターリング近似n!および(np)!、 p!をそのままにしておくなど)
私自身のベンチマークで最も速く合理的な近似は、Apache Commons Mathsライブラリで使用される近似です:http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/special/Gamma.html #logGamma(double)
同僚と私は、概算ではなく正確な計算を使用して、それを打ち負かすことができるかどうかを確認しようとしました。2〜3倍遅い1つを除いて、すべてのアプローチが惨めに失敗しました(多くの注文が遅くなりました)。最高のパフォーマンスを発揮するアプローチは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
末尾再帰を使用してパスカルの三角形を実装しようとするさまざまな試みでの本当に悪いアプローチ。
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;
}
同じものを探していましたが見つかりませんでした。そのため、最終結果が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;
}
質問の表記法を理解している場合は、nCpだけでなく、実際にはnC1、nC2、... nC(n-1)のすべてが必要です。これが正しければ、次の関係を利用して、これをかなり簡単にすることができます。
このアプローチを実装する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>ceiling(n / 2)までのすべての係数が必要な場合は、対称性を使用して、ceiling(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
時間計算量:O(分母)空間計算量: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));
}
}