@codechef - SERSUM@ Series Sum

[toc]

---

@description@

給定 N 個整數 a1,a2,...,aN，定義函數 f 和 g： f(x,k) = (x + a1)^k + (x + a2)^k +···+ (x + aN)^k g(t,k) = f(0,k) + f(1,k) +···+ f(t,k) 給定整數 T 和 K，對於每一個 0 ∼ K 之間的 i，請計算 g(T,i) 對 10^9 + 7 取模的結果。

輸入格式 輸入的第一行包含三個整數 N,K,T。第二行包含 N 個整數 a1,a2,...,aN。

輸出格式 輸出一行，包含 K + 1 個整數，表明 g(T,0),g(T,1),...,g(T,K) 對 10^9 + 7 取模的結果。

數據範圍與子任務 • 1 ≤ N ≤ 10^5 • 1 ≤ K ≤ 5·10^4 • 1 ≤ T ≤ 10^18 • 0 ≤ ai < 10^9 + 7

樣例輸入 2 3 4 0 1 樣例輸出 10 25 85 325

@solution@

@part - 1@

先寫式子： $$g(T, k) = \sum_{i=0}^{T}f(i, k) = \sum_{i=0}^{T}\sum_{j=1}^{N}(i + a_j)^k = \sum_{i=0}^{T}\sum_{j=1}^{N}\sum_{p=0}^{k}C_k^pi^pa_j^{k-p}$$

經過將組合數拆成階乘形式，並做適當地變形，能夠獲得： $$g(T, k) = k!\sum_{p=0}^{k}((\frac{\sum_{i=0}^{T}i^p}{p!})(\frac{\sum_{j=1}^{N}a_j^{k-p}}{(k-p)!}))$$

令 $P(x) = \frac{\sum_{i=0}^{T}i^x}{x!}, Q(x) = \frac{\sum_{i=1}^{N}a_i^x}{x!}$，則有 $g(T, k) = k!*\sum_{p=0}^{k}P(p)*Q(k-p)$。 顯然是一個卷積，可使用 fft（可是這裏由於模數的緣由，要使用神奇的 mtt）。

咱們接下來的問題轉爲求解 P 和 Q。

@part - 2@

先考慮 P，發現它是一個很裸的天然數冪和，可使用伯努利數求解。 因此如下介紹使用伯努利數求解天然數冪和的過程，熟悉的人能夠直接跳過。 雖然之前寫過一點，不過感受不大詳細。

記 $S(n, k) = \sum_{i=0}^{n}i^k$。先考慮一種使用裂項相消求解天然數冪和的方法： 首先由二項式定理獲得 $(n+1)^{k+1} - n^{k+1} = \sum_{i=0}^{k}C_{k+1}^{i}*n^i$。 因而： $$\sum_{i=0}^{n}((i+1)^{k+1} - i^{k+1}) = \sum_{i=0}^{n}\sum_{j=0}^{k}C_{k+1}^{j}*i^j$$ $$(n+1)^{k+1} = \sum_{j=0}^{k}C_{k+1}^{j}S(n, j)$$ $$\frac{(n+1)^{k+1}}{(k+1)!} = \sum_{j=0}^{k}\frac{S(n, j)}{j!}\frac{1}{(k+1-j)!}$$

這個式子感受很是卷積，並且感受很是指數型生成函數，因此咱們就往這個方向思考。

令 $F_n(x) = \sum_{i=0}\frac{S(n, k)}{k!}x^k$，即 S(n, k) 的指數型生成函數。 則有：$e^{(n+1)x} - 1 = (e^x - 1)*F_n$，因而就有 $F_n = \frac{e^{(n+1)x} - 1}{e^x - 1}$。 由於 $e^x - 1$ 的常數項爲 0，使用牛頓迭代找不出逆元，因此咱們能夠找 $\frac{e^x - 1}{x}$ 的逆元（即伯努利數的生成函數）。

因而： $$F_n = \frac{e^{(n+1)x} - 1}{x}*\frac{x}{e^x - 1} = (\sum_{i=0}\frac{(n+1)^{i+1}}{(i+1)!}x^i)(\frac{1}{\sum_{i=0}\frac{1}{(i+1)!}*x^i})$$

寫一個多項式求逆便可。同樣，由於模數緣由要使用 mtt。

@part - 3@

考慮怎麼求解 Q。這應該算是本題最關鍵的部分。

咱們令 R(x) 表示從序列 a 中不重複地選出 x 個數相乘的結果之和。或者說，R 是多項式 $\prod_{i=1}^{N}(a_i + 1)$ 的係數。 可使用分治 fft 簡單地求出 R。

經過容斥能夠求出 Q(x) 的遞推式： $$Q(x) = \sum_{i=1}^{x-1}(-1)^{i-1}*Q(x-i)*R(i) + (-1)^{x-1}*R(x)*x$$

