FFT/NTT卷积神级副本

因为对卷积还不是很清楚，所以开始疯狂的刷题之路，会慢慢更新的。
前面没有来源标注的题目题号洛谷上面都有。

$FFT,NTT$是卷积运算中常见的优化$\to O(nlogn)$
因为文字量极大，很有可能会有细节错误，欢迎指出，谢谢o(=·ω·=)m

传送至只有小怪的村庄——请开始你的逆天之路

A：P1919

【模板】A*B Problem升级版（FFT快速傅里叶）

将输入的$a,b$每一位拆成对应的多项式系数

手玩一下普通乘法的计算法则，发现从左到右对应多项式次数从$0,n-1$依次递增

两数相乘是符合卷积形式，下标是相加的

所以将这个串反转，从左到右从$0$开始编号

然后套$FFT$跑出乘积结果，再转回系数表达式

此时某些位置上的值可能很大，这就涉及到进位的问题

从左到右扫一遍即可完成进位

最后倒着输出便是$AC$

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define maxn 4000100
struct complex {
	double x, i;
	complex(){}
	complex( double X, double I ) {
		x = X, i = I;
	}
}A[maxn], B[maxn];

double pi = acos( -1.0 );

complex operator + ( complex a, complex b ) {
	return complex( a.x + b.x, a.i + b.i );
}

complex operator - ( complex a, complex b ) {
	return complex( a.x - b.x, a.i - b.i );
}

complex operator * ( complex a, complex b ) {
	return complex( a.x * b.x - a.i * b.i, a.x * b.i + a.i * b.x );
}

char a[maxn], b[maxn];
int r[maxn], ans[maxn];
int len;

void FFT( complex *v, int opt ) {
	for( int i = 0;i < len;i ++ )
		if( i < r[i] ) swap( v[i], v[r[i]] );
	for( int i = 1;i < len;i <<= 1 ) {
		complex omega( cos( pi / i ), opt * sin( pi / i ) );
		for( int j = 0;j < len;j += ( i << 1 ) ) {
			complex w( 1, 0 );
			for( int k = 0;k < i;k ++, w = w * omega ) {
				complex x = v[j + k], y = v[j + k + i] * w;
				v[j + k] = x + y;
				v[j + k + i] = x - y;
			}
		}
	}
}

int main() {
	scanf( "%s %s", a, b );
	int n = strlen( a ), m = strlen( b );
	for( int i = 0;i < n;i ++ ) A[n - i - 1] = complex( a[i] - '0', 0 );
	for( int i = 0;i < m;i ++ ) B[m - i - 1] = complex( b[i] - '0', 0 );
	len = 1; int l = 0;
	while( len < ( n + m ) ) len <<= 1, l ++;
	for( int i = 0;i < len;i ++ )
		r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
	FFT( A, 1 );
	FFT( B, 1 );
	for( int i = 0;i < len;i ++ )
		A[i] = A[i] * B[i];	
	FFT( A, -1 );
	for( int i = 0;i < len;i ++ )
		ans[i] = ( A[i].x / len + 0.5 );
	for( int i = 0;i < len;i ++ )
		if( ans[i] > 9 ) ans[i + 1] += ans[i] / 10, ans[i] %= 10;
	while( ! ans[len - 1] ) len --;
	for( int i = len - 1;~ i;i -- )
		printf( "%d", ans[i] );
	return 0;
}

---

B：P4157

[SCOI2006]整数划分

这道题显然大整数是可做的

但是都学到现在了，一般我是不怎么打大整数的，而且现在可能也打不出来

有个结论，尽可能得分$3$，最后若是$4/2$，就不分了

这个结论不需要证吧，考场上肯定是找规律找出来的

可以套在$FFT$上，不就是多项式只有个常数项罢了，照样可以算

$x$个$3$求乘积，就对$FFT$进行快速幂，处理进位问题即可

---

刷怪升级——转战玄灵大陆

C：P6300

悔改

第一次见到这种应用，记在小本本上！！ ✍

想了一会儿，列了个$n^3$的$dp$转移，只能搞$5$分，✧*｡٩(^㉨)و✧｡ 果断放弃

然后，木头拼接起来长度相加，突然想到了生成函数

发现可以用长度为$i$的木头数量作为$x^i$的系数，两端拼在一起就是生成函数的二次幂

火速开敲，马上测样例，(⊙_⊙;)？$39$怎么肥事！！

仔细想了一下，ε=(･д･｀*)ﾊｧ…两个长度为$1$的木头和一个长度为$8$的木头拼接了，算成了两根新木头，但是按题目来说必须有一个长度$1$木头不用

<。)#)))≦我陷入了沉思…… 🤔

我应该取两种拼接木头个数更少的一方，但是这样子怎么写呢？？

---

到此为止，我走向了题解： 一类取min类型的卷积到普通乘法卷积的转化

${\sum }_{i+j=k}min(to{t}_i,to{t}_j)={\sum }_{d=1}{\sum }_{i+j=k}[a_i\ge d][a_j\ge d]$

也就是说，枚举木头个数，每次都卷积

$F_d(x)=\sum [a_i\ge d]x^i$

$an{s}_k=\frac{1}{2}{\sum }_{d=1}[x^k]F_d(x{)}^2$

将$tot$数组离散化直接暴力卷积

据题解说：本质不同的$F_d$是$O(\sqrt{n})$，最坏情况为$1+2+...=O(D^2)$

复杂度$O(m\sqrt{n}logm)$常数不大

现在才明白为什么有的题目明明没有取模，也有题解代码用的$NTT$

因为$NTT$不用敲复数那么多玩意儿

只需要取一个非常大的模数使得答案永远不可能超过即可

妙(・。・)！

#include <cmath>
#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
#define maxn 400100
struct complex {
	double x, i;
	complex(){}
	complex( double X, double I ) { 
		x = X, i = I;
	}
}A[maxn];

const double pi = acos( -1.0 );

complex operator + ( complex a, complex b ) {
	return complex( a.x + b.x, a.i + b.i );
}

complex operator - ( complex a, complex b ) {
	return complex( a.x - b.x, a.i - b.i );
}

complex operator * ( complex a, complex b ) {
	return complex( a.x * b.x - a.i * b.i, a.x * b.i + a.i * b.x );
}

int n, m, len;
int tot[maxn], r[maxn], ans[maxn], tmp[maxn];

void FFT( complex *v, int opt ) {
	for( int i = 0;i < len;i ++ )
		if( i < r[i] ) swap( v[i], v[r[i]] );
	for( int i = 1;i < len;i <<= 1 ) {
		complex omega( cos( pi / i ), opt * sin( pi / i ) );
		for( int j = 0;j < len;j += ( i << 1 ) ) {
			complex w( 1, 0 );
			for( int k = 0;k < i;k ++, w = w * omega ) {
				complex x = v[j + k], y = v[j + k + i] * w;
				v[j + k] = x + y;
				v[j + k + i] = x - y;
			}
		}
	}
}

int main() {
	scanf( "%d %d", &n, &m );
	for( int i = 1, x;i <= n;i ++ ) {
		scanf( "%d", &x );
		tot[x] ++;
	}
	len = 1; int l = 0;
	while( len <= ( m << 1 ) ) len <<= 1, l ++;
	for( int i = 0;i < len;i ++ )
		r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
	for( int i = 1;i <= m;i ++ ) tmp[i] = tot[i];
	sort( tmp + 1, tmp + m + 1 );
	int T = unique( tmp + 1, tmp + m + 1 ) - tmp - 1;
	for( int i = 1;i <= T;i ++ ) {
		for( int j = 0;j < len;j ++ )
			A[j] = complex( tot[j] >= tmp[i], 0 );
		FFT( A, 1 );
		for( int j = 0;j < len;j ++ )
			A[j] = A[j] * A[j];
		FFT( A, -1 );
		for( int j = 0;j < len;j ++ )
			ans[j] += ( tmp[i] - tmp[i - 1] ) * int( A[j].x / len + 0.5 );
	}
	int maxx = 0, ansLen;
	for( int i = 0;i < len;i ++ )
		if( ans[i] / 2 > maxx ) maxx = ans[i] / 2, ansLen = i;
	printf( "%d %d\n", maxx, ansLen );
	return 0;
}

---

D：P3763

[TJOI2017]DNA

曾经放在后缀数组板块的题又被拉出来处决了，不，应该是来处决我了🙅

字符集大小只有 { A , T , G , C } \{A,T,G,C\} {A,T,G,C}，考虑分开枚举处理

---

等于枚举字符的设为$1$，反之为$0$

以样例$‘C’$为例

$S_0:ATCGCCCTA$

$S:CTTCA$

$S_0:001011100$

$S:10010$

假设$S_0$的长度为$n$，$S$的长度为$m$

对于$S_0[k,k+m-1]$与$S[0,m-1]$相同字符的个数则表示${\sum }_{i=0}^{m-1}{S}_0[i+k]\ast S[i]$

令$T[m-i-1]=S[i]$，代替$S$，式子可化为${\sum }_{i=0}^{m-1}{S}_0[i+k]\ast T[m-i-1]$

发现此时是可以$FFT$卷积的，将下标看作多项式$x^i$

${\sum }_{i+j=m+k-1}{S}_0[i]\ast T[j]$

对于四种字符，$FFT$做四遍，加在一起就是子串与碱基序列相同字符的个数了

最后判断相差是否能控制在$3$以内，就知道这个子串碱基能否被改成吃藕碱基了

卷积后的$x^{m+k-1}$项才代表原碱基序列中$[k,k+m-1]$与吃藕碱基的匹配，要转化一下

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define maxn 400100 
struct complex {
	double x, i;
	complex(){}
	complex( double X, double I ) {
		x = X, i = I;
	}
}A[maxn], B[maxn];

const double pi = acos( -1.0 );

complex operator + ( complex a, complex b ) {
	return complex( a.x + b.x, a.i + b.i );
}

complex operator - ( complex a, complex b ) {
	return complex( a.x - b.x, a.i - b.i );
}

complex operator * ( complex a, complex b ) {
	return complex( a.x * b.x - a.i * b.i, a.x * b.i + a.i * b.x );
}

int len;
int r[maxn];

void FFT( complex *v, int opt ) {
	for( int i = 0;i < len;i ++ )
		if( i < r[i] ) swap( v[i], v[r[i]] );
	for( int i = 1;i < len;i <<= 1 ) {
		complex omega( cos( pi / i ), opt * sin( pi / i ) );
		for( int j = 0;j < len;j += ( i << 1 ) ) {
			complex w( 1, 0 );
			for( int k = 0;k < i;k ++, w = w * omega ) {
				complex x = v[j + k], y = v[j + k + i] * w;
				v[j + k] = x + y;
				v[j + k + i] = x - y;
			}
		}
	}
}

int T;
char S0[maxn], S[maxn], ch[4] = { 'A', 'G', 'T', 'C' };
int tot[maxn];

int main() {
	scanf( "%d", &T );
	while( T -- ) {
		scanf( "%s %s", S0, S );
		memset( tot, 0, sizeof( tot ) );
		int n = strlen( S0 ), m = strlen( S );
		len = 1; int l = 0;
		while( len < ( n + m ) ) len <<= 1, l ++;
		for( int i = 0;i < len;i ++ )
			r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
		for( int j = 0;j < 4;j ++ ) {
			for( int i = 0;i < n;i ++ )
				A[i] = complex( S0[i] == ch[j], 0 );
			for( int i = 0;i < m;i ++ )
				B[m - i - 1] = complex( S[i] == ch[j], 0 );
			FFT( A, 1 );
			FFT( B, 1 );
			for( int i = 0;i < len;i ++ )
				A[i] = A[i] * B[i];
			FFT( A, -1 );
			for( int i = m - 1;i <= m + n - 1;i ++ ) 
				tot[i - m + 1] += int( A[i].x / len + 0.5 );
			for( int i = 0;i < len;i ++ )
				A[i] = B[i] = complex( 0, 0 );
		}
		int ans = 0;
		for( int i = 0;i < n - m + 1;i ++ ) //tot[i]:S0[i,i+m-1]和S[0,m-1]相同字符的个数 
			if( m - tot[i] <= 3 ) ans ++;
		printf( "%d\n", ans );
	}
	return 0;
}

---

E：P3321

[SDOI2015]序列统计

虽然以前做过一次，但是肯定是迷迷糊糊的，现在重新再做一遍，代码就找以前的博客吧哈哈哈哈哈,ԾㅂԾ,

相乘目前来说是超越已有工具的，考虑相乘转化为相加$\to$ 指数相加$\to$生成函数

但是好像没有哪个数学定理告诉我相乘转化为指数相加后取模答案不变的

卡在这里了(ˉ▽ˉ；)…

再考虑设$d{p}_{i,j}$表示选了$i$个数相乘的结果取模$m$等于$j$的方案数，则有很简单的转移

$f[i+1][j\times k\mathrm{m}\mathrm{o}\mathrm{d}m]={\sum }_{k\in S}f[i][j]$

好像也转换不了形式使其长得像$FFT$啊凸(艹皿艹 )

