luogu P5488 差分与前缀和 FFT

又是一道FFT 好题。

首先来看一看求前缀和。

求一次前缀和就先当于卷上一个系数全为1的多项式,即(displaystyle sum_{i=0}^{infin}x^i)(想一想,为什么),这个东西就等于 (displaystyle frac{1}{1-x}),简单证明一下。

[S=1+x+x^2+...\ xS=x+x^2+x^3...\ 上边的式子减去下边的式子得到\ S-xS=1\ S=frac{1}{1-x} ]

是不是感觉天衣无缝但又十分扯淡?没错,这只有在(-1<x<1)时才成立。但生成函数是形式幂级数,我们不用关心x的具体取值。

回到我们刚才的问题,因为卷积具有结合律,所以我们要求的就是(displaystylefrac{1}{(1-x)^k})

结论:这个式子的n次系数是 (C_{n+k-1}^{k-1})(C是组合数)。

证明:回想一下(displaystyle sum_{i=0}^{infin}x^i)的每一次相乘的含义,可知(displaystyle (sum_{i=0}^{infin}x^i)^k)中n次系数的含义就是经过k次组成n的方案数,我们可以将n看成是n个小球,k看成是k个盒子,因为组成n的每个 “1”是一样的,每个多项式是不一样的,所以球相同,盒子不同,方案数就是(C_{n+k-1}^{k-1})

因为k很大,所以我们需要先让k对1004535809取模。至于为什么是可以的,可以回想一下P5245 多项式快速幂 这道题,考虑在模的意义下,对(k imes lnF(x))做多项式指数函数,显然是可以的。

但是取模后的k依然很大,我们将 (C_n^m)拆开(displaystyle C_n^m=frac{n!}{m!(n-m)!}=frac{prod_{i=n-m+1}^{n}i}{prod_{i=1}^{m}i}), 又因为(C_n^0=1)然后我们就能递推组合数啦。具体递推过程请参考代码。

接下来看一看求差分。

其实差分一次就相当于卷上(1-x)(想一想,为什么),所以我们要求的就是((1-x)^k).

求法1: 是不是感觉这个式子似曾相识?,没错,我们其实只用对刚刚求(displaystylefrac{1}{(1-x)^k})求个逆就行了,

求法2:根据二项式定理(displaystyle (1-x)^k=sum_{i=0}^kC_{k}^{i}(-x)^i=sum_{i=0}^kC_{k}^{i}(-1)^ix^i).然后就能求啦。

1004535809的原根是3.

#include<bits/stdc++.h>
#define LL long long
using namespace std;
int n, opt, len;
LL k;
const int N = 400010, mod = 1004535809, G = 3, Ginv = (mod + 1) / 3;
int r[N];
LL a[N], b[N], c[N], Y[N];
char ch[2501];
int read() 
{
	int x = 0; int f = 0; char c = getchar();
	while (!isdigit(c)) f |= c == '-', c = getchar();
	while (isdigit(c)) x = (x << 3) + (x << 1) + (c ^ 48), c = getchar();
	return f ? -x : x;
}
LL ksm(LL a, LL b, LL mod) 
{
	LL res = 1;
	for (; b; b >>= 1, a = a * a % mod)
		if (b & 1)res = res * a % mod;
	return res;
}
void NTT(LL *A, int lim, int opt) 
{
	for (int i = 0; i < lim; ++i)
		r[i] = (r[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0);
	for (int i = 0; i < lim; ++i)
		if (i < r[i])swap(A[i], A[r[i]]);
	int len;
	LL wn, w, x, y;
	for (int mid = 1; mid < lim; mid <<= 1) 
	{
		len = mid << 1;
		wn = ksm(opt == 1 ? G : Ginv, (mod - 1) / len, mod);
		for (int j = 0; j < lim; j += len) 
		{
			w = 1;
			for (int k = j; k < j + mid; ++k, w = w * wn % mod) 
			{
				x = A[k]; y = A[k + mid] * w % mod;
				A[k] = (x + y) % mod;
				A[k + mid] = (x - y + mod) % mod;
			}
		}
	}
	if (opt == 1)return;
	int ni = ksm(lim, mod - 2, mod);
	for (int i = 0; i < lim; ++i)A[i] = A[i] * ni % mod;
}
void MUL(LL *A, int n, LL *B, int m) 
{
	int lim = 1;
	while (lim < (n + m))lim <<= 1;
	NTT(A, lim, 1); NTT(B, lim, 1);
	for (int i = 0; i < lim; ++i)A[i] = A[i] * B[i] % mod;
	NTT(A, lim, -1);
}
void INV(int siz, LL *A, LL *B) 
{
	if (siz == 1) 
	{
		B[0] = ksm(A[0], mod - 2, mod);
		return;
	}
	INV((siz + 1) >> 1, A, B);
	int lim = 1;
	while (lim < (siz << 1))lim <<= 1;
	for (int i = 0; i < siz; ++i)c[i] = A[i];
	for (int i = siz; i < lim; ++i)c[i] = 0;
	NTT(c, lim, 1); NTT(B, lim, 1);
	for (int i = 0; i < lim; ++i)B[i] = B[i] * (2 - c[i] * B[i] % mod + mod) % mod;
	NTT(B, lim, -1);
	for (int i = siz; i < lim; ++i)B[i] = 0;
}
int main() 
{
	cin >> n; scanf("%s", ch + 1); cin >> opt;
	len = strlen(ch + 1);
	for (int i = 1; i <= len; ++i)k = k * 10 + ch[i] - '0', k %= mod;
	for (int i = 1; i <= n; ++i)a[i] = read();
	b[0] = 1;
	for (int i = 1; i <= n; ++i)b[i] = b[i - 1] * (i + k - 1) % mod * ksm(i, mod - 2, mod) % mod;
	if (opt == 0) 
	{
		MUL(a, n, b, n);
		for (int i = 1; i <= n; ++i)printf("%lld ", a[i]);
	} 
	else 
	{
		INV(n + 1, b, Y);
		MUL(a, n, Y, n);
		for (int i = 1; i <= n; ++i)printf("%lld ", a[i]);
	}
	return 0;
}
原文地址:https://www.cnblogs.com/wljss/p/12020771.html