垃圾ACMer的暑假训练220714

垃圾ACMer的暑假训练220714

Codeforces

Codeforces Round #624 (Div. 3) B. WeirdSort

原题指路:https://codeforces.com/contest/1311/problem/B

题意 ($2\mathrm{s}$)

有$t(1\le t\le 100)$组测试数据.给定一个下标从$1$开始的、长度为$n(1<n\le 100)$的整数序列$a[]$,元素范围$[1,100]$.现给定一个由$m(1\le m<n)$相异的整数$1\le p_i<n$组成的集合$p[]$,若$p_i\in p[]$,则可交换$a[p_i]$和$a[p_{i+1}]$.若能将$a[]$非降序排列,输出"YES";否则输出"NO".

思路I

交换相邻元素类似于冒排,则可在冒排交换前检查下标是否在$p[]$中.注意有交换下标限制后每一趟排序可能不能将最大数放到序列末尾,则原本的两层for循环可能难控制终止下标的位置,不妨改用while循环.

注意有交换下标限制的冒排可能无法将序列非降序排列,则跳出while循环后检查序列是否非降序,即检查是否存在相邻的逆序对,或直接用is_sorted()函数.

每组测试样例的时间复杂度$O(n^2)$.

代码I

int main() {
	CaseT{
		int n, m; cin >> n >> m;
		vi a(n + 1);
		for (int i = 1; i <= n; i++) cin >> a[i];
		set<int> pos;
		while (m--) {
			int x; cin >> x;
			pos.insert(x);
		}

		while (true) {
			bool flag = false;
			for (int i = 1; i < n; i++) {
				if (pos.count(i) && a[i] > a[i + 1]) {
						swap(a[i], a[i + 1]);
						flag = true;
				}
			}
			if (!flag) break;
		}

		if (is_sorted(all(a))) cout << "YES" << endl;
		else cout << "NO" << endl;
	}
}

---

思路II

注意到可交换的下标构成若干条线段,且交换只能在同一条线段内进行,不能在不同线段之间进行.可用双指针找到每条线段的起止点,只对$a[]$在起止点间的部分排序,最后检查$a[]$是否非降序排列.

每组测试样例的时间复杂度$O(n\log n)$.

代码II

int main() {
	CaseT{
		int n, m; cin >> n >> m;
		vi a(n + 1);
		for (int i = 1; i <= n; i++) cin >> a[i];
		set<int> pos;
		while (m--) {
			int x; cin >> x;
			pos.insert(x);
		}

		for (int l = 1; l <= n; l++) {  // 枚举线段左端点
			if (!pos.count(l)) continue;

			int r = l;  // 线段右端点
			while (r <= n && pos.count(r)) r++;
			sort(a.begin() + l, a.begin() + r + 1);  // 注意左闭右开
			l = r;
		}

		if (is_sorted(all(a))) cout << "YES" << endl;
		else cout << "NO" << endl;
	}
}

---

---

FFT

多项式的表示法:

①系数表示法: f ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a n x n ⇔ f ( x ) = { a 0 , a 1 , ⋯   , a n } f(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n\Leftrightarrow f(x)=\{a_0,a_1,\cdots,a_n\} f(x)=a0​+a1​x+a2​x2+⋯+an​xn⇔f(x)={a0​,a1​,⋯,an​}.

②点值表示法:用多项式函数上的$(n+1)$个点唯一确定该函数.

​ 设$\{\begin{array}{l}f(x_0)=y_0=a_0+a_1{x}_0+a_2{x}_0^2+\cdots +a_n{x}_0^n\\ f(x_1)=y_1=a_0+a_1{x}_1+a_2{x}_1^2+\cdots +a_n{x}_1^n\\ \cdots \\ f(x_n)=y_n=a_0+a_1{x}_n+a_2{x}_n^2+\cdots +a_n{x}_n^n\end{array}$,则 f ( x ) = { ( x 0 , y 0 ) , ( x 1 , y 1 ) , ⋯   , ( x n , y n ) } f(x)=\{(x_0,y_0),(x_1,y_1),\cdots,(x_n,y_n)\} f(x)={(x0​,y0​),(x1​,y1​),⋯,(xn​,yn​)}.

DFT:多项式由系数表示法转化为点值表示法.

IDFT:多项式由点值表示法转化为系数表示法.

FFT:通过取特殊的$x$加速DFT和IDFT.

对多项式乘法$F(x)=f(x)\cdot g(x)$,利用DFT从多项式的系数表示法转为点值表示法,将点值相乘后,利用IDFT还原为系数表示法即可.

设DFT过程对两多项式选取的$x$序列相同,即 { f ( x ) = { ( x 0 , f ( x 0 ) ) , ⋯   , ( x n , f ( x n ) ) } g ( x ) = { ( x 0 , g ( x 0 ) ) , ⋯   , ( x n , g ( x n ) ) } \begin{cases}f(x)=\{(x_0,f(x_0)),\cdots,(x_n,f(x_n))\} \\ g(x)=\{(x_0,g(x_0)),\cdots,(x_n,g(x_n))\} \end{cases} {f(x)={(x0​,f(x0​)),⋯,(xn​,f(xn​))}g(x)={(x0​,g(x0​)),⋯,(xn​,g(xn​))}​,

​ 则 F ( x ) = { ( x 0 , f ( x 0 ) g ( x 0 ) ) , ⋯   , ( x n , f ( x n ) g ( x n ) ) } F(x)=\{(x_0,f(x_0)g(x_0)),\cdots,(x_n,f(x_n)g(x_n))\} F(x)={(x0​,f(x0​)g(x0​)),⋯,(xn​,f(xn​)g(xn​))}.

求系数的过程若用Gauss消元,求乘积的每一项$x_i$的系数时间复杂度$O(n^2)$.

注意到$n$次单位根${\omega }_n={\mathrm{e}}^{\frac{2\pi i}{n}}=\cos \frac{2\pi }{n}+i\sin \frac{2\pi }{n}$有性质:①${\omega }_n^n=1$;②${\omega }_n^k={\omega }_{2n}^{2k}$;③${\omega }_{2n}^{k+n}=-{\omega }_{2n}^k$.

FFT的基本思想是分治.就DFT而言,它分治地求$x={\omega }_n^k$时$f(x)$的值,并将多项式分为奇次项和偶次项处理.

对长度为$2^m(m\in {\mathbb{N}}^{\ast })$、最高次为$(2^m-1)$的多项式,代入${\omega }_n^0,{\omega }_n^1,\cdots ,{\omega }_n^{n-1}(n=2^m)$共$2^m$个不同值.

以$8$项$7$次多项式$f(x)=a_0+a_{1}x+\cdots +a_7{x}^7$为例,将其按奇次项、偶次项分为两组:

​ $f(x)=(a_0+a_2{x}^2+a_4{x}^4+a_6{x}^6)+x(a_1+a_3{x}^2+a_5{x}^4+a_7{x}^6)$.

令$G(x)=a_0+a_{1}x+a_4{x}^2+a_6{x}^3,H(x)=a_1+a_{3}x+a_5{x}^2+a_7{x}^3$,则$f(x)=G(x^2)+x\cdot H(x^2)$.

注意到$\mathrm{D}\mathrm{F}\mathrm{T}[f({\omega }_n^k)]=\mathrm{D}\mathrm{F}\mathrm{T}[G(({\omega }_n^k{)}^2)]+{\omega }_n^k\cdot \mathrm{D}\mathrm{F}\mathrm{T}[H(({\omega }_n^k{)}^2)]$

​ $=\mathrm{D}\mathrm{F}\mathrm{T}[G({\omega }_n^{2k})]+{\omega }_n^k\cdot \mathrm{D}\mathrm{F}\mathrm{T}[H({\omega }_n^{2k})]=\mathrm{D}\mathrm{F}\mathrm{T}[G({\omega }_{n/2}^k)]+{\omega }_n^k\cdot \mathrm{D}\mathrm{F}\mathrm{T}[H({\omega }_{n/2}^k)]$.

​ $\mathrm{D}\mathrm{F}\mathrm{T}[f({\omega }_n^{k+n/2})]=\mathrm{D}\mathrm{F}\mathrm{T}[G({\omega }_n^{2k+n})]+{\omega }_n^{k+n/2}\cdot \mathrm{D}\mathrm{F}\mathrm{T}[H({\omega }_n^{2k+n})]$

​ $=\mathrm{D}\mathrm{F}\mathrm{T}[G({\omega }_n^{2k})]-{\omega }_n^k\cdot \mathrm{D}\mathrm{F}\mathrm{T}[H({\omega }_n^{2k})]=\mathrm{D}\mathrm{F}\mathrm{T}[G({\omega }_{n/2}^k)]-{\omega }_n^k\cdot \mathrm{D}\mathrm{F}\mathrm{T}[H({\omega }_{n/2}^k)]$.

则求出$\mathrm{D}\mathrm{F}\mathrm{T}[G({\omega }_{n/2}^k)]$和$\mathrm{D}\mathrm{F}\mathrm{T}[H({\omega }_{n/2}^k)]$后,可求出$\mathrm{D}\mathrm{F}\mathrm{T}[F({\omega }_n^k)]$和$\mathrm{D}\mathrm{F}\mathrm{T}[f({\omega }_n^{k+n/2})]$,

​ 对$G$和$H$分别递归DFT即可.总时间复杂度$O(n\log n)$.

因分治DFT只能处理长度是$2$的幂次的多项式,否则分治时左右不一样长,右边取不到系数,故在第一次DFT前需将多项式长度补成$2$的幂次,即高次系数补$0$.

[朴素递归FFT]

const cp I(0, 1);  // 虚数单位i
const int MAXN = 1 << 20;
cp tmp[MAXN];  // DFT中暂时存放系数的数组

void DFT(cp* f, int n, int op) {  // op=1时为DFT,op=-1时为IDFT
	if (n == 1) return;

	for (int i = 0; i < n; i++) tmp[i] = f[i];  // 系数拷贝到tmp[]中

	for (int i = 0; i < n; i++) {
		if (i & 1) f[i / 2 + n / 2] = tmp[i];  // 奇次项放右边
		else f[i / 2] = tmp[i];  // 偶次项放左边
	}

	cp* g = f, * h = f + n / 2;  // 偶次项、奇次项
	DFT(g, n / 2, op), DFT(h, n / 2, op);  // 递归DFT求g和h的系数

	cp cur(1, 0);  // 用于记录单位根的乘方
	cp omega(cos(2 * pi / n), sin(2 * pi * op / n));  // n次单位根,op=-1时倒数
	// cp omega = exp(I * (2 * pi / n * op));  // Euler公式的写法

	for (int k = 0; k < n / 2; k++) {
		tmp[k] = g[k] + cur * h[k];  // 偶次项
		tmp[k + n / 2] = g[k] - cur * h[k];  // 奇次项
		cur *= omega;
	}

	for (int i = 0; i < n; i++) f[i] = tmp[i];  // 系数拷贝回f[]中
}

用分治优化.注意到上述过程中每一步都将整个多项式的奇次项和偶次项系数分开,直至分到只剩下一个系数,但递归过程需更多内存,可"模仿递归"地将系数从原数组中拆分出来,再倍增地合并计算出来的值.

以$8$项$7$次多项式$f(x)=a_0+a_{1}x+\cdots +a_7{x}^7$为例:

①初始序列: { x 0 , x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 } \{x_0,x_1,x_2,x_3,x_4,x_5,x_6,x_7\} {x0​,x1​,x2​,x3​,x4​,x5​,x6​,x7​}.

②一次二分后: { x 0 , x 2 , x 4 , x 6 } , { x 1 , x 3 , x 5 , x 7 } \{x_0,x_2,x_4,x_6\},\{x_1,x_3,x_5,x_7\} {x0​,x2​,x4​,x6​},{x1​,x3​,x5​,x7​}.

③两次二分后: { x 0 , x 4 } , { x 2 , x 6 } , { x 1 , x 5 } , { x 3 , x 7 } \{x_0,x_4\},\{x_2,x_6\},\{x_1,x_5\},\{x_3,x_7\} {x0​,x4​},{x2​,x6​},{x1​,x5​},{x3​,x7​}.

④三次二分后: { x 0 } , { x 4 } , { x 2 } , { x 6 } , { x 1 } , { x 5 } , { x 3 } , { x 7 } \{x_0\},\{x_4\},\{x_2\},\{x_6\},\{x_1\},\{x_5\},\{x_3\},\{x_7\} {x0​},{x4​},{x2​},{x6​},{x1​},{x5​},{x3​},{x7​}.

$000$	$001$	$010$	$011$	$100$	$101$	$110$	$111$
$a_0$	$a_1$	$a_2$	$a_3$	$a_4$	$a_5$	$a_6$	$a_7$
$a_0$	$a_4$	$a_2$	$a_6$	$a_1$	$a_5$	$a_3$	$a_7$
$000$	$100$	$010$	$110$	$001$	$101$	$011$	$111$

规律:最后一次二分得到的排列是原序列每一项的下标的二进制表示镜像的结果,这称为位逆序置换(Bit-reversal Permutation).

[朴素位逆序置换] 时间复杂度:$O(n\log n)$

void BRP(cp* y, int len) {  // 位逆序置换,len为2的幂次,O(nlog n)
	for (int i = 1, j = len / 2; i < len - 1; i++) {  // i做正常的+1
		if (i < j) swap(y[i], y[j]);  // 交换下标互为镜像的项的系数

		// j做镜像的+1
		int k = len / 2;
		while (j >= k) {
			j -= k;
			k >>= 1;
		}
		if (j < k) j += k;
	}
}

事实上BRP可$O(n)$递推求.设$len=2^k,R(x)$为长度为$k$的二进制数$x$镜像后的数,现需求$R(0),\cdots ,R(n-1)$.

初始条件$R(0)=0$.

注意到求$R(x)$时已求得$R\left(\lfloor \frac{x}{2}\rfloor \right)$,将$x$右移一位后取反,再右移一位即得$x$除LSB外镜像的结果.

考察$x$的LSB镜像的结果:

​ ①若个位为$0$,则镜像所得数的MSB为$0$.

​ ②若个位为$1$,则镜像所的数的MSB为$1$,即加上$\frac{len}{2}=2^{k-1}$.

综上:$R(x)=\lfloor \frac{R\left(\lfloor \frac{x}{2}\rfloor \right)}{2}\rfloor +(x\mathrm{m}\mathrm{o}\mathrm{d}2)\cdot \frac{len}{2}$.

[递推位逆序置换] 时间复杂度:$O(n)$

int rev[MAXN];  // rev[x]表示二进制数x镜像后的结果
void BRP(cp* y, int len) {  // 位逆序置换,len为2的幂次,O(n)
	for (int i = 0; i < len; i++) {
		rev[i] = rev[i >> 1] >> 1;
		if (i & 1) rev[i] |= (len >> 1);  // 若i是奇数,则翻转为len/2
	}

	for (int i = 0; i < len; i++)  // 交换下标互为镜像的项的系数
		if (i < rev[i]) swap(y[i], y[rev[i]]);
}

IDFT将多项式的点值表示法转化为系数表示法,而DFT本身是线性变换,可将多项式视为向量,左乘一个矩阵得到变换后的向量,模拟将单位根代入多项式的过程,即$\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ ⋮\\ y_{n-1}\end{array}\right]\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& {\omega }_n^1& {\omega }_n^2& \cdots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& \cdots & {\omega }_n^{2(n-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& \cdots & {\omega }_n^{(n-1{)}^2}\end{array}\right]=\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_{n-1}\end{array}\right]$.

现已知$(y_0,y_1,\cdots ,y_{n-1}{)}^{\prime }$,只需左乘$\mathrm{R}\mathrm{H}\mathrm{S}$的大矩阵$M$的逆即得$(a_0,a_1,\cdots ,a_{n-1}{)}^{\prime }$.其中$M$的逆是每一项取倒数再除以$n$.

注意到$\frac{1}{{\omega }_k}={\omega }_k^{-1}={\mathrm{e}}^{-\frac{2\pi i}{k}}$,只需将单位根取为${\mathrm{e}}^{-\frac{2\pi i}{k}}$,计算结果即变为原来的倒数.故可在函数中传参数$op$,取单位根为${\mathrm{e}}^{op\cdot \frac{2\pi i}{k}}$,则$op=1$时为DFT,$op=-1$时为IDFT.

利用单位根的周期性可理解DFT与IDFT的关系.

设原多项式$f(x)=\sum _{i=0}^{n-1}{a}_i{x}^i$.已知$y_i=f({\omega }_n^i)(0\le i\le n-1)$,求$a_i$.

构造多项式$A(x)=\sum _{i=0}^{n-1}{y}_i{x}^i$,即将 { y 0 , y 1 , ⋯   , y n − 1 } \{y_0,y_1,\cdots,y_{n-1}\} {y0​,y1​,⋯,yn−1​}视为$A$的系数表示法.

[法I] 设$b_i={\omega }_n^{-i}$,则$A(x)$在$x=b_0,b_1,\cdots ,b_{n-1}$处的点值表示法为 { ( b 0 , A ( b 0 ) ) , ( b 1 , A ( b 1 ) ) , ⋯   , ( b n − 1 , A ( b n − 1 ) ) } \{(b_0,A(b_0)),(b_1,A(b_1)),\cdots,(b_{n-1},A(b_{n-1}))\} {(b0​,A(b0​)),(b1​,A(b1​)),⋯,(bn−1​,A(bn−1​))}.

而$A(b_k)=\sum _{i=0}^{n-1}f({\omega }_n^i){\omega }_n^{-ik}=\sum _{i=0}^{n-1}{\omega }_n^{-ik}\sum _{j=0}^{n-1}{a}_j({\omega }_n^i{)}^j=\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j{\omega }_n^{i(j-k)}=\sum _{j=0}^{n-1}{a}_j\sum _{i=0}^{n-1}({\omega }^{j-k}{)}^i$,

令$S({\omega }_n^t)=\sum _{i=0}^{n-1}({\omega }_n^t{)}^i$,则$t\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}n)$时,$S({\omega }_n^t)=n$;

​ $t\not\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}n)$时,${\omega }_n^{t}S({\omega }_n^t)=\sum _{i=1}^n({\omega }_n^t{)}^i$,错位相减得$S({\omega }_n^t)=\frac{({\omega }_n^t{)}^n-({\omega }_n^t{)}^0}{{\omega }_n^t-1}=0$,