---

尽管做过，但是还是没想出来啊/(ㄒoㄒ)/，去看一下之前写的题解

b（￣▽￣）d　题目没说清楚$m$是个质数！！尽管说了我也不太想得到

指数的相加可以用$m$的原根来做底数

$dp$转移也稍有不同

$f[i+i][j]={\sum }_{j_1\times j_2\%m=j}f[i][j_1]\times f[i][j_2]$

但是大致想了$75\%$左右出来，已经相比于之前完全$0\%$好多了

---

F：P5641

【CSGRound2】开拓者的卓识、

(>ρ<")推半天推的是个什么玩意儿柿子噢

---

这个柿子从上往下推似乎推不出什么玩意儿，推了有一会儿确实不会处理

考虑反过来从下往上推

令$an{s}_r=su{m}_{k,1,r}$，考虑每个$a_i$的贡献，假设$a_i$出现了$c_i$次

显然答案表示为$an{s}_r={\sum }_{i=1}^n{a}_i\cdot c_i$

$c_i$还有一种含义理解：位置$i$在$[1,r]$中被$k-1$重区间覆盖的方案数

（减一是因为最大的区间已经是被确定为$[1,r]$

等价于选择$k-1$个区间$[l_1,r_1],[l_2,r_2]...[l_{k-1},r_{k-1}]$

使得$1\le l_i,r_i\le r,[l_i,r_i]\subseteq [l_{i+1},r_{i+1}],i\in [l_i,r_i]$的方案数

---

或者…༼ ∗ •̀ (oo) •́ ∗ ༽可以枚举理解

$su{m}_{1,l,r}$：就是$a_l$到$a_r$的和

$su{m}_{2,l,r}$：可以看作选择一个区间$[l_1,r_1]\subseteq [l,r]$，然后把$a_{l_1}$到$a_{r_1}$的和再加到答案里面

$su{m}_{3,l,r}$：选择一个区间$[l_2,r_2]\subseteq [l_1,r_1]$，再对其做一遍$su{m}_{2,l_2,r_2}$，相当于再选一个区间$[l_2,r_2]\subseteq [l_1,r_1]$，将$a_{l_2}$到$a_{r_2}$的和再累计到答案中

………….

$su{m}_{k,l,r}$：就是选择$k-1$个区间

---

继续等价于$l_i$在$[1,i]$中选，$r_i$在$[i,r]$中选，一共在$[1,i]$中选$k-1$个（可重位置，在$[i,r]$中选$k-1$个（可重位置，的方案数

令$T_n^m$表示在$n$个盒子中选择$m$个盒子的方案数

则$c_i$方案数即为$T_i^{k-1}\times T_{r-i+1}^{k-1}$

$T_n^m$也是非常好求的，就是数学中经典的隔板法

用$m-1$个板将盒子分成$m$份，每份合并看作是最后的一个超级盒子，那么一共就选择了$m$个超级盒子，符合$T$的含义

所以就可以用组合数来计算，此时就是不重位置的挑选了，$T_n^m=C_{n+m-1}^m$

那么$c_i=C_{(i)+(k-1)-1}^{k-1}\times C_{(r-i+1)+(k-1)-1}^{k-1}=C_{i+k-2}^{k-1}\times C_{r-i+k-1}^{k-1}$

$an{s}_r={\sum }_{i=1}^n{a}_i{C}_{i+k-2}^{k-1}\times C_{r-i+k-1}^{k-1}$

令$g_i=C_{i+k-1}^{k-1},f_i=C_{i+k-2}^{k-1}$

则$an{s}_r={\sum }_{i=1}^n{g}_{r-i}\ast f_i$

就可以卷积了，有模数选择$NTT$

由于$k$非常的大，预处理组合数显然不可做，那就使用组合数的递推式即可

巧妙发现$f_i=a_i\times g_{i-1}$，$f_i$也可递推获得

#include <cstdio>
#include <iostream>
using namespace std;
#define mod 998244353
#define int long long
#define maxn 400100
int n, k, len;
int a[maxn], r[maxn], f[maxn], g[maxn];

int qkpow( int x, int y ) {
	int ans = 1;
	while( y ) {
		if( y & 1 ) ans = ans * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return ans;
}

void NTT( int *v, int opt ) {
	for( int i = 0;i < len;i ++ )
		if( i < r[i] ) swap( v[i], v[r[i]] );
	for( int i = 1;i < len;i <<= 1 ) {
		int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ) );
		for( int j = 0;j < len;j += ( i << 1 ) )
			for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {
				int x = v[j + k], y = v[j + k + i] * w % mod;
				v[j + k] = ( x + y ) % mod;
				v[j + k + i] = ( x - y + mod ) % mod;
			}
	}
	if( opt == -1 ) {
		int inv = qkpow( len, mod - 2 );
		for( int i = 0;i < len;i ++ )
			v[i] = v[i] * inv % mod;
	}
}

signed main() {
	scanf( "%lld %lld", &n, &k );
	for( int i = 1;i <= n;i ++ )
		scanf( "%lld", &a[i] );
	g[0] = 1, g[1] = k % mod, f[1] = a[1];
	for( int i = 2;i <= n;i ++ ) {
		g[i] = g[i - 1] * ( i + k - 1 ) % mod * qkpow( i, mod - 2 ) % mod;
		f[i] = a[i] * g[i - 1] % mod;
	}
	len = 1; int l = 0;
	while( len <= ( n << 1 ) ) len <<= 1, l ++;
	for( int i = 0;i < len;i ++ )
		r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
	NTT( f, 1 );
	NTT( g, 1 );
	for( int i = 0;i < len;i ++ )
		f[i] = f[i] * g[i] % mod;
	NTT( f, -1 );
	for( int i = 1;i <= n;i ++ )
		printf( "%lld ", f[i] );
	return 0;
}

---

G：P4986

逃离

题意大概是这个亚子。走的范围是一个环，红线代表$C(x)\ast t$，绿线代表$A(x)\ast t$，蓝线代表$B(x)\ast t$

非常小学鸡化容易操作，将蓝绿线各自平移成一条线段，然后与圆的半径构成直角三角形，要从一个点一起出去，则斜边就是$ hdxrie $的走法，直角边就是$ Althen $的走法

根据直角三角形法则，有$(A(x)\ast t{)}^2+(B(x)\ast t{)}^2=(C(x)\ast t{)}^2$

化简得$A^2(x)+B^2(x)=C^2(x)$

定义一个新函数$f(x)=C^2(x)-A^2(x)-B^2(x)$

最后的即是求函数$f(x)$的零点，当然必须满足$x\in [L,R]$

而求一个函数的根近似解，套用一阶牛顿迭代

$x=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$

---

随便找一个点$x_0$，然后做这个点的切线，斜率就是对点$x_0$求导（$f’(x_0)$）

从这个切线的根（与$x$轴的交点）出发，做一根垂线，和曲线相交于$x$点

不停重复上述过程，直到逼近原函数的零点附近

稍微解释一下牛顿迭代的公式

对于一条直线解析式$y=kx+b$，那么关于$x_0$的切线即为$f(x_0)=f^{\prime }(x_0)x_0+b$

根据初中知识知道，切线斜率还可以写成$\frac{y-y_0}{x-x_0}$，那么用切线的零点和$x_0$重写切线解析式$f^{\prime }(x_0)=\frac{0-f(x_0)}{x-x_0}\Rightarrow x=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$

最后就是设置一个牛顿迭代次数$ti$，暴算去逼近零点（如果有的话

当然这过程中要和$L,R$比大小保证是答案是在范围内的

#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
#define maxn 400100
#define eps 1e-10
struct complex {
	double x, i;
	complex(){}
	complex( double X, double I ) {
		x = X, i = I;
	}
}A[maxn];

const double pi = acos( -1.0 );

complex operator + ( complex a, complex b ) {
	return complex( a.x + b.x, a.i + b.i );
}

complex operator - ( complex a, complex b ) {
	return complex( a.x - b.x, a.i - b.i );
}

complex operator * ( complex a, complex b ) {
	return complex( a.x * b.x - a.i * b.i, a.x * b.i + a.i * b.x );
}

int len;
int r[maxn];

void FFT( complex *v, int opt ) {
	for( int i = 0;i < len;i ++ )
		if( i < r[i] ) swap( v[i], v[r[i]] );
	for( int i = 1;i < len;i <<= 1 ) {
		complex omega( cos( pi / i ), opt * sin( pi / i ) );
		for( int j = 0;j < len;j += ( i << 1 ) ) {
			complex w( 1, 0 );
			for( int k = 0;k < i;k ++, w = w * omega ) {
				complex x = v[j + k], y = v[j + k + i] * w;
				v[j + k] = x + y, v[j + k + i] = x - y;
			}
		}
	}
}

int La, Lb, Lc, n;
double L, R;
int a[maxn], b[maxn], c[maxn], f[maxn], g[maxn];

void work( int &n, int *ans ) {
	len = 1; int l = 0;
	while( len <= ( n << 1 ) ) len <<= 1, l ++;
	for( int i = 0;i < len;i ++ )
		r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
	for( int i = 0;i < len;i ++ )
		A[i] = complex( ans[i], 0 );
	FFT( A, 1 );
	for( int i = 0;i < len;i ++ )
		A[i] = A[i] * A[i];
	FFT( A, -1 );
	for( int i = 0;i < len;i ++ )
		ans[i] = int( A[i].x / len + 0.5 );
	n = len;
}

double calc_f( double x ) {
	double ans = 0, mi = 1;
	for( int i = 0;i <= n;i ++ )
		ans += f[i] * mi, mi *= x;
	return ans;
}

double calc_g( double x ) {
	double ans = 0, mi = 1;
	for( int i = 0;i < n;i ++ )
		ans += g[i] * mi, mi *= x;
	return ans;
}

int main() {
	scanf( "%d %d %d %lf %lf", &La, &Lb, &Lc, &L, &R );
	for( int i = 0;i <= La;i ++ )
		scanf( "%d", &a[i] );
	for( int i = 0;i <= Lb;i ++ )
		scanf( "%d", &b[i] );
	for( int i = 0;i <= Lc;i ++ )
		scanf( "%d", &c[i] );
	work( La, a );
	work( Lb, b );
	work( Lc, c );
	n = max( La, max( Lb, Lc ) );
	for( int i = 0;i <= n;i ++ )
		f[i] = c[i] - a[i] - b[i];
	for( int i = 0;i < n;i ++ )
		g[i] = f[i + 1] * ( i + 1 );
	int ti = 50;
	double x = ( L + R ) / 2.0;
	while( ti ) {
		double tmp = calc_f( x );
		if( fabs( tmp ) < eps ) break;
		x = x - tmp / calc_g( x );
		x = max( x, L ), x = min( x, R );
		ti --;
	}
	if( ti ) printf( "%.8f\n", x );
	else printf( "Inconsistent!" );
	return 0;
}

---

H：P4721——获得知识药剂一瓶——分治$FFT$

【模板】分治 FFT

对于蒟蒻来说还是算一个新知识，嘚学Y(>_<、)Y

分治FFT是一种基于CDQ分治的算法，主要用于在$O(nlo{g}^{2}n)$的时间复杂度计算已知$g_1$到$g_{n-1}$，$f_0=1$，求$f_i={\sum }_{j=1}^i{f}_{i-j}\cdot g_j$

考虑分治$l,r,mid=\lfloor \frac{l+r}{2}\rfloor$，计算左半段$[l,mid]$对于右半段$[mid+1,r]$的贡献

设$x\in [mid+1,r]$，那么对$f_x$的贡献转移$w_x$即为$w_x={\sum }_{i=l}^{mid}{f}_i{g}_{x-i}$

不妨将推导范围扩大一点，先保证$f_i=0,i\in [mid+1,r]$

$w_x={\sum }_{i=l}^{x-1}{f}_i{g}_{x-i}={\sum }_{i=0}^{x-l-1}{f}_{i+l}\cdot g_{x-l-i}$

令$A_i=f_{i+l},B_i=g_{i+1}$，代回上式$w_x={\sum }_{i=0}^{x-l-1}{A}_i{B}_{x-l-i-1}$

发现此时的$A,B$可以进行$FFT$卷积

每次操作的长度为$r-l+1$，总长度$nlogn$，一共$logn$层，所以时间复杂度大概在$O(nlo{g}^{2}n)$

#include <cstdio>
#include <iostream>
using namespace std;
#define mod 998244353
#define int long long
#define maxn 400100
int r[maxn], f[maxn], g[maxn], A[maxn], B[maxn];

int qkpow( int x, int y ) {
	int ans = 1;
	while( y ) {
		if( y & 1 ) ans = ans * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return ans;
}

void NTT( int *v, int opt, int len ) {
	for( int i = 0;i < len;i ++ )
		if( i < r[i] ) swap( v[i], v[r[i]] );
	for( int i = 1;i < len;i <<= 1 ) {
		int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ) );
		for( int j = 0;j < len;j += ( i << 1 ) )
			for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {
				int x = v[j + k], y = v[j + k + i] * w % mod;
				v[j + k] = ( x + y ) % mod;
				v[j + k + i] = ( x - y + mod ) % mod;
			}
	}
	if( opt == -1 ) {
		int inv = qkpow( len, mod - 2 );
		for( int i = 0;i < len;i ++ )
			v[i] = v[i] * inv % mod;
	}
}

