51nod 1172 Partial Sums V2 任意模FFT

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)

Solution

---

之前模拟赛见过类似的结论。。

我们考虑一个位置前缀和k次后对后面每一个位置的贡献，记s[i,j]表示往后第i位，第k次前缀和时这一位的贡献，可以发现s[i,j]=f[i+j-1,i]，这里f是杨辉三角形。感性地讲就是斜过来的杨辉三角形

然后就好办了，显然有$an{s}_n=\sum _{i=0}^n\left(\genfrac{}{}{0ex}{}{k-1-i}{k-1}\right)\cdot a_{n-i}$
这是一个卷积形式，我们可以用NTT来做。由于模数的问题，我们可以考虑CRT合并也可以拆系数

记$\sqrt{MOD}=M$，那么所有系数都可以表示成$k\cdot M+b$，我们分别从原本的两个多项式中抽出四个多项式，然后利用乘法分配律就可以了

非常诡异的事情是我的代码和beginend大爷对于同一组数据答案相同，但是交上去只有我wa了。。
拆系数FFT需要很高的精度，并且要预处理单位复根

Code

---

#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <math.h>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)
#define drp(i,st,ed) for (int i=st;i>=ed;--i)
#define fill(x,t) memset(x,t,sizeof(x))

typedef long long LL;
const long double pi=acos((long double)-1);
const int N=400005;

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

int rev[N],n,k;

LL ny[N];

void FFT(com *a,int n,int f) {
	for (int i=0;i<n;++i) if (i<rev[i]) std:: swap(a[i],a[rev[i]]);
	static com w[N];
	w[0]=(com) {1,0};
	for (int i=1;i<n;i<<=1) {
		com wn=(com) {cos(pi/i),f*sin(pi/i)};
		for (int j=i;j>-1;j-=2) w[j]=w[j/2];
		for (int j=1;j<i;j+=2) w[j]=w[j-1]*wn;
		for (int j=0;j<n;j+=(i<<1)) {
			for (int k=0;k<i;++k) {
				com u=a[j+k],v=a[j+k+i]*w[k];
				a[j+k]=u+v; a[j+k+i]=u-v;
			}
		}
	}
	if (f==-1) for (int i=0;i<n;++i) a[i]=a[i]/n;
}

void solve(LL *a,LL *b,int M) {
	static com a1[N],a2[N],b1[N],b2[N],w1[N],w2[N],w3[N];
	for (int i=0;i<n;++i) {
		a1[i]=(com) {(int)a[i]/M,0}; b1[i]=(com) {a[i]%M,0};
		a2[i]=(com) {(int)b[i]/M,0}; b2[i]=(com) {b[i]%M,0};
	}
	int len=1,lg=0; for (;len<=n*2;) lg++,len<<=1;
	for (int i=0;i<len;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
	FFT(a1,len,1); FFT(b1,len,1);
	FFT(a2,len,1); FFT(b2,len,1);
	for (int i=0;i<len;++i) {
		w1[i]=a1[i]*a2[i];
		w2[i]=b1[i]*a2[i]+b2[i]*a1[i];
		w3[i]=b1[i]*b2[i];
	}
	FFT(w1,len,-1);
	FFT(w2,len,-1);
	FFT(w3,len,-1);
	const int MOD=1000000007;
	for (int i=0;i<n;++i) {
		LL ans=(LL)(w1[i].x+0.5)%MOD*M*M%MOD;
		ans=(ans+(LL)(w2[i].x+0.5)%MOD*M%MOD)%MOD;
		ans=(ans+(LL)(w3[i].x+0.5)%MOD)%MOD;
		ans=(ans%MOD+MOD)%MOD;
		printf("%lld\n", ans);
		// printf("%.2Lf %.2Lf\n", w1[i].x,w1[i].y);
		// printf("%.2Lf %.2Lf\n", w2[i].x,w2[i].y);
		// printf("%.2Lf %.2Lf\n", w3[i].x,w3[i].y);
	}
}

int main(void) {
	freopen("data.in","r",stdin);
	freopen("myp.out","w",stdout);
	static LL a[N],c[N];
	const int MOD=1000000007;
	ny[0]=ny[1]=1;
	rep(i,2,N-1) ny[i]=1LL*(MOD-MOD/i)*ny[MOD%i]%MOD;
	rep(i,2,N-1) ny[i]=ny[i-1]*ny[i]%MOD;
	scanf("%d%d",&n,&k);
	c[0]=1; c[1]=k;
	rep(i,2,N-1) c[i]=c[i-1]*(k-1+i)%MOD;
	rep(i,2,N-1) c[i]=c[i]*ny[i]%MOD;
	rep(i,0,n-1) scanf("%lld",&a[i]);
	int M=sqrt(MOD);
	solve(a,c,M);
	return 0;
}