​ 即$S({\omega }_n^t)=\{\begin{array}{l}n,a=0\\ 0,a\ne 0\end{array}$,进而$A(b_k)=\sum _{j=0}^{n-1}{a}_{j}S({\omega }_n^{j-k})=a_k\cdot n$,

进而$A(x)$的点值表示法为 { ( b 0 , a 0 ⋅ n ) , ( b 1 , a 1 ⋅ n ) , ⋯   , ( b n − 1 , a n − 1 ⋅ n ) } \{(b_0,a_0\cdot n),(b_1,a_1\cdot n),\cdots,(b_{n-1},a_{n-1}\cdot n)\} {(b0​,a0​⋅n),(b1​,a1​⋅n),⋯,(bn−1​,an−1​⋅n)}.

做法:取单位根为其倒数,对 { y 0 , y 1 , ⋯   , y n − 1 } \{y_0,y_1,\cdots,y_{n-1}\} {y0​,y1​,⋯,yn−1​}跑一遍FFT,再除以$n$即得$f(x)$的系数表示法.

[迭代FFT写法I] 时间复杂度:$O(n\log n)$

void BRP(cp* y, int len) {  // 位逆序置换,len为2的幂次,O(nlog n)
	for (int i = 1, j = len / 2; i < len - 1; i++) {  // i做正常的+1
		if (i < j) swap(y[i], y[j]);  // 交换下标互为镜像的项的系数

		// j做镜像的+1
		int k = len / 2;
		while (j >= k) {
			j -= k;
			k >>= 1;
		}
		if (j < k) j += k;
	}
}

void FFT(cp* y, int len, int op) {  // len为2的幂次,op=1时为DFT,op=-1时为IDFT
	BRP(y, len);  // DFT前做一遍BRP,将系数位置调整为递归后的位置

	for (int l = 2; l <= len; l <<= 1) {  // 模拟合并过程,l为当前区间长度
		cp omega(cos(2 * pi / l), sin(op * 2 * pi / l));  // 单位根

		for (int j = 0; j < len; j += l) {
			cp cur(1, 0);  // 记录单位根的乘方
			for (int k = j; k < j + l / 2; k++) {
				cp p = y[k], q = cur * y[k + l / 2];
				y[k] = p + q, y[k + l / 2] = p - q;  // 偶次项、奇次项

				cur *= omega;
			}
		}
	}

	if (op == -1)  // IDFT还需将系数除以len
		for (int i = 0; i < len; i++) y[i] = { y[i].real() / len , y[i].imag() };
}

[法II] 将${\omega }_n^i$代入$A(x)$,类似于法I推得$A({\omega }_n^k)=\sum _{j=0}^{n-1}{a}_{j}S({\omega }_n^{j+k})$,其中$S({\omega }_n^{j+k})=\{\begin{array}{l}n,i+k\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}n)\\ 0,i+k\not\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}n)\end{array}$.

