数学算法汇总

数学算法介绍

主要介绍一些数学相关的算法。

多项式

• DFT:离散傅里叶变换，方式有两种， FFT 和 NTT。
• IDFT:逆离散傅里叶变换。

FFT 多项式乘法

最开始学这玩意的时候感觉非常迷，后面数学水平上去了其实也不难。

多项式有两种表示方法，一种是系数法，另一种是点值法，总所周知 $n$ 个不同点唯一确定一个 $n-1$ 次多项式。

原理：多项式的两个系数表达式相乘是 $O(n^2)$ 的，但是其点值表达式相乘却是 $O(n)$ 的，所以考虑将系数表达式转成点值表达式然后相乘。

事实上，点值表达式和系数表达式的互相转化，如果点值取特殊点，可以做到 $O(n\log n)$，即使是任意点，即多项式多点求值和多项式多点插值也可以做到 $O(n \log^2n)$

设最终多项式次数为 $n-1$，我们进行多项式乘法时选择的点值叫单位根，即 $x^n=1$ 在复数域上的所有根。

这玩意有一些性质，不过我们得先把次数变为 $n=2^k$ 形式。

无法想象发明这个东西的人是怎么想到的，可能这就是被记在历史书上的人的水平。

以下内容如果将坐标系视为极坐标系会更好理解

$$\begin{align} W_n^i &= -W_n^{i+\frac{n}{2}}\\ W_{\frac{n}{2}}^{i} &= W_{n}^{2i}\\ \end{align}$$

考虑这样一个问题，对于一个多项式 $a_0+a_1x+a_2x^2 \cdots a_{n-1}x^{n-1}$ ，我们需要同时求出它在 $W_n^{0},W_n^{1}\cdots W_n^{n-1}$ 处的取值。发现由于第一个性质，貌似可以偷个懒，因为后 $\frac{n}{2}$ 个 $W_i$ 就是前 $\frac{n}{2}$ 个数的相反数。

相反数的性质，奇数次幂的符号要改变。考虑对系数按奇偶性分类，式子变成了这个样子。

$(a_0+a_2x^2+\cdots +a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots + a_{n-1}x^{n-2})$

然后考虑前后两个部分，都是一个多项式，需要对他们各自求 $x=W_\frac{2}{n}^0,W_\frac{2}{n}^1\cdots W_\frac{2}{n}^{\frac{n}{2}-1}$ 处的取值，本质上是求 $x^2=W_{\frac{2}{n}}^0,W_{\frac{2}{n}}^2\cdots W_{\frac{2}{n}}^{n-2}$ 处的取值，结合第二个性质，woc，就是两个子问题，解决之后就可以 $O(n)$ 得到原问题的解，边界显然是 $n=1$。

复杂度 $T(n)=2T(\frac{n}{2}) + O(n) = O(n\log n)$

void FFT(int now,com *a,int op)
{
	if(now==0)
	return;
	com a1[1<<now],a2[1<<now];
	for(i=0;i<1<<now;i+=2)
	{
		a1[i/2]=a[i];
		a2[i/2]=a[i+1];
	}
	FFT(now-1,a1,op);
	FFT(now-1,a2,op);
	com w0=(com{cos(2.0*Pi/(1<<now)),op*sin(2.0*Pi/(1<<now))}),w=(com){1,0};
	for(i=0;i<1<<now-1;i++,w=w*w0)
	{
		a[i]=a1[i]+w*a2[i];
		a[i+(1<<(now-1))]=a1[i]-w*a2[i];
	}
}

但是如果递归的话，常数会比较拉跨。因为递归必然需要复制数组重新弄成一个下标 $1-n$ 的问题，无论用什么办法解决，你的高速缓存都会表示意见很大，所以考虑迭代写法。

本质上递归是一层一层合并了两个数组，那么能不能直接模拟这个合并的过程呢，答案是可以的。

观察发现本质上是将下标二进制 reverse 之后逐层合并的，我们也这么做就行。

求 reverse 可以 $O(n)$，如下（如果你不了解运算顺序，请老老实实打括号）

这个原理很简单，不看最后一位，其它位先 reverse，然后处理一下最后一位就行。

for(int i=1;i<1<<lim;i++)res[i]=res[i>>1]>>1|(i&1)<<1-lim;

下面是迭代写法代码，本质是模拟了递归合并的过程。