void mul( int *A, int *B, int n, int l ) {
	for( int i = 0;i < n;i ++ ) 
		r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
	NTT( A, 1, n );
	NTT( B, 1, n );
	for( int i = 0;i < n;i ++ ) A[i] = A[i] * B[i] % mod;
	NTT( A, -1, n );
}

void CDQ( int L, int R ) {
	if( L == R ) return;
	int mid = ( L + R ) >> 1;
	CDQ( L, mid );
	int len = 1, l = 0;
	while( len <= R - L + 1 ) len <<= 1, l ++;
	for( int i = L;i <= mid;i ++ ) A[i - L] = f[i];
	for( int i = 1;i <= R - L;i ++ ) B[i - 1] = g[i];
	mul( A, B, len, l );
	for( int i = mid + 1;i <= R;i ++ )
		f[i] = ( f[i] + A[i - L - 1] ) % mod;
	for( int i = 0;i < len;i ++ ) A[i] = B[i] = 0;
	CDQ( mid + 1, R );
}

signed main() {
	int n;
	scanf( "%lld", &n );
	for( int i = 1;i < n;i ++ )
		scanf( "%lld", &g[i] );
	f[0] = 1;
	CDQ( 0, n - 1 );
	for( int i = 0;i < n;i ++ )
		printf( "%lld ", f[i] );
	return 0;
}

---

I：P3338

[ZJOI2014]力

$$F_j=\sum _{i=1}^{j-1}\frac{q_i\times q_j}{(i-j{)}^2}-\sum _{i=j+1}^n\frac{q_i\times q_j}{(i-j{)}^2},E_j=\frac{F_j}{q_j}\Rightarrow E_j=\sum _{i=1}^{j-1}\frac{q_i}{(i-j{)}^2}-\sum _{i=j+1}^n\frac{q_i}{(i-j{)}^2}$$

把$i=j$加进去，其实是不影响的，代码又不会真的去除，我们只一心想要卷积

${\sum }_{i=1}^j\frac{q_i}{(i-j{)}^2}-{\sum }_{i=j}^n\frac{q_i}{(i-j{)}^2}$

设$f_i=q_i,g_i=\frac{1}{i^2}$

${\sum }_{i=1}^j{f}_i{g}_{j-i}-{\sum }_{i=j}^n{f}_i{g}_{i-j}$

令$f_0=g_0=0$

${\sum }_{i=0}^j{f}_i\ast g_{j-i}-{\sum }_{i=j}^n{f}_i{g}_{i-j}$

左边已经是一个卷积的形式了，接下来只考虑右边即可

一般卷积形式都会把求和符号的下标换成从$0$开始

并且用两个数组的下标去凑成一个求和上限的恒定值

${\sum }_{i=j}^n{f}_i{g}_{i-j}={\sum }_{i=0}^{n-j}{f}_{i+j}{g}_i$

设$h_i=f_{n-i}$

${\sum }_{i=0}^{n-j}{h}_{n-i-j}\ast g_i$

两边都变成了完美的卷积形式(๑′ᴗ‵๑)

${\sum }_{i=0}^{j}f[i]\ast g[j-i]-{\sum }_{i=0}^{n-j}h[n-i-j]\ast g[i]$

设多项式$F(x)={\sum }_{i=0}^n{f}_i,G(x)={\sum }_{i=0}^n{g}_i,H(x)={\sum }_{i=0}^n{h}_i$

令$L(x)=F(x)\ast G(x),R(x)=G(x)\ast H(x)$

则$an{s}_i=l_i-r_{n-i}$，$l_i,r_{n-i}$分别表示多项式对应项的系数

#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
#define maxn 400100
struct complex {
	double x, i;
	complex(){}
	complex( double X, double I ) {
		x = X, i = I;
	}
}F[maxn], G[maxn], H[maxn];

const double pi = acos( -1.0 );

complex operator + ( complex a, complex b ) {
	return complex( a.x + b.x, a.i + b.i );
}

complex operator - ( complex a, complex b ) {
	return complex( a.x - b.x, a.i - b.i );
}

complex operator * ( complex a, complex b ) {
	return complex( a.x * b.x - a.i * b.i, a.x * b.i + a.i * b.x );
}

int len;
int r[maxn];

void FFT( complex *v, int opt ) {
	for( int i = 0;i < len;i ++ )
		if( i < r[i] ) swap( v[i], v[r[i]] );
	for( int i = 1;i < len;i <<= 1 ) {
		complex omega( cos( pi / i ), opt * sin( pi / i ) );
		for( int j = 0;j < len;j += ( i << 1 ) ) {
			complex w( 1, 0 );
			for( int k = 0;k < i;k ++, w = w * omega ) {
				complex x = v[j + k], y = v[j + k + i] * w;
				v[j + k] = x + y;
				v[j + k + i] = x - y;
			}
		}
	}
}

int n;
int main() {
	scanf( "%d", &n );
	for( int i = 1;i <= n;i ++ ) {
		scanf( "%lf", &F[i].x ), F[i].i = 0;
		H[n - i] = F[i];
		G[i] = complex( 1.0 / i / i, 0 );
	}
	len = 1; int l = 0;
	while( len <= ( n << 1 ) ) len <<= 1, l ++;
	for( int i = 0;i < len;i ++ )
		r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
	FFT( F, 1 );
	FFT( G, 1 );
	FFT( H, 1 );
	for( int i = 0;i < len;i ++ ) F[i] = F[i] * G[i], H[i] = H[i] * G[i];
	FFT( F, -1 );
	FFT( H, -1 );
	for( int i = 1;i <= n;i ++ )
		printf( "%.3f\n", F[i].x / len - H[n - i].x / len );
	return 0;
}

---

J：P4173

残缺的字符串

读完题，迅速反应类似的本系列的D题DNA

考虑分成$26$个字母独立考虑（包括$\ast$，所在位置一直为$1$即可

将当前枚举的字母所在位置赋为$1$，其余为$0$

对于$B[k,k+m-1]$和$A[0,m-1]$的相同的字符长度个数表示为${\sum }_{i=0}^{m-1}B[k+i]A[i]$

令$T_{m-i-1}=A_i$

${\sum }_{i=0}^{m-1}{B}_{k+i}\ast T_{m-i-1}$

同样的卷积

${\sum }_{i+j=m-k-1}{B}_i\ast T_j$

最后$to{t}_i=[m-i-1]B\ast T$，将$27$个相同字符个数加起来，判断是否长度达到$m$即可

麻溜敲完，很不幸——$30$，剩下全是$TLE$，╮(╯-╰)╭好吧…

你以为是再次中场休息？？还是太年轻了骚年；又是一个新的好算法呢

---

