HDU 5519 Kykneion asma （2015 ICPC 沈阳 K）状压dp+容斥

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=5519

题意

给你n(<=15000)，再给你0~4这五个数字的可用数量a[i](<=30000)，你需要用这些数字构造长度为n的序列，不能有前导零，求合法的方案数。

分析

网上有多种解法，主要如下：

①生成函数FFT（窝不会）：https://blog.csdn.net/Quack_quack/article/details/50748753?utm_source=blogxgwz4

②另一种做法是三维数组状压dp+容斥：https://blog.csdn.net/u014610830/article/details/49816237

③还有一种非常神奇的状压dp+容斥，只用了两维数组（然鹅窝没看懂，如果有人看懂了欢迎给我讲讲）：https://www.cnblogs.com/myx12345/p/10033413.html

窝只是题解的搬运工

窝只看懂了第二种：

设f(n,a0,a1,a2,a3,a4)为用a0~a4个对应的数字构造长度为n的序列的总数（注意这里可以有前导零）

那么我们要求的就是f(n,a0,a1,a2,a3,a4)-f(n-1,a0-1,a1,a2,a3,a4)（a0=0时特判）

对于f，我们直接求会很麻烦，考虑求出不合法的方案数，然后用5^n-不合法的方案数（记为g）。

不合法，即至少有一个数字超出了限制。那么g就可以用容斥求出来了：

g(n,a0,a1,a2,a3,a4)=g（一个数字超出限制）- g（两个数字超出限制）+ g（三个数字超出限制）- g（四个数字超出限制）+ g（五个数字超出限制）

因此设置dp数组时，要设置三维，dp[i][state][j]表示构造i位，超出限制的数字状态为state，最终j位超出限制的方案数（j影响着我们后面的容斥）

转移有两种：

①dp[i][state][j]=dp[i-1][state][j]*(5-j-cnt[state])，注意最终超出范围但状态里还未超出的数字不能选。

②dp[i][state][j]=sum(dp[i-a[k]-1][state^(1<<k)][j]*C(i-1,a[i])，k=0,1,...,4

相当于枚举这一位的数字，然后直接让它超出限制，只需要从前i-1位中挑选a[i]位插入该数，然后其他数相对位置不变。

最后，令tmp=5^n，再减去不合法的方案数即为答案。 不合法的方案数，对超出数字的奇偶个数进行容斥即可。

代码：

#include<bits/stdc++.h>
#define ll long long
#define inf 0x3f3f3f3f
#define mst(head,x,n) memset(head+1,x,n*sizeof(head[0]))
#define rep(i,a,b) for(int i=(a);i<=(b);i++)
#define dep(i,a,b) for(int i=(a);i>=(b);i--)
using namespace std;
const int maxn=4e4+5;
//const double pi=acos(-1.0);
//const double eps=1e-9;
const ll mo=1e9+7;
int n,m,q,k,flag,x,f,y,p;
long long ni[maxn];
long long a[maxn];
ll c[maxn],ans,tmp;
ll dp[15005][35][6];
ll cnt[maxn];
long long zuhe(int x,int y){  //组合数
    return a[x]*ni[y]%mo*ni[x-y]%mo;
}
long long calc(long long x,long long y){
    long long z=1;
    while (y){
        if (y&1)(z*=x)%=mo;
        (x*=x)%=mo,y/=2;
    }
    return z;
}
template <typename T>
inline void read(T &X){
    X=0;int w=0; char ch=0;
    while(!isdigit(ch)) {w|=ch=='-';ch=getchar();}
    while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();
    if(w) X=-X;
}
ll solve(){
    memset(dp,0,sizeof(dp));
    int mm=31;
    rep(i,1,5) dp[0][0][i]=1;
    rep(i,1,n){
        rep(j,1,5){
            rep(st,0,mm) if(cnt[st]<=j){
                dp[i][st][j]=(dp[i][st][j]+dp[i-1][st][j]*(5-j+cnt[st])%mo)%mo;
                    rep(k,0,4)if(st&(1<<k)&&c[k]<i){
                        dp[i][st][j]=(dp[i][st][j]+dp[i-c[k]-1][st^(1<<k)][j]*zuhe(i-1,c[k])%mo)%mo;
                    }
            }
        }
    }
    ll tmp=calc(5LL,n);
    rep(i,1,mm){
        if(cnt[i]&1){
            tmp=(tmp-dp[n][i][cnt[i]]+mo)%mo;
        }
        else {
            tmp=(tmp+dp[n][i][cnt[i]]+mo)%mo;
        }
    }
    return tmp;
}
int main(){

#ifdef ONLINE_JUDGE
#else
    freopen("D:/Temp/in.txt", "r", stdin);
#endif
    a[0]=1;
    for (int i=1;i<maxn;i++)a[i]=a[i-1]*i%mo;
    ni[maxn-1]=calc(a[maxn-1],mo-2);
    for (int i=maxn-2;i>=0;i--)
    ni[i]=ni[i+1]*(i+1)%mo;
    rep(i,0,maxn-2) cnt[i]=cnt[i>>1]+(i&1);
    int T,cas=1;
    read(T);
    while(T--)
    {
        read(n);
        rep(i,0,4) read(c[i]);
        if(!c[0]){
            ans=solve();
        }
        else {
            ans=solve();
            c[0]--;n--;
            ans=(ans-solve()+mo)%mo;
        }
        printf("Case #%d: %lld\n",cas++,ans);
    }
    return 0;
}