洛谷P4389 付公主的背包 [生成函数,NTT]

传送门


同样是回过头来发现不会做了,要加深一下记忆。


思路

只要听说过生成函数的人相信第一眼都可以想到生成函数。

所以我们要求

\[ans=\prod \sum_n x^{nV}=\prod \frac{1}{1-x^V} \]

也就是\(\prod (1-x^V)\)

但这玩意好像还是不会做,怎么办呢?

按照套路,可以先\(\ln\)一下,加起来,再\(\exp\)回去。

所以现在要求

\[\sum \ln(1-x^V) \]

……

……

……

不会。

不会怎么办?

打表找规律!

经过打表,可以发现\(\ln(1-x^V)=\sum \frac{-1}{n}x^{nV}\)

然后开个桶记录每种体积有多少,\(n\log n\)加一下系数,就做完了。


代码

#include<bits/stdc++.h>
namespace my_std{
    using namespace std;
    #define pii pair<int,int>
    #define fir first
    #define sec second
    #define MP make_pair
    #define rep(i,x,y) for (int i=(x);i<=(y);i++)
    #define drep(i,x,y) for (int i=(x);i>=(y);i--)
    #define go(x) for (int i=head[x];i;i=edge[i].nxt)
    #define sz 404040
    #define mod 998244353
    typedef long long ll;
    template<typename T>
    inline void read(T& t)
    {
        t=0;char f=0,ch=getchar();
        double d=0.1;
        while(ch>'9'||ch<'0') f|=(ch=='-'),ch=getchar();
        while(ch<='9'&&ch>='0') t=t*10+ch-48,ch=getchar();
        if(ch=='.')
        {
            ch=getchar();
            while(ch<='9'&&ch>='0') t+=d*(ch^48),d*=0.1,ch=getchar();
        }
        t=(f?-t:t);
    }
    template<typename T,typename... Args>
    inline void read(T& t,Args&... args){read(t); read(args...);}
    void file()
    {
        #ifndef ONLINE_JUDGE
        freopen("a.txt","r",stdin);
        #endif
    }
//	inline ll mul(ll a,ll b){ll d=(ll)(a*(double)b/mod+0.5);ll ret=a*b-d*mod;if (ret<0) ret+=mod;return ret;}
}
using namespace my_std;

ll ksm(ll x,int y)
{
    ll ret=1;
    for (;y;y>>=1,x=x*x%mod) if (y&1) ret=ret*x%mod;
    return ret;
}
ll inv(ll x){return ksm(x,mod-2);}

int limit,r[sz];
void NTT_init(int n)
{
    limit=1;int l=-1;
    while (limit<=n+n) limit<<=1,++l;
    rep(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<l);
}
void NTT(ll *a,int type)
{
    rep(i,0,limit-1) if (i<r[i]) swap(a[i],a[r[i]]);
    rep(i,0,limit-1) a[i]%=mod;
    for (int mid=1;mid<limit;mid<<=1)
    {
        ll Wn=ksm(3,(mod-1)/mid>>1);if (type==-1) Wn=inv(Wn);
        for (int j=0,len=mid<<1;j<limit;j+=len)
        {
            ll w=1;
            for (int k=0;k<mid;k++,w=w*Wn%mod)
            {
                ll x=a[j+k],y=a[j+k+mid]*w;
                a[j+k]=(x+y)%mod;a[j+k+mid]=(1ll*mod*mod-y+x)%mod;
            }
        }
    }
    if (type==1) return;
    ll I=inv(limit);
    rep(i,0,limit-1) a[i]=a[i]*I%mod;
}
ll tmp1[sz],tmp2[sz],tmp3[sz],tmp4[sz];
void PolyInv(ll *a,ll *f,int n) // f=a^{-1}
{
    if (n==1) return (void)(f[0]=inv(a[0]));
    int mid=(n+1)>>1;
    PolyInv(a,f,mid);
    NTT_init(n);
    rep(i,0,mid-1) tmp1[i]=f[i];rep(i,mid,limit-1) tmp1[i]=0;
    rep(i,0,n-1) tmp2[i]=a[i];rep(i,n,limit-1) tmp2[i]=0;
    NTT(tmp1,1);NTT(tmp2,1);
    rep(i,0,limit-1) tmp1[i]=tmp1[i]*(mod+2-tmp1[i]*tmp2[i]%mod)%mod;
    NTT(tmp1,-1);
    rep(i,0,n-1) f[i]=tmp1[i];
    rep(i,0,limit-1) tmp1[i]=tmp2[i]=0;
}
void Derivative(ll *a,ll *b,int n){rep(i,0,n-2) b[i]=a[i+1]*(i+1)%mod;b[n-1]=0;}
void Integrate(ll *a,ll *b,int n){drep(i,n-1,1) b[i]=a[i-1]*inv(i)%mod;b[0]=0;}
void PolyLn(ll *a,ll *f,int n) // f=ln a
{
    NTT_init(n);
    PolyInv(a,tmp3,n);Derivative(a,tmp4,n);
    NTT(tmp3,1);NTT(tmp4,1);
    rep(i,0,limit-1) tmp1[i]=tmp3[i]*tmp4[i]%mod;
    NTT(tmp1,-1);
    Integrate(tmp1,f,n);
    rep(i,0,limit-1) tmp1[i]=tmp3[i]=tmp4[i]=0;
}
void PolyExp(ll *a,ll *f,int n)
{
    if (n==1) return (void)(f[0]=1);
    int mid=(n+1)>>1;
    PolyExp(a,f,mid);
    rep(i,mid,n-1) f[i]=0;
    PolyLn(f,tmp2,n);
    rep(i,0,n-1) tmp1[i]=f[i];
    rep(i,0,n-1) tmp2[i]=(a[i]-tmp2[i]+mod)%mod;
    ++tmp2[0];
    NTT_init(n);
    NTT(tmp1,1);NTT(tmp2,1);
    rep(i,0,limit-1) tmp1[i]=tmp1[i]*tmp2[i]%mod;
    NTT(tmp1,-1);
    rep(i,0,n-1) f[i]=tmp1[i];
    rep(i,0,limit-1) tmp1[i]=tmp2[i]=0;
}

int n,m;
int cnt[sz];
ll f[sz],ans[sz];

int main()
{
    file();
    read(n,m);
    int x;
    rep(i,1,n) read(x),++cnt[x];
    rep(i,1,m)
        if (cnt[i])
            for (int j=i;j<=m;j+=i)
                (f[j]+=cnt[i]*inv(j/i)%mod)%=mod;
    PolyExp(f,ans,m+1);
    rep(i,1,m) printf("%lld\n",ans[i]);
}
posted @ 2019-03-13 17:06  p_b_p_b  阅读(260)  评论(0编辑  收藏  举报