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

10 篇文章 0 订阅
9 篇文章 0 订阅

【题意】
给出一串数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;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值