n、r、mが非常に大きな数の場合、プログラミングでnCr % m (つまり、「(n choose r)係数m」)を計算する方法は?
4 に答える
1
まず第一に、「非常に大きい」とはどういう意味かによります。問題が単純に中間値が標準の 64 ビット整数に対して大きくなりすぎることである場合は、gmpやBigIntegerなどを使用できます。
関連する数値が非常に大きくなり、メモリや忍耐力が不足し、中間値を完全な精度で計算できなくなる可能性があります。この場合、最善の策は、最初に各階乗の素因数分解を決定し、これらの因数分解を使用して二項式の因数分解を決定し、次に各中間ステップでモジュラスを使用してこの因数分解を乗算することです。
n までのすべての素数のリストが必要になります。
以下は疑似コードです。を使用してintいますが、使用している多数のライブラリに置き換える必要があります。
int factorial_prime_power(int f, int p) {
// When p is prime, returns k such that p^k divides f!
int k = 0;
while (f > p) {
f = f / p;
k = k + f;
}
return k;
}
int binomial_prime_power(int n, int r, int p) {
// when p is prime, returns k such that p^k divides nCr
return factorial_prime_power(n,p) - factorial_prime_power(r,p) - factorial_prime_power(n-r,p);
}
int powmod(int p, int k, int m) {
// quickly calculates p^k mod m
int res = 1;
int q = p;
int j = k;
while (j > 0) {
// invariant: p^k is congruent to res * q^j
if (j is odd) {
res = (res * q) % m;
j = (j-1)/2;
} else {
j = j / 2;
}
q = (q * q) % m;
}
return res;
}
int big_binomial(int n, int r, int m) {
if (n < r or r < 0) {
return 0;
}
int res = 1;
for(p in all primes from 2 to n) {
k = binomial_prime_power(n,r,p);
res = (res * powmod(p,k,m)) % m;
}
return res;
}
于 2012-11-11T15:25:38.380 に答える
1
の素因数がそれぞれ大きすぎmないという条件で、整数 を法とする二項係数の計算の複雑さを軽減する手法があります。あまりにも大きな力でm分割します。m
最初のステップは因数分解ですm。
m = ∏ (p_i ^ e_i)
次に、これらの素数ベキのそれぞれを法とする二項係数を計算し、その結果を中国剰余定理と組み合わせます。
素数 ( e_i == 1) を法とする二項係数の計算は、一般的な場合よりも少し簡単に計算できます。たとえば、この回答ですが、そこに含めることもできます。
素数pとについて、n >= 0定義しましょう
F(p, n) = ∏ k = p^(n/p) * (n/p)!
1<=k<=n
p | k
と
G(p, n) = ∏ k
1<=k<=n
gcd(k,n)=1
次に、
n! = F(p, n) * G(p, n)
(n/p)!に現れるに対して同じ分割を繰り返し使用しF(p, n)、
m
n! = p^K * ∏ G(p, n/(p^j))
j=0
どこでp^m <= n < p^(m+1)。のすべての因数はG(p, x)と互いに素p^eであるため、二項係数の分母の対応する因数は を法として逆変換できます。また、 を法としてp^e計算する効率的な方法が見つかれば、を法として二項係数を効率的に計算する方法が得られます。G(p, x)p^ep^e
二項係数については、次のようになります。
n! / (r! * (n-r)!) = p^M * (∏ G(p, (n/p^j)) * [ ∏ G(p, r/(p^j)) * ∏ G(p, (n-r)/(p^j)) ]^(-1)
しましょうH(p, e, n) = G(p, n) % (p^e)。重要な点は、すべての数の互いに素であって次p^eを超えない積p^eは非常に単純であるということです。とでない限り、-1モジュロに合同であり、その場合は 1 に合同です。p^ep = 2e > 2
そう
H(p, e, n) ≡ (-1)^(n/(p^e)) * H(p, e, n % (p^e)) (mod p^e)
(最初の因数が 1 の場合は と をp = 2除きます)、 を計算するだけでよいので、結果を調べることができます。e > 2H(p, e, k)0 <= k < p^e
コード：
// invert k modulo p, k and p are supposed coprime
unsigned long long invertMod(unsigned long long k, unsigned long long p) {
unsigned long long q, pn = 1, po = 0, r = p, s = k;
unsigned odd = 1;
do {
q = r/s;
q = pn*q + po;
po = pn;
pn = q;
q = r%s;
r = s;
s = q;
odd ^= 1;
}while(pn < p);
return odd ? p-po : po;
}
// Calculate the binomial coefficient (n choose k) modulo (prime^exponent)
// requires prime to be a prime, exponent > 0, and 0 <= k <= n,
// furthermore supposes prime^exponent < 2^32, otherwise intermediate
// computations could have mathematical results out of range.
// If k or (n-k) is small, a direct computation would be more efficient.
unsigned long long binmod(unsigned long long prime, unsigned exponent,
unsigned long long n, unsigned long long k) {
// The modulus, prime^exponent
unsigned long long ppow = 1;
// We suppose exponent is small, so that exponentiation by repeated
// squaring wouldn't gain much.
for(unsigned i = 0; i < exponent; ++i) {
ppow *= prime;
}
// array of remainders of products
unsigned long long *remainders = malloc(ppow * sizeof *remainders);
if (!remainders) {
fprintf(stderr, "Allocation failure\n");
exit(EXIT_FAILURE);
}
for(unsigned long long i = 1; i < ppow; ++i) {
remainders[i] = i;
}
for(unsigned long long i = 0; i < ppow; i += prime) {
remainders[i] = 1;
}
for(unsigned long long i = 2; i < ppow; ++i) {
remainders[i] *= remainders[i-1];
remainders[i] %= ppow;
}
// Now to business.
unsigned long long pmult = 0, ntemp = n, ktemp = k, mtemp = n-k,
numer = 1, denom = 1, q, r, f;
if (prime == 2 && exponent > 2) {
f = 0;
} else {
f = 1;
}
while(ntemp) {
r = ntemp % ppow;
q = ntemp / ppow;
numer *= remainders[r];
numer %= ppow;
if (q & f) {
numer = ppow - numer;
}
ntemp /= prime;
pmult += ntemp;
}
while(ktemp) {
r = ktemp % ppow;
q = ktemp / ppow;
denom *= remainders[r];
denom %= ppow;
if (q & f) {
denom = ppow - denom;
}
ktemp /= prime;
pmult -= ktemp;
}
while(mtemp) {
r = mtemp % ppow;
q = mtemp / ppow;
denom *= remainders[r];
denom %= ppow;
if (q & f) {
denom = ppow - denom;
}
mtemp /= prime;
pmult -= mtemp;
}
// free memory before returning, we don't use it anymore
free(remainders);
if (pmult >= exponent) {
return 0;
}
while(pmult > 0) {
numer = (numer * prime) % ppow;
--pmult;
}
return (numer * invertMod(denom, ppow)) % ppow;
}
これは段階的に計算n choose k modulo p^eしO(p^e + log n)ます。
于 2012-11-08T20:10:47.957 に答える
1
nCr のようなものを大量に計算し、最適化される可能性のあるライブラリ関数を検索したいと思います。その場合、 GNU Scientific Libraryを見ることができます。
于 2012-11-08T06:53:28.127 に答える
0
チャンクでそれを行います。
例: n!/(r!(nr)!) = 1*2*...*n/(1*2*...r*1*2*...*(nr))=1/ 1*1 * 2/2*2 * 3/3*3 * ... *... .
各チャンクは簡単に計算できるため、計算でオーバーフローを避ける必要があります。チャンクを計算し、現在の結果にそれを掛けます。
また、その前に n と r をモジュロする価値があります。
于 2012-11-08T07:01:02.147 に答える