http://codeforces.com/gym/100633/problem/J
其实这个解法不难学的，不需要太多的数学。但是证明的话，我可能给不了严格的证明。可以看看这篇文章
http://www.cnblogs.com/jianglangcaijin/p/3446839.html 膜拜
#include <cstdio> #include <cstdlib> #include <cstring> #include <cmath> #include <algorithm> #include <assert.h> #define IOS ios::sync_with_stdio(false) using namespace std; #define inf (0x3f3f3f3f) typedef long long int LL; #include <iostream> #include <sstream> #include <vector> #include <set> #include <map> #include <queue> #include <string> #include <bitset> const int maxn = 1e6 + 20; LL MOD[maxn], r[maxn]; LL hasPrime(LL val, LL pi) { //求解val!含有多少个pi LL ans = 0; while (val) { ans += val / pi; val /= pi; } return ans; } LL quick_pow(LL a, LL b, LL MOD) { //求解 a^b%MOD的值 LL base = a % MOD; LL ans = 1; //相乘，所以这里是1 while (b) { if (b & 1) { ans = (ans * base) % MOD; //如果这里是很大的数据，就要用quick_mul } base = (base * base) % MOD; //notice。注意这里,每次的base是自己base倍 b >>= 1; } return ans; } LL exgcd(LL a, LL mod, LL &x, LL &y) { //求解a关于mod的逆元 ★：当且仅当a和mod互质才有用 if (mod == 0) { x = 1; y = 0; return a;//保留基本功能，返回最大公约数 } LL g = exgcd(mod, a % mod, x, y); LL t = x; //这里根据mod==0 return回来后， x = y; //x,y是最新的值x2,y2,改变一下，这样赋值就是为了x1=y2 y = t - (a / mod) * y; // y1=x2(变成了t)-[a/mod]y2; return g; //保留基本功能，返回最大公约数 } LL get_inv(LL a, LL MOD) { //求逆元。记得要a和MOD互质才有逆元的 LL x, y; //求a关于MOD的逆元，就是得到的k值是a*k%MOD==1 LL GCD = exgcd(a, MOD, x, y); if (GCD == 1) //互质才有逆元可说 return (x % MOD + MOD) % MOD; //防止是负数 else { assert(false); return -1;//不存在 } } LL factorialMod(LL n, LL pi, LL cnt) { //求解n! % pi^cnt if (!n) return 1; LL piPow = quick_pow(pi, cnt, 7e18), temp = 1; LL y = n / piPow;//分成y段，不要写在上面，piPow变量还没定义出来。 for (LL i = 1; i <= piPow; ++i) { //每piPow为一段，然后同余piPow if (i % pi == 0) continue; //pi的倍数早已算出 temp = temp * i % piPow; } //1 * 2 * 4 * 5 * 7 * 8和10 * 11 * 13 * 14 * 16 * 17模9的结果是一样的 LL ans = quick_pow(temp, y, piPow); //分成了y段然后同余piPow for (LL i = y * piPow + 1; i <= n; ++i) { //剩下的数字要暴力，例如19 if (i % pi == 0) continue; //pi的倍数早已算出 ans = ans * (i % piPow) % piPow; //取模两次，i会爆LL } return ans * factorialMod(n / pi, pi, cnt) % piPow; } LL CRT(LL r[], LL mod[], int n) { // X % mod[i] = r[i] LL M = 1; LL ans = 0; for (int i = 1; i <= n; ++i) { M *= mod[i]; assert(M > 0); } for (int i = 1; i <= n; ++i) { LL MI = M / mod[i]; //排除这个数 ans += r[i] * (MI * get_inv(MI, mod[i])); //使得MI * get_inv(MI, mod[i]) % mod[i] = 1 ans %= M; } if (ans < 0) ans += M; return ans; } LL calc(LL n, LL m, LL pi, LL cnt) { //求解C(n, m) % pi^cnt LL piPow = quick_pow(pi, cnt, 3e18); LL hasA = hasPrime(n, pi); LL hasB = hasPrime(n - m, pi); LL hasC = hasPrime(m, pi); LL ans = quick_pow(pi, hasA - hasB - hasC, piPow); hasA = factorialMod(n, pi, cnt); hasB = factorialMod(n - m, pi, cnt); hasC = factorialMod(m, pi, cnt); return ans * hasA % piPow * get_inv(hasB, piPow) % piPow * get_inv(hasC, piPow) % piPow; } LL exLucas(LL n, LL m, LL p) { //扩展lucas定理 if (n <= m) return 1 % p; int lenMod = 0; for (LL i = 2; i * i <= p; ++i) { if (p % i == 0) { //i是p的质因子 int cnt = 0; while (p % i == 0) { cnt++; p /= i; } ++lenMod; MOD[lenMod] = quick_pow(i, cnt, 7e18); r[lenMod] = calc(n, m, i, cnt); } } if (p > 1) { ++lenMod; MOD[lenMod] = p; r[lenMod] = calc(n, m, p, 1); } return CRT(r, MOD, lenMod); } void work() { LL n, m, p; cin >> n >> m >> p; cout << exLucas(n, m, p) << endl; } int main() { #ifdef local freopen("data.txt", "r", stdin); // freopen("data.txt", "w", stdout); #endif work(); return 0; }
组合数取模的所有应该都做完了吧，不会还有什么奇淫技巧吧。。如果有，求读者留言呀，Orz
1、组合数Cnm 防溢出公式 1、如果，(n&m)==m 那么C(n,m)为奇数，否则为偶数 2、求解C(n,0)、C(n,1)……C(n,n)有多少个奇数：ans=1<<(n二进制中1的个数) 3、Cn¬¬¬¬¬¬¬¬¬¬¬¬ – 1m - 1 + Cn – 1m = Cnm 4、求组合数的时候可能会溢出，这个时候我们可以边乘边除，来防止溢出。 因为Cnm = = 那么把上面的1用i来循环，从右到左计算即可 组合数是很大的，C(100,50)也会爆ULL，这个只能求些小的数，例如C(1e8,4)也不会爆。 LL C(LL n, LL m) { if (n < m) return 0; //防止sb地在循环 if (n == m) return 1; //C(0,0)也是1的 LL ans = 1; LL mx = max(n - m, m); //这个也是必要的。能约就约最大 LL mi = n - mx; for (int i = 1; i <= mi; ++i) { ans = ans * (mx + i) / i; } return ans; } Lucas 定理 解决很大的组合数问题 时间O（logp（n）*p）用在%很小的数比较有用。 求解C(n,m)%p 其中：n, m, p (1 <= m <= n <= 10^9, p是质数且p <= 1e5) 当MOD的数真的很小，MOD = 110119的话，可以预处理阶乘，这样快很多。 LL C(LL n, LL m, LL MOD) { if (n < m) return 0; //防止sb地在循环，在lucas的时候 if (n == m) return 1 % MOD; LL ans1 = 1; LL ans2 = 1; LL mx = max(n - m, m); //这个也是必要的。能约就约最大的那个 LL mi = n - mx; for (int i = 1; i <= mi; ++i) { ans1 = ans1 * (mx + i) %MOD; ans2 = ans2 * i % MOD; } return (ans1 * quick_pow(ans2, MOD - 2, MOD) % MOD); //这里放到最后进行,不然会很慢 } LL Lucas(LL n, LL m, LL MOD) { LL ans = 1; while (n && m && ans) { ans = ans * C(n % MOD, m % MOD, MOD) % MOD; n /= MOD; m /= MOD; } return ans; } NEFU 628 求解：C(n,m)%p的值。n, m and p (1 <= n, m, p <= 10^5)。 ★并且p有可能是合数 思路：设X = C(n, m) % p，那么X肯定可以分解成p1a * p2b …. * pzz这样的东西，那么把最后每个质因子剩下的个数算出来，进行快速幂对p取模即可。这里只进行了乘法，无须判断是否有逆元。 快速幂那里没有进行求逆元操作。 LL calc(LL n, int p) { LL ans = 0; while (n) { ans += n / p; // 2的倍数贡献一个2，然后4的倍数继续贡献一个。 n /= p; } return ans; } LL solve(LL n, LL m, LL p) { LL ans = 1; for (int i = 1; i <= total && prime[i] <= n ; i++) { LL t = calc(n, prime[i]); //calc是算出n!中有多少个prime[i]这个因子。 t -= calc(n - m, prime[i]); t -= calc(m, prime[i]); // t最小也是0 if (t) { // 就是当t不是0的时候 ans *= quick_pow(prime[i], t, p); ans %= p; } } return ans; } 有时候，对于p很少的情况p<=1e4，然后我们数据大T<=10000，这样，我们可以预处理出fac[i][j]表示 (j的阶乘)%prime[i]的值。inv[i][j]表示 (j!)关于prime[i]的逆元。然后O(1)处理。注意的是这个公式的话，fac[1][2]以及后面那些fac[1][3]…..都是0的，因为很简单，你如果阶乘中有数字>=prime[i]，那么%prime[i]后结果都是0。但是这样的后果就是C(6,2)%2等于0了。所以这里的组合数要用Lucas辅助来求得。(只能用Lucas，Lucas能避免这个情况) int fac[maxn][maxn]; // fac[i][j]表示(j!)%prime[i]的值 j<prime[i],如果j==prime[i]，后面的都是0 int inv[maxn][maxn]; // inv[i][j]表示 (j!)对prime[i]求逆元 void init() { for (int i = 1; i <= total; i++) { fac[i][0] = 1; inv[i][0] = 1; // (0!)=1 for (int j = 1; j < prime[i]; j++) { //等于prime[i]的话，%后是0了，没用 fac[i][j] = (j * fac[i][j - 1]) % prime[i]; inv[i][j] = quick_pow(fac[i][j], prime[i] - 2, prime[i]); } } return ; } int C(int n, int m, int MOD) { if (m > n) return 0; if (m == n) return 1; int pos = book[MOD]; //book[prime[i]]=i;表明这个素数下标是几多 return (fac[pos][n] * (inv[pos][n - m] * inv[pos][m] % MOD)) % MOD; } 这里想得到C(6, 2)%2的话。要调用Lucas(6, 2, 2);