Codechef SERSUM Series Sum

Link
运用Newton二项式定理可以得到(f(x,k)=sumlimits_{i=0}^k{kchoose i}x^isumlimits_{j=0}^na_j^{j-i})。
因此我们有(g(t,k)=sumlimits_{i=0}^tsumlimits_{j=0}^k{kchoose j}i^jsumlimits_{l=1}^na_l^{k-j}).
考虑化成卷积形式(frac{g(t,k)}{k!}=sumlimits_{i+j=k}frac{sumlimits_{x=0}^tx^i}{i!}frac{sumlimits_{x=1}^na_x^j}{j!})。
设(G(x)=sumlimits_{i=0}^kfrac{g(t,i)}{i!}x^i,P(x)=sumlimits_{i=0}^kfrac{sumlimits_{j=0}^tj^i}{i!}x^i,Q(x)=sumlimits_{i=0}^kfrac{sumlimits_{j=1}^na_j^i}{i!}x^i)，那么我们有(G(x)=P(x)Q(x))。
(Q(x))怎么算之前已经讲过了。Link
(P(x))实际上也可以用相同的办法算，不过这里我们考虑用Bernoulli数算。
我们知道([x^k]P(x)=frac{sumlimits_{i=0}^ti^k}{k!}=frac1{(k+1)!}sumlimits_{i=0}^k{k+1choose i}B_i(t+1)^{k+1-i}=sumlimits_{i+j=k}frac{B_i}{i!}frac{(t+1)^{j+1}}{(k-j+1)!})。
那么我们现在要考虑的就是怎么算(B(x)=sumlimits_{i=0}^kfrac{B_i}{i!}x^i)。
我们知道(sumlimits_{n=0}^{+infty}frac{B_n}{n!}x^n=frac x{e^x-1}=frac{x}{sumlimits_{n=0}^{+infty}frac{x^n}{n!}-1}=frac1{sumlimits_{n=0}^{+infty}frac{x^n}{(n+1)!}})，然后多项式求逆就好了。