故$A({\omega }_n^k)=a_{n-k}\cdot n$.这表面:对 { y 0 , y 1 , ⋯   , y n − 1 } \{y_0,y_1,\cdots,y_{n-1}\} {y0​,y1​,⋯,yn−1​}跑一遍DFT后再除以$n$,反转后$(n-1)$个元素即得$f(x)$的系数表示法.

[迭代FFT写法II] 时间复杂度:$O(n\log n)$

void BRP(cp* y, int len) {  // 位逆序置换,len为2的幂次,O(nlog n)
	for (int i = 1, j = len / 2; i < len - 1; i++) {  // i做正常的+1
		if (i < j) swap(y[i], y[j]);  // 交换下标互为镜像的项的系数

		// j做镜像的+1
		int k = len / 2;
		while (j >= k) {
			j -= k;
			k >>= 1;
		}
		if (j < k) j += k;
	}
}

void FFT(cp* y, int len, int op) {  // len为2的幂次,op=1时为DFT,op=-1时为IDFT
	BRP(y, len);  // DFT前做一遍BRP,将系数位置调整为递归后的位置

	for (int l = 2; l <= len; l <<= 1) {  // 模拟合并过程,l为当前区间长度
		cp omega(cos(2 * pi / l), sin(op * 2 * pi / l));  // 单位根

		for (int j = 0; j < len; j += l) {
			cp cur(1, 0);  // 记录单位根的乘方
			for (int k = j; k < j + l / 2; k++) {
				cp p = y[k], q = cur * y[k + l / 2];
				y[k] = p + q, y[k + l / 2] = p - q;  // 偶次项、奇次项

				cur *= omega;
			}
		}
	}

	if (op == -1) {  // IDFT还需将系数除以len
		reverse(y + 1, y + len);  // 反转后(n-1)个元素
		for (int i = 0; i < len; i++) y[i] = { y[i].real() / len , y[i].imag() };
	}
}

