51nod 1172 Partial Sums V2 任意模数FFT

题意

给出一个数组A，经过一次处理，生成一个数组S，数组S中的每个值相当于数组A的累加，比如：A = {1 3 5 6} => S = {1 4 9 15}。如果对生成的数组S再进行一次累加操作，{1 4 9 15} => {1 5 14 29}，现在给出数组A，问进行K次操作后的结果。(输出结果 Mod 10^9 + 7)
2 <= n <= 50000, 0 <= k <= 10^9, 0 <= a[i] <= 10^9

分析

首先有个结论就是$ans[n]=\sum_{i=0}^nC_{i+k-1}^i*a[n-i]$。
这条式子是怎么来的呢？我的理解就是，把每一次前缀和后的数组放在一起，初始数组在第0行，前缀和放在第一行，如此类推。那么a[i]对a[n]贡献的系数就相当于每次可以从(x,y)走到(x+1,y+l)（l>=0），从(0,i)走到(k,n)的不同方案数。根据隔板法不难得到系数就是$C_{n-i+k-1}^{n-i}$。
上面的式子是一个卷积的形式，可以用FFT来优化。但注意到这里的模数并不是NTT模数，那么就要用到任意模数FFT。
具体来讲就是取$M=\sqrt P$，设$f(x)=k(x)*M+r(x)$，那么有
$f_1(x)*f_2(x)=(k_1(x)*M+r_1(x))*(k_2(x)*M+r_2(x))$
$=M^2*k_1(x)*k_2(x)+M*(k_1(x)*r_2(x)+k_2(x)*r_1(x))+r_1(x)*r_2(x)$
总共需要做7次FFT且卷积后元素的数量级为O(nP)，可以用long double来进行FFT而不是高精度。
一开始打完后被卡了精度，后面把pi从define变成const就过了。。。
第一次手打复数类型，发现比自带的要快将近一倍。。。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<cmath>

using namespace std;

typedef long long LL;

const int MOD=1000000007;
const int N=50005;
const long double pi=acos((long double)-1);

struct com
{
    long double x,y;
    com operator + (const com &d) {return (com){x+d.x,y+d.y};}
    com operator - (const com &d) {return (com){x-d.x,y-d.y};}
    com operator * (const com &d) {return (com){x*d.x-y*d.y,x*d.y+d.x*y};}
    com operator / (const long double &d) {return (com){x/d,y/d};}
};

int n,M,ny[N],c[N],a[N],k,L,rev[N*4];
com s1[N*4],s2[N*4],s3[N*4],k1[N*4],r1[N*4],k2[N*4],r2[N*4];

int read()
{
    int x=0,f=1;char ch=getchar();
    while (ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}

void init()
{
    n=read();k=read();
    for (int i=0;i<n;i++) a[i]=read();
    ny[0]=ny[1]=c[0]=1;c[1]=k;
    for (int i=2;i<n;i++) ny[i]=(LL)(MOD-MOD/i)*ny[MOD%i]%MOD,c[i]=(LL)c[i-1]*(k+i-1)%MOD;
    for (int i=2;i<n;i++) ny[i]=(LL)ny[i]*ny[i-1]%MOD,c[i]=(LL)c[i]*ny[i]%MOD;
}

void fft(com *a,int f)
{
    for (int i=0;i<L;i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
    for (int i=1;i<L;i<<=1)
    {
        com wn=(com){cos(pi/i),f*sin(pi/i)};
        for (int j=0;j<L;j+=(i<<1))
        {
            com w=(com){1,0};
            for (int k=0;k<i;k++)
            {
                com u=a[j+k],v=a[j+k+i]*w;
                a[j+k]=u+v;a[j+k+i]=u-v;
                w=w*wn;
            }
        }
    }
    if (f==-1) for (int i=0;i<L;i++) a[i]=a[i]/L;
}

int main()
{
    init();M=sqrt(MOD);
    int lg=0;
    for (L=1;L<=n*2;L<<=1,lg++);
    for (int i=0;i<L;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
    for (int i=0;i<n;i++) k1[i]=(com){a[i]/M,0.0},r1[i]=(com){a[i]%M,0.0},k2[i]=(com){c[i]/M,0.0},r2[i]=(com){c[i]%M,0.0};
    fft(k1,1);fft(r1,1);fft(k2,1);fft(r2,1);
    for (int i=0;i<L;i++) s1[i]=k1[i]*k2[i],s2[i]=k1[i]*r2[i]+k2[i]*r1[i],s3[i]=r1[i]*r2[i];
    fft(s1,-1);fft(s2,-1);fft(s3,-1);
    for (int i=0;i<n;i++)
    {
        int x1=(LL)(s1[i].x+0.5)%MOD*M*M%MOD,x2=(LL)(s2[i].x+0.5)%MOD*M%MOD,x3=(LL)(s3[i].x+0.5)%MOD;
        int ans=x1+x2;ans-=ans>=MOD?MOD:0;ans+=x3;ans-=ans>=MOD?MOD:0;
        printf("%d\n",ans);
    }
    return 0;
}