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 + ;
LL MOD[maxn], r[maxn];
LL hasPrime(LL val, LL pi) { //求解val!含有多少个pi
LL ans = ;
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
while (b) {
if (b & ) {
ans = (ans * base) % MOD; //如果这里是很大的数据，就要用quick_mul
}
base = (base * base) % MOD; //notice。注意这里,每次的base是自己base倍
b >>= ;
}
return ans;
}
LL exgcd(LL a, LL mod, LL &x, LL &y) {
//求解a关于mod的逆元 ★：当且仅当a和mod互质才有用
if (mod == ) {
x = ;
y = ;
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 == ) //互质才有逆元可说
return (x % MOD + MOD) % MOD; //防止是负数
else {
assert(false);
return -;//不存在
}
} LL factorialMod(LL n, LL pi, LL cnt) { //求解n! % pi^cnt
if (!n) return ;
LL piPow = quick_pow(pi, cnt, 7e18), temp = ;
LL y = n / piPow;//分成y段，不要写在上面，piPow变量还没定义出来。
for (LL i = ; i <= piPow; ++i) { //每piPow为一段，然后同余piPow
if (i % pi == ) 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 + ; i <= n; ++i) { //剩下的数字要暴力，例如19
if (i % pi == ) 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 = ;
LL ans = ;
for (int i = ; i <= n; ++i) {
M *= mod[i];
assert(M > );
}
for (int i = ; 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 < ) 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 % p;
int lenMod = ;
for (LL i = ; i * i <= p; ++i) {
if (p % i == ) { //i是p的质因子
int cnt = ;
while (p % i == ) {
cnt++;
p /= i;
}
++lenMod;
MOD[lenMod] = quick_pow(i, cnt, 7e18);
r[lenMod] = calc(n, m, i, cnt);
}
}
if (p > ) {
++lenMod;
MOD[lenMod] = p;
r[lenMod] = calc(n, m, p, );
}
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 ;
}
组合数取模的所有应该都做完了吧，不会还有什么奇淫技巧吧。。如果有，求读者留言呀，Orz
、组合数Cnm 防溢出公式
、如果，(n&m)==m 那么C(n,m)为奇数，否则为偶数
、求解C(n,)、C(n,)……C(n,n)有多少个奇数：ans=<<(n二进制中1的个数)
、Cn¬¬¬¬¬¬¬¬¬¬¬¬ – 1m - + Cn – 1m = Cnm
、求组合数的时候可能会溢出，这个时候我们可以边乘边除，来防止溢出。
因为Cnm = = 那么把上面的1用i来循环，从右到左计算即可
组合数是很大的，C(,)也会爆ULL，这个只能求些小的数，例如C(1e8,)也不会爆。
LL C(LL n, LL m) {
if (n < m) return ; //防止sb地在循环
if (n == m) return ; //C(0,0)也是1的
LL ans = ;
LL mx = max(n - m, m); //这个也是必要的。能约就约最大
LL mi = n - mx;
for (int i = ; i <= mi; ++i) {
ans = ans * (mx + i) / i;
}
return ans;
} Lucas 定理 解决很大的组合数问题 时间O（logp（n）*p）用在%很小的数比较有用。
求解C(n,m)%p 其中：n, m, p ( <= m <= n <= ^, p是质数且p <= 1e5)
当MOD的数真的很小，MOD = 110119的话，可以预处理阶乘，这样快很多。
LL C(LL n, LL m, LL MOD) {
if (n < m) return ; //防止sb地在循环，在lucas的时候
if (n == m) return % MOD;
LL ans1 = ;
LL ans2 = ;
LL mx = max(n - m, m); //这个也是必要的。能约就约最大的那个
LL mi = n - mx;
for (int i = ; i <= mi; ++i) {
ans1 = ans1 * (mx + i) %MOD;
ans2 = ans2 * i % MOD;
}
return (ans1 * quick_pow(ans2, MOD - , MOD) % MOD); //这里放到最后进行,不然会很慢
} LL Lucas(LL n, LL m, LL MOD) {
LL ans = ;
while (n && m && ans) {
ans = ans * C(n % MOD, m % MOD, MOD) % MOD;
n /= MOD;
m /= MOD;
}
return ans;
} NEFU
求解：C(n,m)%p的值。n, m and p ( <= n, m, p <= ^)。 ★并且p有可能是合数
思路：设X = C(n, m) % p，那么X肯定可以分解成p1a * p2b …. * pzz这样的东西，那么把最后每个质因子剩下的个数算出来，进行快速幂对p取模即可。这里只进行了乘法，无须判断是否有逆元。 快速幂那里没有进行求逆元操作。
LL calc(LL n, int p) {
LL ans = ;
while (n) {
ans += n / p; // 2的倍数贡献一个2，然后4的倍数继续贡献一个。
n /= p;
}
return ans;
}
LL solve(LL n, LL m, LL p) {
LL ans = ;
for (int i = ; 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<=，这样，我们可以预处理出fac[i][j]表示 (j的阶乘)%prime[i]的值。inv[i][j]表示 (j!)关于prime[i]的逆元。然后O()处理。注意的是这个公式的话，fac[][]以及后面那些fac[][]…..都是0的，因为很简单，你如果阶乘中有数字>=prime[i]，那么%prime[i]后结果都是0。但是这样的后果就是C(,)%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 = ; i <= total; i++) {
fac[i][] = ;
inv[i][] = ; // (0!)=1
for (int j = ; j < prime[i]; j++) { //等于prime[i]的话，%后是0了，没用
fac[i][j] = (j * fac[i][j - ]) % prime[i];
inv[i][j] = quick_pow(fac[i][j], prime[i] - , prime[i]);
}
}
return ;
}
int C(int n, int m, int MOD) {
if (m > n) return ;
if (m == n) return ;
int pos = book[MOD]; //book[prime[i]]=i;表明这个素数下标是几多
return (fac[pos][n] * (inv[pos][n - m] * inv[pos][m] % MOD)) % MOD;
}
这里想得到C(, )%2的话。要调用Lucas(, , );