---

---

31.1.2 A*B Problem 升级版

题意 ($1.5\mathrm{s}$)

输入整数$a,b(1\le a,b\le 10^{1\mathrm{e}6})$,输出$a\times b$.

思路I

将整数视为多项式当$x=10$时的值,做FFT.

代码I

struct Complex {
	double x, y;

	Complex(double _x = 0.0, double _y = 0.0) :x(_x), y(_y) {}

	Complex operator+(const Complex& b)const { return Complex(x + b.x, y + b.y); }
	Complex operator-(const Complex& b)const { return Complex(x - b.x, y - b.y); }
	Complex operator*(const Complex& b)const { return Complex(x * b.x - y * b.y, x * b.y + y * b.x); }
};

void BRP(Complex* y, int len) {  // 位逆序置换,len为2的幂次,O(nlog n)
	for (int i = 1, j = len / 2; i < len - 1; i++) {  // i做正常的+1
		if (i < j) swap(y[i], y[j]);  // 交换下标互为镜像的项的系数

		// j做镜像的+1
		int k = len / 2;
		while (j >= k) {
			j -= k;
			k >>= 1;
		}
		if (j < k) j += k;
	}
}

void FFT(Complex* y, int len, int op) {  // len为2的幂次,op=1时为DFT,op=-1时为IDFT
	BRP(y, len);  // DFT前做一遍BRP,将系数位置调整为递归后的位置

	for (int l = 2; l <= len; l <<= 1) {  // 模拟合并过程,l为当前区间长度
		Complex omega(cos(2 * pi / l), sin(op * 2 * pi / l));  // 单位根

		for (int j = 0; j < len; j += l) {
			Complex cur(1, 0);  // 记录单位根的乘方
			for (int k = j; k < j + l / 2; k++) {
				Complex p = y[k], q = cur * y[k + l / 2];
				y[k] = p + q, y[k + l / 2] = p - q;  // 偶次项、奇次项

				cur = cur * omega;
			}
		}
	}

	if (op == -1)  // IDFT还需将系数除以len
		for (int i = 0; i < len; i++) y[i].x /= len;
}

const int MAXN = 2e6 + 5;
char str1[MAXN], str2[MAXN];
Complex a[MAXN << 1], b[MAXN << 1];
int ans[MAXN << 1];

