BZOJ1042 HAOI2008硬币购物（任意模数NTT+多项式求逆+生成函数/容斥原理+动态规划）

　　第一眼生成函数。四个等比数列形式的多项式相乘，可以化成四个分式。其中分母部分是固定的，可以多项式求逆预处理出来。而分子部分由于项数很少，询问时2^4算一下贡献就好了。这个思路比较直观。只是常数巨大，以及需要敲一发类似任意模数ntt的东西来避免爆精度。成功以这种做法拿下luogu倒数rank1，至于bzoj不指望能过了。

#include<iostream> 
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<iomanip>
using namespace std;
int read()
{
    int x=0,f=1;char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
    while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
#define N 550000
#define T 100000
#define P1 998244353
#define P2 1004535809
int r[N],c1,c2,c3,c4,tot,d1,d2,d3,d4,s,t;
int a[N],b[N],c[N],e[2][N];
long long f[N];
void inc(int &x,int P){x++;if (x>=P) x-=P;}
void dec(int &x,int P){x--;if (x<0) x+=P;}
int ksm(int a,int k,int P)
{
    if (a==0) return 0;
    if (k==0) return 1;
    int tmp=ksm(a,k>>1,P);
    if (k&1) return 1ll*tmp*tmp%P*a%P;
    else return 1ll*tmp*tmp%P;
}
long long ksc(long long a,long long b,long long P)
{
    long long t=a*b-(long long)((long double)a*b/P+0.5)*P;
    return t<0?t+P:t;
}
void DFT(int n,int *a,int p,int P)
{
    for (int i=0;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]);
    for (register int i=2;i<=n;i<<=1)
    {
        int wn=ksm(p,(P-1)/i,P);
        for (register int j=0;j<n;j+=i)
        {
            int w=1;
            for (register int k=j;k<j+(i>>1);k++,w=1ll*w*wn%P)
            {
                int x=a[k],y=1ll*w*a[k+(i>>1)]%P;
                a[k]=(x+y)%P,a[k+(i>>1)]=(x-y+P)%P;
            }
        }
    }
}
void mul(int n,int *a,int *b,int P,int inv3)
{
    DFT(n,a,3,P),DFT(n,b,3,P);
    for (int i=0;i<n;i++) a[i]=1ll*a[i]*(P+2-1ll*a[i]*b[i]%P)%P;
    DFT(n,a,inv3,P);
    int inv=ksm(n,P-2,P);
    for (int i=0;i<n;i++) a[i]=1ll*a[i]*inv%P;
}
void solve(int P,int inv3,int op)
{
    memset(a,0,sizeof(a));
    memset(b,0,sizeof(b));
    memset(c,0,sizeof(c));
    if (c1+c2+c3+c4<=T) inc(a[c1+c2+c3+c4],P);
    if (c1+c2+c3<=T) dec(a[c1+c2+c3],P);
    if (c1+c2+c4<=T) dec(a[c1+c2+c4],P);
    if (c4+c2+c3<=T) dec(a[c4+c2+c3],P);
    if (c1+c4+c3<=T) dec(a[c1+c4+c3],P);
    if (c1+c2<=T) inc(a[c1+c2],P);
    if (c1+c3<=T) inc(a[c1+c3],P);
    if (c1+c4<=T) inc(a[c1+c4],P);
    if (c2+c3<=T) inc(a[c2+c3],P);
    if (c4+c2<=T) inc(a[c4+c2],P);
    if (c3+c4<=T) inc(a[c3+c4],P);
    dec(a[c1],P);dec(a[c2],P);dec(a[c3],P);dec(a[c4],P);
    inc(a[0],P);
    t=1;b[0]=1;
    while (t<=T)
    {
        t<<=1;
        for (int i=0;i<t;i++) c[i]=a[i];
        for (int i=0;i<(t<<1);i++) r[i]=(r[i>>1]>>1)|(i&1)*t;
        mul(t<<1,b,c,P,inv3);
        for (int i=t;i<(t<<1);i++) b[i]=0;
    }
    memcpy(e[op],b,sizeof(e[op]));
}
void crt()
{
    long long P=1ll*P1*P2,inv1=ksm(P2%P1,P1-2,P1),inv2=ksm(P1%P2,P2-2,P2);
    for (int i=0;i<=T;i++)
    f[i]=(ksc(1ll*e[0][i]*P2%P,inv1,P)+ksc(1ll*e[1][i]*P1%P,inv2,P))%P;
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("bzoj1042.in","r",stdin);
    freopen("bzoj1042.out","w",stdout);
    const char LL[]="%I64d\n";
#else
    const char LL[]="%lld\n";
#endif
    c1=read(),c2=read(),c3=read(),c4=read(),tot=read();
    solve(P1,332748118,0);
    solve(P2,334845270,1);
    crt();
    while (tot--)
    {
        d1=read(),d2=read(),d3=read(),d4=read(),s=read();
        d1=min(1ll*s+1,1ll*(d1+1)*c1);
        d2=min(1ll*s+1,1ll*(d2+1)*c2);
        d3=min(1ll*s+1,1ll*(d3+1)*c3);
        d4=min(1ll*s+1,1ll*(d4+1)*c4);
        long long ans=f[s];
        if (d1+d2+d3+d4<=s) ans+=f[s-(d1+d2+d3+d4)];
        if (d1+d2+d3<=s) ans-=f[s-(d1+d2+d3)];
        if (d1+d2+d4<=s) ans-=f[s-(d1+d2+d4)];
        if (d4+d2+d3<=s) ans-=f[s-(d4+d2+d3)];
        if (d1+d4+d3<=s) ans-=f[s-(d1+d4+d3)];
        if (d1+d2<=s) ans+=f[s-(d1+d2)];
        if (d1+d3<=s) ans+=f[s-(d1+d3)];
        if (d1+d4<=s) ans+=f[s-(d1+d4)];
        if (d2+d3<=s) ans+=f[s-(d2+d3)];
        if (d4+d2<=s) ans+=f[s-(d4+d2)];
        if (d3+d4<=s) ans+=f[s-(d3+d4)];
        if (d1<=s) ans-=f[s-d1];
        if (d2<=s) ans-=f[s-d2];
        if (d3<=s) ans-=f[s-d3];
        if (d4<=s) ans-=f[s-d4];
        printf(LL,ans);
    }
    return 0;
}

