51nod 1172 Partial Sums V2 任意模FFT

Partial Sums V2

题目连接：

https://www.51nod.com/onlineJudge/questionCode.html#!problemId=1172

Description

给出一个数组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)

 Input

第1行，2个数N和K，中间用空格分隔，N表示数组的长度，K表示处理的次数（2 <= n <= 50000, 0 <= k <= 10^9, 0 <= a[i] <= 10^9)

Output

共N行，每行一个数，对应经过K次处理后的结果。输出结果 Mod 10^9 + 7。

Sample Input

4 2
1
3
5
6

Sample Output

1
5
14
29

题意：

题解：

考虑每个位置的贡献，所求即$ans[n]=\sum_{i+j=n}^{}C(k-1+i,i)\cdot a[j]$ ，其中 $C(n,i) $表示从 n 个物品里选出 i 个物品的方案数。
上述卷积形式可以利用FFT或NTT加速运算，难点在于答案模 109+7，设其为 P ，这意味这参与卷积的每个元素值为不大于 (P−1) 的非负整数，卷积后的值为不大于 (n+1)⋅(P−1)2=O(nP2) 的值，这个上界可以达到 5⋅1022 。
一种想法是利用3种模意义下的NTT结果进行合并，得到这个可能达到 5⋅1022 的数，再取模，一共需要9次NTT，涉及高精度。
还有一种想法就是降低卷积前的元素值上界，令 M 为满足 P≤M2 的最小正整数，将多项式 f 拆成 f1+f2⋅M 的形式，那么每个元素的值上界是 O(M)=O($\sqrt{p}$) 的，卷积后的值上界是 O(nP) 的，用long double来做FFT差不多可以接受，而卷积可以写为 f∗g=(f1+f2⋅M)∗(g1+g2⋅M)=(f1∗f2)+(f1∗g2+f2∗g1)⋅M+(f2∗g2)⋅M2 ，一共需要7次FFT，不涉及高精度。

 

很邪，51nod上c++会wa，c++11会过..

以前只在大牛的博客上看到过第一种做法，但没有实用性，第二种方法很强。

昨天刚好看到，拉出来封个版


----------


2016-8-18 更新
看到另一种黑科技优化，大概是fft中用了相位的部分，只用4次fft，膜

代码   

//#include <bits/stdc++.h>
#include <stdio.h>
#include <iostream>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include <limits.h>
#include <algorithm>
#include <queue>
#include <vector>
#include <set>
#include <map>
#include <stack>
#include <bitset>
#include <string>
#include <time.h>
using namespace std;
long double esp=1e-11;
//#pragma comment(linker, "/STACK:1024000000,1024000000")
#define fi first
#define se second
#define all(a) (a).begin(),(a).end()
#define cle(a) while(!a.empty())a.pop()
#define mem(p,c) memset(p,c,sizeof(p))
#define mp(A, B) make_pair(A, B)
#define pb push_back
#define lson l , m , rt << 1
#define rson m + 1 , r , rt << 1 | 1
typedef long long int LL;
const long double PI = acos((long double)-1);
const LL INF=0x3f3f3f3fll;
const int MOD =1000000007ll;
const int maxn=140100;
struct complex
{
    long double r,i;
    complex(long double _r = 0.0,long double _i = 0.0){r = _r; i = _i;}
    complex operator +(const complex &b){return complex(r+b.r,i+b.i);}
    complex operator -(const complex &b){return complex(r-b.r,i-b.i);}
    complex operator *(const complex &b){return complex(r*b.r-i*b.i,r*b.i+i*b.r);}
};
complex conj(complex a){return complex(a.r,-a.i);}
void change(complex y[],int len)
{
    int i,j,k;
    for(i = 1, j = len/2;i < len-1; i++)
    {
        if(i < j)swap(y[i],y[j]);
        k = len/2;
        while( j >= k)
        {
            j -= k;
            k /= 2;
        }
        if(j < k) j += k;
    }
}
void FFT(complex y[],int len,int on)  //len=2^k
{
    change(y,len);
    for(int h = 2; h <= len; h <<= 1)
    {
        complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
        for(int j = 0;j < len;j+=h)
        {
            complex w(1,0);
            for(int k = j;k < j+h/2;k++)
            {
                complex u = y[k];
                complex t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1)
        for(int i = 0;i < len;i++)
            y[i].r /= len;
}
int callen(int len1,int len2){
    int len=1;
    while(len < (len1<<1) || len < (len2<<1))len<<=1;
    return len;
}
LL fftans[maxn];
complex A[maxn],B[maxn],dft[4][maxn],dt[4];
int td[4];
void fft(LL* y1,int len1,LL* y2,int len2,LL mod){
    int len=callen(len1,len2);
    for(int x=0;x<len1;x++)A[x]=complex(y1[x]&32767,y1[x]>>15);
    for(int x=len1;x<len;x++)A[x]=complex(0,0);
    for(int x=0;x<len2;x++)B[x]=complex(y2[x]&32767,y2[x]>>15);
    for(int x=len2;x<len;x++)B[x]=complex(0,0);
    FFT(A,len,1);FFT(B,len,1);
    int j;
    for(int x=0;x<len;x++)
    {
        j=(len-x)&(len-1);
        dt[0]=(A[x]+conj(A[j]))*complex(0.5,0);
        dt[1]=(A[x]-conj(A[j]))*complex(0,-0.5);
        dt[2]=(B[x]+conj(B[j]))*complex(0.5,0);
        dt[3]=(B[x]-conj(B[j]))*complex(0,-0.5);
        dft[0][j]=dt[0]*dt[2];
        dft[1][j]=dt[0]*dt[3];
        dft[2][j]=dt[1]*dt[2];
        dft[3][j]=dt[1]*dt[3];
    }
    for(int x=0;x<len;x++)
    {
        A[x]=dft[0][x]+dft[1][x]*complex(0,1);
        B[x]=dft[2][x]+dft[3][x]*complex(0,1);
    }
    FFT(A,len,1);FFT(B,len,1);
    for(int x=0;x<len;x++)
    {
        td[0]=(LL)(A[x].r/len+0.5)%mod;
        td[1]=(LL)(A[x].i/len+0.5)%mod;
        td[2]=(LL)(B[x].r/len+0.5)%mod;
        td[3]=(LL)(B[x].i/len+0.5)%mod;
        fftans[x]=(td[0]+((LL)(td[1]+td[2])<<15)+((LL)td[3]<<30))%mod;
    }
}
LL inv[50001];
LL a[maxn],b[maxn];

int main()
{
    //freopen("in.txt", "r", stdin);
    //freopen("inlay.in", "r", stdin);
    //freopen("out.txt", "w", stdout);      //%I64d
    //vector<int>::iterator iter;
    //for(int x=1;x<=n;x++)
    //for(int y=1;y<=n;y++)
    //scanf("%d",&a);
    //printf("%d\n",ans);
    int n;
    LL k;
    scanf("%d%lld",&n,&k);
    inv[1]=1;
    for(int x=2;x<n;x++)inv[x]=(MOD-MOD/x)*inv[MOD%x]%MOD;
    for(int x=0;x<n;x++)
        scanf("%lld",&a[x]);
    LL tem=1;
    b[0]=1;
    for(int x=1;x<n;x++)
    {
        tem=tem*(x+k-1)%MOD*inv[x]%MOD;
        b[x]=tem;
    }
    fft(a,n,b,n,MOD);
    for(int x=0;x<n;x++)
        printf("%lld\n",(fftans[x]+MOD)%MOD);
    return 0;
}

2016-08-18