HDU 6061 推导 NTT

复函数,递归代入,可以得到最终的式子为$f(x-sum_{i=1}^{m}{a_i})$,且$f(x) = sum_{i = 0}^{n}{c_ix^i}$,求最终各个x项的系数。

设$S=sum_{i=1}^{m}{a_i}$

先二项式展开

egin{eqnarray*} f(x-S)&=&sum_{i=0}^{n}{c_i{(x-S)}^i} ewline &=&sum_{i=0}^{n}{ c^isum_{j=0}^{i}{ inom{j}{i}x^{j}(-S)^{i-j} } }end{eqnarray*}

交换求和符号

egin{eqnarray*}&=&sum_{j=0}^{n}{ x^{j} sum_{i=0}^{n-j}{c_{i+j}inom{j}{i+j}(-S)^{i}} } ewline &=&sum_{j=0}^{n}{ x^{j} sum_{i=0}^{n-j}{c_{i+j}frac{(i+j)!}{{i!}{j!}}(-S)^{i}} } ewline &=&sum_{j=0}^{n}{ frac{x^{j}}{j!} sum_{i=0}^{n-j}{c_{i+j}(i+j)! frac{(-S)^{i}}{i!} } }end{eqnarray*}

因为只要系数,也就是说对于每个j,得到对应的${j!}b_j$,求后一个求和符号里的值就行了

注意到卷积公式$sum^{n}_{i=0}{f(i)h(n-i)}= f(n)*h(n)$,而且模数998244353是费马素数

故设$k = n - j, A[i] = frac{(-S)^{i}}{i!},B[i]=c_{k-i-i}(k-i-j)! $,则$C[n - j] =sum_{i=0}^{j}A[i]B[n-j-i]=sum_{i=0}^{n-j} c_{i+j}(i+j)! frac{(-S)^{i}}{i!} $

对A和B使用NTT加速卷积的计算,最后结果在A[i]上,记得反转为A[n-i]。推公式化卷积式,套模板。

#include <bits/stdc++.h>
#define LL long long
using namespace std;

const LL N = 1e5 + 10;
const LL mod = 998244353;
const int g = 3;
const int maxlen = 1 << 18;


LL wn[maxlen], fac[N], inv[N];

LL fpow(LL a, LL n)
{
	LL ans = 1;
	while(n)
	{
		if(n & 1)
			ans = (ans * a % mod + mod) %mod;
		a = (a * a + mod) % mod;
		n >>= 1;
	}
	return ans;
}

void init()
{
	wn[0] = 1;
	wn[1] = fpow(g, ((mod - 1)>>18));
	for(int i = 2; i < maxlen; i++)
		wn[i] = wn[i - 1] * wn[1] % mod;
	fac[1] = inv[1] = 1;
	fac[0] = inv[0] = 1;
	for(int i = 2; i < N; i++)
	{
		fac[i] = fac[i - 1] * i % mod;
		inv[i] = (mod - mod / i) * inv[mod % i] % mod;
	}
	for(int i = 1; i < N; i++)
		(inv[i] *= inv[i - 1]) %= mod;
}

void rader(LL f[], int len)
{
	for(int i = 1, j = len >> 1; i < len - 1; i++)
	{
		if(i < j) swap(f[i], f[j]);
		int k = len >> 1;
		while(j >= k)
		{
			j -= k;
			k >>= 1;
		}
		if(j < k) j += k;
	}
}
void ntt(LL f[], int len, int on)
{
	/*for(int i = 0, j = 0; i < len;i++)
	{
		if(i > j) swap(f[i], f[j]);
		for(int l = len >> 1; (j^=l) < l; l>>=1);
	}*/
	rader(f, len);
	for(int i = 1, d = 1; d < len; i++, d <<= 1)
	{
		//LL wnn = fpow(g, (mod-1)/(d<<1));
		for(int j = 0; j < len; j += (d << 1))
		{
			//LL w = 1;
			for(int k = 0; k < d; k++)
			{
				LL t = wn[(maxlen >> i) * k] * f[j + k + d] % mod;
				//LL t = w*f[j+k+d]%mod;
				//w = w*wnn % mod;
				f[j + k + d] = ((f[j + k] - t) % mod + mod) % mod;
				f[j + k] = ((f[j + k] + t) % mod + mod) % mod;
			}
		}
	}
	if(on == -1)
	{
		reverse(f + 1, f + len);
		LL inv2 = fpow(len, mod - 2);
		for(int i = 0; i < len; i++)
			f[i] = f[i] % mod * inv2 % mod;
	}
}

void work(LL a[], LL b[], int len)
{
	ntt(a, len, 1);
	ntt(b, len, 1);
	for(int i = 0; i < len; i++)
		a[i] = (a[i] * b[i] % mod + mod) % mod;
	ntt(a, len, -1);
}

LL A[maxlen], B[maxlen], Suma, c[N];
int n, m; 
int main()
{
	init();
	while(~scanf("%d", &n))
	{
		for(int i = 0; i <= n; i++) scanf("%lld", c + i);
		scanf("%d", &m);
		Suma = 0;
		LL t;
		for(int i = 0; i < m; i++)
			scanf("%lld", &t), Suma -= t;
		Suma = (Suma + mod) % mod;
		while(Suma < 0)
			Suma += mod;
		if(Suma == 0)
		{
			for(int i = 0; i <= n; i++) 
				printf("%lld ", c[i]);
			printf("
");
			continue;
		}
		//getLength
		int len = 1;
		while(len <= 2*n) len <<= 1;
		//
		LL ae = 1;
		for(int i = 0; i < len; i++)
		{
			if(i <= n)
			{
				B[i] = ae * inv[i] % mod;
				A[i] =(fac[n - i] * c[n - i]) % mod;
			}
			else A[i] = B[i] = 0;
			ae = (ae * Suma) % mod;
			//cout << A[i] << "~"<< B[i] << endl;
		}
		work(A, B, len);
		for(int i = 0; i <= n; i++)
			A[i] = A[i] * inv[n - i] % mod;
		for(int i = 0; i <= n; i++)
			printf("%lld ", A[n - i]);
		printf("
");
	}
}


原文地址:https://www.cnblogs.com/Yumesenya/p/7470832.html