怎麼理解呢？ 首先： $$Q(x-1)*R(1) = \sum_{i=1}^{N}a_i^x + \sum_{i=1}^{N}\sum_{j=1,i\not =j}^{N}a_i^{x-1}*a_j$$

前一半是咱們須要的，後一半是咱們想要容斥掉的。 而後： $$Q(x-2)*R(2) = \sum_{i=1}^{N}\sum_{j=1,i\not =j}^{N}a_i^{x-1}*a_j + \sum_{i=1}^{N}\sum_{j=1,i\not =j}^{N}\sum_{k=1,i\not =k,j\not =k}^{N}a_i^{x-2}a_ja_k$$

同上面同樣，前一半是咱們想要的，後一半是咱們想要容斥掉的。 因而容斥到最後，全部 ai 的指數都爲 1 時（不難發現此時項數爲 x），能夠發現一組 (ai, aj, ak, ...) 的貢獻被重複計算了 x 次（在每一項都被計算了一次）。 因而在最後消掉這些貢獻便可，即 $(-1)^{x-1}*R(x)*x$ 一項的含義。

怎麼想出來的呢？我怎麼知道，這是人類的智慧。

而後就是套路通常的生成函數了。 令 $G(x) = \sum_{i=1}(-1)^{i-1}*R(i)*x^i, H(x) = \sum_{i=1}(-1)^{i-1}R(i)ix^i$。 則 $Q = QG + H$，因而 $Q = \frac{H}{1-G}$。

要考慮 Q 的邊界狀況，即 $a_1^0 + a_2^0 + ... + a_N^0$ 的狀況。