void FFT(com a[],com b[],int op)
{
	for(i=0;i<1<<maxn;i++)
    if(r[i]>i)swap(a[r[i]],a[i]),swap(b[r[i]],b[i]);
	for(i=1;i<=maxn;i++){
		for(j=0;j<1<<maxn;j+=1<<i){
			com w0=(com){cos(2.0*Pi/(1<<i)),op*sin(2.0*Pi/(1<<i))},w=(com){1,0};
			for(k=0;k<1<<i-1;k++,w=w*w0){
				com x=a[j+k],y=a[j+k+(1<<i-1)]*w;
				a[j+k]=x+y;
				a[j+k+(1<<i-1)]=x-y;
                
				x=b[j+k],y=b[j+k+(1<<i-1)]*w;
				b[j+k]=x+y;
				b[j+k+(1<<i-1)]=x-y;
			}
		}
	}
}

搞完了系数转点值，接下来是点值转系数。

前人告诉我们只需要将单位根改为 $W_n^{0},W_n^{-1}\cdots W_n^{-n+1}$，再做一遍系数转点值的过程就可以得到系数，但是系数会变成原来的 $n$ 倍，除掉就行。

给出简要证明，$i$ 次项的系数为 $a_i$，转一次点值

之后变为 $b_i$，再做一次变成 $c_i$

$$\begin{eqnarray} c_x &=& \sum\limits_{i=0}^{n-1}b_iW_n^{-ix} \\ &=& \sum\limits_{i=0}^{n-1} W_n^{-ix}\sum\limits_{j=0}^{n-1}a_jW_n^{ij}\\ &=& \sum\limits_{i=0}^{n-1} \sum\limits_{j=0}^{n-1}a_jW_n^{i(j-x)}\\ \end{eqnarray}$$

对于 $j=x$，贡献显然为 $\sum\limits_{i=0}^{n-1}a_xW_n^{i\times0} = na_x$

对于 $j\neq x$ 贡献为 $a_j \sum\limits_{i=0}^{n-1}(W_n^{j-x})^{i}$

对这个式子的求和用等比数列求和公式有贡献为 $\dfrac{W_n^{n(j-x)}-1}{W_n^{j-x}-1}$

显然分子为 $0$，分母不为 $0$，所以贡献是 $0$，所以结果就是 $na_x$

搞定。

NTT 多项式乘法

FFT 多项式乘法是由缺陷的，由于浮点数精度和运算速度问题，FFT 可能并不能很好的解决一些问题，所以引入了 NTT，NTT 从有限整数域中找到了这样一组具有同样优秀性质的 $W_n$，即 $g$，也就是原根。

原根的内容可以参考数学证明总结中的介绍。

注意，和 FFT 一样，NTT 也需要严格的按照 $2^k$ 取次数，因为我们利用了 $W_n^{2i} = W_{\frac{n}{2}}^{i}$ 这一重要性质

所以能取出较大的 $2^k$ 作为阶的质数才可以作为 NTT 的模数，常见的 NTT 模数是 $998244353=2^{23}\times 7\times 17 +1$ ，我们可以取它的原根 $g=3$ 作为基本单位根带入，实际上如果要找到一个应用于 $n$ 的单位根 $W_n$，需要取 $W_n=g^{\frac{p-1}{n}}$。这样它就满足了我们在 FFT 证明中用到的一切性质。

然后照着 FFT 打一遍就行，只是基本运算这些换为模 $p$ 意义下的运算就行。

void NTT(int *a,int type)
{
	for(i=0;i<1<<s;i++)
	if(rk[i]>i)swap(a[rk[i]],a[i]);
	for(int len=1;len<=s;len++)
	{
		int w=1,wn=quick(g,mod-1>>len);
		if(type==-1)wn=quick(wn,mod-2);
		for(j=0;j+(1<<len)<=1<<s;j+=1<<len,w=1)
		{
			for(k=j;k<j+(1<<len-1);k++,w=1ll*w*wn%mod)
			{
				int x=a[k],y=a[k+(1<<len-1)];
				a[k]=(x+1ll*w*y%mod)%mod,a[k+(1<<len-1)]=(x-1ll*w*y%mod)%mod;
			}
		}
	}
}

其它多项式乘法和一些优化

FFT 三次变两次

