相关知识
时间域上的函数f(t)经过傅里叶变换（Fourier Transform）变成频率域上的F(w)，也就是用一些不同频率正弦曲线的加 权叠加得到时间域上的信号。
傅里叶逆变换是将频率域上的F(w)变成时间域上的函数f(t)，一般称(f(t))为原函数，称(F(w))为象函数。原函数和象函数构成一个傅里叶变换对。
离散傅里叶变换（DFT）是傅里叶变换在时域和频域都呈现离散的形式。
当然也就有离散傅里叶逆变换（IDFT）：
快速傅里叶变换FFT
快速傅里叶变换（FFT）利用分治的思想，将 DFT 分解成更小的 DFT进行计算，可以在(O(nlog n))时间复杂度内完成 DFT 的算法。可用来加速卷积和多项式乘法。
多项式乘法
对于(C(x)=A(x)B(x)=sum_{i=0}^{n-1}sum_{j=0}^{i}a_jb_{i-j}x^i)，（A的最高项次数加上B的最高项次数为n-1）直接计算是(O(n^2))的。我们将它们变成点值表达式
那么就可以(O(n))计算出C(x)。
将其转为点值表达式需要带入n个不同的x，利用霍纳法则，带入一个x需要(O(n))时间：
从点值表达式转换为系数表达式，用插值法，将n个不同的值带入A，唯一确定了一个系数表达式（证明可见《算法导论》）。
利用 FFT 可以在(O(nlog n)) 内完成系数到点值，以及从点值到系数表达式。
n次单位复数根
n次单位复数根(w)是满足(w^n=1)的复数。 主n次单位根为(w_n=e^{frac {2pi } ni})，其它所有n次单位复数根都是(w_n)的幂次。
我们知道复数的指数形式的定义是：
所以n个复数根均匀分布在以复平面原点为圆心的单位圆上。
单位根的性质
-
消去引理 (w^{dk}_{dn}=w^k_n)，带入定义即可证明。推论:(w^{n/2}_n=w_2=-1)
-
折半引理((w^{k+n/2}_n)^2=(w_n^k)^2)，展开即可证明。推论:(w^{k+n/2}_n=-w^{k}_n)
为了方便分治，我们将n补到2的次幂（多项式后面加0）。将n个n次单位根带入A(x)的过程就是对系数向量进行 DFT 。
考虑按指数奇偶分类:
那么
当(k< n/2)时，由消去引理
也就是相当于对两个长度为n/2的多项式(A_{even}(x)),(A_{odd}(x))进行 DFT。
指数不小于(n/2)的部分与(k+n/2)一一对应，由折半引理
所以这个过程是可以递归的，且(T(n)=2T(n/2)+O(n))。
递归的FFT（c++代码）
typedef complex<double> CD;
const double pi = acos(-1);
CD tmp[N],epsilon[N];
void init_epsilon(int n){
for(int i = 0; i < n; ++i){
epsilon[i] = CD(cos(2.0 * pi * i / n), sin(2.0 * pi * i / n));
arti_epsilon[i] = conj(epsilon[i]);
}
}
void recursive_fft(int n, CD* A,int offset, int step, CD* w){
if(n==1)return;
int m=n>>1;
recursive_fft(m,A,offset,step<<1,w);
recursive_fft(m,A,offset+step,step<<1,w);
for(int k=0;k<m;++k){
int pos=2*step*k;
tmp[k] =A[pos+offset]+w[k*step]*A[pos+offset+step];
tmp[k+m]=A[pos+offset]-w[k*step]*A[pos+offset+step];
}
for(int i=0;i<n;++i)
A[i*step+offset]=tmp[i];
}
迭代实现
但是递归需要比较大的空间，如何实现迭代的写法？
观察递归的过程，第一步：
0（000）2（010）4（100）6（110），1（001）3（011）5（101）7（111）
第二步：
0 （000）4（100），2（010） 6（110），1（001）5（101）3（011）7（111）
将下标的二进制翻转过来就是：
000，001，010，011，100，101，110，111
对应了0，1，2，3，...
用reverse(i)表示i的二进制翻转后的数。就相当于二进制从高位到低位的+1，是从左往右遇到第一个0，就改为1，左边的1改为0。
int reverse(int x){
for(int l=1<<bit_length;(x^=l)<l;l>>=l);
return x;
}
我们得到递归最后一步的数组：
0，4，2，6，1，5，3，7
就可以从下到上迭代了。
void bit_reverse(CD* A,int n){
for(int i=0,j=0;i<n;++i){
if(i>j)swap(A[i],A[j]);
for(int l=n>>1;(j^=l)<l;l>>=1);
}
}
void fft(CD* A, int n, CD* w){
bit_reverse(A,n);
for(int i=2;i<=n;i<<=1)//自下到上，i为每一层的步长，或者说子问题长度
for(int j=0,m=i>>1;j<n;j+=i)//j为偏移量，或者说这一层的每一个子问题的起点
for(int k=0;k<m;++k){//k=0..i/2，计算出子问题的第k个和第k+i/2个的值
CD b=w[n/i*k]*A[j+m+k];
A[j+m+k]=A[j+k]-b;
A[j+k]+=b;
}
}
离散傅里叶逆变换（IDFT）
将点值转回系数表达式，相当于解方程组(vec{y}=V_nvec{a})，其中(vec {a})是系数向量。(V_n)是如下范德蒙德矩阵：
那么(V_n^{-1}vec{y}=vec{a})。
可令([V_n^{-1}]_{kj}=w^{-kj}_n/n)，只需证明(V_n^{-1}V_n=I_n)：
所以只要用(w_n^{-1})替换(w_n)，并将结果每个元素除以n，就可以计算出(DDF_n^{-1})。
FFT解决高精度乘法，c++代码
注意FFT因为用了cos和sin，以及是浮点数计算，所以会有精度误差，故加了0.5。
#include <bits/stdc++.h>
using namespace std;
#define rep(i,l,r) for(int i=l;i<r;++i)
#define per(i,l,r) for(int i=r-1;i>=l;--i)
#define SZ(x) ((int)(x).size())
typedef double dd;
typedef complex<dd> CD;
const dd PI=acos(-1.0);
const int L=18,N=1<<L;
CD eps[N],inv_eps[N],f[N],g[N];
void init_eps(int p){
rep(i,0,p)eps[i]=CD(cos(PI*i*2/p),sin(PI*i*2/p)),inv_eps[i]=conj(eps[i]);
}
void fft(CD p[], int n, CD w[]){
for(int i=0,j=0;i<n;++i){
if(i>j)swap(p[i],p[j]);
for(int l=n>>1;(j^=l)<l;l>>=1);
}
for(int i=2;i<=n;i<<=1)
for(int j=0,m=i>>1;j<n;j+=i)
rep(k,0,m){
CD b=w[n/i*k]*p[j+m+k];
p[j+m+k]=p[j+k]-b;
p[j+k]+=b;
}
}
int ans[N];
int main(){
string a,b;
cin>>a>>b;
int n=max(SZ(a),SZ(b)),p=1;
while(p<n)p<<=1;p<<=1;
rep(i,0,p)f[i]=g[i]=0;
n=0;per(i,0,SZ(a))f[n++]=a[i]-'0';
n=0;per(i,0,SZ(b))g[n++]=b[i]-'0';
init_eps(p);
fft(f,p,eps);fft(g,p,eps);
rep(i,0,p)f[i]*=g[i];
fft(f,p,inv_eps);
int t=0;
rep(i,0,p){
ans[i]=t+(f[i].real()+0.5)/p;
if(ans[i]>9){t=ans[i]/10;ans[i]%=10;}
else t=0;
}
bool flag=0;
per(i,0,p)if(ans[i]||flag){
printf("%d",ans[i]);flag=1;
}
if(flag==0)puts("0");
return 0;
}
数论变换（NTT）
为了解决FFT产生的精度误差，就有了一种在模意义下的方法，运算过程只有整数。它就是NTT（Number-Theoretic Transform）。
原根
设m是正整数，a是整数，则使得同余式(a^r equiv 1 mod m) 成立的最小正整数 (r) 叫做 (a) 对于模 (m) 的指数。
若a模m的指数等于φ(m)（欧拉函数），则称a为模m的一个原根。
原根的性质
- 对于质数p来说原根 (g) 满足 (g^0,g^1,g^2,…,g^{p-1}) 构成模 (p) 的简化剩余系。
我们令(p=ccdot 2^k+1) ，那么对于2的幂n有 (n|(p-1))。设
考虑FFT中，需要用到单位根的以下性质：
-
(w_n^k(0le k < n))互不相同，保证点值表示的合法；
-
(omega_{n} ^ {2k} = omega_{n/2} ^k)，用于分治；
-
(omega_n ^ { k + frac{n}{2} } = -omega_n ^ k)，用于分治；
-
仅当(k eq 0)时，(sum_{j=0}^{n-1}(w_n^k)^j=0)，用于逆变换。
原根是否也有以上性质呢？
令(p=ccdot 2^k+1)，取(n=2^m)，则(g_n^k(0le k<n)=g^{frac {(p-1)k}="" n}="g0,g{ccdot" 2{k-m}},g{2ccdot="" 2{k-m}}…,g{(p-1){ccdot="" 2^{k-m}}})<="" span="">，互不相等。</n)=g^{frac>
这里有个小知识点，若$a_1,a_2,..a_p$ 是模p的完全剩余系，那么$aa_1+b,aa_2+b,..aa_p+b$也是模p的完全剩余系，证明：若$aa_j+baequiv a_i+b mod p$则$a_i equiv a_j mod p$，矛盾。
(g_n^{2k}=g^{frac {2k(p-1)} n}=g^{frac{k(p-1)}{n/2}}=g_{n/2}^k)
(g_n^n=g^{p-1}equiv 1mod pRightarrow g^{frac {p-1}2}equivsqrt 1 mod p)，因为(g^k)互不相同，所以(g^{frac {p-1}2}equiv -1mod p)，于是(g^{k+frac n 2}_n=g_n^kcdot g_n^{frac n 2}=g_n^kcdot g^{frac {p-1} 2}=-g_n^k)。
- 当(k eq 0)时，(sum_{j=0}^{n-1}(g_n^k)^j=frac{1-(g_n^k)^n}{1-g_n^k})，由3得，(g_n^n=1)，所以(1-(g_n^k)^n=1-(g_n^n)^k=0)。否则等于n。
因此在FFT中用原根代替单位根就是NTT了。如果要求模数任意，只要再用 中国剩余定理（CRT）合并即可。
NTT的c++代码
仍然是上面高精度乘法那题
#include <bits/stdc++.h>
using namespace std;
#define rep(i,l,r) for(int i=l;i<r;++i)
#define per(i,l,r) for(int i=r-1;i>=l;--i)
#define SZ(x) ((int)(x).size())
typedef long long LL;
const int L=18,N=1<<L;
const LL C = 479;
const LL P = (C << 21) + 1;
const LL G = 3;
LL qpow(LL a, LL b, LL m){
LL ans = 1;
for(a%=m;b;b>>=1,a=a*a%m)if(b&1)ans=ans*a%m;
return ans;
}
LL eps[N],inv_eps[N],f[N],g[N];
void init_eps(int n){
LL t=(P-1)/n, invG=qpow(G,P-2,P);
rep(i,0,n) eps[i]=qpow(G,t*i,P),inv_eps[i]=qpow(invG,t*i,P);
}
void fft(LL p[], int n, LL w[]){
for(int i=0,j=0;i<n;++i){
if(i>j)swap(p[i],p[j]);
for(int l=n>>1;(j^=l)<l;l>>=1);
}
for(int i=2;i<=n;i<<=1)
for(int j=0,m=i>>1;j<n;j+=i)
rep(k,0,m){
LL b=w[n/i*k]*p[j+m+k]%P;
p[j+m+k]=(p[j+k]-b+P)%P;
p[j+k]=(p[j+k]+b)%P;
}
}
LL ans[N];
int main(){
string a,b;
cin>>a>>b;
int n=max(SZ(a),SZ(b)),p=1;
while(p<n)p<<=1;p<<=1;
rep(i,0,p)f[i]=g[i]=0;
n=0;per(i,0,SZ(a))f[n++]=a[i]-'0';
n=0;per(i,0,SZ(b))g[n++]=b[i]-'0';
init_eps(p);
fft(f,p,eps);fft(g,p,eps);
rep(i,0,p)f[i]=f[i]*g[i]%P;
fft(f,p,inv_eps);
int t=0;
LL invp=qpow(p,P-2,P);
rep(i,0,p){
ans[i]=(t+f[i]*invp%P)%P;
if(ans[i]>9){t=ans[i]/10;ans[i]%=10;}
else t=0;
}
bool flag=0;
per(i,0,p)if(ans[i]||flag){
printf("%lld",ans[i]);flag=1;
}
if(flag==0)puts("0");
return 0;
}