@accepted code@

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<iostream>
using namespace std;
typedef long long ll;
typedef long double ld;
struct complex{
	ld r, i;
	complex(ld _r=0, ld _i=0):r(_r), i(_i) {}
	friend complex operator + (complex a, complex b) {
		return complex(a.r + b.r, a.i + b.i);
	}
	friend complex operator - (complex a, complex b) {
		return complex(a.r - b.r, a.i - b.i);
	}
	friend complex operator * (complex a, complex b) {
		return complex(a.r*b.r - a.i*b.i, a.i*b.r + b.i*a.r);
	}
	friend complex operator / (complex a, ld k) {
		return complex(a.r / k, a.i / k);
	}
	friend complex conj(complex a) {
		return complex(a.r, -a.i);
	}
};
typedef complex comp;
const int MAXN = 200000;
const int MOD = int(1E9) + 7;
const int SQ = sqrt(MOD);
const ld PI = acos(-1);
int pow_mod(int b, int p) {
	int ret = 1;
	while( p ) {
		if( p & 1 ) ret = 1LL*ret*b%MOD;
		b = 1LL*b*b%MOD;
		p >>= 1;
	}
	return ret;
}
struct poly{
	comp w[20];
	poly() {
		for(int i=0;i<20;i++)
			w[i] = comp(cos(2*PI/(1<<i)), sin(2*PI/(1<<i)));
	}
	void clear(comp *A, int l, int r) {
		for(int i=l;i<r;i++)
			A[i] = comp(0, 0);
	}
	void copy(comp *A, comp *B, int n) {
		for(int i=0;i<n;i++)
			A[i] = B[i];
	}
	void debug(comp *A, int n) {
		for(int i=0;i<n;i++)
			cout << A[i].r << " " << A[i].i << endl;
		puts("");
	}
	void debug(int *A, int n) {
		for(int i=0;i<n;i++)
			cout << A[i] << endl;
		puts("");
	}
	int length(int n) {
		int len; for(len = 1; len < n; len <<= 1) ;
		return len;
	}
	void fft(comp *A, int n, int type) {
		for(int i=0,j=0;i<n;i++) {
			if( i < j ) swap(A[i], A[j]);
			for(int k=(n>>1);(j^=k)<k;k>>=1);
		}
		for(int i=1;(1<<i)<=n;i++) {
			int s = (1<<i), t = (s>>1);
			comp u = comp(w[i].r, type*w[i].i);
			for(int j=0;j<n;j+=s) {
				comp p = complex(1, 0);
				for(int k=0;k<t;k++,p=p*u) {
					comp x = A[j+k], y = p*A[j+k+t];
					A[j+k] = x + y, A[j+k+t] = x - y;
				}
			}
		}
		if( type == -1 ) {
			for(int i=0;i<n;i++)
				A[i] = A[i] / n;
		}
	}
	comp tmp2[MAXN + 5], tmp3[MAXN + 5];
	void poly_mul(int *A, int *B, int *C, int n, int m, int k) {
		int len = length(n + m - 1);
		clear(tmp2, 0, len), clear(tmp3, 0, len);
		for(int i=0;i<n;i++) tmp2[i] = comp(A[i]/SQ, A[i]%SQ);
		for(int i=0;i<m;i++) tmp3[i] = comp(B[i]/SQ, B[i]%SQ);
		fft(tmp2, len, 1), fft(tmp3, len, 1);
		for(int i=0;i<=len/2;i++) {
			comp x0 = tmp2[i], y0 = tmp2[(len-i)%len];
			comp x1 = tmp3[i], y1 = tmp3[(len-i)%len];
			comp a0 = (x0 + conj(y0))/2, a1 = (conj(y0) - x0)/2*comp(0, 1);
			comp b0 = (x1 + conj(y1))/2, b1 = (conj(y1) - x1)/2*comp(0, 1);
			tmp2[i] = a0*b0 + comp(0, 1)*a1*b1, tmp3[i] = a0*b1 + a1*b0;
			a0 = (y0 + conj(x0))/2, a1 = (conj(x0) - y0)/2*comp(0, 1);
			b0 = (y1 + conj(x1))/2, b1 = (conj(x1) - y1)/2*comp(0, 1);
			tmp2[(len-i)%len] = a0*b0 + comp(0, 1)*a1*b1, tmp3[(len-i)%len] = a0*b1 + a1*b0;
		}
		fft(tmp2, len, -1), fft(tmp3, len, -1);
		for(int i=0;i<len;i++)
			C[i] = (ll(tmp2[i].r+0.5)%MOD*SQ%MOD*SQ%MOD + ll(tmp3[i].r+0.5)%MOD*SQ%MOD + ll(tmp2[i].i+0.5)%MOD)%MOD;
	}
	int tmp1[MAXN + 5];
	void poly_inv(int *A, int *B, int n) {
		if( n == 1 ) {
			B[0] = pow_mod(A[0], MOD - 2);
			return ;
		}
		poly_inv(A, B, (n + 1) >> 1);
		poly_mul(A, B, tmp1, n, (n+1)>>1, n);
		for(int i=0;i<n;i++)
			tmp1[i] = (MOD - tmp1[i])%MOD;
		tmp1[0] = (tmp1[0] + 2)%MOD;
		poly_mul(tmp1, B, B, n, (n+1)>>1, n);
	}
}oper;
int fct[MAXN + 5], ifct[MAXN + 5];
void init() {
	fct[0] = 1;
	for(int i=1;i<=MAXN;i++)
		fct[i] = 1LL*fct[i-1]*i%MOD;
	ifct[MAXN] = pow_mod(fct[MAXN], MOD-2);
	for(int i=MAXN-1;i>=0;i--)
		ifct[i] = 1LL*ifct[i+1]*(i+1)%MOD;
}
int a[MAXN + 5], N, K; ll T;
int f[MAXN + 5], B[MAXN + 5], tmp[MAXN + 5];
void func1() {
	for(int i=0;i<K;i++)
		tmp[i] = ifct[i+1];
	oper.poly_inv(tmp, B, K);
	for(int i=0,p=T;i<K;i++,p=1LL*p*T%MOD)
		tmp[i] = 1LL*p*ifct[i+1]%MOD;
	oper.poly_mul(tmp, B, f, K, K, K);
}
int g[MAXN + 5], h[MAXN + 5], hl[20][MAXN + 5], hr[20][MAXN + 5];
void func3(int le, int ri, int dep, int *A) {
	if( le + 1 == ri ) {
		A[0] = 1, A[1] = a[le];
		return ;
	}
	int mid = (le + ri) >> 1;
	func3(le, mid, dep + 1, &hl[dep][le]);
	func3(mid, ri, dep + 1, &hr[dep][le]);
	oper.poly_mul(&hl[dep][le], &hr[dep][le], A, mid-le+1, ri-mid+1, ri-le+1);
}
void func2() {
	func3(0, N, 0, h);
	for(int i=0,f=1;i<K;i++,f=1LL*f*(MOD-1)%MOD)
		h[i] = 1LL*f*h[i]%MOD;
	oper.poly_inv(h, g, K);
	for(int i=0;i<K;i++)
		h[i] = 1LL*(MOD-1)*h[i]%MOD*i%MOD;
	oper.poly_mul(h, g, g, K, K, K);
	g[0] = N;
	for(int i=0;i<K;i++)
		g[i] = 1LL*g[i]*ifct[i]%MOD;
}
int ans[MAXN + 5];
int main() {
	init();
	scanf("%d%d%lld", &N, &K, &T), K++, T++, T %= MOD;
	for(int i=0;i<N;i++)
		scanf("%d", &a[i]);
	func1(); func2(); oper.poly_mul(f, g, ans, K, K, K);
	for(int i=0;i<K;i++)
		printf("%lld\n", 1LL*fct[i]*ans[i]%MOD);
}

@details@

由於 fft 時咱們須要的是 2 的冪的長度，若是不另外建 temp 來存儲而是直接在原數組上存的話，須要注意兩個數組之間不會互相影響。 因此我也不知道爲何我不建個temp就行了