LOJ#2320. 「清华集训 2017」生成树计数 prufer序列+生成函数+NTT

好神仙的多项式啊！

需要用到：prufer序列，生成函数，多项式求逆，多项式取 ln，多项式 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=200005;
const int mod=998244353;              
int A[N],B[N],w[2][N],mem[N*100],*ptr=mem,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)
    {            
        int lim=1;
        while(lim<=l) lim<<=1;     
        poly b(lim);
        getinv(a,b.a,lim,len);  
        b.get_mod(l);                                  
        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);
        b.get_mod(l);     
        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);    
        b.get_mod(l); 
        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;
    }                
}po,_A,Ai,Bi;                      
poly solve(int l,int r,int *ar)
{
    if(l==r)
    {  
        poly c;   
        c.fix(2); 
        c.a[0]=1; 
        c.a[1]=mod-ar[l];     
        return c;
    }
    int mid=(l+r)>>1;  
    return solve(l,mid,ar)*solve(mid+1,r,ar); 
}
#define MAX 30002             
int n,m;   
int al[N],fac[N],inv[N];           
void Initialize() 
{
    int i,j;       
    inv[0]=fac[0]=1;  
    for(i=1;i<=MAX;++i) fac[i]=(ll)fac[i-1]*i%mod,inv[i]=INV(fac[i]);
}
void calc_sum() 
{
    int i,j;                    
    _A=solve(1,n,al);              
    _A=_A.dao()*_A.Inv(MAX);                           
    poly tmp(2);  
    tmp.a[0]=0,tmp.a[1]=1;            
    _A=_A*tmp; 
    for(i=0;i<MAX;++i) _A.a[i]=mod-_A.a[i];   
    _A.a[0]=(ll)(_A.a[0]+n)%mod;    
}                                 
int main()
{
    int i,j;  
    // IO::setIO("input");   
    scanf("%d%d",&n,&m);  
    for(i=1;i<=n;++i) scanf("%d",&al[i]);                                  
    Initialize(),calc_sum();                           
    Ai.fix(n),Bi.fix(n);    
    for(i=0;i<n;++i)  
    {
        Ai.a[i]=(ll)qpow(i+1,m)*inv[i]%mod;    
        Bi.a[i]=(ll)qpow(i+1,2*m)*inv[i]%mod;   
    }
    Bi=Bi*Ai.Inv(n);      
    Ai=Ai.ln(n);         
    for(i=0;i<n;++i)  
    {

        Ai.a[i]=(ll)Ai.a[i]*_A.a[i]%mod;       
        Bi.a[i]=(ll)Bi.a[i]*_A.a[i]%mod;      
    }
    Bi=Bi*Ai.exp(n);               
    int ans=Bi.a[n-2];     
    for(i=1;i<=n;++i) ans=(ll)ans*al[i]%mod;    
    ans=(ll)ans*fac[n-2]%mod;    
    printf("%d
",ans);     
    return 0;   
}