——$luogu$第一篇题解真的很良心٩(๑`н´๑)۶

普通的单模式串匹配

长度为$m$的模式串$A$，长度为$n$的文本串$B$，求所有位置，满足$B$串从第$x$个字符开始的连续$m$个字符与$A$完全相同

Step1：定义匹配函数

定义字符串的匹配函数为$C(x,y)=A(x)-B(y)$

在这里形式化的定义“匹配”，若$C(x,y)=0$，则$A$的第$x$个字符与$B$的第$y$个字符相匹配

Step2：定义完全匹配函数

定义$D(x)={\sum }_{i=0}^{m-1}C(i,x-m+i+1)$，若$D(x)=0$，则$A$与$B$以$x$结尾的连续$m$个字符完全匹配

但到这里，发现因为带了求和符号，所以只要两串所含字符一样，势必会有$D(x)=0$

这就与我们的初衷背道而驰了，而导致这种情况发生的原因——$C(x,y)$自带$\pm$符号

考虑如何去掉这个符号的影响，使得只要两个串的某一位置不匹配，完全匹配函数就一定不为$0$

非常初中数学老师教的——加绝对值，但是在解不等式的时候，老师教过绝对值也不好弄，处理成——平方

此时有，$D(x)={\sum }_{i=0}^{m-1}C(i,x-m+i+1{)}^2={\sum }_{i=0}^{m-1}(A(i)-B(x-m+i+1){)}^2$

$Step3:$暴力整式子，观察结构，不谋手段の优化

先展开看看

$D(x)={\sum }_{i=0}^{m-1}(A^2(i)-2A(i)B(x-m+i+1)+B^2(x-m+i+1))$

很套路的尝试反转$A$串，定义$A^{\prime }(i)=A(m-i-1)$

$D(x)={\sum }_{i=0}^{m-1}(A^{\prime 2}(m-i-1)-2A^{\prime }(m-i-1)B(x-m+i+1)+B^2(x-m+i+1))$

一般只想求和符号里面放一个东西，求和里面的加减乘除考虑各自拆开

$D(x)={\sum }_{i=0}^{m-1}{A}^{\prime }(m-i-1{)}^2+{\sum }_{i=0}^{m-1}B(x-m+i+1{)}^2-2{\sum }_{i=0}^{m-1}{A}^{\prime }(m-i-1)B(x-m+i+1)$

第一项，可以$O(m)$预处理，与$B$无瓜，$sovled$

第二项，可以$O(n)$预处理前缀和，$solved$

第三项，看到ta，应该是满心欢喜的，๑乛◡乛๑卷积小宝贝，改写成$2{\sum }_{i+j=x}{A}^{\prime }(i)\ast B(j)$

最后整理一下，设$T={\sum }_{i=0}^{m-1}{A}^{\prime }(i{)}^2,pre(x)={\sum }_{i=0}^{x}B(i{)}^2,g(x)={\sum }_{i+j=x}{A}^{\prime }(i)\ast B(j)$

进而有，$D(x)=T+pre(x)-pre(x-m)-2g(x)$

$FFT$优化算$g(x)$，最后$O(n)$算$D(x)$

通配符的单模式串匹配

考虑依葫芦画瓢，一样的步骤

Step1：定义匹配函数

设通配符的数值为$0$，$0$乘任何数为$0$，既然想把与通配符的匹配算出来的$C(x,y)$也变成$0$，不放把数值乘进去

$C(x,y)=(A(x)-B(y){)}^2\cdot A(x)\cdot B(y)$

Step2：完全匹配函数

$D(x)={\sum }_{i=0}^{m-1}(A(i)-B(x-m+i+1){)}^2\cdot A(i)\cdot B(x-m+i+1)$

Step3：暴力化式子

$D(x)={\sum }_{i=0}^{m-1}A(i{)}^{3}B(x-m+i+1)+{\sum }_{i=0}^{m-1}B(x-m+i+1{)}^{3}A(i)-2{\sum }_{i=0}^{m-1}A(i{)}^{2}B(x-m+i+1{)}^2$

定义$A^{\prime }(i)=A(m-i-1)$

$D(x)={\sum }_{i=0}^{m-1}{A}^{\prime }(m-i-1{)}^{3}B(x-m+i+1)+{\sum }_{i=0}^{m-1}B(x-m+i+1{)}^3{A}^{\prime }(m-i-1)-2{\sum }_{i=0}^{m-1}{A}^{\prime }(x-m-1{)}^{2}B(x-m+i+1{)}^2$

更妙的地方出现了！！发现所有小括号加起来都是$x$！！！卷积宝贝作法了ԅ(¯㉨¯ԅ)

$D(x)={\sum }_{i+j=x}{A}^{\prime }(i{)}^{3}B(j)+{\sum }_{i+j=x}B(j{)}^3{A}^{\prime }(i)-2{\sum }_{i+j=x}{A}^{\prime }(i{)}^{2}B(j{)}^2$

三次多项式乘法即可！！！（〜^㉨)〜(^▽)

常数稍稍大了点，但这不重要，吸氧能过的都不叫卡常p(#￣▽￣#)o

#include <cmath>
#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
#define maxn 1200100
struct complex {
	double x, i;
	complex(){}
	complex( double X, double I ) {
		x = X, i = I;
	}
}a[maxn], b[maxn], D[maxn];

const double pi = acos( -1.0 );

complex operator + ( complex a, complex b ) {
	return complex( a.x + b.x, a.i + b.i );
}

complex operator - ( complex a, complex b ) {
	return complex( a.x - b.x, a.i - b.i );
}

complex operator * ( complex a, complex b ) {
	return complex( a.x * b.x - a.i * b.i, a.x * b.i + a.i * b.x );
}

int len;
int r[maxn];

void FFT( complex *v, int opt ) {
	for( int i = 0;i < len;i ++ )
		if( i < r[i] ) swap( v[i], v[r[i]] );
	for( int i = 1;i < len;i <<= 1 ) {
		complex omega( cos( pi / i ), opt * sin( pi / i ) );
		for( int j = 0;j < len;j += ( i << 1 ) ) {
			complex w( 1, 0 );
			for( int k = 0;k < i;k ++, w = w * omega ) {
				complex x = v[j + k], y = v[j + k + i] * w;
				v[j + k] = x + y;
				v[j + k + i] = x - y;
			}
		}
	}
}

int m, n, cnt;
char s1[maxn], s2[maxn];
int ans[maxn], A[maxn], B[maxn];

int main() {
	scanf( "%d %d %s %s", &m, &n, s1, s2 );
	len = 1; int l = 0;
	while( len < m + n ) len <<= 1, l ++;
	for( int i = 0;i < len;i ++ )
		r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
	for( int i = 0;i < m;i ++ ) 
		A[i] = ( s1[i] == '*' ) ? 0 : s1[i] - 'a' + 1;
	for( int i = 0;i < n;i ++ )
		B[i] = ( s2[i] == '*' ) ? 0 : s2[i] - 'a' + 1;
	reverse( A, A + m );
	for( int i = 0;i < len;i ++ ) {
		a[i] = complex( A[i] * A[i] * A[i], 0 );
		b[i] = complex( B[i], 0 );
	}
	FFT( a, 1 ), FFT( b, 1 );
	for( int i = 0;i < len;i ++ )
		D[i] = D[i] + a[i] * b[i];
	for( int i = 0;i < len;i ++ ) {
		a[i] = complex( A[i], 0 );
		b[i] = complex( B[i] * B[i] * B[i], 0 );
	}
	FFT( a, 1 ), FFT( b, 1 );
	for( int i = 0;i < len;i ++ )
		D[i] = D[i] + a[i] * b[i];
	for( int i = 0;i < len;i ++ ) {
		a[i] = complex( A[i] * A[i], 0 );
		b[i] = complex( B[i] * B[i], 0 );
	}
	FFT( a, 1 ), FFT( b, 1 );
	for( int i = 0;i < len;i ++ )
		D[i] = D[i] - a[i] * b[i] * complex( 2, 0 );
	FFT( D, -1 );
	for( int i = m - 1;i < n;i ++ )
		if( fabs( D[i].x / len ) < 0.5 ) ans[++ cnt] = i - m + 2;
    //本来D[i]算的是i-m+1 但此题要求下标从1开始 所以要多加1
	printf( "%d\n", cnt );
	for( int i = 1;i <= cnt;i ++ )
	 	printf( "%d ", ans[i] );
	return 0;
}

---

K：P5748

集合划分计数

题意：求$n\le 1e5$的贝尔数，$bell(n)$即将$n$个有标号的球划分若干集合的方案数

设$f_{i,j}$表示把$i$个数划分为$j$个集合的方案数

那么就有两种转移，要么$i$自成一个集合，要么加入$j$个当中的任意一个集合

$f_{i,j}=f_{i-1,j-1}+f_{i-1,j}\times j$

$ans={\sum }_{i=1}^n{f}_{n,i}$

$O(n^2)$简单的暴力$DP$转移

设$f_i$表示把$n$个数划分为$i$个集合的方案数

(╯▽╰ ) 卡住了——(￣△￣；)

---

法一：生成函数

指数型生成函数用来解决排列问题，普通型生成函数用来解决组合问题

设一个非空集合的指数型生成函数$F(x)$，答案的指数型生成函数为$G(x)$

枚举划分的集合，集合间不区分，有$G(x)={\sum }_{i=0}^{\infty }\frac{F^i(x)}{i!}=e^{F(x)}$

一个非空集合的指数型生成函数，显然为$e^x-1$（${\sum }_{i=0}^{\infty }\frac{x^i}{i!}$减去一个空集$f_0=1$

多项式$exp$即可，题目求的是有标号的方案数，最后乘以$n!$即可

法二：分治$FFT$

枚举第一个元素所在集合大小，列出状态转移方程

$f_i={\sum }_{j=1}^i{C}_{i-1,j-1}\times f_{i-j}={\sum }_{j=1}^i\frac{(i-1)!}{(j-1)!(i-j)!}\times f_{i-j}$

把带$i,j$的项分开

$\frac{f_i}{(i-1)!}={\sum }_{j=1}^i\frac{1}{(j-1)!}\times \frac{f_{i-j}}{(i-j)!}$

令$F_i=\frac{1}{i!},G_i=\frac{f_i}{i!}$，卷积！！

$\frac{f_i}{(i-1)!}={\sum }_{j=0}^{i-1}{F}_j\ast G_{i-j-1}={\sum }_{j+k=i-1}{F}_j\ast G_k$

$F$是已知多项式，预处理即可，$G$是一个只与$f_j,j<i$有关的多项式，最后答案别忘记乘$(i-1)!$

这不是左边对右边造成贡献的好套路吗？？——中场休息的分治$FFT$登场！

#include <cstdio>
#include <iostream>
using namespace std;
#define int long long
#define mod 998244353
#define maxn 400100
#define n 100000
int r[maxn], fac[maxn], inv[maxn], f[maxn], g[maxn], A[maxn], B[maxn];

int qkpow( int x, int y ) {
	int ans = 1;
	while( y ) {
		if( y & 1 ) ans = ans * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return ans;
}

void NTT( int *v, int opt, int len ) {
	for( int i = 0;i < len;i ++ )
		if( i < r[i] ) swap( v[i], v[r[i]] );
	for( int i = 1;i < len;i <<= 1 ) {
		int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ) );
		for( int j = 0;j < len;j += ( i << 1 ) )
			for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {
				int x = v[j + k], y = v[j + k + i] * w % mod;
				v[j + k] = ( x + y ) % mod;
				v[j + k + i] = ( x - y + mod ) % mod;	
			}
	}
	if( opt == -1 ) {
		int Inv = qkpow( len, mod - 2 );
		for( int i = 0;i < len;i ++ )
			v[i] = v[i] * Inv % mod;
	}
}

void init() {
	fac[0] = inv[0] = 1;
	for( int i = 1;i <= n;i ++ ) fac[i] = fac[i - 1] * i % mod;
	inv[n] = qkpow( fac[n], mod - 2 );
	for( int i = n - 1;i;i -- ) inv[i] = inv[i + 1] * ( i + 1 ) % mod;
}

void CDQ( int L, int R ) {
	if( L == R ) { if( L ) f[L] = f[L] * fac[L - 1] % mod; return; }
	int mid = ( L + R ) >> 1;
	CDQ( L, mid );
	int len = 1, l = 0;
	while( len <= R - L + 1 ) len <<= 1, l ++;
	for( int i = 0;i < len;i ++ ) 
		r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
	for( int i = L;i <= mid;i ++ ) A[i - L] = f[i] * inv[i] % mod;
	for( int i = 0;i < len;i ++ ) B[i] = inv[i];
	NTT( A, 1, len );
	NTT( B, 1, len );
	for( int i = 0;i < len;i ++ ) A[i] = A[i] * B[i] % mod;
	NTT( A, -1, len );
	for( int i = mid + 1;i <= R;i ++ ) f[i] = ( f[i] + A[i - 1 - L] ) % mod;
	for( int i = 0;i < len;i ++ ) A[i] = B[i] = 0;
	CDQ( mid + 1, R );
}

signed main() {
	init();
	f[0] = 1;
	CDQ( 0, n );
	int T;
	scanf( "%lld", &T );
	while( T -- ) {
		int x;
		scanf( "%lld", &x );
		printf( "%lld\n", f[x] );
	}
	return 0;
}

---

L：P3702

[SDOI2017]序列计数

至少有一个为质数(;¬_¬)俺讨厌至少。。。

转化为满足要求的序列总数减去一个质数都没有的满足要求的序列数

一个质数都没有的序列数，不妨把$m$以内的所有质数全都丢掉

设$f_{i,j}$：前$i$个数的总和取模$p$为$j$的方案数（一个质数都没有，枚举第$i$个选的数

$f_{i,j}={\sum }_{(j_1+j_2)\%p=j}{f}_{i-1,j_2}$，$j_1$为合数

下面带取模，好熟悉，$E$题的序列统计⊙﹏⊙!，修改转移方程

$f_{i<<1,j}={\sum }_{(j_1+j_2)\%p=j}{f}_{i,j_1}\ast f_{i,j_2}$

突然想起，$p$好像没保证是质数。。。淦((٩(//̀Д/́/)۶))没有原根，指数上的取模就不能进行了

不！٩(๑`н´๑)۶我是sb，这比序列统计还要简单，下面的条件是加又不是乘，干嘛甩成指数做，多此一举！！

