[BZOJ3129][Sdoi2013]方程（容斥原理+扩展lucas）

题目描述

传送门

题解

这题刚开始sb了
首先所有的强制选1个
n1+1~n2的限制就再强制选ai-1个
剩余的需要容斥一下，答案=至少0个超过限制的-至少1个超过限制的+至少2个超过限制的…
$2^{n1}$枚举那些限制强制超过，强制超过的就是要再强制选ai个
然后最后剩下了一个m，x1+x2+…+xn=m，可以利用插板法得出答案极为$C_{m+n-1}^{n-1}$
扩展lucas好久不打又快忘干净了…

段子：
（由于我事先已经知道了这题要用到扩展lucas）
我：诶诶诶ATP你看这题又没有说什么质因子小于$10^5$什么的为啥能用扩展lucas啊？
ATP：诶呀你不用管它反正它不会不让你过哒！
我：…

代码

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
#define LL long long
#define N 1000005

int T,Mod,n,n1,n2,m;
int a[N],p[N],prime[N];
LL ans;

LL fast_pow(LL a,int p,int Mod)
{
    LL ans=1;
    for (;p;p>>=1,a=a*a%Mod)
        if (p&1)
            ans=ans*a%Mod;
    return ans;
}
LL exgcd(LL a,LL b,LL &x,LL &y)
{
    if (!b) x=1,y=0;
    else exgcd(b,a%b,y,x),y-=a/b*x;
}
LL inv(LL a,LL b)
{
    if (!a) return 0;
    LL x=0,y=0;
    exgcd(a,b,x,y);
    x=(x%b+b)%b;
    if (!x) x+=b;
    return x;
}
LL Mul(int n,int pi,int pk)
{
    if (!n) return 1;
    LL ans=1LL;
    if (n/pk)
    {
        for (int i=2;i<=pk;++i)
            if (i%pi) ans=ans*(LL)i%pk;
        ans=fast_pow(ans,n/pk,pk);
    }
    for (int i=2;i<=n%pk;++i)
        if (i%pi) ans=ans*(LL)i%pk;
    return ans*Mul(n/pi,pi,pk)%pk;
}
LL C(int n,int m,int pi,int pk)
{
    if (m>n) return 0;
    LL a=Mul(n,pi,pk),b=Mul(m,pi,pk),c=Mul(n-m,pi,pk);
    int k=0;
    for (int i=n;i;i/=pi) k+=i/pi;
    for (int i=m;i;i/=pi) k-=i/pi;
    for (int i=n-m;i;i/=pi) k-=i/pi;
    LL ans=a*inv(b,pk)%pk*inv(c,pk)%pk*fast_pow(pi,k,pk)%pk;
    return ans*(Mod/pk)%Mod*inv(Mod/pk,pk)%Mod;
}
LL exlucas(LL n,LL m)
{
    LL ans=0;
    for (int i=1;i<=p[0];++i)
        ans=(ans+C(n,m,p[i],prime[i]))%Mod;
    return ans;
}
LL solve(int state)
{
    int opt=0,now=m;
    for (int i=0;i<n1;++i)
        if ((state>>i)&1)
            ++opt,now-=a[i+1];
    if (now<0) return 0;
    if (opt&1) opt=-1;
    else opt=1;
    return exlucas(now+n-1,n-1)*opt;
}
int main()
{
    scanf("%d%d",&T,&Mod);
    int x=Mod;
    for (int i=2;i*i<=x;++i)
        if (x%i==0)
        {
            p[++p[0]]=i;prime[p[0]]=1;
            while (x%i==0) prime[p[0]]*=i,x/=i;
        }
    if (x>1) p[++p[0]]=x,prime[p[0]]=x;
    while (T--)
    {
        scanf("%d%d%d%d",&n,&n1,&n2,&m);
        m-=n;
        for (int i=1;i<=n1+n2;++i) scanf("%d",&a[i]);
        for (int i=n1+1;i<=n1+n2;++i) m-=a[i]-1;
        if (m<0)
        {
            puts("0");
            continue;
        }
        ans=0;
        for (int i=0;i<1<<n1;++i)
            ans=(ans+solve(i))%Mod;
        ans=(ans%Mod+Mod)%Mod;
        printf("%lld\n",ans);
    }
}