定义一个操作：
[ A_k(i)=sum_{j=0}^{i}A_{k-1}(j)
]
已知(A_0(x))，求(A_k(x) mod 10^9+7)。A的次数界为(nleq 50,000，kleq 10^9)。
解析：
(A_k(i))只会对(jgeq i)的(A_m(i))产生贡献，不妨假设(A_0(i)=[i=0])，来看一看这个多项式变换k次后的(A_k(i))是多少。
[egin{align}
k=0 & & A= &{1,0,0,0,0,0cdots}
onumber\
k=1 & & A= &{1,1,1,1,1,1cdots}
onumber\
k=2 & & A= &{1,2,3,4,5,6cdots}
onumber\
k=3 & & A= &{1,3,6,10,15,21cdots}
onumber\
end{align}
]
可以发现，(A_0(0))对(A_k(j))的贡献倍率为(C_{k+j-1}^{k-1})。
同理，(A_0(i))对(A_k(j))的贡献倍率就是(C_{k+(j-i)-1}^{k-1})。
因此：
[ A_k(j)=sum_{i=0}^{j}{k+(j-i)-1choose k-1}A_0(i)\
令B(i)={i+k-1choose k-1}，则\
A_k(j)=sum_{i=0}^{j}B(j-i)A_0(i)\
]
只要构建出B数组，就可以卷积求出(A_k(j))
在(mod 10^9+7)意义下递推B数组。
[ B(0)=1\
B(i)=frac{B(i-1)*(i+k-1)}{i}
]
复杂度(O(nlog n))。若是线性预处理1~n的逆元的话复杂度为(O(n))。
[ps:线性求[1,n]在素数p模意义下的逆元:\inv(1)=1，inv(i)=(p-p/i)cdot inv(p mod i) mod p
]
由于模数较大，直接FFT每一项的值数量级为(10^9*10^9*n=10^{23})，double 和long double 的精度都不足以满足。
一般有两种姿势。
第一种是多模数NTT，一般最后合并时要用到高精度中国剩余定理或者一些特殊技巧，而且NTT次数较多，效率较低。
另一种是拆系数FFT，即将两个多项式A，B拆成四个。
[ 定义A_p(i)和A_q(i): A(i)=A_p(i)+A_q(i)cdot Base\
即: A_q(i)=leftlfloorfrac{A(i)}{Base}
ight
floor，A_p(i)=A(i)mod Base\
B同理\
Base 一般取sqrt{mod}最优
]
这样的话
[egin{align}
C(i)= & sum_{j=0}^{i}A(j)B(i-j)
onumber\
= & sum_{j=0}^{i}(A_p(j)+Basecdot A_q(j))(B_p(i-j)+Basecdot B_q(i-j))
onumber\
= & sum_{j=0}^{i}A_p(j)B_p(i-j)+Basesum_{j=0}^{i}A_p(j)B_q(i-j)
onumber\
& +Basesum_{j=0}^{i}A_q(j)B_p(i-j)+Base^2sum_{j=0}^{i} A_q(j)B_q(i-j)
onumber\
end{align}
]
所以变成了4遍卷积。每一项的阈值数量级为(10^4)，卷积结果的阈值数量级为(10^{13})double精度是可以承受的。
做完卷积之后乘以对应的(Base)再加起来就是(C(i))了。
复杂度(O(nlog n))，常数略大。
#include<stdio.h>
#include<math.h>
#include<iostream>
using namespace std;
typedef long long LL;
typedef double LDB;
const int mod=1000000007;
const LDB PI=2*acos(-1.0);
inline int Mod1(int x){return x>=mod?x-mod:x;}
const int Block=31623;
int n,K,N,L;
int A[50010],B[50010],Inv[50010];
struct Complex{
LDB a,b;
Complex(LDB x=0,LDB y=0){a=x,b=y;}
Complex operator + (const Complex& x)const{return Complex(a+x.a,b+x.b);}
Complex operator * (const Complex& x)const{return Complex(a*x.a-b*x.b,a*x.b+b*x.a);}
Complex operator - (const Complex& x)const{return Complex(a-x.a,b-x.b);}
void operator *= (const LDB& x){a*=x,b*=x;}
void print(){
printf("(%.3f,%.3f) ",(double)a,(double)b);
}
}s1[131082],s2[131082],t1[131082],t2[131082],wn[2][131082];
int Rev[131082];
void FFT(Complex* f,int Ty){
int i,j,l,c,t,k;
for(i=0;i<N;++i) if(Rev[i]<i) swap(f[Rev[i]],f[i]);
Complex w,x;
for(l=2,c=1,t=(N>>1);l<=N;c=l,l<<=1,t>>=1){
for(i=0;i<N;i+=l){
for(j=0,k=0;j<c;++j,k+=t){
x=f[i+j+c]*wn[Ty][k];
f[i+j+c]=f[i+j]-x;
f[i+j]=f[i+j]+x;
}
}
}
if(Ty){
LDB inv=(LDB)1/N;
for(int i=0;i<N;++i) f[i]*=inv;
}
}
void Work(){
for(int i=0;i<n;++i){
s1[i].a=A[i]/Block;
s2[i].a=A[i]%Block;
t1[i].a=B[i]/Block;
t2[i].a=B[i]%Block;
}
for(N=1,L=0;N<n;N<<=1,++L); N<<=1;
for(int i=0;i<N;++i){
Rev[i]=((Rev[i>>1]>>1)|((i&1)<<L));
wn[0][i]=Complex(cos(PI*i/N),sin(PI*i/N));
wn[1][i]=Complex(cos(-PI*i/N),sin(-PI*i/N));
}
FFT(s1,0); FFT(s2,0);
FFT(t1,0); FFT(t2,0);
Complex a,b,c,d;
for(int i=0;i<N;++i){
a=s1[i],b=s2[i],c=t1[i],d=t2[i];
s1[i]=a*c;
s2[i]=b*c;
t1[i]=a*d;
t2[i]=b*d;
}
FFT(s1,1); FFT(s2,1);
FFT(t1,1); FFT(t2,1);
int q,w,e,r,ans;
for(int i=0;i<n;++i){
ans=0;
q=(LL)(s1[i].a+0.5)%mod;
w=(LL)(s2[i].a+0.5)%mod;
e=(LL)(t1[i].a+0.5)%mod;
r=(LL)(t2[i].a+0.5)%mod;
ans=Mod1(ans+r);
ans=Mod1(ans+(LL)Block*Mod1(w+e)%mod);
ans=Mod1(ans+(LL)Block*Block*q%mod);
printf("%d
",ans);
}
}
int main(){
scanf("%d%d",&n,&K);
for(int i=0;i<n;++i) scanf("%d",&A[i]);
if(K==0){
for(int i=0;i<n;++i) printf("%d
",A[i]);
return 0;
}
--K;
Inv[1]=1; for(int i=2;i<n;++i) Inv[i]=(LL)(mod-mod/i)*Inv[mod%i]%mod;
B[0]=1; for(int i=1;i<n;++i) B[i]=(LL)B[i-1]*Inv[i]%mod*(K+i)%mod;
Work();
getchar(); getchar();
return 0;
}