已经是个卷积形式了。。。难道三模数$NTT$？？害怕码力．<{=．．．．（嘎_嘎嘎~）

---

有矩阵快速幂+容斥的，容斥实在对于现在的我要求太高了

暴力卷积快速幂就可以跑过去(•̀⌄•́乁，这算不算搞出这道题了p(#￣▽￣#)o

#include <cstdio>
#define int long long
#define mod 20170408
#define maxn 20000005
#define maxp 400
int n, m, p, cnt;
int f[maxp], g[maxp], c[maxp], F[maxp], G[maxp], prime[maxn];
bool vis[maxn];

void mul( int *a, int *b ) {
	for( int i = 0;i < p;i ++ )
		for( int j = 0;j < p;j ++ )
			c[i + j] = ( c[i + j] + a[i] * b[j] % mod ) % mod;
	for( int i = 0;i < p;i ++ )
		a[i] = ( c[i] + c[i + p] ) % mod, c[i] = c[i + p] = 0;
}

signed main() {
	scanf( "%lld %lld %lld", &n, &m, &p );
	f[1] = g[1] = F[0] = G[0] = 1;
	for( int i = 2;i <= m;i ++ ) {
		f[i % p] ++;
		if( ! vis[i] ) prime[++ cnt] = i;
		else g[i % p] ++;
		for( int j = 1;j <= cnt && i * prime[j] <= m;j ++ ) {
			vis[i * prime[j]] = 1;
			if( i % prime[j] == 0 ) break;
		}
	}
	while( n ) {
		if( n & 1 ) mul( F, f ), mul( G, g );
		mul( f, f ), mul( g, g );
		n >>= 1;
	}
	printf( "%lld\n", ( F[0] - G[0] + mod ) % mod );
	return 0;
}

这题的难度放错位置了，就当中场休息好了

M：P5667

拉格朗日插值2

套拉格朗日差值公式$f_{m+x}={\sum }_{i=0}^{n}f(i){\prod }_{i\ne j}\frac{m+x-j}{i-j}$

${\prod }_{i\ne j}m+x-j\to (m+x)(m+x-1)+...(m+x-n)/(m+x-i)=\frac{(m+x)!}{(m+x-n-1)!(m+x-i)}$

对于枚举的$i$，有$[0,i)$小于$i$，相减为正，而对于$(i,n]$的共$n-i$个数而言，分母为负

所以乘积的符号便由$n-i$的奇偶决定，$(-1{)}^{n-i}$

${\prod }_{j<i}i-j=(i-0)(i-1)(i-2)...(2)(1)=i!$

${\prod }_{i<j}i-j=(n-i)(n-i+1)...(2)(1)=(n-i)!$（不考虑正负，上面已经判断了

$\Rightarrow f_{m+x}={\sum }_{i=0}^n(f(i)\times \frac{(m+x)!}{(-1{)}^{n-i}i!(n-i)!(m+x-n-1)!(m+x-i)})$

把只与$x$有关的式子提出来，$(-1{)}^{n-i}$只带来符号改变，不妨放在分子上（分母太肥了(;¬_¬)

$\Rightarrow f_{m+x}=\frac{(m+x)!}{(m+x-n-1)!}\times {\sum }_{i=0}^n(\frac{(-1{)}^{n-i}f(i)}{i!(n-i)!}\ast \frac{1}{m+x-i})$

发现是个卷积形式，但是$i$的取值不是$0\le i\le x$，而达到了$n$，与一般的有所区别(・。・)

但是没关系啊ԅ(¯㉨¯ԅ)，我们卷积最喜欢的不就是重新定义新函数改变原多项式顺序吗？只要最后能很快速地反映所求答案位置

老子把整个序列硬往后挪$n$位，第$x+n$位计算的答案存的就是原来$f_{m+x}$

令$A_i=\frac{(-1{)}^{n-i}f(i)}{i!(n-i)!},B_i=\frac{1}{m-n+i}$

$f_{x+m}=\frac{(x+m)!}{(m+x-n-1)!}{\sum }_{i=0}^n{A}_i\ast B_{n+x-i}=\frac{(x+m)!}{(m+x-n-1)!}{\sum }_{j+k=n+x}{A}_j\ast B_k$

令$i=x+n$（$x\in [0,n]\Rightarrow i\in [n,2n]$

$=\frac{(i+m-n)!}{(i+m-n-n-1)!}{\sum }_{j+k=i}{A}_j\ast B_k$

预处理阶乘$fa{c}_i=i!$，$i$逆元$in{v}_i$，$i!$阶乘逆元$invfa{c}_i$，差位的阶乘$MS\_fa{c}_i={\prod }_{j=0}^{i-1}(m-n+j)$，$m-n+i$差位逆元$MS\_in{v}_{i+1}$，$MS\_fa{c}_i$差位阶乘逆元$MS\_facin{v}_i$，则有$\frac{1}{(i+m-n)!}=MS\_fa{c}_{i+1},(i+m-n-1-n)!=MS\_facin{v}_{i-n}$

最后bb这么多，是为了解释代码的一些地方，避免不必要的时间投入看代码，FFT就是这种玩下标的东西，经常容易玩懵，找准哪个位置是自己原来求解数组的对应位置是关键

#include <cstdio>
#include <iostream>
using namespace std;
#define int long long
#define mod 998244353
#define maxn 1000000
int n, m, len;
int f[maxn], r[maxn], A[maxn], B[maxn];
int inv[maxn], facinv[maxn], fac[maxn];
int MS_fac[maxn], MS_inv[maxn], MS_facinv[maxn];

int qkpow( int x, int y ) {
	int ans = 1;
	while( y ) {
		if( y & 1 ) ans = ans * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return ans;
}

void NTT( int *v, int opt ) {
	for( int i = 0;i < len;i ++ )
		if( i < r[i] ) swap( v[i], v[r[i]] );
	for( int i = 1;i < len;i <<= 1 ) {
		int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ) );
		for( int j = 0;j < len;j += ( i << 1 ) )
			for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {
				int x = v[j + k], y = v[j + k + i] * w % mod;
				v[j + k] = ( x + y ) % mod;
				v[j + k + i] = ( x - y + mod ) % mod;
			}
	}
	if( opt == -1 ) {
		int Inv = qkpow( len, mod - 2 );
		for( int i = 0;i < len;i ++ )
			v[i] = v[i] * Inv % mod;
	}
}

void init( int n, int m ) {
	fac[0] = MS_fac[0] = MS_inv[0] = inv[0] = 1;
	for( int i = 1;i <= ( n << 1 | 1 );i ++ ) {
		fac[i] = fac[i - 1] * i % mod;
		MS_fac[i] = MS_fac[i - 1] * ( m - n + i - 1 ) % mod;
	}
	facinv[n << 1 | 1] = qkpow( fac[n << 1 | 1], mod - 2 );
	MS_facinv[n << 1 | 1] = qkpow( MS_fac[n << 1 | 1], mod - 2 );
	for( int i = ( n << 1 );~ i;i -- ) {
		facinv[i] = facinv[i + 1] * ( i + 1 ) % mod;
		MS_facinv[i] = MS_facinv[i + 1] * ( m - n + i ) % mod;
	}
	for( int i = ( n << 1 | 1 );i;i -- ) {
		inv[i] = facinv[i] * fac[i - 1] % mod;
		MS_inv[i] = MS_facinv[i] * MS_fac[i - 1] % mod;
	}
}

signed main() {
	scanf( "%lld %lld", &n, &m );
	init( n, m );
	for( int i = 0;i <= n;i ++ )
		scanf( "%lld", &f[i] );
	len = 1; int l = 0;
	while( len <= n * 3 ) len <<= 1, l ++;
	for( int i = 0;i < len;i ++ )
		r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
	for( int i = 0;i <= n;i ++ ) {
		A[i] = f[i] * facinv[i] % mod * facinv[n - i] % mod;
		if( ( n - i ) & 1 ) A[i] = mod - A[i];
	}
	for( int i = 0;i <= ( n << 1 );i ++ )
		B[i] = MS_inv[i + 1];
	NTT( A, 1 ), NTT( B, 1 );
	for( int i = 0;i < len;i ++ )
		A[i] = A[i] * B[i] % mod;
	NTT( A, -1 );
	for( int i = n;i <= ( n << 1 );i ++ )
		printf( "%lld ", MS_fac[i + 1] * MS_facinv[i - n] % mod * A[i] % mod );
	return 0;
}

N：P4245——拾得知识药丸一颗——三模数$NTT$

【模板】任意模数多项式乘法

别称：三模数$NTT$ 嘿嘿嘿——才发现这种三模数$NTT$就是$MTT$，那没事了V●ᴥ●V

当然这题有很多好的解法，但是鄙人认为还是应该掌握基础的不要脑子的sb基础做法

$\{\begin{array}{l}x\equiv x_1(mod{M}_1)\\ x\equiv x_2(mod{M}_2)\\ x\equiv x_3(mod{M}_3)\end{array}$

$x=x_1+k_1{M}_1=x_2+k_2{M}_2\Rightarrow x_1+k_1{M}_1\equiv x_2(\mathrm{m}\mathrm{o}\mathrm{d}{M}_2)\Rightarrow k_1=\frac{x_2-x_1}{M_1}(\mathrm{m}\mathrm{o}\mathrm{d}{M}_2)$

$x\equiv x_1+k_1{M}_1(\mathrm{m}\mathrm{o}\mathrm{d}{M}_1{M}_2)$，令$t=x_1+k_1{M}_1$

$t+k^{\prime }{M}_1{M}_2=x_3+k_3{M}_3\Rightarrow t+k^{\prime }{M}_1{M}_2\equiv x_3(\mathrm{m}\mathrm{o}\mathrm{d}{M}_3)\Rightarrow k^{\prime }=\frac{x_3-t}{M_1{M}_2}(\mathrm{m}\mathrm{o}\mathrm{d}{M}_3)$

$\because x<M_1{M}_2{M}_3\Rightarrow x=k^{\prime }{M}_1{M}_2+t$

中国剩余定理及扩展剩余定理证明

#include <cstdio>
#include <iostream>
using namespace std;
#define int long long
#define maxn 400100
int len, n, m, p;
int r[maxn], ans[maxn], B[maxn], F[maxn], G[maxn];
int Mod[5] = { 0, 998244353, 469762049, 1004535809 };

int qkpow( int x, int y, int mod ) {
	x %= mod;
	int ret = 1;
	while( y ) {
		if( y & 1 ) ret = ret * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return ret;
}

struct poly {
	int A[maxn], mod;
	void NTT( int *v, int opt ) {
		for( int i = 0;i < len;i ++ )
			if( i < r[i] ) swap( v[i], v[r[i]] );
		for( int i = 1;i < len;i <<= 1 ) {
			int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ), mod );
			for( int j = 0;j < len;j += ( i << 1 ) )
				for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {
					int x = v[j + k], y = v[j + k + i] * w % mod;
					v[j + k] = ( x + y ) % mod;
					v[j + k + i] = ( x - y + mod ) % mod;
				}
		}
		if( opt == -1 ) {
			int inv = qkpow( len, mod - 2, mod );
			for( int i = 0;i < len;i ++ )
				v[i] = v[i] * inv % mod;
		}
	}
}ntt[5];

int qkmul( int x, int y, int mod ) {
	x %= mod;
	int ret = 0;
	while( y ) {
		if( y & 1 ) ret = ( ret + x ) % mod;
		x = ( x + x ) % mod;
		y >>= 1;
	}
	return ret;
}

void CRT() {
	int mod_12 = Mod[1] * Mod[2];
	int inv1 = qkpow( Mod[1], Mod[2] - 2, Mod[2] );
	int inv2 = qkpow( Mod[2], Mod[1] - 2, Mod[1] );
	int inv_12 = qkpow( mod_12, Mod[3] - 2, Mod[3] );
	for( int i = 0;i < len;i ++ ) {
		int x1 = ntt[1].A[i];
		int x2 = ntt[2].A[i];
		int x3 = ntt[3].A[i];
		int t = ( qkmul( x1 * Mod[2], inv2, mod_12 ) + qkmul( x2 * Mod[1], inv1, mod_12 ) ) % mod_12;
		int k = ( x3 - t % Mod[3] + Mod[3] ) % Mod[3] * inv_12 % Mod[3];
		ans[i] = ( k % p * ( mod_12 % p ) % p + t % p ) % p;
	}
}

signed main() {
	scanf( "%lld %lld %lld", &n, &m, &p );
	for( int i = 0;i <= n;i ++ ) scanf( "%lld", &F[i] );
	for( int i = 0;i <= m;i ++ ) scanf( "%lld", &G[i] );
	len = 1;int l = 0;
	while( len <= n + m ) len <<= 1, l ++;
	for( int i = 0;i < len;i ++ )
		r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
	for( int k = 1;k <= 3;k ++ ) {
		ntt[k].mod = Mod[k];
		for( int i = 0;i <= n;i ++ ) ntt[k].A[i] = F[i];
		for( int i = 0;i <= m;i ++ ) B[i] = G[i];
		for( int i = m + 1;i < len;i ++ ) B[i] = 0;
		ntt[k].NTT( ntt[k].A, 1 );
		ntt[k].NTT( B, 1 );
		for( int i = 0;i < len;i ++ ) 
			ntt[k].A[i] = ntt[k].A[i] * B[i] % Mod[k];
		ntt[k].NTT( ntt[k].A, -1 );
	}
	CRT();
	for( int i = 0;i <= n + m;i ++ )
		printf( "%lld ", ans[i] );
	return 0;
}

---

O：P5488

差分与前缀和

一.前缀和

$s_i={\sum }_{j=0}^i{a}_j$，考虑前缀和的式子与卷积式子$s_i={\sum }_{j+k=i}{a}_j\ast b_k$

发现前缀和相当于是乘了一个系数全是$1$的多项式${\sum }_{i=0}^{\infty }1x_i$

由等比数列公式可得该多项式的闭形式$\frac{1}{1-x}$

$k$阶前缀和，就相当于卷了$k$次，卷积有结合律，所以最后相当于原来的数组卷积$\frac{1}{(1-x{)}^k}$

这个式子的第$n$项系数为$C_{n+k-1}^n=C_{n+k-1}^{k-1}$，递推即可

至于为什么是这个系数，可以打表验证，其实发现前缀和是跟杨辉三角有着某种联系的，意会也可

二.差分

差分相当于是乘了一个系数全是$-1$的多项式${\sum }_{i=0}^{\infty }-1x^i$，闭形式为$1-x$

最后相当于卷积$(1-x{)}^k={\sum }_{i=0}^{\infty }(-1{)}^i{C}_k^i{x}^i$

递推即可$C_k^i=C_k^{i-1}\times \frac{k-i+1}{i}$

#include <cstdio>
#include <iostream>
using namespace std;
#define mod 1004535809
#define int long long
#define maxn 400100
int len;
int r[maxn], inv[maxn], A[maxn], B[maxn];

int read() {
	int x = 0;char s = getchar();
	while( s < '0' || s > '9' ) s = getchar();
	while( '0' <= s && s <= '9' ) {
		x = ( ( x << 1 ) + ( x << 3 ) + ( s - 48 ) ) % mod;
		s = getchar();
	}
	return x;
}

int qkpow( int x, int y ) {
	int ans = 1;
	while( y ) {
		if( y & 1 ) ans = ans * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return ans;
}

void NTT( int *v, int opt ) {
	for( int i = 0;i < len;i ++ )
		if( i < r[i] ) swap( v[i], v[r[i]] );
	for( int i = 1;i < len;i <<= 1 ) {
		int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ) );
		for( int j = 0;j < len;j += ( i << 1 ) ) 
			for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {
				int x = v[j + k], y = v[j + k + i] * w % mod;
				v[j + k] = ( x + y ) % mod;
				v[j + k + i] = ( x - y + mod ) % mod;
			}
	}
	if( opt == -1 ) {
		int Inv = qkpow( len, mod - 2 );
		for( int i = 0;i < len;i ++ )
			v[i] = v[i] * Inv % mod;
	}
}

signed main() {
	int n, type, k;
	scanf( "%lld", &n );
	k = read();
	scanf( "%lld", &type );
	for( int i = 0;i < n;i ++ ) scanf( "%lld", &A[i] );
	inv[1] = 1;
	for( int i = 2;i < n;i ++ ) inv[i] = ( mod - mod / i ) * inv[mod % i] % mod;
	B[0] = 1;
	if( ! type ) {
		for( int i = 1;i < n;i ++ )
			B[i] = B[i - 1] * ( i + k - 1 ) % mod * inv[i] % mod;
	}
	else {
		for( int i = 1;i < n;i ++ )
			B[i] = ( B[i - 1] * ( k - i + 1 + mod ) % mod ) % mod * inv[i] % mod;
		for( int i = 1;i < n;i ++ )
			if( i & 1 ) B[i] = mod - B[i];
	}
	len = 1;int l = 0;
	while( len < ( n << 1 ) ) len <<= 1, l ++;
	for( int i = 0;i < len;i ++ )
		r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
	NTT( A, 1 ), NTT( B, 1 );
	for( int i = 0;i < len;i ++ ) A[i] = A[i] * B[i] % mod;
	NTT( A, -1 );
	for( int i = 0;i < n;i ++ )	
		printf( "%lld ", A[i] );
	return 0;
}

---

P：P4199

万径人踪灭

问：不连续的关于某根对称轴对称的回文子序列个数

非常难做，正面硬刚刚不动 (-ロ-)(-ロ-)

考虑转化问题：关于某根对称轴对称的回文子序列个数-回文子串个数

而对于一个序列求其回文子串个数，直接运用$manacher$即可求解

现在考虑怎么求关于某根对称轴对称的回文子序列个数

对称轴有两种情况，假设左右对称两个位置字符相等个数为$k$（不算对称轴

①以某个位置为轴

答案为$2^{k+1}-1$，那些位置以及对称轴选与不选再减去全都不选的情况

②以某两个位置的间隙为轴

答案为$2^k-1$

考虑每个位置，$f_i$：以$i$为对称轴的两个位置字符相等的个数

$f_i={\sum }_{j=0}^i[s_j,s_{i\times 2-j}]$

非常可以，这很卷积，套用$D$题$DNA$的套路，分字符自己卷自己，再把个数加起来

注意：其实卷积后求出的$f_i$表示以$\frac{i}{2}$为对称轴（除不尽就是以位置间隙为轴的意思

位置字符相等个数$\times 2$

所以真正的个数应该为$\frac{f_i}{2}$

但是如果$\frac{i}{2}$除得尽，代表一个位置的话，卷积中的$g_{\frac{i}{2}}\ast g_{\frac{i}{2}}$便只算了一次，

不妨用$\lfloor \frac{f_i}{2}\rfloor +[2|i]$来特别加上即可，统计答案就是${\sum }_{i=0}^{2n-2}{2}^{\lfloor \frac{f_i}{2}\rfloor +[2|i]}-1$

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define maxn 400100
#define int long long
#define MOD 1000000007
#define mod 998244353
int len;
int rev[maxn], A[maxn], B[maxn], f[maxn];

int qkpow( int x, int y, int Mod ) {
	int ans = 1;
	while( y ) {
		if( y & 1 ) ans = ans * x % Mod;
		x = x * x % Mod;
		y >>= 1;
	}
	return ans;
}

void NTT( int *v, int opt ) {
	for( int i = 0;i < len;i ++ )
		if( i < rev[i] ) swap( v[i], v[rev[i]] );
	for( int i = 1;i < len;i <<= 1 ) {
		int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ), mod );
		for( int j = 0;j < len;j += ( i << 1 ) )
			for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {
				int x = v[j + k], y = v[j + k + i] * w % mod;
				v[j + k] = ( x + y ) % mod;
				v[j + k + i] = ( x - y + mod ) % mod;
			}
	}
	if( opt == -1 ) {
		int inv = qkpow( len, mod - 2, mod );
		for( int i = 0;i < len;i ++ )
			v[i] = v[i] * inv % mod;
	}
}