int main() {
	cin >> str1 >> str2;
	int len1 = strlen(str1), len2 = strlen(str2);

	int len = 1;
	while (len < (len1 << 1) || len < (len2 << 1)) len <<= 1;  // 多项式项数补成2的幂次

	// 倒序存多项式系数,长度不足补0
	for (int i = 0; i < len1; i++) a[i] = Complex(str1[len1 - 1 - i] & 15, 0);
	for (int i = len1; i < len; i++) a[i] = Complex(0, 0);
	for (int i = 0; i < len2; i++) b[i] = Complex(str2[len2 - 1 - i] & 15, 0);
	for (int i = len2; i < len; i++) b[i] = Complex(0, 0);

	FFT(a, len, 1), FFT(b, len, 1);  // DFT
	for (int i = 0; i < len; i++) a[i] = a[i] * b[i];  // 点值相乘
	FFT(a, len, -1);  // IDFT

	for (int i = 0; i < len; i++) ans[i] = int(a[i].x + 0.5);  // +0,5防止double的精度误差
	for (int i = 0; i < len; i++) {  // 整理为十进制表示
		ans[i + 1] += ans[i] / 10;
		ans[i] %= 10;
	}

	len = len1 + len2 - 1;
	while (len && !ans[len]) len--;  // 去除前导零

	for (int i = len; i >= 0; i--) cout << ans[i];  // 倒序输出答案的数码
}

---

思路II

注意到$(a+bi{)}^2=(a^2-b^2)+2abi$,可将两多项式的值作为同一复数的实部和虚部,计算该复数的平方后,虚部除以$2$即得结果.

代码II:三次转两次优化

struct Complex {
	double x, y;

	Complex(double _x = 0.0, double _y = 0.0) :x(_x), y(_y) {}

	Complex operator+(const Complex& b)const { return Complex(x + b.x, y + b.y); }
	Complex operator-(const Complex& b)const { return Complex(x - b.x, y - b.y); }
	Complex operator*(const Complex& b)const { return Complex(x * b.x - y * b.y, x * b.y + y * b.x); }
};

void BRP(Complex* y, int len) {  // 位逆序置换,len为2的幂次,O(nlog n)
	for (int i = 1, j = len / 2; i < len - 1; i++) {  // i做正常的+1
		if (i < j) swap(y[i], y[j]);  // 交换下标互为镜像的项的系数

		// j做镜像的+1
		int k = len / 2;
		while (j >= k) {
			j -= k;
			k >>= 1;
		}
		if (j < k) j += k;
	}
}

void FFT(Complex* y, int len, int op) {  // len为2的幂次,op=1时为DFT,op=-1时为IDFT
	BRP(y, len);  // DFT前做一遍BRP,将系数位置调整为递归后的位置

	for (int l = 2; l <= len; l <<= 1) {  // 模拟合并过程,l为当前区间长度
		Complex omega(cos(2 * pi / l), sin(op * 2 * pi / l));  // 单位根

		for (int j = 0; j < len; j += l) {
			Complex cur(1, 0);  // 记录单位根的乘方
			for (int k = j; k < j + l / 2; k++) {
				Complex p = y[k], q = cur * y[k + l / 2];
				y[k] = p + q, y[k + l / 2] = p - q;  // 偶次项、奇次项

				cur = cur * omega;
			}
		}
	}

	if (op == -1)  // IDFT还需将系数除以len
		for (int i = 0; i < len; i++) y[i].x /= len;
}

const int MAXN = 2e6 + 5;
char str1[MAXN], str2[MAXN];
Complex a[MAXN << 1];
int ans[MAXN << 1];

int main() {
	cin >> str1 >> str2;
	int len1 = strlen(str1), len2 = strlen(str2);
	
	int len = 1;
	while (len < (len1 << 1) || len < (len2 << 1)) len <<= 1;  // 多项式项数补成2的幂次

	// 倒序存多项式系数
	for (int i = 0; i < min(len1, len2); i++)
		a[i] = Complex(str1[len1 - 1 - i] & 15, str2[len2 - 1 - i] & 15);
	
	// 长度不足补0
	if (len1 >= len2)
		for (int i = len2; i < len1; i++) a[i] = Complex(str1[len1 - 1 - i] & 15, 0);
	else
		for (int i = len1; i < len2; i++) a[i] = Complex(0, str2[len2 - 1 - i] & 15);
	for (int i = max(len1, len2); i < len; i++) a[i] = Complex(0, 0);

	FFT(a, len, 1);  // DFT
	for (int i = 0; i < len; i++) a[i] = a[i] * a[i];
	FFT(a, len, -1);  // IDFT

	double k = 0.5 / len;  // 用于计算虚部
	for (int i = 0; i < len; i++) ans[i] = int(a[i].y * k + 0.5);  // +0,5防止double的精度误差

	for (int i = 0; i < len; i++) {  // 整理为十进制表示
		ans[i + 1] += ans[i] / 10;
		ans[i] %= 10;
	}

	len = len1 + len2 - 1;
	while (len && !ans[len]) len--;  // 去除前导零

	for (int i = len; i >= 0; i--) cout << ans[i];  // 倒序输出答案的数码
}

---

---

31.1.2 多项式乘法

题意 ($2\mathrm{s}$)

给定$n$次多项式$f(x)$和$m$次多项式$g(x)$,求它们的卷积.

第一行输入整数$n,m(1\le n,m\le 1\mathrm{e}6)$.第二行输入$(n+1)$个整数,从低到高表示$f(x)$的系数.第三行输入$(m+1)$个整数,从低到高表示$g(x)$的系数.数据保证系数再$[0,9]$范围内.

输出$(n+m-1)$个数字,从低到高表示$f(x)\cdot g(x)$的系数.

代码

struct Complex {
	double x, y;