　　还有一种更优秀的做法。考虑如果硬币没有个数限制的话，就是一个完全背包。添加限制可以想到容斥。我们枚举有哪几种硬币超过了个数限制，就可以容斥斥斥容容容斥把多重背包转化成完全背包了。

 

#include<iostream> 
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<iomanip>
using namespace std;
int read()
{
    int x=0,f=1;char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
    while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
#define N 100010
#define ll long long
int c[4],t,d[4],s;
ll f[N],ans;
int calc(int k,int x){if (k<x) return 0;else return f[k-x];}
void dfs(int k,int sum,ll tot)
{
    if (tot>s) return;
    if (k==4) {ans+=((sum&1)?-1:1)*f[s-tot];return;} 
    dfs(k+1,sum+1,tot+1ll*(d[k]+1)*c[k]);
    dfs(k+1,sum,tot);
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("bzoj1042.in","r",stdin);
    freopen("bzoj1042.out","w",stdout);
    const char LL[]="%I64d\n";
#else
    const char LL[]="%lld\n";
#endif
    for (int i=0;i<4;i++) c[i]=read();
    t=read();
    f[0]=1;
    for (int i=0;i<4;i++)
        for (int j=c[i];j<=N-10;j++)
        f[j]+=f[j-c[i]];
    while (t--)
    {
        for (int i=0;i<4;i++) d[i]=read();
        s=read();
        ans=0;
        dfs(0,0,0);
        printf(LL,ans);
    }
    return 0;
}

 

 　　仔细考虑一下会发现两个做法本质上其实是一样的。分子部分所乘的多项式就是一个容斥的过程，而求逆所得的结果就是完全背包。