int r[maxn], p[maxn];
char s[maxn];

int manacher() {
	int n = ( strlen( s ) << 1 | 1 ), idx = 0;
	for( int i = 0;i < n;i ++ )
		p[i] = ( i & 1 ) ? s[idx ++] : '#';
	int R = -1, C = -1;
	for( int i = 0;i < n;i ++ ) {
		r[i] = ( R > i ) ? min( R - i, r[C * 2 - i] ) : 1;
		while( i + r[i] < n && ~ ( i - r[i] ) )
			if( p[i + r[i]] == p[i - r[i]] ) r[i] ++;
			else break;
		if( i + r[i] > R ) R = i + r[i], C = i;
	}
	int ans = 0;
	for( int i = 0;i < n;i ++ )
		ans += ( r[i] >> 1 );
	return ans % MOD;
}

signed main() {
	scanf( "%s", s );
	int n = strlen( s );
	len = 1;int l = 0;
	while( len < ( n << 1 ) ) len <<= 1, l ++;
	for( int i = 0;i < len;i ++ )
		rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
	for( int i = 0;i < n;i ++ )
		if( s[i] == 'a' ) A[i] = 1;
	for( int i = 0;i < n;i ++ )
		if( s[i] == 'b' ) B[i] = 1;
	NTT( A, 1 ), NTT( B, 1 );
	for( int i = 0;i < len;i ++ )
		f[i] = ( A[i] * A[i] % mod + B[i] * B[i] % mod ) % mod;
	NTT( f, -1 );
	int ans = 0;
	for( int i = 0, k = 1;i <= n * 2 - 2;i ++, k ^= 1 )
		ans = ( ans + qkpow( 2, ( f[i] >> 1 ) + k, MOD ) - 1 + MOD ) % MOD;
	printf( "%lld\n", ( ans - manacher() + MOD ) % MOD );
	return 0;
}

---

Q：P4061

[HEOI2016/TJOI2016]求和

$f_n={\sum }_{i=0}^n{\sum }_{j=0}^{i}S(i,j)\times 2^j\times j!$

$\because S(i,j)=0[i<j]$

$\therefore f_n={\sum }_{i=0}^n{\sum }_{j=0}^{n}S(i,j)\times 2^j\times j!$

$S(n,m)=\frac{1}{m!}{\sum }_{i=0}^m(-1{)}^i{C}_m^i(m-i{)}^n={\sum }_{i=0}^n\frac{(-1{)}^i}{i!}\frac{(m-i{)}^n}{(m-i)!}$

$f_n={\sum }_{j=0}^n{2}^j\times j!{\sum }_{i=0}^n{\sum }_{k=0}^j\frac{(-1{)}^k}{k!}\cdot \frac{(j-k{)}^i}{(j-k)!}={\sum }_{j=0}^n{2}^j\times j!{\sum }_{k=0}^j\frac{(-1{)}^k}{k!}\cdot \frac{{\sum }_{i=0}^n(j-k{)}^i}{(j-k)!}$

（换一下变量名更好看V●ᴥ●V

$f_n={\sum }_{i=0}^n{2}^i\times i!{\sum }_{j=0}^i\frac{(-1{)}^j}{j!}\cdot \frac{{\sum }_{k=0}^n(i-j{)}^k}{(i-j)!}$

令$g_i=\frac{(-1{)}^i}{i!},h_i=\frac{{\sum }_{k=0}^n{i}^k}{i!}=\frac{i^{n+1}-1}{(i-1)i!}$（等比数列公式不可能还有人呢不会吧嘿嘿嘿qwq

$f_n={\sum }_{i=0}^n{2}^i\times i!{\sum }_{j=0}^n{g}_j\ast h_{i-j}$（特别的，$g_0=0,g_1=n+1$

$NTT$冲鸭！！

#include <cstdio>
#include <iostream>
using namespace std;
#define int long long
#define mod 998244353
#define maxn 400100
int n, m, len;
int r[maxn], f[maxn], g[maxn];
int fac[maxn], inv[maxn], facinv[maxn];

int qkpow( int x, int y ) {
	int ans = 1;
	while( y ) {
		if( y & 1 ) ans = ans * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return ans;
}

void NTT( int *v, int opt ) {
	for( int i = 0;i < len;i ++ )
		if( i < r[i] ) swap( v[i], v[r[i]] );
	for( int i = 1;i < len;i <<= 1 ) {
		int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ) );
		for( int j = 0;j < len;j += ( i << 1 ) )
			for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {
				int x = v[j + k], y = v[j + k + i] * w % mod;
				v[j + k] = ( x + y ) % mod;
				v[j + k + i] = ( x - y + mod ) % mod;
			}
	} 
	if( opt == -1 ) {
		int Inv = qkpow( len, mod - 2 );
		for( int i = 0;i < len;i ++ )
			v[i] = v[i] * Inv % mod;
	}
}

void init( int n ) {
	fac[0] = inv[0] = inv[1] = facinv[0] = 1;
	for( int i = 1;i <= n;i ++ )
		fac[i] = fac[i - 1] * i % mod;
	facinv[n] = qkpow( fac[n], mod - 2 );
	for( int i = n - 1;i;i -- )
		facinv[i] = facinv[i + 1] * ( i + 1 ) % mod;
	for( int i = 2;i <= n;i ++ )
		inv[i] = ( mod - mod / i ) * inv[mod % i] % mod;
}

signed main() {
	scanf( "%lld", &n );
	init( n );
	for( int i = 0;i <= n;i ++ ) {
		g[i] = facinv[i];
		if( i & 1 ) g[i] = mod - g[i];
	}
	f[0] = 1, f[1] = n + 1;
	for( int i = 2;i <= n;i ++ )
		f[i] = ( qkpow( i, n + 1 ) - 1 ) * inv[i - 1] % mod * facinv[i] % mod;
	len = 1; int l = 0;
	while( len <= ( n << 1 ) ) len <<= 1, l ++;
	for( int i = 0;i < len;i ++ )
		r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
	NTT( g, 1 ), NTT( f, 1 ); 
		for( int i = 0;i < len;i ++ )
			f[i] = f[i] * g[i] % mod;
	NTT( f, -1 );
	int ans = 0;
	for( int i = 0;i <= n;i ++ )
		ans = ( ans + qkpow( 2, i ) * fac[i] % mod * f[i] % mod ) % mod;
	printf( "%lld\n", ans );
	return 0;
}

---

缘遇逆天秘籍——金牌打野

搭配生成函数修炼事半功倍哦♪(´ε｀)

R：P4841

[集训队作业2013]城市规划

求$n$个点的简单有标号无向连通图数目

法一：分治$FFT$

设$f_i$表示$i$个点的有标号无向连通图数，考虑列状态转移方程

呃—— (-ロ-) 好像转移不动

再设$g_i$表示$i$个点的图的数量，有$\frac{i\times (i-1)}{2}$条边，每条可选可不选，很容易就得出$g_i=2^{\frac{i\times (i-1)}{2}}$

$g_n={\sum }_{i=1}^n{C}_{n-1}^{i-1}{f}_i{g}_{n-i}$

枚举$1$号点所在连通图的大小，联通的方案数即为$f_i$，有标号所以从剩下的点中选$i-1$个，乘$C_{n-1}^{i-1}$，其余点分成什么样的 图都可以，方案数为$g_{n-i}$

$g_n={\sum }_{i=1}^n{C}_{n-1}^{i-1}{f}_i{g}_{n-i}=C_{n-1}^{n-1}{f}_n{g}_0+{\sum }_{i=1}^{n-1}{C}_{n-1}^{i-1}{f}_i{g}_{n-i}=f_n+{\sum }_{i=1}^{n-1}\frac{(n-1)!}{(i-1)!(n-i)!}{f}_i{g}_{n-i}$

$\Rightarrow \frac{f_n}{(n-1)!}=\frac{g_n}{(n-1)!}-{\sum }_{i=1}^{n-1}\frac{f_i}{(i-1)!}\frac{g_{n-i}}{(n-i)!}$

令$F_i=\frac{f_i}{(i-1)!},G_i=\frac{g_i}{i!}$，则$F_n=G_n\times n-{\sum }_{i=1}^{n-1}{F}_i\ast G_{n-i}$

分治$FFT$卷积宝贝又出现了mua!(*╯3╰)

#include <cstdio>
#include <iostream>
using namespace std;
#define mod 1004535809
#define int long long
#define maxn 520100
int len;
int r[maxn], fac[maxn], inv[maxn], f[maxn], g[maxn], F[maxn], G[maxn];

