快速數論變換 NTT

前言：
在學習NTT之前，應當先掌握FFT（快速傅立葉變換）的基本知識，並能動手完成代碼實現。如果有時間（心情）我會寫一篇FFT的算法介紹。

在FFT中起到相當重要的作用的就是那個主n次單位根\(w_n=e^{\frac{2i\pi}{n}}\)，一切的一切都圍繞這個神奇的復數展開。但是復數的運算設計到兩個double變量，浮點數的計算不僅很慢，而且經常出現精度誤差。

為此我們找到了NTT作為代替，其中NTT本身是在模一個費馬素數\(P\)的意義下展開的，我們設\(G\)為\(P\)的原根，會發現\(G^{\frac{P-1}{n}}\)具有和\(w_n=e^{\frac{2i\pi}{n}}\)相似的性質，可以代替\(w_n\)作為單位根。這個就不寫詳細證明了，反正也沒什么要緊，知道這個性質的存在就好。

所以，我們現在用\(G^{\frac{P-1}{n}}\)代替原來的\(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