[caioj1456][FFT][拆系数板子]累加

来源:互联网 发布:java弹出消息框 编辑:程序博客网 时间:2024/06/10 22:31

【题意】
给出一串数a[i],定义数组b经过累加变换a得到。
如a{1,2,3,4}
则第一个b是{1,3,6,10}
此后的b则为累加变换b得到
如b{1,3,6,10}
则第二个b{1,4,10,20}
我们想知道第k次变换后的b数组
【输入】
第1行,2个数N和K,中间用空格分隔,N表示数组的长度,K表示处理的次数(2 <= n <= 50000, 0 <= k <= 10^9, 0 <= a[i] <= 10^9)
【输出】
共N行,每行一个数,对应经过K次处理后的结果。输出结果 Mod 10^9 + 7。
【样例输入】
4 1
1
2
3
4
【样例输出】
1
3
6
10
【题解】
引入了新的内容,拆系数FFT(膜myy)
本题最关键的一点就是在于模数,因为模数过大,FFT可能会炸精度。所以需要拆掉系数
设模数为M,则每个数都可以表示为K*√M+B
设两个多项式分别为A(x)和B(x),其系数序列分别为ai和bi
则可以拆出来4个序列k[ai],k[bi],b[ai],b[bi]
分别对4个序列进行FFT操作然后合并起来就达到了要求啦
注意拆系数FFT的精度要求较高,需要long double!!!(坑点)
之后我们来看题
这里写图片描述
上图是第几次操作后ai对后面数的贡献
可以发现系数是一个斜着的杨辉三角,那么可以用组合数快速求出来
这里写图片描述
除法使用逆元来求
那么求出c数组,将给出的数列和c数组拆系数后做卷积运算即可

#include<cstdio>#include<cstring>#include<cstdlib>#include<algorithm>#include<cmath>using namespace std;#define double long double//***卡了我一天的地方,这精度得用long double!!!! typedef long long LL;const int N=50010;const LL mod=(1e9+7);const LL m=sqrt(mod)+5;const double PI=acos(-1.0);const int MAXN=(1<<17);struct Complex{    double r,i;//real imag    Complex() {}    Complex(double _r,double _i){r=_r;i=_i;}    friend Complex operator + (const Complex &x,const Complex &y){return Complex(x.r+y.r,x.i+y.i);}    friend Complex operator - (const Complex &x,const Complex &y){return Complex(x.r-y.r,x.i-y.i);}    friend Complex operator * (const Complex &x,const Complex &y){return Complex(x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r);}    friend Complex operator * (const Complex &x,const double t){return Complex(x.r*t,x.i*t);}};int R[MAXN],L;Complex x1[MAXN],x2[MAXN],x3[MAXN],x4[MAXN];Complex x5[MAXN],x6[MAXN],x7[MAXN],x8[MAXN];void fft(Complex *y,int len,int on){    for(int i=0;i<len;i++)if(i<R[i])swap(y[i],y[R[i]]);    for(int i=1;i<len;i<<=1)    {        Complex wn(cos(PI/i),sin(on*PI/i));        for(int j=0;j<len;j+=(i<<1))        {            Complex w(1,0);            for(int k=0;k<i;k++,w=w*wn)            {                Complex u=y[j+k];                Complex v=w*y[j+k+i];                y[j+k]=u+v;                y[j+k+i]=u-v;            }        }    }    if(on==-1)    {        double t=1.0/len;        for(int i=0;i<len;i++)y[i].r*=t;//保护精度    }}LL pow_mod(LL a,LL b)//a^b%mod{    LL ans=1;a%=mod;    while(b)    {        if(b%2==1)ans=ans*a%mod;        a=a*a%mod;b/=2;    }    return ans;}LL p[550000];int n,k,a[51000];int len;LL inv(int x){return pow_mod(x,mod-2);}void pre(){   p[0]=1;   for(int i=1;i<len;i++)p[i]=p[i-1]*((i+k-1)%mod)%mod*inv(i)%mod;   //p[i]=p[i-1]*(i+k-1)/i   //inv(i)是i的逆元 乘i的逆元实质上就是除i }int main(){    scanf("%d%d",&n,&k);    for(len=1;len<=n*2;len<<=1)L++;    for(int i=0;i<len;i++)R[i]=(R[i>>1]>>1)|(i&1)<<(L-1);    pre();//组合数    for(int i=1;i<=n;i++)scanf("%d",&a[i]);    for(int i=0;i<=n;i++)    {        x1[i].r=a[i]/m;        x2[i].r=a[i]%m;        //a[i]=x1[i].r*m+x2.r        x3[i].r=p[i]%mod/m;        x4[i].r=p[i]%mod%m;        //p[i]=x3[i].r*m+x4.r    }    //如上拆掉系数,分别做fft    fft(x1,len,1);fft(x2,len,1);fft(x3,len,1);fft(x4,len,1);    //卷积写为a*p=(x1*m+x2)*(x3*m+x4)=(x1*x3*m*m)+(x1*m*x4)+(x2*x3*m)+(x2*x4)     //先把他们乘起来,m最后乘进去就行    for(int i=0;i<len;i++)    {        x5[i]=x1[i]*x3[i];        x6[i]=x1[i]*x4[i];        x7[i]=x2[i]*x3[i];        x8[i]=x2[i]*x4[i];    }    fft(x5,len,-1);fft(x6,len,-1);fft(x7,len,-1);fft(x8,len,-1);//IFFT    for(int i=1;i<=n;i++)    {        LL ans=0;        ans+=(long long)(x5[i].r+0.5)*(m*m%mod)%mod;        ans+=((long long)(x6[i].r+0.5)+(long long)(x7[i].r+0.5))%mod*m%mod;        ans+=(long long)(x8[i].r+0.5);        //套卷积式即可,为了防止爆精度我们多模几个 mod        printf("%d\n",ans%mod);    }    return 0;}
原创粉丝点击