int qkpow( int x, int y ) {
	int ans = 1;
	while( y ) {
		if( y & 1 ) ans = ans * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return ans;
}

void NTT( int *v, int opt, int len ) {
	for( int i = 0;i < len;i ++ )
		if( i < r[i] ) swap( v[i], v[r[i]] );
	for( int i = 1;i < len;i <<= 1 ) {
		int omega = qkpow( opt == 1 ? 3 : mod / 3 + 1, ( mod - 1 ) / ( i << 1 ) );
		for( int j = 0;j < len;j += ( i << 1 ) )
			for( int k = 0, w = 1;k < i;k ++, w = w * omega % mod ) {
				int x = v[j + k], y = v[j + k + i] * w % mod;
				v[j + k] = ( x + y ) % mod;
				v[j + k + i] = ( x - y + mod ) % mod;
			}
	}
	if( opt == -1 ) {
		int Inv = qkpow( len, mod - 2 );
		for( int i = 0;i < len;i ++ )
			v[i] = v[i] * Inv % mod;
	}
}

void solve( int L, int R ) {
	if( L == R ) return;
	int mid = ( L + R ) >> 1;
	solve( L, mid );
	int len = 1, l = 0;
	while( len <= R - L + 1 ) len <<= 1, l ++;
	for( int i = 0;i < len;i ++ )
		r[i] = ( r[i >> 1] >> 1 ) | ( ( i & 1 ) << ( l - 1 ) );
	for( int i = L;i <= mid;i ++ ) F[i - L] = f[i];
	for( int i = 1;i <= R - L;i ++ ) G[i - 1] = g[i];
	NTT( F, 1, len ), NTT( G, 1, len );
	for( int i = 0;i < len;i ++ ) F[i] = F[i] * G[i] % mod;
	NTT( F, -1, len );
	for( int i = mid + 1;i <= R;i ++ ) f[i] = ( f[i] - F[i - L - 1] + mod ) % mod;
	for( int i = 0;i < len;i ++ ) F[i] = G[i] = 0;
	solve( mid + 1, R );
}

void init( int n ) {
	fac[0] = inv[0] = 1;
	for( int i = 1;i <= n;i ++ ) fac[i] = fac[i - 1] * i % mod;
	inv[n] = qkpow( fac[n], mod - 2 );
	for( int i = n - 1;i;i -- ) inv[i] = inv[i + 1] * ( i + 1 ) % mod;
	for( int i = 1;i <= n;i ++ )
		g[i] = qkpow( 2, i * ( i - 1 ) / 2 ) * inv[i] % mod, f[i] = g[i] * i % mod;
}

signed main() {
	int n;
	scanf( "%lld", &n );
	init( maxn - 1 );
	solve( 0, n );
	printf( "%lld\n", f[n] * fac[n - 1] % mod );
	return 0;
}

法二：生成函数+多项式求逆

$g_n={\sum }_{i=1}^n{C}_{n-1}^{i-1}{f}_i{g}_{n-i}$

换一种推式子的方式

$g_n={\sum }_{i=1}^n\frac{(n-1)!}{(i-1)!(n-i)!}{f}_i{g}_{n-i}\Rightarrow \frac{g_n}{(n-1)!}={\sum }_{i=1}^n\frac{f_i}{(i-1)!}\frac{g_{n-i}}{(n-i)!}$

令$G^{\prime },F,G$分别为它们的生成函数，则$G^{\prime }(x)={\sum }_{i=0}^{\infty }{g}_{i+1}\frac{x^i}{i!},F(x)={\sum }_{i=0}^{\infty }{f}_{i+1}\frac{x^i}{i!},G(x)={\sum }_{i=0}^{\infty }{g}_i\frac{x^i}{i!}$

$G^{\prime }=F\ast G\Rightarrow F=\frac{G^{\prime }}{G}$

多项式求逆即可

#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
#define int long long
#define mod 1004535809
#define maxn 600000
int n, len, l;
int r[maxn], a[maxn], b[maxn], c[maxn], d[maxn], fac[maxn];

int qkpow( int x, int y ) {
    int ans = 1;
    while( y ) {
        if( y & 1 ) ans = ans * x % mod;
        x = x * x % mod;
        y >>= 1;
    }
    return ans;
}

void NTT( int *c, int f ) {
    for( int i = 0;i < len;i ++ ) 
        if( i < r[i] ) swap( c[i], c[r[i]] );
    for( int i = 1;i < len;i <<= 1 ) {
        int omega = qkpow( f == 1 ? 3 : 334845270, ( mod - 1 ) / ( i << 1 ) );
        for( int j = 0;j < len;j += ( i << 1 ) ) {
            int w = 1;
            for( int k = 0;k < i;k ++, w = w * omega % mod ) {
                int x = c[j + k], y = c[j + k + i] * w % mod;
                c[j + k] = ( x + y ) % mod;
                c[j + k + i] = ( x - y + mod ) % mod;
            }
        }
    }
    if( f == -1 ) {
        int inv = qkpow( len, mod - 2 );
        for( int i = 0;i < len;i ++ ) c[i] = c[i] * inv % mod;
    }
}

void polyinv( int deg, int *b ) {
    if( deg == 1 ) { b[0] = 1; return; }
    polyinv( deg + 1 >> 1, b );
    for( len = 1, l = 0;len < ( deg << 1 ) - 1;len <<= 1, l ++ );
    for( int i = 0;i < len;i ++ ) r[i] = r[i >> 1] >> 1 | ( ( i & 1 ) << l - 1 );
    for( int i = 0;i < deg;i ++ ) c[i] = a[i];
    for( int i = deg;i < n;i ++ ) c[i] = 0;
    NTT( c, 1 ), NTT( b, 1 );
    for( int i = 0;i < len;i ++ ) b[i] = ( 2 - b[i] * c[i] % mod + mod ) % mod * b[i] % mod;
    NTT( b, -1 );
    for( int i = deg;i < n;i ++ ) b[i] = 0;
}

signed main() {
    scanf( "%lld", &n ); fac[0] = 1, n ++;
    for( int i = 1;i <= n;i ++ ) fac[i] = fac[i - 1] * i % mod;
    a[0] = 1;
    for( int i = 1;i < n;i ++ ) {
        a[i] = qkpow( 2, i * ( i - 1 ) / 2 ) * qkpow( fac[i], mod - 2 ) % mod;
        d[i] = qkpow( 2, i * ( i - 1 ) / 2 ) * qkpow( fac[i - 1], mod - 2 ) % mod;
    }
    polyinv( n, b );
    for( len = 1, l = 0;len < ( n << 1 );len <<= 1, l ++ );
    for( int i = 0;i < len;i ++ ) r[i] = r[i >> 1] >> 1 | ( ( i & 1 ) << l - 1 );
    NTT( d, 1 ), NTT( b, 1 );
    for( int i = 0;i < len;i ++ ) d[i] = d[i] * b[i] % mod;
    NTT( d, -1 );
    printf( "%lld\n", d[n - 1] * fac[n - 2] % mod );
    return 0;
}

法三：生成函数+多项式$ln$

设$f_i$表示$i$个点的无标号的无向连通图数量，$g_i$表示无标号的图数量，令$F(x),G(x)$为别为其的$EGF$，枚举连通块个数

$G(x)={\sum }_{i=0}^{\infty }\frac{F(x{)}^i}{i!}=e^{F(x)}\Rightarrow F(x)=lnG(x)$

最后$[x^n]$系数即为答案，因为题目是有标号的连通图数，所以要乘$n!$，同理$g_i$是无标号的数，$G(x)$的系数应为$\frac{2^{\frac{i\times (i-1)}{2}}}{i!}$

---

S：[UVA4671]K-neighbor substrings

读完题目发现这是 D 的弱化版。同样的思路，枚举字符设为 $1$ 其余为 $0$，卷积求出该字符匹配数量。

所有字符匹配数量求和，然后判断是否 $\ge m-K$ 即可。

唯一的不同就是，这个是求 本质不同的子串 数量，首先想到的就是哈希。

但是很离谱，单模双模哈希怎么样都会挂，unsigned long long 自然溢出就可以通过。🙄

#include <map> 
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define ull unsigned long long
#define maxn 800005
map < ull, bool > mp;
int len, l;
int r[maxn], ans[maxn];
ull Hash[maxn], mi[maxn];
char a[maxn], b[maxn];

struct complex {
    double x, i;
    complex(){}
    complex( double X, double I ) { x = X, i = I; }
}f[maxn], g[maxn], h[maxn];
double pi = acos( -1.0 );
complex operator + ( complex u, complex v ) {
    return complex( u.x + v.x, u.i + v.i );
}
complex operator - ( complex u, complex v ) {
    return complex( u.x - v.x, u.i - v.i );
}
complex operator * ( complex u, complex v ) {
    return complex( u.x * v.x - u.i * v.i, u.i * v.x + u.x * v.i );
}

void FFT( complex *c, int opt ) {
    for( int i = 0;i < len;i ++ ) if( i < r[i] ) swap( c[i], c[r[i]] );
    for( int i = 1;i < len;i <<= 1 ) {
        complex omega( cos( pi / i ), opt * sin( pi / i ) );
        for( int j = 0;j < len;j += ( i << 1 ) ) {
            complex w( 1, 0 );
            for( int k = 0;k < i;k ++, w = w * omega ) {
                complex x = c[j + k], y = c[j + k + i] * w;
                c[j + k] = x + y, c[j + k + i] = x - y;
            }
        }
    }
}

bool check( int l, int r ) {
    ull num = Hash[r] - Hash[l - 1] * mi[r - l + 1];
    if( mp[num] ) return 0;
    else { mp[num] = 1; return 1; }
}

int  main() {
    mi[0] = 1;
    for( int i = 1;i < maxn;i ++ ) mi[i] = mi[i - 1] * 131;
    int K, Case = 0;
    while( scanf( "%d", &K ) and ~ K ) {
        scanf( "%s %s", a, b );
        int n = strlen( a ), m = strlen( b );
        for( len = 1, l = 0;len < ( n + m );len <<= 1, l ++ ) 
            len <<= 1, l ++;
        for( int i = 0;i < len;i ++ ) 
            r[i] = r[i >> 1] >> 1 | ( i & 1 ) << l - 1; 
        for( char ch = 'a';ch <= 'b';ch ++ ) {
            for( int i = 0;i < len;i ++ ) {
                f[i] = complex( ch == a[i], 0 );
                h[i] = complex( ch == b[i], 0 );
                g[i] = complex( 0, 0 );
            }
            for( int i = 0;i < m;i ++ ) g[i] = h[m - i - 1];
            FFT( f, 1 ), FFT( g, 1 );
            for( int i = 0;i < len;i ++ ) f[i] = f[i] * g[i];
            FFT( f, -1 );
            for( int i = m - 1;i < len;i ++ ) 
                ans[i - m + 1] += (int)( f[i].x / len + 0.5 );
        }
        Hash[0] = a[0] - 'a';
        for( int i = 1;i < n;i ++ ) 
            Hash[i] = Hash[i - 1] * 131 + (a[i] - 'a');
        int ret = 0;
        for( int i = 0;i <= n - m;i ++ ) 
            if( m - ans[i] <= K and check( i, i + m - 1 ) ) 
                ret ++;
        printf( "Case %d: %d\n", ++ Case, ret );
        for( int i = 0;i < len;i ++ ) ans[i] = 0;
        mp.clear();
    }
    return 0;
}

---

T：SP8372 TSUM - Triple Sums

洛谷链接

$A_i+A_j+A_k=S$ 已经属于是条件反射了，将值当成下标，直接卷积就行 $an{s}_{A_i+A_j+A_k=S}$ 。

但是这里有两个问题，一个是忽略了 $i<j<k$ 的顺序要求，在最后答案 $/6$ 即可。另一个是可能算重，因为 $A$ 并没有互不相同。

算重就容斥好了。$B_{2A_i},C_{3A_i}$。

$B,C$ 同样考虑顺序问题，所以最后的 $an{s}_i=a_i^3-3a_i{b}_i+2c_i$。再 $\text{IDFT}$ 回去。

注意 $A$ 有负的可能，直接整体平移 $20000$ 即可，这些都是小细节。

#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
#define int long long
#define maxn 500000
int len, l;
int r[maxn];

struct complex {
    double x, i;
    complex(){}
    complex( double X, double I ) { x = X; i = I; }
}a[maxn], b[maxn], c[maxn], ans[maxn];
double pi = acos( -1.0 );
complex operator + ( complex u, complex v ) {
    return complex( u.x + v.x, u.i + v.i );
}
complex operator - ( complex u, complex v ) {
    return complex( u.x - v.x, u.i - v.i );
}
complex operator * ( complex u, complex v ) {
    return complex( u.x * v.x - u.i * v.i, u.x * v.i + u.i * v.x );
}

void FFT( complex *c, int opt ) {
    for( int i = 0;i < len;i ++ ) if( i < r[i] ) swap( c[i], c[r[i]] );
    for( int i = 1;i < len;i <<= 1 ) {
        complex omega( cos( pi / i ), opt * sin( pi / i ) );
        for( int j = 0;j < len;j += ( i << 1 ) ) {
            complex w( 1, 0 );
            for( int k = 0;k < i;k ++, w = w * omega ) {
                complex x = c[j + k], y = w * c[j + k + i];
                c[j + k] = x + y, c[j + k + i] = x - y;
            }
        }
    }
}

signed main() {
    int n; scanf( "%lld", &n );
    for( int i = 1, x;i <= n;i ++ ) {
        scanf( "%lld", &x ); 
        x += 20000;
        a[x].x ++;
        b[x * 2].x ++;
        c[x * 3].x ++;
    }
    for( len = 1, l = 0;len <= 150000;len <<= 1, l ++ );
    for( int i = 0;i < len;i ++ ) r[i] = r[i >> 1] >> 1 | (i & 1) << l - 1;
    FFT( a, 1 ), FFT( b, 1 ), FFT( c, 1 );
    for( int i = 0;i < len;i ++ ) 
        ans[i] = a[i] * a[i] * a[i] - complex( 3, 0 ) * a[i] * b[i] + complex( 2, 0 ) * c[i];
    FFT( ans, -1 );
    for( int i = 0;i < len;i ++ ) {
        int x = (int)( ans[i].x / len + 0.5 );
        x /= 6;
        if( x > 0 ) printf( "%lld : %lld\n", i - 60000, x );
    }
    return 0;
}

---

U：[CodeChef - COUNTARI]Arithmetic Progressions

分块和 $FFT$。

选择一个适当的块大小 $block$，然后枚举每个块。

将合法的三元组 $(A_i,A_j,A_k)$ 分为两种情况考虑。

• 在块内。
  枚举块内的两个点，一个做 $A_i$，一个做 $A_j/A_k$。
  如果第三个是后面块的，就需要知道后面未操作块中是该值的个数。
  如果第三个是前面块或者本块的，因为本块内元素的枚举有序所以不会算重，就需要知道前面已操作块和本块已枚举 $A_i$ 的数是该值的个数。
  分别用 $suf,pre,cnt$ 来记录后面未操作块，前面已操作块和本块已操作数的信息。
• 在块外。
  $A_j$ 在本块内，$A_i$ 在已操作块中，$A_k$ 在未操作块中。
  这就是很常见的卷积了。这也是为什么上一种情况需要单独记录本块已操作数量，而不是直接与 $pre$ 合并。

最后每个块统计答案的时候，顺便将 $cnt\to pre$ 中。

选择好合适的块大小，假设值最大值为 $m$，时间复杂度为 $O(\frac{n^2}{block}+block\cdot m\log m)$。

#include <cmath>
#include <cstdio>
#include <iostream> 
using namespace std;
#define maxn 200000
int len, l, n, block, Max;
long long ans;
int r[maxn], a[maxn], pre[maxn], cnt[maxn], suf[maxn];

struct complex {
    double x, i;
    complex(){}
    complex( double X, double I ) { x = X, i = I; }
}f[maxn], g[maxn];
double pi = acos( -1.0 );
complex operator + ( complex u, complex v ) {
    return complex( u.x + v.x, u.i + v.i );
}
complex operator - ( complex u, complex v ) {
    return complex( u.x - v.x, u.i - v.i );
}
complex operator * ( complex u, complex v ) {
    return complex( u.x * v.x - u.i * v.i, u.x * v.i + u.i * v.x );
}

void FFT( complex *c, int opt ) {
    for( int i = 0;i < len;i ++ ) if( i < r[i] ) swap( c[i], c[r[i]] );
    for( int i = 1;i < len;i <<= 1 ) {
        complex omega( cos( pi / i ), opt * sin( pi / i ) );
        for( int j = 0;j < len;j += (i << 1) ) {
            complex w( 1, 0 );
            for( int k = 0;k < i;k ++, w = w * omega ) {
                complex x = c[j + k], y = c[j + k + i] * w;
                c[j + k] = x + y, c[j + k + i] = x - y;
            }
        }
    }
    if( opt == -1 ) for( int i = 0;i < len;i ++ ) c[i].x /= len;
}

int main() {
    scanf( "%d", &n );
    block = sqrt( n );
    for( int i = 1;i <= n;i ++ ) {
        scanf( "%d", &a[i] );
        Max = max( Max, a[i] );
        suf[a[i]] ++;
    }
    for( len = 1;len <= (Max << 1);len <<= 1, l ++ );
    for( int i = 0;i < len;i ++ ) r[i] = r[i >> 1] >> 1 | (i & 1) << l - 1;
    for( int l = 1;l <= n;l += block ) {
        int r = min( l + block - 1, n ), k;
        for( int i = l;i <= r;i ++ ) suf[a[i]] --;
        for( int i = l;i <= r;i ++ ) {
            for( int j = i + 1;j <= r;j ++ ) {
                k = (a[i] << 1) - a[j];
                if( 1 <= k and k <= Max ) ans += cnt[k] + pre[k];
                k = (a[j] << 1) - a[i];
                if( 1 <= k and k <= Max ) ans += suf[k];
            }
            cnt[a[i]] ++;
        }
        for( int i = 0;i <= Max;i ++ ) {
            f[i] = complex( pre[i], 0 );
            g[i] = complex( suf[i], 0 );
        }
        for( int i = Max + 1;i < len;i ++) {
        	f[i] = complex( 0, 0 );
        	g[i] = complex( 0, 0 );
        }
        FFT( f, 1 ), FFT( g, 1 );
        for( int i = 0;i < len;i ++ ) f[i] = f[i] * g[i];	
        FFT( f, -1 );
        for( int i = l;i <= r;i ++ ) {
            ans += (long long)( f[a[i] << 1].x + 0.5 );
            pre[a[i]] ++, cnt[a[i]] --;
        }
    }
    printf( "%lld\n", ans );
    return 0;
}

---

V：CodeForces528D Fuzzy Search

洛谷链接

做了很多个这种 $\text{{A,T,G,C}}$ 的字符 $\text{FFT/NTT}$ 的题后，我发现 这种多个字符之间的适配问题，一般都是拆开单个单个字符做的。

考虑最基础的适配是要求完全相同，即 $k=0$ 的严苛条件，我们套路地采取的是翻转一个串，$g_i^{\prime }=g_{m-i-1}$，使得卷积下标和为定值，分字符卷积后求和，判断是否和数组长度相同。

这是解决字符适配问题的通解基础，所有题目都只是在这个的基础上进行系列的改动罢了。

对于一个字符 $ch$。

设 $f_i=1/0$ 表示 $S_i$ 偏差 $k$ 以内是否有字符 $ch$ 出现，即 $f_i=[{\exists }_{j\in [i-k+1,i+k-1]}{S}_j=ch]$。

设 $g_i=1/0$ 表示 $T_i$ 是否为字符 $ch$，即 $g_i=[T_i=ch]$。

同样地，翻转 $g$，$g_i^{\prime }=g_{m-i-1}$，然后和 $f_i$ 卷积后为 $h_i$。

则，如果 $h_i\ne 0$ 也就意味着有 $h_i$ 个 $T_a=ch$ 在 $k$ 偏差内有对应的 $S_b=ch$。

注意：通过题意可以知道有可能 $T_{a_1},T_{a_2}$ 对应的是同一个 $b$。

多个字符累和 $h_i$ 后，对于 $h_i=m$ 的 $i$ 位置就意味着整个 $T$ 都适配成功。

#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
#define maxn 1000000
int r[maxn];
int len, l;

struct complex {
    double x, i;
    complex(){}
    complex( double X, double I ) { x = X, i = I; }
}f[maxn], g[maxn];
double pi = acos( -1.0 );
complex operator + ( complex u, complex v ) {
    return complex( u.x + v.x, u.i + v.i );
}
complex operator - ( complex u, complex v ) {
    return complex( u.x - v.x, u.i - v.i );
}
complex operator * ( complex u, complex v ) {
    return complex( u.x * v.x - u.i * v.i, u.x * v.i + u.i * v.x );
}

void FFT( complex *c, int opt ) {
    for( int i = 0;i < len;i ++ ) if( i < r[i] ) swap( c[i], c[r[i]] );
    for( int i = 1;i < len;i <<= 1 ) {
        complex omega( cos( pi / i ), sin( pi / i ) * opt );
        for( int j = 0;j < len;j += ( i << 1 ) ) {
            complex w( 1, 0 );
            for( int k = 0;k < i;k ++, w = w * omega ) {
                complex x = c[j + k], y = c[j + k + i] * w;
                c[j + k] = x + y, c[j + k + i] = x - y;
            }
        }
    }
    if( opt == -1 ) for( int i = 0;i < len;i ++ ) c[i].x /= len;
}

#define inf 0x3f3f3f3f
int n, m, k;
int cnt[maxn];
char s[maxn], t[maxn];

void work( char ch ) {
    for( int i = 0;i < len;i ++ ) f[i] = g[i] = complex( 0, 0 );
    int lst = -inf;
    for( int i = 0;i < n;i ++ ) {
        if( s[i] == ch ) lst = i;
        if( i - lst <= k ) f[i] = complex( 1, 0 );
    }
    lst = inf;
    for( int i = n - 1;~ i;i -- ) {
        if( s[i] == ch ) lst = i;
        if( lst - i <= k ) f[i] = complex( 1, 0 );
    }
    for( int i = 0;i < m;i ++ )
        if( t[i] == ch ) g[m - i - 1] = complex( 1, 0 );
    FFT( f, 1 ), FFT( g, 1 );
    for( int i = 0;i < len;i ++ ) f[i] = f[i] * g[i];
    FFT( f, -1 );
    for( int i = 0;i < len;i ++ ) cnt[i] += int( f[i].x + 0.5 );
}

signed main() {
    scanf( "%d %d %d %s %s", &n, &m, &k, s, t );
    for( len = 1;len <= n + m;len <<= 1, l ++ );
    for( int i = 0;i < len;i ++ ) r[i] = r[i >> 1] >> 1 | ( i & 1 ) << l - 1;
    work( 'A' ), work( 'T' ), work( 'G' ), work( 'C' );
    int ans = 0;
    for( int i = 0;i < len;i ++ ) if( cnt[i] == m ) ans ++;
    printf( "%d\n", ans );
    return 0;
}

---

W：CodeForces 954I Yet Another String Matching Problem

洛谷链接

两个字符串比对，通过一些操作使得相应位是一样的。这种“比对”一般都是转化为图论，然后与连通块有关。

基础模型可以认为是两个序列，连边 $A_i\to B_i$ 求连通块个数。

先只考虑两个串的距离如何计算？

连边 $S_i\to T_i$ （忽略原本就相同的位置）后，会形成若干的连通块。

一个联通块中只需要操作 $|连通块|-1$ 次就能将整个连通块变为相同的字符。

多个连通块之间是有向边的存在，是一个有向无环图，最后会归中于某个字符。

所以两个串的距离即为 $(\sum |连通块|)-1$。

具体实现可以用并查集维护字符集 { a , b , c , d , e , f } \{a,b,c,d,e,f\} {a,b,c,d,e,f}。每次合并原来不同的两个字符集的时候，答案 $+1$ 即可。

现在考虑有多个 $S$ 的子集 $S^{\prime }$ 在与 $T$ 适配。

设 $an{s}_i:$ 考虑 $S[i-|T|+1,i]$ 与 $T$ 的适配距离答案。

同样地，分字符计算贡献再累和。

枚举 $S^{\prime }$ 的字符 $x$，和 $T$ 的字符 $y$ 进行适配，$x\ne y$。

我们需要判断是否有 $S_i^{\prime }=x,T_i=y$ 的位置，如果存在就意味着是需要操作的。

就执行并查集合并，答案 $+1$，如果之前的某些字符配对组合已经让 $x,y$ 在同一个并查集，就不要加这个贡献。

每个 $i$ 不同结尾，适配距离也不同，对应的并查集也不同。所以每个位置都要自开一个独立的并查集。

本题的距离适配完全是适配的基础模型，翻转 $T$ 再卷积为 $f$，判断 $f_i$ 是否为 $1$，在 $f_i$ 的并查集里面进行合并即可。

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
#define maxn 1000000
int len, l, n, m;
int r[maxn], ans[maxn];
char s[maxn], t[maxn];

struct complex {
    double x, i;
    complex(){}
    complex( double X, double I ) { x = X; i = I; }
}f[maxn], g[maxn];
double pi = acos( -1.0 );
complex operator + ( complex u, complex v ) {
    return complex( u.x + v.x, u.i + v.i );
}
complex operator - ( complex u, complex v ) {
    return complex( u.x - v.x, u.i - v.i );
}
complex operator * ( complex u, complex v ) {
    return complex( u.x * v.x - u.i * v.i, u.x * v.i + v.x * u.i );
}

void FFT( complex *c, int opt ) {
    for( int i = 0;i < len;i ++ ) if( i < r[i] ) swap( c[i], c[r[i]] );
    for( int i = 1;i < len;i <<= 1 ) {
        complex omega( cos( pi / i ), opt * sin( pi / i ) );
        for( int j = 0;j < len;j += ( i << 1 ) ) {
            complex w( 1, 0 );
            for( int k = 0;k < i;k ++, w = w * omega ) {
                complex x = c[j + k], y = c[j + k + i] * w;
                c[j + k] = x + y, c[j + k + i] = x - y;
            }
        }
    }
    if( opt == -1 ) for( int i = 0;i < len;i ++ ) c[i].x /= len;
}

struct DSU {
    int f[10];
    void init() { for( int i = 0;i < 6;i ++ ) f[i] = i; }
    int find( int x ) { return f[x] == x ? x : f[x] = find( f[x] ); }
    void merge( int x, int y ) { x = find( x ), y = find( y ), f[y] = x; }
}dsu[maxn];

void work( int x, int y ) {
    for( int i = 0;i < len;i ++ ) f[i] = g[i] = complex( 0, 0 );
    for( int i = 0;i < n;i ++ ) f[i].x = ( s[i] - 'a' == x );
    for( int i = 0;i < m;i ++ ) g[m - i - 1].x = ( t[i] - 'a' == y );
    FFT( f, 1 ), FFT( g, 1 );
    for( int i = 0;i < len;i ++ ) f[i] = f[i] * g[i];
    FFT( f, -1 );
    for( int i = 0;i < n;i ++ ) {
        int c = int( f[i].x + 0.5 );
        if( ! c ) continue;
        if( dsu[i].find( x ) ^ dsu[i].find( y ) ) ans[i] ++, dsu[i].merge( x, y );
    }
}

int main() {
    scanf( "%s %s", s, t );
    n = strlen( s ), m = strlen( t );
    for( len = 1;len <= n + m;len <<= 1, l ++ );
    for( int i = 0;i < len;i ++ ) r[i] = r[i >> 1] >> 1 | ( i & 1 ) << l - 1;
    for( int i = 0;i < n;i ++ ) for( int j = 0;j < 6;j ++ ) dsu[i].init();
    for( int i = 0;i < 6;i ++ ) for( int j = 0;j < 6;j ++ ) if( i ^ j ) work( i, j );
    for( int i = m - 1;i < n;i ++ ) printf( "%d\n", ans[i] );
    return 0;
}