@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就行了