多项式算法2：NTT（快速数论变换）

前言

算法简介
上次的FFT大量使用浮点数，有精度不好控制的问题，使用可以保证精度的方法又容易超时。这里NTT是在取模的一个前提下找到一个新的且为整数的单位根（在取模意义下），我们记其为$G$，那么用这个根去替代原来的${\omega }_n$来进行运算。因为所有运算就会变成整数运算，所以精度误差大大降低。使用NTT也方便按照题目要求进行取模运算，缺点就是只能进行以整数为系数的多项式运算，同时数组长度和每一项的值（如果题目原本不要求取模的话）都有限制。

前置知识

阶
若正整数$a$，$p$互质且满足$p>1$，则：
对于$a^n\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$最小的$n$，我们称之为$a$模$p$的阶，记作${\delta }_p(a)$。
对于$i\in [0,{\delta }_p(a))$，所有的$a^i\mathrm{m}\mathrm{o}\mathrm{d}p$结果互不相同。
我们假设存在两个整数$j$，$k$满足$$a^j\equiv a^k(\mathrm{m}\mathrm{o}\mathrm{d}p),0\le j<k<{\delta }_p(a)$$则有$a^{k-j}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，与阶的定义矛盾。
原根
设$n$是正整数，$g$是整数，若$g$模$n$的阶${\delta }_n(g)=\phi (n)$且$\gcd (g,n)=1$，则称$g$为$n$的一个原根。
若$\gcd (g,n)=1$且$n>0$，则$g$为$n$的一个原根的充要条件为： S = { g 1 , g 2 , g 3 , ⋯   , g φ ( n ) } S=\{g^1,g^2,g^3,\cdots,g^{\varphi(n)} \} S={g1,g2,g3,⋯,gφ(n)}为$n$的一组简化剩余系（即一组包含所有模$n$后与$n$互质且取模后的数两两不相同的数）。
由上文可知若$g$为原根，则$S$中任意两个元素在模$n$的意义下不同余。
又因为$\gcd (g,n)=1$，故而$S$为$n$的一组简化剩余系。
FFT
FFT详解

正文

关于原根和模数的选择
由于NTT大量二分数组，我们希望选一些形如$p=q\times 2^k+1$的质数$p$，其中$q$为奇自然数，$k$为整数。
我们有以下的原根可以选。

模数	原根	原根的逆元	最大长度
$469762049=7\times 2^{26}+1$	3	156587350	$2^{26}$
$998244353=119\times 2^{23}+1$	3	332748118	$2^{23}$
$2281701377=17\times 2^{27}+1$	3	760567126	$2^{27}$
$1004535809=479\times 2^{21}+1$	3	334845270	$2^{21}$

