出longlong 的组合数对 素数求余

 

 

原文：http://hi.baidu.com/yu_zhuoran/blog/item/4993c6f3e24c7d59342acc0e.html

 

 

1) C(m, n) % p (p为素数, m, n <= 10^8)

算法伪代码:

1) 预处理 f[x] = (x!) mod p. (0 <= x < p)

2) 求m!中把p因子去掉，余下的数乘起来 mod p.

函数 G(m)求解 2)问题:
伪代码:

G(m)
1) if (0 == m) return 1;
2) m = k*p + t;
3) G[m] = f[p-1] ^ k * f[t] * G[m / k] mod p;

(其中 f[p-1]^k mod p 二分做）

3)
C(m, n) = m! / (n! * (m-n)!);

A = m! mod p;
B = n!*(m-n)! mod p;
C = A / B mod p;

在这里有俩种算法:
1) 求B的逆元; ( O(n^2) 枚举, 或者.. )
2) 扩展欧几里德求 BC = A(mod p) (= 是同余的意思)
最小的 C 就是答案; 这里的 A, B 转换成 对 2)问题
求解.

 

#include <iostream> #include <map> #include <set> #include <string> #include <algorithm> #include <math.h> #include <queue> using namespace std; template <class T> void out(T x, int n){ for(int i = 0; i < n; ++i) cout << x[i] << ' '; cout << endl; } template <class T> void out(T x, int n, int m){ for(int i = 0; i < n; ++i) out(x[i], m); cout << endl; } #define OUT(x) (cout << #x << " = " << x << endl) #define FOR(i, a, b) for(int i = (int)a; i < (int)b; ++i) #define REP(i, b) FOR(i, 0, b) #define FORD(i, a, b) for(int i = (int)a; i >= (int)b; --i) #define N 10007 #define MOD(x, y) ((x) % (y) + y) % (y) long long f[N+1]; long long ext_gcd(long long a, long long b, long long &x, long long &y) { if (0 == b) { x = 1; y = 0; return a; } else { long long tx, ty, d; d = ext_gcd(b, MOD(a, b), tx, ty); x = ty; y = tx - a / b * ty; return d; } } long long modular_linear(long long a, long long b, long long n) { long long tx, ty, d, x0, i; d = ext_gcd(a, n, tx, ty); if (0 == b % d) { x0 = MOD(tx * b / d, n); return x0; } else { return 1; } } void init(){ f[0] = 1; FOR (i, 1, N) { f[i] = (f[i-1] * i) % N; } } // m! 有多少p因子 long long calP(long long m, long long p){ long long ret = 0; while (m) { ret += m / p; m /= p; } return ret; } long long square(long long x){ return x * x; } long long bigmod(long long b, long long p, long long m) { if (0 == p) return 1; else if (0 == p % 2) return square(bigmod(b, p/2, m)) % m; else return ((b % m) * bigmod(b, p-1, m)) % m; } // 计算m!-mod p long long calM(long long m, long long p){ if (0 == m) return 1; long long ret = 0, k, t; k = m / p; t = m % p; ret = (bigmod(f[p-1], k, N) * f[t]) % p; ret *= calM(m / p, p); return ret % p; } long long calCmnModP(long long m, long long n, long long p){ long long A, B, u, v; u = calP(m, p); v = calP(n, p) + calP(m-n, p); if (u > v) return 0; else { A = calM(m, p); B = (calM(n, p) * calM(m-n, p)) % p; return modular_linear(B, A, p); } } int main(){ int t, n, m; scanf("%d", &t); init(); while (t--) { scanf("%d %d", &m, &n); printf("%d/n", calCmnModP(m, n, N)); } return 0; }