FFT 三次变两次，把 $b$ 扔到 $a$ 的虚部去，变成了 $A(x) = (a_0+b_0i) + (a_1+b_1i)x+ \cdots + (a_{n-1}+b_{n-1}i)x^{n-1}$

然后求 $A^2(x)$，得到的系数表达式的虚部就是 $2ab$。

证明显然，会比 NTT 略快。

任意模数 NTT

NTT 解决的是模意义下的乘法问题，当然，如果确保值域不超出模数，那么模意义下的结果就是正确结果。

如果模数不是 NTT 模数，就需要写任意模数 NTT。

实现方式有两种，其一是拆系数 FFT，其二是多模数 NTT 后用 CRT 合并。

多模数 NTT

一般选取 998244353 1004535809 469726049 作为 NTT 模数，它们的原根都是 3。分别进行 NTT 后，可以得到在模它们意义下的结果，用中国剩余定理合并，得到一个模约 $5\times 10^{26}$ 的数的结果，一般来说任意的模数不会超过 $10^9$，所以乘起来的结果不超过 $10^{24}$，直接用合并的结果就行。

提一下实现细节。

一种方式是写个 Num 类，里面放三个 int，重写一下加减乘除，正常做。

另一种方式是写个 Poly 类，封装个乘法，用 for 做三次，都比较阳间。

用 CRT 时记得 __int128。或者用一个科技

一共三个方程，考虑这样合并。

其中 $mod_1^{-1}$ 表示 $mod_1$ 在 $mod_2$ 意义下的逆元。

$$x = x_1 \pmod {mod1} \\ x = x_2 \pmod {mod2} \\ k_1 mod_1 + x_1 = x_2 \pmod{mod2}\\ x = (x_2-x_1)mod_1^{-1}mod_1 + x_1 \pmod {mod_2*mod_1}$$

发现这样第二次合并的时候可以直接在模需要的模数 $p$ 下进行计算。

注意，务必考虑好取模的顺序，计算第二个时不能先对 p 取模，因为可能结果还需要先模上 $mod_1mod_2mod_3$，

怎样避免 __int128 ?

把 $(x_2-x_1)mod_1^{-1}$ 对 $mod_2$ 取模即可，注意处理正负号，保证任意时刻结果在 [0,对应模数) 内。

拆系数 FFT

FFT 精度不高，一般来说 double 跑个 $10^{14}$ 问题不大（最后除掉 $n$ 以后的结果，也就是乘出来的多项式的系数，举个例子就是如果跑 $n=10^4$，那需要保证 $a_i,b_i\leq 10^5$），但是跑 $10^{24}$ 次方， FFT 表示我做不到，用 long double 也不行。

所以考虑拆系数，把每个数写成 $a_1\times2^{15} + a_2$ 的形式。

然后乘法就变成了 $a_1b_12^{30}+(a_1b_2+a_2b_1)2^{15}+a_2b_2$。

设新的四个多项式为 $A_1,A_2,B_1,B_2$

一共 8 次 FFT，但可以优化，考虑复多项式

$$P=B_1+B_2i\\ T_1=P\times A_1 \\ T_2=P\times A_2$$

发现需要的系数在 $T_1,T_2$ 的实部和虚部，对这俩做 IDFT，收工。

一共 5 次 DFT。

也可以这样:

先考虑

$$P=A_1+A_2i\\ P'=A_1-A_2i\\ Q=B_1+B_2i$$

再考虑：

$$T_1=PQ=(A_1B_1-A_2B_2)+(A_1B_2+A_2B_1)i\\ T_2=P'Q=(A_1B_1+A_2B_2)+(A_1B_2-A_2B_1)i$$

最后：

$$T_1+T_2=2(A_1B_1+A_1B_2i)\\ T_1-T_2=-2(A_2B_2-A_2B_1i)$$

弄出 $T_1,T_2$ 的系数表示法就行。其实复系数多项式乘法和普通系数多项式乘法没有什么区别。

这种 5 次的思路还有其它方式，不一一介绍了。

听说有种方法可以把两次实数域多项式的 DFT 和 IDFT 整成一次复数域的，但我觉得没啥必要会。

DFT 思路是构造一个多项式 $P=A+Bi$，对其做 DFT，然后线性求 $Q=A-Bi$ 的 DFT 结果。

IDFT 思路不会。

更推荐第一种写法，不会出现精度问题，如果效率不够再考虑改成第二种