DFT
通过观察，我们可以发现原根有一些和复数单位根很像的性质：
我们令$n$为大于$1$的$2$的幂，$p$为素数且$n|(p-1)$，$g$为$p$的一个原根。
我们设$$g_n=g^{\frac{p-1}{n}}$$所以$$g_n^n=g^{n\times \frac{p-1}{n}}=g^{p-1}$$$$g_n^{\frac{n}{2}}=g^{\frac{p-1}{2}}$$$$g_{an}^{ak}=g^{\frac{ak(p-1)}{an}}=g^{\frac{k(p-1)}{n}}=g_n^k$$通过观察不难得出：$$g_n^n\equiv g^{p-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$$又因为$$(g_n^{\frac{n}{2}}{)}^2\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$$，且同属于一个简化剩余系，不能相等，
故$$g_n^{\frac{n}{2}}\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}p)$$同时我们有$$(g_n^{k+\frac{n}{2}}{)}^2\equiv (g_n^k{)}^2\times g_n^n\equiv (g_n^k{)}^2(\mathrm{m}\mathrm{o}\mathrm{d}p)$$又因为同属于一个简化剩余系，不能相等，可得$$g_n^{k+\frac{n}{2}}\equiv -g_n^k(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
总而言之，单位根在FFT中用到的性质原根都有。
$$g_n^n\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
$$g_{an}^{ak}\equiv g_n^k(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
$$g_n^{\frac{n}{2}}\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
$$g_n^{k+\frac{n}{2}}\equiv -g_n^k(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
所以在DFT的过程中带入$g_n^k$和$g_n^{k+\frac{n}{2}}$实质上和带入${\omega }_n^k$和${\omega }_n^{k+\frac{n}{2}}$一样。
IDFT
和FFT一样，用矩阵来辅助思考：
$$\left[\begin{array}{ccccc}{g}_n^0& g_n^0& g_n^0& \cdots & g_n^0\\ g_n^0& g_n^1& g_n^2& \cdots & g_n^{n-1}\\ g_n^0& g_n^2& g_n^4& \cdots & g_n^{2n-2}\\ g_n^0& g_n^3& g_n^6& \cdots & g_n^{3n-3}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ g_n^0& g_n^{n-1}& g_n^{2n-2}& \cdots & g_n^{(n-1{)}^2}\end{array}\right]\ast \left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ ⋮\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}{F}_0\\ F_1\\ F_2\\ F_3\\ ⋮\\ F_{n-1}\end{array}\right]$$
和FFT不太一样，乘以单位根的共轭复数的过程将变成乘以原根在模$p$意义下的逆元，最后再乘以长度$n$模$p$下的逆元即可。
证明如下（以下运算均在模$p$的条件下进行）：
$$g_n^k=g^{(\frac{p-1}{n})k}$$
如上矩阵，$$F(k)=\sum _{j=0}^{n-1}a(j)(g_n^i{)}^{jk}$$
逆变换为：
$$\frac{1}{n}\sum _{k=0}^{n-1}F(k)(g_n^i{)}^{-jk}=\frac{1}{n}\sum _{k=0}^{n-1}(\sum _{l=0}^{n-1}a(l)(g_n^i{)}^{lk})(g_n^i{)}^{-jk}=\frac{1}{n}\sum _{k=0}^{n-1}(\sum _{l=0}^{n-1}a(l)(g_n^i{)}^{(l-j)k})$$
式子$$=\frac{1}{n}\sum _{k=0}^{n-1}(a(0)(g_n^i{)}^{(0-j)k}+a(1)(g_n^i{)}^{(1-j)k}+\cdots +a(n-1)(g_n^i{)}^{((n-1)-j)k})$$
展开$k$的求和式得
$$\begin{gathered} =\frac{1}{n}[a(0)((g_n^i{)}^{(0-j)0}+(g_n^i{)}^{(0-j)1}+\cdots +(g_n^i{)}^{(0-j)(n-1)})+ \\ a(1)((g_n^i{)}^{(1-j)0}+(g_n^i{)}^{(1-j)1}+\cdots +(g_n^i{)}^{(1-j)(n-1)})+ \\ \cdots + \\ a(n-1)((g_n^i{)}^{((n-1)-j)0}+(g_n^i{)}^{((n-1)-j)1}+\cdots +(g_n^i{)}^{((n-1)-j)(n-1)})] \end{gathered}$$
由等比数列的求和公式得：
$$a(x)((g_n^i{)}^{(x-j)0}+(g_n^i{)}^{(x-j)1}+\cdots +(g_n^i{)}^{(x-j)(n-1)})=a(x)\frac{1-((g_n^i{)}^{n(x-j)}{)}^n}{1-(g_n^i{)}^{(x-j)}}$$当$j\ne x$，可得：$$a(x)\frac{1-((g_n^i{)}^{(x-j)}{)}^n}{1-(g_n^i{)}^{(x-j)}}=a(x)\frac{1-(g^i{)}^{(x-j)}}{1-(g_n^i{)}^{(x-j)}}=a(x)\frac{1-1^{(x-j)}}{1-(g_n^i{)}^{(x-j)}}=0$$
当$j=x$，可得：$$a(x)((g_n^i{)}^{(x-j)0}+(g_n^i{)}^{(x-j)1}+\cdots +(g_n^i{)}^{(x-j)(n-1)})=na(x)$$
因而：$$\frac{1}{n}\sum _{k=0}^{n-1}F(k)(g_n^i{)}^{-jk}=\frac{1}{n}\sum _{k=0}^{n-1}(\sum _{l=0}^{n-1}a(l)(g_n^i{)}^{lk})(g_n^i{)}^{-jk}=na(x)$$
该矩阵的逆矩阵得证。
模板题
贴上DFT与IDFT结合起来的代码：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#define ll long long
using namespace std;
const int N = 1 << 22;
const int g = 3 , gi = 332748118 , mod = 998244353;//不加const时间翻四倍 
ll qw( ll a , ll b ) {
	ll ans = 1;
	while ( b ) {
		if( b & 1 ) {
			ans = ans * a % mod;
		}
		a = a * a % mod;
		b >>= 1;
	}
	return ans;
}
int rev[N];
int n , m , l;
ll a[N] , b[N];
void pre( int bit ) {
	for ( int i = 0 ; i < ( 1 << bit ) ; ++i ) {
		rev[i] = (rev[i>>1]>>1)|((i&1)<<(bit - 1));
	}
}
void NTT( ll *F , int len , int on ) {
	for ( int i = 0 ; i < len ; ++i ) {//把递归的底层交换好
		if ( i < rev[i] ) {
			swap( F[i] , F[rev[i]] );
		}
	}
	for ( int i = 2 ; i <= len ; i <<= 1 ) {//枚举步长，从递归的下面往上走 
		ll gn = qw( on ? g : gi , ( mod - 1 ) / ( i ) );
		for ( int j = 0 ; j <= len - 1 ; j += i ) {//走一遍步长 
			ll gg = 1;
			for ( int k = j ; k < j + i / 2 ; ++k ) {//枚举每块区间内的每一个元素
				ll u = F[k];
				ll t = gg * F[k + i / 2] % mod;
				F[k] = (u + t) % mod;
				F[k + i / 2] = ( u - t  + mod ) % mod;
				gg = gg * gn % mod;
			}
		}
	}
	return;
}

int main () {
	int len = 0;
	scanf("%d%d",&n,&m);
	for ( int i = 0 ; i <= n ; ++i ) {
		scanf("%lld",a + i);
	}
	for ( int i = 0 ; i <= m ; ++i ) {
		scanf("%lld",b + i);
	}
	len = n + m;
	int tim = 1;
	l = 0;
	while( tim <= len ) {
		tim <<= 1;
		l++;
	}
	len = tim;
	pre(l);
	NTT( a , len , 1 );
	NTT( b , len , 1 );
	for ( int i = 0 ; i <= len - 1 ; ++i ) {
		a[i] = a[i] * b[i] % mod;
	}
	NTT( a , len , 0 );
	ll inv = qw( len , mod - 2 );
	for ( int i = 0 ; i <= n + m ; ++i ) {
		printf("%lld ",a[i] * inv % mod);
	}
	return 0;
}

任意模数的坑以后再填吧QWQ。