组合数取模小结

最新推荐文章于 2020-05-27 23:06:03 发布

原创最新推荐文章于 2020-05-27 23:06:03 发布 · 388 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#组合数学

ACM--组合数学专栏收录该内容

27 篇文章

订阅专栏

本文介绍了如何使用中国剩余定理（CRT）与Lucas定理解决组合数问题，并提供了四个具体的编程实例，涵盖不同的CRT应用及优化技巧。

非常好的学习博客：https://blog.youkuaiyun.com/skywalkert/article/details/52553048

下图转至：https://blog.youkuaiyun.com/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);
}