	Complex(double _x = 0.0, double _y = 0.0) :x(_x), y(_y) {}

	Complex operator+(const Complex& b)const { return Complex(x + b.x, y + b.y); }
	Complex operator-(const Complex& b)const { return Complex(x - b.x, y - b.y); }
	Complex operator*(const Complex& b)const { return Complex(x * b.x - y * b.y, x * b.y + y * b.x); }
};

void BRP(Complex* y, int len) {  // 位逆序置换,len为2的幂次,O(nlog n)
	for (int i = 1, j = len / 2; i < len - 1; i++) {  // i做正常的+1
		if (i < j) swap(y[i], y[j]);  // 交换下标互为镜像的项的系数

		// j做镜像的+1
		int k = len / 2;
		while (j >= k) {
			j -= k;
			k >>= 1;
		}
		if (j < k) j += k;
	}
}

void FFT(Complex* y, int len, int op) {  // len为2的幂次,op=1时为DFT,op=-1时为IDFT
	BRP(y, len);  // DFT前做一遍BRP,将系数位置调整为递归后的位置

	for (int l = 2; l <= len; l <<= 1) {  // 模拟合并过程,l为当前区间长度
		Complex omega(cos(2 * pi / l), sin(op * 2 * pi / l));  // 单位根

		for (int j = 0; j < len; j += l) {
			Complex cur(1, 0);  // 记录单位根的乘方
			for (int k = j; k < j + l / 2; k++, cur = cur * omega) {
				Complex p = y[k], q = cur * y[k + l / 2];
				y[k] = p + q, y[k + l / 2] = p - q;  // 偶次项、奇次项
			}
		}
	}

	if (op == -1)  // IDFT还需将系数除以len
		for (int i = 0; i < len; i++) y[i].x /= len;
}

const int MAXN = 1e6 + 5;
int a[MAXN], b[MAXN];
Complex A[MAXN << 1];
int ans[MAXN << 1];

int main() {
	int n, m; cin >> n >> m;
	n++, m++;  // 注意系数分别有(n+1)、(m+1)个
	for (int i = 0; i < n; i++) cin >> a[i];
	for (int i = 0; i < m; i++) cin >> b[i];

	int len = 1;
	while (len < (n << 1) || len < (m << 1)) len <<= 1;  // 多项式项数补成2的幂次

	// 倒序存多项式系数
	for (int i = 0; i < min(n, m); i++) A[i] = Complex(a[i], b[i]);
	
	// 长度不足补0
	if (n >= m)
		for (int i = m; i < n; i++) A[i] = Complex(a[i], 0);
	else
		for (int i = n; i < m; i++) A[i] = Complex(0, b[i]);
	for (int i = max(n, m); i < len; i++) A[i] = Complex(0, 0);

	FFT(A, len, 1);  // DFT
	for (int i = 0; i < len; i++) A[i] = A[i] * A[i];
	FFT(A, len, -1);  // IDFT

	double k = 0.5 / len;  // 用于计算虚部
	for (int i = 0; i < len; i++) ans[i] = int(A[i].y * k + 0.5);  // +0,5防止double的精度误差

	for (int i = 0; i < n + m - 1; i++) cout << ans[i] << ' ';  // 倒序输出答案的数码
}

---

---

NTT

FFT涉及到浮点数运算,常数较大、精度较低.NTT是FFT在数论基础上的实现,用于解决多项式乘法带模数的问题,受模数限制.

NTT通过将DFT化为$F=\mathbb{Z}/p$,其中$p$为素数.这是一个有限域,只要整数$n$可被$(p-1)$整除,就存在本原$n$次方根,不妨设$p=qn+1(n=2^m)$,原根$gs.t.g^{qn}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$.$g_n=g^q(\mathrm{m}\mathrm{o}\mathrm{d}p)$与${\omega }_n$有类似的性质,如$g_n^n\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p),g_n^{n/2}\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}p)$.NTT中的$n$可比FFT中的$n$大,不妨将NTT中的$n$记作$N$,则将$\frac{qN}{n}$视为上述的$q$即可避免大小问题.此处$g^{qn}$的地位相当于FFT中的${\mathrm{e}}^{2\pi n}$.

迭代到长度$l$时,$g_l=g^{\frac{p-1}{l}}$,类似于${\omega }_n=g_l=g_N^{\frac{N}{l}}=g_N^{\frac{p-1}{l}}$.

常取①$p=998244353=7\times 17\times 2^{23}+1$,其原根$g=3$;

​ ②$p=1004535809=479\times 2^{21}+1$,其原根$g=3$.

---

31.2.1 A*B Problem 升级版

题意 ($1.5\mathrm{s}$)

输入整数$a,b(1\le a,b\le 10^{1\mathrm{e}6})$,输出$a\times b$.

思路

将整数视为多项式当$x=10$时的值,做NTT.

代码

const int MAXN = 3e6 + 5;
const int MOD = 998244353;

int rev[MAXN];  // rev[x]表示二进制数x镜像后的结果