#include<cmath>
#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
using ld=double;
const int N=131073,P=1000000007;const ld pi=2*acos(-1);
namespace IO
{
    char ibuf[(1<<21)+1],obuf[(1<<21)+1],st[15],*iS,*iT,*oS=obuf,*oT=obuf+(1<<21);
    char Get(){return (iS==iT? (iT=(iS=ibuf)+fread(ibuf,1,(1<<21)+1,stdin),(iS==iT? EOF:*iS++)):*iS++);}
    void Flush(){fwrite(obuf,1,oS-obuf,stdout),oS=obuf;}
    void Put(char x){*oS++=x;if(oS==oT)Flush();}
    int read(){int x=0,c=Get();while(!isdigit(c))c=Get();while(isdigit(c))x=x*10+c-48,c=Get();return x;}
    int raed(){int x=0,c=Get();while(!isdigit(c))c=Get();while(isdigit(c))x=(x*10ll+c-48)%P,c=Get();return x;}
    void write(int x){int top=0;if(!x)Put('0');while(x)st[++top]=(x%10)+48,x/=10;while(top)Put(st[top--]);Put(' ');}
}
using IO::read;
struct complex{ld x,y;complex(ld a=0,ld b=0){x=a,y=b;}}w[N];
complex operator+(complex a,complex b){return {a.x+b.x,a.y+b.y};}
complex operator-(complex a,complex b){return {a.x-b.x,a.y-b.y};}
complex operator*(complex a,ld x){return {a.x*x,a.y*x};}
complex operator/(complex a,ld x){return {a.x/x,a.y/x};}
complex operator*(complex a,complex b){return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
complex conj(complex a){return {a.x,-a.y};}
int inc(int a,int b){return a+=b-P,a+(a>>31&P);}
int dec(int a,int b){return a-=b,a+(a>>31&P);}
int mul(int a,int b){return 1ll*a*b%P;}
int pow(int a,int k){int r=1;for(;k;k>>=1,a=mul(a,a))if(k&1)r=mul(a,r);return r;}
int n,k,m,len,pw[N],fac[N],inv[N],ifac[N],rev[N];
int getlen(int n){return 1<<(32-__builtin_clz(n));}
void init(int n)
{
    int lim=1<<(len=32-__builtin_clz(n)),p=lim>>1;ld x;
    w[lim>>1]={1,0},fac[0]=ifac[0]=inv[0]=fac[1]=ifac[1]=inv[1]=pw[0]=1,pw[1]=m;
    for(int i=1;i<lim;++i) rev[i]=(rev[i>>1]>>1)|(i&1? p:0);
    for(int i=0;i<p;++i) x=pi*i/lim,w[p+i]={cos(x),sin(x)};
    for(int i=p-1;i;--i) w[i]=w[i<<1];
    for(int i=2;i<=lim;++i) fac[i]=mul(fac[i-1],i),ifac[i]=mul(ifac[i-1],inv[i]=mul(inv[P%i],P-P/i)),pw[i]=mul(pw[i-1],m);
}
void FFT(complex*a,int lim,int f)
{
    static complex x;
    if(!~f) std::reverse(a+1,a+lim);
    for(int i=0,x=len-__builtin_ctz(lim);i<lim;++i) if(i<rev[i]>>x) std::swap(a[i],a[rev[i]>>x]);
    for(int i=1;i<lim;i<<=1) for(int j=0,d=i<<1;j<lim;j+=d) for(int k=0;k<i;++k) x=a[i+j+k]*w[i+k],a[i+j+k]=a[j+k]-x,a[j+k]=a[j+k]+x;
    if(!~f) for(int i=0;i<lim;++i) a[i]=a[i]/lim;
}
void DFT(complex*a,complex*b,int lim)
{
    static complex t[2][N],x,y;
    for(int i=0;i<lim;++i) t[0][i]=a[i]+b[i]*complex{0,1};
    FFT(t[0],lim,1),memcpy(t[1],t[0],lim<<4),std::reverse(t[1]+1,t[1]+lim);
    for(int i=0;i<lim;++i) x=t[0][i],y=conj(t[1][i]),a[i]=(x+y)/2,b[i]=(x-y)*complex{0,-1}/2;
}
void IDFT(complex*a,complex*b,int lim)
{
    for(int i=0;i<lim;++i) a[i]=a[i]+b[i]*complex{0,1};
    FFT(a,lim,-1);
    for(int i=0;i<lim;++i) b[i]={a[i].y,0},a[i].y=0;
}
void MTT(int*a,int*b,int*c,int n,int m)
{
    static complex t[4][N],A,B,C,D;static int r[4]={1<<30,1<<15,1<<15,1},s[N];
    int lim=getlen(n+m-2);
    if(n<=50||m<=50)
    {
	memset(s,0,(n+m)<<2);
	for(int i=0;i<n;++i) for(int j=0;j<m;++j) s[i+j]=inc(s[i+j],mul(a[i],b[j]));
	memcpy(c,s,(n+m)<<2);
	return ;
    }
    for(int i=0;i<4;++i) memset(t[i],0,lim<<4);
    for(int i=0;i<n;++i) t[0][i].x=a[i]>>15,t[1][i].x=a[i]&32767;
    for(int i=0;i<m;++i) t[2][i].x=b[i]>>15,t[3][i].x=b[i]&32767;
    DFT(t[0],t[1],lim),DFT(t[2],t[3],lim);
    for(int i=0;i<lim;++i) A=t[0][i]*t[2][i],B=t[0][i]*t[3][i],C=t[1][i]*t[2][i],D=t[1][i]*t[3][i],t[0][i]=A,t[1][i]=B,t[2][i]=C,t[3][i]=D;
    IDFT(t[0],t[1],lim),IDFT(t[2],t[3],lim),memset(c,0,(n+m)<<2);
    for(int i=0;i<4;++i) for(int j=0;j<n+m;++j) c[j]=inc(c[j],mul(r[i],(long long)(t[i][j].x+0.5)%P));
}
void Inv(int*a,int*b,int deg)
{
    if(deg==1) return b[0]=pow(a[0],P-2),void();
    static int t[N];
    Inv(a,b,(deg+1)>>1),MTT(a,b,t,deg,deg),MTT(b,t,t,deg,deg);
    for(int i=0;i<deg;++i) b[i]=dec(inc(b[i],b[i]),t[i]);
}
void Der(int*a,int*b,int deg){for(int i=1;i<deg;++i)b[i-1]=mul(a[i],i);b[deg-1]=0;}
void Int(int*a,int*b,int deg){for(int i=1;i<deg;++i)b[i]=mul(a[i-1],inv[i]);b[0]=0;}
void Ln(int*a,int*b,int deg)
{
    static int t[N];
    Inv(a,t,deg),Der(a,b,deg),MTT(b,t,t,deg,deg),Int(t,b,deg);
}
int a[N],f[17][N],t[2][N];
void solve(int l,int r,int d)
{
    if(l==r) return f[d][0]=1,f[d][1]=dec(0,a[l]),void();
    int mid=(l+r)>>1;
    solve(l,mid,d),solve(mid+1,r,d+1),MTT(f[d],f[d+1],f[d],mid-l+2,r-mid+1);
}
int main()
{
    n=read(),k=read(),m=IO::raed()+1,init(std::max(k*2,n));
    for(int i=1;i<=n;++i) a[i]=read();
    for(int i=0;i<=k;++i) t[0][i]=ifac[i+1];
    Inv(t[0],t[1],k+1); 
    for(int i=0;i<=k;++i) t[0][i]=mul(pw[i+1],ifac[i+1]);
    MTT(t[0],t[1],t[1],k+1,k+1),solve(1,n,0),Ln(f[0],t[0],k+1),t[0][0]=n;  
    for(int i=1;i<=k;++i) t[0][i]=mul(ifac[i],dec(0,mul(i,t[0][i])));
    MTT(t[1],t[0],t[0],k+1,k+1);
    for(int i=0;i<=k;++i) IO::write(mul(fac[i],t[0][i]));
    IO::Flush();
}