luoguP4389 付公主的背包 生成函数+多项式exp

标签：生成函数，多项式 exp

题意：一共有 $n$ 个物品，每件物品体积为 $v_{i}$，数量为 $10^5$，求：装满一个容量为 $m$ 的背包有多少种不同方案.

比较友好的生成函数，并不是很难推导（似乎有些套路QAQ）

对于每一件物品列生成函数 $F_{i}(x)=\sum_{i=0}^{\infty} [i\%v=0]x^i$，则 $Ans=\prod_{i=1}^{n}F_{i}(x)$

对于 $F_{i}(x)$，有 $F_{i}(x)=\sum_{i=0}^{\infty} (x^v)^i=\frac{1}{1-x^v}$

所有有 $Ans=\prod_{i=1}^{n}\frac{1}{1-x^{v_{i}}}$

我们发现这个东西很难处理，因为这是 $n$ 个多项式相乘.

那么我们这个时候九有一个常用的套路：对多项式取 $\ln$，将乘法 $\Rightarrow $ 加法（因为多项式加法比较简单，而乘法涉及卷积）

$\ln Ans=\ln \prod_{i=1}^{n}\frac{1}{1-x^{v_{i}}}=\sum_{i=1}^{n}\ln (\frac{1}{1-x^{v_{i}}})$

然后又有一个套路：$\ln(1+x)=\sum_{i=1}^{\infty} (-1)^{i+1}\frac{x^i}{i}$ ，所以我们试着将 $\ln$ 中的分数提出.

$\Rightarrow \ln \frac{1}{1-x^{v_{i}}}=-\ln (1-x^{v_{i}})$

$\Rightarrow \ln Ans=-\sum_{i=1}^{n}\ln (1-x^{v_{i}})$

$\Rightarrow \ln(Ans)=\sum_{j=1}^{n}\sum_{i=1}^{\infty}\frac{x^{v_{j}i}}{i}$

利用调和级数来记录每个物品的贡献，然后最后再求一个多项式 $\exp$ 就好了.

code:

#include <cmath> 
#include <cstring> 
#include <algorithm> 
#include <cstdio> 
#include <string>   
#define ll long long
#define ull unsigned long long    
using namespace std; 
namespace IO 
{ 
    char buf[100000],*p1,*p2;
    #define nc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++)
    int rd()
    {
        int x=0; char s=nc();
        while(s<'0') s=nc();
        while(s>='0') x=(((x<<2)+x)<<1)+s-'0',s=nc();
        return x;
    }        
    void print(int x) {if(x>=10) print(x/10);putchar(x%10+'0');}
    void setIO(string s) 
    {
        string in=s+".in"; 
        string out=s+".out";   
        freopen(in.c_str(),"r",stdin);  
        // freopen(out.c_str(),"w",stdout); 
    }    
};     
const int G=3; 
const int N=2000005;  
const int mod=998244353;                  
int A[N],B[N],w[2][N],mem[N*100],*ptr=mem,inv[N],tmpa[N],tmpb[N],aa[N],bb[N];         
inline int qpow(int x,int y)
{
    int tmp=1;    
    for(;y;y>>=1,x=(ll)x*x%mod)     if(y&1) tmp=(ll)tmp*x%mod; 
    return tmp;   
}     
inline int INV(int a) { return qpow(a,mod-2); }       
inline void ntt_init(int len)
{
    int i,j,k,mid,x,y;     
    w[1][0]=w[0][0]=1,x=qpow(3,(mod-1)/len),y=qpow(x,mod-2);
    for (i=1;i<len;++i) w[0][i]=(ll)w[0][i-1]*x%mod,w[1][i]=(ll)w[1][i-1]*y%mod;        
}
void NTT(int *a,int len,int flag)
{
    int i,j,k,mid,x,y;              
    for(i=k=0;i<len;++i)
    {
        if(i>k)    swap(a[i],a[k]); 
        for(j=len>>1;(k^=j)<j;j>>=1); 
    }  
    for(mid=1;mid<len;mid<<=1)           
        for(i=0;i<len;i+=mid<<1)
            for(j=0;j<mid;++j)         
            {
                x=a[i+j], y=(ll)w[flag==-1][len/(mid<<1)*j]*a[i+j+mid]%mod; 
                a[i+j]=(x+y)%mod; 
                a[i+j+mid]=(x-y+mod)%mod;  
            }  
    if(flag==-1) 
    {
        int rev=INV(len);  
        for(i=0;i<len;++i)    a[i]=(ll)a[i]*rev%mod;  
    }
}             
inline void getinv(int *a,int *b,int len,int la)          
{
    if(len==1) { b[0]=INV(a[0]);   return; }
    getinv(a,b,len>>1,la);   
    int l=len<<1,i;  
    memset(A,0,l*sizeof(A[0]));             
    memset(B,0,l*sizeof(A[0]));   
    memcpy(A,a,min(la,len)*sizeof(a[0]));                                                     
    memcpy(B,b,len*sizeof(b[0]));            
    ntt_init(l);  
    NTT(A,l,1),NTT(B,l,1);     
    for(i=0;i<l;++i)  A[i]=((ll)2-(ll)A[i]*B[i]%mod+mod)*B[i]%mod;
    NTT(A,l,-1);                                
    memcpy(b,A,len<<2);         
}        
void get_dao(int *a,int *b,int len) 
{
    for(int i=1;i<len;++i) b[i-1]=(ll)i*a[i]%mod;    
    b[len-1]=0;  
}                    
void get_jifen(int *a,int *b,int len) 
{
    for(int i=1;i<len;++i) b[i]=(ll)INV(i)*a[i-1]%mod;   
    b[0]=0;  
}
void get_ln(int *a,int *b,int len,int la) 
{
    int l=len<<1,i;   
    memset(tmpa,0,l<<2); 
    memset(tmpb,0,l<<2);   
    get_dao(a,tmpa,min(len,la));    
    getinv(a,tmpb,len,la);  
    ntt_init(l);          
    NTT(tmpa,l,1),NTT(tmpb,l,1);  
    for(i=0;i<l;++i) tmpa[i]=(ll)tmpa[i]*tmpb[i]%mod;    
    NTT(tmpa,l,-1);   
    get_jifen(tmpa,b,len);      
}       
void get_exp(int *a,int *b,int len,int la) 
{
    if(len==1) { b[0]=1; return; }                           
    int l=len<<1,i; 
    get_exp(a,b,len>>1,la);               
    for(i=0;i<l;++i)  aa[i]=bb[i]=0;  
    for(i=0;i<(len>>1);++i) aa[i]=b[i];            
    get_ln(b,bb,len,len>>1);                                               
    for(i=0;i<len;++i) bb[i]=(ll)(mod-bb[i]+(i>la?0:a[i]))%mod;                               
    bb[0]=(bb[0]+1)%mod;   
    ntt_init(l);   
    NTT(aa,l,1),NTT(bb,l,1);   
    for(i=0;i<l;++i) aa[i]=(ll)aa[i]*bb[i]%mod;     
    NTT(aa,l,-1);      
    for(i=0;i<len;++i)  b[i]=aa[i];     
}
struct poly
{
    int len,*a;   
    poly(){}      
    poly(int l) {len=l,a=ptr,ptr+=l; }           
    inline void rev() { reverse(a,a+len); }      
    inline void fix(int l) {len=l,a=ptr,ptr+=l;}  
    inline void get_mod(int l) { for(int i=l;i<len;++i) a[i]=0;  len=l;  }
    inline poly dao()
    {       
        poly re(len-1);  
        for(int i=1;i<len;++i)  re.a[i-1]=(ll)i*a[i]%mod;        
        return re;   
    }   
    inline poly jifen() 
    {
        poly c;   
        c.fix(len+1);       
        c.a[0]=0;
        for(int i=1;i<=len;++i) c.a[i]=(ll)a[i-1]*INV(i)%mod;               
        return c;  
    }          
    inline poly Inv(int l)
    { 
        poly b(l);                      
        getinv(a,b.a,l,len);                                 
        return b;                       
    }                  
    inline poly ln(int l) 
    {
        int lim=1; 
        while(lim<=l) lim<<=1;                     
        poly b(lim);    
        get_ln(a,b.a,lim,len); 
        return b;    
    }                         
    inline poly exp(int l) 
    {   
        int lim=1;  
        while(lim<=l) lim<<=1;   
        poly b(lim);  
        get_exp(a,b.a,lim,len);        
        return b;  
    }                               
    inline poly operator*(const poly &b) const
    {
        poly c(len+b.len-1);  
        if(c.len<=500)
        {        
            for(int i=0;i<len;++i)  
                if(a[i])   for(int j=0;j<b.len;++j)  c.a[i+j]=(c.a[i+j]+(ll)(a[i])*b.a[j])%mod;     
            return c;
        }
        int n=1;   
        while(n<(len+b.len)) n<<=1;
        memset(A,0,n<<2); 
        memset(B,0,n<<2);  
        memcpy(A,a,len<<2);                            
        memcpy(B,b.a,b.len<<2);                 
        ntt_init(n);       
        NTT(A,n,1), NTT(B,n,1);    
        for(int i=0;i<n;++i) A[i]=(ll)A[i]*B[i]%mod;  
        NTT(A,n,-1);  
        memcpy(c.a,A,c.len<<2); 
        return c;      
    }   
    poly operator+(const poly &b) const
    {
        poly c(max(len,b.len));   
        for(int i=0;i<c.len;++i)  c.a[i]=((i<len?a[i]:0)+(i<b.len?b.a[i]:0))%mod;  
        return c;   
    }
    poly operator-(const poly &b) const
    {   
        poly c(len);      
        for(int i=0;i<len;++i)  
        {
            if(i>=b.len)   c.a[i]=a[i]; 
            else c.a[i]=(a[i]-b.a[i]+mod)%mod;   
        }
        return c; 
    }
    poly operator/(poly u)
    { 
        int n=len,m=u.len,l=1; 
        while(l<(n-m+1)) l<<=1;                          
        rev(),u.rev();           
        poly v=u.Inv(l);   
        v.get_mod(n-m+1);       
        poly re=(*this)*v;  
        rev(),u.rev();   
        re.get_mod(n-m+1);        
        re.rev(); 
        return re;  
    }     
    poly operator%(poly u)
    {     
        poly re=(*this)-u*(*this/u);       
        re.get_mod(u.len-1);      
        return re;   
    }                    
}p[N<<2],pr;     
// 插值部分 
namespace Inter
{ 
    int xx[N],yy[N];             
    #define lson now<<1 
    #define rson now<<1|1          
    inline void pushup(int l,int r,int now)
    {
        int mid=(l+r)>>1;     
        if(r>mid)   p[now]=p[lson]*p[rson];
        else p[now]=p[lson];  
    }
    void build(int l,int r,int now,int *pp)
    {
        if(l==r)
        {    
            p[now].fix(2); 
            p[now].a[0]=mod-pp[l]; 
            p[now].a[1]=1;  
            return;
        } 
        int mid=(l+r)>>1;  
        if(l<=mid)  build(l,mid,lson,pp);    
        if(r>mid)   build(mid+1,r,rson,pp);         
        p[now]=p[lson]*p[rson];  
    }   
    // 多项式多点求值 
    void get_val(int l,int r,int now,poly b,int *pp,int *t)
    {
        if(b.len<=500)    
        {  
            for(int i=l;i<=r;++i)
            {
                ull s=0;            
                for(int j=b.len-1;j>=0;--j)    
                {
                    s=((ull)s*pp[i]+b.a[j])%mod; 
                    if(!(j&7))   s%=mod;      
                }
                t[i]=s%mod;  
            }
            return; 
        }
        int mid=(l+r)>>1;    
        if(l<=mid)   get_val(l,mid,lson,b%p[lson],pp,t); 
        if(r>mid)    get_val(mid+1,r,rson,b%p[rson],pp,t);    
    }  
    // 多项式快速插值   
    poly solve_polate(int l,int r,int now,int *t)
    {
        if(l==r)
        {
            poly re(1);  
            re.a[0]=t[l];  
            return re;  
        }
        int mid=(l+r)>>1;   
        poly L,R; 
        L=solve_polate(l,mid,lson,t);  
        R=solve_polate(mid+1,r,rson,t);  
        return L*p[rson]+R*p[lson];          
    }         
};           
int bu[N];  
int main() 
{
    // IO::setIO("input");        
    int i,j,n,m; 
    scanf("%d%d",&n,&m); 
    for(i=1;i<=n;++i)
    {
        int x; 
        scanf("%d",&x),++bu[x]; 
    }                 
    poly b(m+1);     
    for(i=1;i<=m;++i) 
    {
        if(bu[i]) 
        {
            for(j=i;j<=m;j+=i) 
            {    
                j/i;     
                b.a[j]=(ll)(b.a[j]+(ll)bu[i]*INV(j/i)%mod)%mod;        
            }
        }
    }
    b=b.exp(m+1);     
    for(i=1;i<=m;++i)  printf("%d\n",b.a[i]);  
    return 0;    
}