快速数论变换 NTT

前言： 在学习NTT之前，应当先掌握FFT（快速傅立叶变换）的基本知识，并能动手完成代码实现。如果有时间（心情）我会写一篇FFT的算法介绍。

在FFT中起到相当重要的作用的就是那个主n次单位根$w_n=e^{frac{2ipi}}$，一切的一切都围绕这个神奇的复数展开。但是复数的运算设计到两个double变量，浮点数的计算不仅很慢，而且经常出现精度误差。

为此我们找到了NTT作为代替，其中NTT本身是在模一个费马素数$P$的意义下展开的，我们设$G$为$P$的原根，会发现$G^{frac}$具有和$w_n=e^{frac{2ipi}}$相似的性质，可以代替$w_n$作为单位根。这个就不写详细证明了，反正也没什么要紧，知道这个性质的存在就好。

所以，我们现在用$G^{frac}$代替原来的$w_n$，继续做原来的分治过程。原本的除法全部改成乘以逆元即可。

在51Nod上1027和1028都可以拿来练习模板，好想Uoj的#37也是高精度乘法的样子。

下面给出我的NTT代码实现。

#include <cstdio>
#include <cstring>

template <class T>
inline void swap(T &a, T &b)
{
	T c;
	c = a;
	a = b;
	b = c;
}

const int siz = 500005;

const int P = 998244353, G = 3;

inline int pow(int a, int b)
{
	int r = 1;
	
	while (b)
	{
		if (b & 1)
			r = 1LL * r * a % P;
		
		b >>= 1, a = 1LL * a * a % P;
	}
	
	return r;
}

inline void calculateNTT(int *s, int n, int f)
{
	{
		int cnt = 0;
		
		static int rev[siz];
		
		while (n >> cnt)++cnt; --cnt;
		
		memset(rev, 0, sizeof rev);
		
		for (int i = 0; i < n; ++i)
		{
			rev[i] |= rev[i >> 1] >> 1;
			rev[i] |= (i & 1) << (cnt - 1);
		}
		
		for (int i = 0; i < n; ++i)if (i < rev[i])swap(s[i], s[rev[i]]);
	}
	
	{
		for (int i = 1; i < n; i <<= 1)
		{
			int wn = pow(G, (P - 1) / (i * 2));
			
			if (f == -1)wn = pow(wn, P - 2);
			
			for (int j = 0; j < n; j += (i << 1))
			{
				int wk = 1;
				
				for (int k = 0; k < i; ++k, wk = 1LL * wk * wn % P)
				{
					int x = s[j + k];
					int y = 1LL * s[i + j + k] * wk % P;
					
					s[j + k] = x + y;
					s[i + j + k] = x - y;
					
					s[j + k] = (s[j + k] % P + P) % P;
					s[i + j + k] = (s[i + j + k] % P + P) % P;
				}
			}
		}
	}
	
	{
		if (f == -1)
		{
			int inv = pow(n, P - 2);
			
			for (int i = 0; i < n; ++i)
				s[i] = 1LL * s[i] * inv % P;
		}
	}
}

signed main(void)
{
	static char sa[siz]; 
	static char sb[siz];
	
	scanf("%s", sa);
	scanf("%s", sb);
	
	static int la, a[siz];
	static int lb, b[siz];
	
	la = strlen(sa);
	lb = strlen(sb);
	
	for (int i = 0; i < la; ++i)a[i] = sa[la - i - 1] - '0';
	for (int i = 0; i < lb; ++i)b[i] = sb[lb - i - 1] - '0';
	
	int len; for (len = 1; len < la || len < lb; len <<= 1);
	
	calculateNTT(a, len << 1, +1);
	calculateNTT(b, len << 1, +1);
	
	
	for (int i = 0; i < len << 1; ++i)a[i] = 1LL * a[i] * b[i] % P;
	
	calculateNTT(a, len << 1, -1); 
	
	for (int i = 0; i < len << 1; ++i)a[i + 1] += a[i] / 10, a[i] = a[i] % 10;
	
	len <<= 1; while (!a[len])--len;
	
	for (int i = len; ~i; --i)printf("%d", a[i]); puts("");
}

@Author: YouSiki