void NTT(int* y, int lim, int op) {  // op=1时为NTT,op=-1时为INTT
	for (int i = 0; i < lim; i++)  // 做一遍BRP,将系数位置调整为递归后的位置
		if (rev[i] < i) swap(y[i], y[rev[i]]);

	for (int l = 2; l <= lim; l <<= 1) {  // 模拟合并过程,l为当前区间长度
		int g = qpow(3, (MOD - 1) / l, MOD);
		for (int i = 0; i < lim; i += l) {
			int cur = 1;
			int k = l >> 1;
			for (int j = 0; j < k; j++, cur = (ll)cur * g % MOD) {
				int tmp = (ll)y[i + j + k] * cur % MOD;
				y[i + j + k] = (y[i + j] - tmp + MOD) % MOD;  // 奇次项
				y[i + j] = ((ll)y[i + j] + tmp) % MOD;  // 偶次项
			}
		}
	}

	if (op == -1) {  // INTT
		reverse(y + 1, y + lim);
		int inv = qpow(lim, MOD - 2, MOD);
		for (int i = 0; i < lim; i++) y[i] = (ll)y[i] * inv % MOD;
	}
}

char a[MAXN], b[MAXN];
int A[MAXN << 1], B[MAXN << 1];
int ans[MAXN << 1];

int main() {
	int lim = 1;  // 多项式的长度

	cin >> a;
	int n = strlen(a);
	for (int i = 0; i < n; i++) A[i] = a[n - i - 1] & 15;  // 倒序存系数,下同
	while (lim < (n << 1)) lim <<= 1;  // 长度补成2的幂次

	cin >> b;
	n = strlen(b);
	for (int i = 0; i < n; i++) B[i] = b[n - i - 1] & 15;
	while (lim < (n << 1)) lim <<= 1;  // 长度补成2的幂次

	for (int i = 0; i < lim; i++)  // 长度不足补0
		rev[i] = (rev[i >> 1] >> 1) + (i & 1) * (lim >> 1);

	NTT(A, lim, 1), NTT(B, lim, 1);  // NTT
	for (int i = 0; i < lim; i++) ans[i] = (ll)A[i] * B[i] % MOD;
	NTT(ans, lim, -1);  // INTT

	int len = 0;  // ans[]的长度
	for (int i = 0; i < lim; i++) {  // 整理成十进制表示
		if (ans[i] >= 10) {
			len = i + 1;
			ans[i + 1] += ans[i] / 10, ans[i] %= 10;
		}
		if (ans[i]) len = max(len, i);
	}
	while (ans[len] >= 10) ans[len + 1] += ans[len] / 10, ans[len++] %= 10;

	for (int i = len; i >= 0; i--) cout << ans[i];
}

---

---

31.2.2 多项式乘法

题意 ($2\mathrm{s}$)

给定$n$次多项式$f(x)$和$m$次多项式$g(x)$,求它们的卷积.

第一行输入整数$n,m(1\le n,m\le 1\mathrm{e}6)$.第二行输入$(n+1)$个整数,从低到高表示$f(x)$的系数.第三行输入$(m+1)$个整数,从低到高表示$g(x)$的系数.数据保证系数再$[0,9]$范围内.

输出$(n+m-1)$个数字,从低到高表示$f(x)\cdot g(x)$的系数.

代码

const int MAXN = 4e6 + 5;
const int MOD = 998244353;

int rev[MAXN];  // rev[x]表示二进制数x镜像后的结果

void NTT(int* y, int lim, int op) {  // op=1时为NTT,op=-1时为INTT
	for (int i = 0; i < lim; i++)  // 做一遍BRP,将系数位置调整为递归后的位置
		if (rev[i] < i) swap(y[i], y[rev[i]]);

	for (int l = 2; l <= lim; l <<= 1) {  // 模拟合并过程,l为当前区间长度
		int g = qpow(3, (MOD - 1) / l, MOD);
		for (int i = 0; i < lim; i += l) {
			int cur = 1;
			int k = l >> 1;
			for (int j = 0; j < k; j++, cur = (ll)cur * g % MOD) {
				int tmp = (ll)y[i + j + k] * cur % MOD;
				y[i + j + k] = (y[i + j] - tmp + MOD) % MOD;  // 奇次项
				y[i + j] = ((ll)y[i + j] + tmp) % MOD;  // 偶次项
			}
		}
	}

	if (op == -1) {  // INTT
		reverse(y + 1, y + lim);
		int inv = qpow(lim, MOD - 2, MOD);
		for (int i = 0; i < lim; i++) y[i] = (ll)y[i] * inv % MOD;
	}
}

int a[MAXN], b[MAXN];
int ans[MAXN];

int main() {
	int lim = 1;  // 多项式的长度

	int n, m; cin >> n >> m;
	n++, m++;  // 注意系数分别有(n+1)、(m+1)个
	for (int i = 0; i < n; i++) cin >> a[i];
	while (lim < (n << 1)) lim <<= 1;  // 长度补成2的幂次

	for (int i = 0; i < m; i++) cin >> b[i];
	while (lim < (m << 1)) lim <<= 1;  // 长度补成2的幂次

	for (int i = 0; i < lim; i++)  // 长度不足补0
		rev[i] = (rev[i >> 1] >> 1) + (i & 1) * (lim >> 1);

	NTT(a, lim, 1), NTT(b, lim, 1);  // NTT
	for (int i = 0; i < lim; i++) ans[i] = (ll)a[i] * b[i] % MOD;
	NTT(ans, lim, -1);  // INTT

	for (int i = 0; i < n + m - 1; i++) cout << ans[i] << ' ';
}

---

---