非常好的学习博客:https://blog.csdn.net/skywalkert/article/details/52553048
下图转至:https://blog.csdn.net/Bfk_zr/article/details/78310541
例题1.hdu5446:lucas+普通CRT
求c(n,m)%p 1<=n,m.p<=1e18 保证p为最多15个不同pi的乘积 pi<=1e5
#include<iostream>
#include<stdio.h>
using namespace std;
typedef long long ll;
//求c(n,m)%p nmp<=1e18 保证p为不同pi的乘积 pi<=1e5
ll mul(ll a, ll b, ll p) //快速乘,计算a*b%p
{
ll ret = 0;
while(b)
{
if(b & 1) ret = (ret + a) % p;
a = (a + a) % p; b >>= 1;
}
return ret;
}
void ex_gcd(ll a, ll b, ll &x, ll &y, ll &d)
{
if (!b){
d = a, x = 1, y = 0;
}
else
{
ex_gcd(b, a % b, y, x, d);
y -= x * (a / b);
}
}
ll inv(ll t, ll p) //如果不存在,返回-1
{
ll d, x, y;
ex_gcd(t, p, x, y, d);
return d == 1 ? (x % p + p) % p : -1;
}
//n个方程:x=a[i](mod m[i]) (0<=i<n) m[i]两两互质
ll china(int n, ll *a, ll *m)
{
ll M = 1, ret = 0;
for(int i = 0; i < n; i ++) M *= m[i];
for(int i = 0; i < n; i ++)
{
ll w = M / m[i];
//ret = (ret + w * inv(w, m[i]) * a[i]) % M;
ret = (ret + mul(w * inv(w, m[i]), a[i], M)) % M;//防爆long long int
}
return (ret + M) % M;
}
ll exp_mod(ll a,ll b,ll p)
{
ll res=1;
while(b!=0)//利用二进制求乘方
{
if(b&1) res=(res*a)%p;
a=(a*a)%p;
b>>=1;
}
return res;
}
ll comb(ll a,ll b,ll p)
{
if(a<b) return 0;
if(a==b) return 1;
if(b>a-b) b=a-b;//组合的性质,n-m和m是相等的
ll ans=1,ca=1,cb=1;
for(ll i=0; i<b; ++i)
{
ca=(ca*(a-i))%p;//这里ca表示n!/(n-m)!
cb=(cb*(b-i))%p;
}
ans= (ca*exp_mod(cb,p-2,p))%p;//这里对m!求逆元即可
return ans;
}
ll lucas(ll n,ll m,ll p)
{
ll ans=1;
while(n&&m&&p)
{
ans=(ans*comb(n%p,m%p,p))%p;
n/=p;
m/=p;
}
return ans;
}
ll mm[100005];//部分解
ll pp[100005];//质数
int main()
{
int T;
cin>>T;
ll n,m,k;
while(T--)
{
scanf("%I64d%I64d%I64d",&n,&m,&k);
ll ans=1;
for(int i=0; i<k; i++)
{
cin>>pp[i];
mm[i]=lucas(n,m,pp[i]);
}
ans=ans*china(k,mm,pp);
cout<<ans<<endl;
}
}
例题2.hdu4773+普通CRT
模数为97*3xxxxx
#include<iostream>
#include<stdio.h>
using namespace std;
typedef long long ll;
//求C(n+len-1,len,p) n<=1e6 p=97*3xxxxx
ll mul(ll a, ll b, ll p) //快速乘,计算a*b%p
{
ll ret = 0;
while(b)
{
if(b & 1) ret = (ret + a) % p;
a = (a + a) % p;
b >>= 1;
}
return ret;
}
void ex_gcd(ll a, ll b, ll &x, ll &y, ll &d)
{
if (!b){
d = a, x = 1, y = 0;
}
else
{
ex_gcd(b, a % b, y, x, d);
y -= x * (a / b);
}
}
ll inv(ll t, ll p) //如果不存在,返回-1
{
ll d, x, y;
ex_gcd(t, p, x, y, d);
return d == 1 ? (x % p + p) % p : -1;
}
//n个方程:x=a[i](mod m[i]) (0<=i<n) m[i]两两互质
ll china(int n, ll *a, ll *m)
{
ll M = 1, ret = 0;
for(int i = 0; i < n; i ++) M *= m[i];
for(int i = 0; i < n; i ++)
{
ll w = M / m[i];
//ret = (ret + w * inv(w, m[i]) * a[i]) % M;
ret = (ret + mul(w * inv(w, m[i]), a[i], M)) % M;//防爆long long int
}
return (ret + M) % M;
}
ll exp_mod(ll a,ll b,ll p)
{
ll res=1;
while(b!=0)//利用二进制求乘方
{
if(b&1) res=(res*a)%p;
a=(a*a)%p;
b>>=1;
}
return res;
}
ll comb(ll a,ll b,ll p)
{
if(a<b) return 0;
if(a==b) return 1;
if(b>a-b) b=a-b;//组合的性质,n-m和m是相等的
ll ans=1,ca=1,cb=1;
for(ll i=0; i<b; ++i)
{
ca=(ca*(a-i))%p;//这里ca表示n!/(n-m)!
cb=(cb*(b-i))%p;
}
ans= (ca*exp_mod(cb,p-2,p))%p;//这里对m!求逆元即可
return ans;
}
ll lucas(ll n,ll m,ll p)
{
ll ans=1;
while(n&&m&&p)
{
ans=(ans*comb(n%p,m%p,p))%p;
n/=p;
m/=p;
}
return ans;
}
ll a[50];ll mm[3];ll pp[3];
int main()
{
ll MOD=364875103;pp[0]=97;pp[1]=3761599;
int T;cin>>T;
ll n,m,k;
int kase=0;
while(T--)
{
scanf("%I64d%I64d%I64d",&n,&m,&k);
a[k]=m;
for(int i=0; i<k; i++) {scanf("%I64d",&a[i]);}
ll ans=1;
for(int i=0; i<k; i++)
{
mm[0]=lucas(a[i+1]-a[i]+n-1,a[i+1]-a[i],pp[0]);
mm[1]=lucas(a[i+1]-a[i]+n-1,a[i+1]-a[i],pp[1]);
ans=ans*china(2,mm,pp);
ans%=MOD;
}
printf("Case #%d: ",++kase);
cout<<ans<<endl;
}
}
例题3.洛谷P4720 扩展CRT+lucas
求c(n,m)%p 1<=n,m<=1e18 p<=1e6的任意数
//C(n,m)%p n,m<=1e18 p<=1e6的任意数
#include<cstdio>
using namespace std;
typedef long long ll;
ll n,m,p,x0,y0;
void gcd(ll a,ll b,ll r){
if(!r){
x0=0;y0=1;return;
}
gcd(b,r,b%r);
ll tmp=x0;
x0=y0;
y0=tmp-a/b*y0;
}
ll inv(ll a,ll mod){
gcd(a,mod,a%mod);
x0=(x0%mod+mod)%mod;
return x0;
}
ll pow(ll a,ll x,ll mod){
if(x==1) return a;
if(x==0) return 1;
ll temp=pow(a,x/2,mod);
temp=temp*temp%mod;
if(x%2) temp=temp*a%mod;
return temp;
}
ll fac(ll a,ll x,ll mod){
if(!a) return 1;
ll res=1;
for(ll i=2;i<=mod;i++)
if(i%x) res=res*i%mod;
res=pow(res,a/mod,mod);
for(ll i=2;i<=a%mod;i++)
if(i%x) res=res*i%mod;
return res*fac(a/x,x,mod)%mod;
}
ll C(ll n,ll m,ll x,ll mod){
ll res=0,t=fac(n,x,mod),temp1=fac(m,x,mod),temp2=fac(n-m,x,mod);
for(ll i=n;i;i/=x) res+=i/x;
for(ll i=m;i;i/=x) res-=i/x;
for(ll i=n-m;i;i/=x) res-=i/x;
return t*inv(temp1,mod)%mod*inv(temp2,mod)%mod*pow(x,res,mod)%mod;
}
ll crt(ll x,ll mod){
ll temp=x*inv(p/mod,mod)%p*(p/mod)%p;
return temp;
}
ll lucas(ll n,ll m){//对每个p^k 进行扩展CRT合并
ll res=0,tmp=p,pk;
for(ll i=2;i*i<=p;i++){
pk=1;
if(tmp%i==0){
while(tmp%i==0){
tmp/=i;pk*=i;
}
res=(res+crt(C(n,m,i,pk),pk))%p;
}
if(tmp==1) break;
}
if(tmp>1) res=(res+crt(C(n,m,tmp,tmp),tmp))%p;
return res;
}
int main(){
scanf("%lld%lld%lld",&n,&m,&p);
printf("%lld",(lucas(n,m)%p+p)%p);
return 0;
}
例题4.BZOJ1485 线性筛预处理最小质因子求C(2n,n)%p
n<=1e6 ,p<=1e9的任意数
//求C(2n,n)/(n+1)%p n<=1e6 p<=1e9的任意数 扩展CRT TLE
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<queue>
#define SF scanf
#define PF printf
using namespace std;
typedef long long LL;
const int MAXN = 2000000;
int Div[MAXN+10], prime[MAXN+10], tot, MOD, n;
int cnt[MAXN+10];
bool vis[MAXN+10];
void init() {
for(int i = 2; i <= MAXN; i++) {
if(!vis[i]) prime[++tot] = i;
for(int j = 1; i * prime[j] <= MAXN; j++) {
Div[i*prime[j]] = prime[j];
vis[i*prime[j]] = true;
if(i % prime[j] == 0) break;
}
}
}
// h(n) = C(2n, n) / (n+1) n<=1e6 P<=1e9 任意数
int pow_mod(int x, int k) {
LL ans = 1;
while(k) {
if(k & 1) ans = 1LL * ans * x % MOD;
x = 1LL * x * x % MOD;
k >>= 1;
}
return ans;
}
int main() {
init();
SF("%d%d", &n, &MOD);
for(int i = 2; i <= n; i++) cnt[i] = -1;//分母2*3*4*......n
for(int i = n+2; i <= 2*n; i++) cnt[i] = 1;//分子n+2*n+3.....*2n
int ans = 1;
for(int i = n*2; i >= 2; i--)
if(!vis[i]) ans = 1LL * ans * pow_mod(i, cnt[i]) % MOD;//是质数
else cnt[Div[i]] += cnt[i], cnt[i/Div[i]] += cnt[i];//不是质数
PF("%d", ans);
}