快速傅里葉變換 & 快速數論變換
[update 3.29.2017]
前言
2月10日初學，記得那時好像是正月十五放假那一天
當時寫了手寫版的筆記
過去近50天差不多忘光了，於是復習一下，具體請看手寫版筆記
參考文獻：picks miskcoo menci 阮一峰
Fast Fourier Transform
單位復數根
虛數 復數
\(i\)，表示逆時針旋轉90度
\(a+bi\)，對應復平面上的向量
復數加法 同向量
復數乘法 “模長相乘，幅角相加”，\((a+bi)*(c+di)=ac-bd+adi+bci\)
共軛復數 實部相等，虛部互為相反數. 單位根的倒數等於共軛復數
歐拉公式 \(e^{iu}=cos(u)+isin(u)\)
單位復數根
n次單位復數根：滿足\(\omega^n=1\)的復數\(\omega, \omega_n^k = e^{\frac{2\pi i}{n}k}\)
主n次單位根 \(\omega_n = e^{\frac{2\pi i}{n}}\)
消去引理，折半引理，求和引理
\(n\)個\(n\)次單位復數根在乘法意義下形成一個群，與\((Z_n,+)\)有相同的結構，因為\(w(n,0)=w(n,n)=1\ \rightarrow\ w(n,j)*w(n,k)=w(n,(j+k) mod n)\)
FFT
離散傅里葉變換DFT
對於多項式\(A(x)=\sum\limits_{j=0}^{n-1}a_jx^j\)，代入n次單位復數根所得到的列向量就是a的離散傅里葉變換
快速傅里葉變換FFT
\(O(nlogn)\)計算離散傅里葉變換
使用分治的思想，按下標奇偶分類，\(A_0(x)\)是偶數項，\(A_1(x)\)是奇數項，則\(A(x)=A_0(x^2)+xA_1(x^2)\)，根據折半引理僅有\(\frac{n}{2}\)次單位復數根組成
\(k < \frac{n}{2},\)
傅里葉逆變換
在單位復數根處插值
矩陣證明略
用\(\omega_n^{-1}\)代替\(\omega_n\)，計算結果每個元素除以\(n\)即可
實現
\(\omega\)可以預處理也可以遞推，預處理精度更高
遞歸結束時每個元素所在的位置就是“二進制翻轉”的位置，可以非遞歸的實現fft
加倍次數界，兩個次數界為n的多項式相乘，次數界為2n-1，加倍到第一個大於等於的2的冪
注意:
- 我傳入的參數是次數界n，最高次數n-1，數組中用0到n-1表示
- 取整用floor向下取整，類型轉換是向0取整
Fast Number-Theoretic Transform
生成子群 & 原根
子群：
\(群(S,\oplus),\ (S',\oplus),\ 滿足S' \subset S，則(S',\oplus)是(S,\oplus)的子群\)
拉格朗日定理:
\(|S'| \mid |S|\)
證明需要用到陪集，得到陪集大小等於子群大小，每個陪集要么不想交要么相等，所有陪集的並是集合S，那么顯然成立。
生成子群
\(a \in S\)的生成子群\(<a>=\{a^{(k)}:\ k\ge 1\}\)，\(a\)是\(<a>\)的生成元
階：
群\(S\)中\(a\)的階是滿足\(a^r=e\)的最小的r,符號\(ord(a)\)
\(ord(a)=|<a>|\)，顯然成立
考慮群\(Z_n^*=\{[a]_n \in Zn:gcd(a,n)=1\},\ |Z_n^*| = \phi(n)\)
階就是滿足\(a^{r} \equiv 1 \pmod n\)的最小的\(r,\ ord(a)=r\)
原根
\(g滿足ord_n(g)=|Z_n^*|=\phi(n)\)，對於質數\(p\)，也就是說\(g^i \mod p, 0\le i <p\)結果互不相同
模n有原根的充要條件 \(n=2,4,p^e,2p^e\)
離散對數
\(g^t \equiv a \pmod n,\ ind_{n,g}(a)=t\)
因為g是原根，所以\(g^t\)每\(\phi(n)\)是一個周期，可以取到\(|Z_n^*|\)的所有元素
對於n是質數時，就是得到\([1,n-1]\)的所有數，就是\([0,n-2]\)到\([1,n-1]\)的映射
離散對數滿足對數的相關性質，如\(ind(ab)\equiv ind(a)+ind(b) \pmod {n-1}\)
求原根
可以證明滿足\(g^{r} \equiv 1 \pmod p\)的最小的r一定是\(p-1\)的約數
對於質數\(p\)，質因子分解\(p-1\)，若\(g^{\frac{p-1}{p_i}} \neq 1 \pmod p\)恆成立，g為p的原根
NTT
對於質數\(p=qn+1,\ n=2^m\)，原根\(g\)，則\(g^{qn} \equiv 1 \pmod p\)
將\(g_n=g^{q} \pmod p\)看做\(w_n\)的等價，滿足\(w_n\)類似的性質，如：
- \(g_n^n \equiv 1 \pmod p,\ g_n^{\frac{n}{2}} \equiv -1 \pmod p\)
這里的n(用N表示吧)可以比原來那個的n(乘法結果的長度擴展到2的冪次后的n)大，只要把\(\frac{qN}{n}\)看做q就行了
常見的\(p=1004535809=479 \cdot 2^{21} + 1,\ g=3,\quad p=998244353= 2 * 17 * 2^{23} + 1,\ g=3 \)
實現
\(g^{qn}\)就是\(e^{2\pi i}\)的等價，迭代到長度\(l\)時，\(g_l=g^{\frac{p-1}{l}}\)
或者\(w_n=g_l=g_n^{\frac{n}{l}}=g^{\frac{p-1}{l}}\)
*** 這里放一個大整數相乘的模板 ```cpp //fft #include
struct meow{
double x, y;
meow(double a=0, double b=0):x(a), y(b){}
};
meow operator +(meow a, meow b) {return meow(a.x+b.x, a.y+b.y);}
meow operator -(meow a, meow b) {return meow(a.x-b.x, a.y-b.y);}
meow operator (meow a, meow b) {return meow(a.xb.x-a.yb.y, a.xb.y+a.y*b.x);}
meow conj(meow a) {return meow(a.x, -a.y);}
typedef meow cd;
struct FastFourierTransform {
int n, rev[N];
cd omega[N], omegaInv[N];
void ini(int lim) {
n=1; int k=0;
while(n<lim) n<<=1, k++;
for(int i=0; i<n; i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(k-1));
for(int k=0; k<n; k++) {
omega[k] = cd(cos(2*PI/n*k), sin(2*PI/n*k));
omegaInv[k] = conj(omega[k]);
}
}
void fft(cd *a, cd *w) {
for(int i=0; i<n; i++) if(i<rev[i]) swap(a[i], a[rev[i]]);
for(int l=2; l<=n; l<<=1) {
int m=l>>1;
for(cd *p=a; p!=a+n; p+=l)
for(int k=0; k<m; k++) {
cd t = w[n/l*k] * p[k+m];
p[k+m]=p[k]-t;
p[k]=p[k]+t;
}
}
}
void dft(cd *a, int flag) {
if(flag==1) fft(a, omega);
else {
fft(a, omegaInv);
for(int i=0; i<n; i++) a[i].x/=n;
}
}
void mul(cd *a, cd *b, int m) {
ini(m);
dft(a, 1); dft(b, 1);
for(int i=0; i<n; i++) a[i]=a[i]*b[i];
dft(a, -1);
}
}f;
int n1, n2, m, c[N];
cd a[N], b[N];
char s1[N], s2[N];
int main() {
freopen("in","r",stdin);
scanf("%s%s",s1,s2);
n1=strlen(s1); n2=strlen(s2);
for(int i=0; i<n1; i++) a[i].x = s1[n1-i-1]-'0';
for(int i=0; i<n2; i++) b[i].x = s2[n2-i-1]-'0';
m=n1+n2-1;
f.mul(a, b, m);
for(int i=0; i<m; i++) c[i]=floor(a[i].x+0.5);
for(int i=0; i<m; i++) c[i+1]+=c[i]/10, c[i]%=10;
if(c[m]) m++;
for(int i=m-1; i>=0; i--) printf("%d",c[i]);
}
```cpp
//ntt
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N=(1<<18)+5, INF=1e9;
const double PI=acos(-1);
inline int read(){
char c=getchar();int x=0,f=1;
while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
return x*f;
}
ll P=1004535809;
ll Pow(ll a, ll b,ll P) {
ll ans=1;
for(; b; b>>=1, a=a*a%P)
if(b&1) ans=ans*a%P;
return ans;
}
struct NumberTheoreticTransform {
int n, rev[N];
ll g;
void ini(int lim) {
g=3;
n=1; int k=0;
while(n<lim) n<<=1, k++;
for(int i=0; i<n; i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(k-1));
}
void dft(ll *a, int flag) {
for(int i=0; i<n; i++) if(i<rev[i]) swap(a[i], a[rev[i]]);
for(int l=2; l<=n; l<<=1) {
int m=l>>1;
ll wn = Pow(g, flag==1 ? (P-1)/l : P-1-(P-1)/l, P);
for(ll *p=a; p!=a+n; p+=l) {
ll w=1;
for(int k=0; k<m; k++) {
ll t = w * p[k+m]%P;
p[k+m]=(p[k]-t+P)%P;
p[k]=(p[k]+t)%P;
w=w*wn%P;
}
}
}
if(flag==-1) {
ll inv=Pow(n, P-2, P);
for(int i=0; i<n; i++) a[i]=a[i]*inv%P;
}
}
void mul(ll *a, ll *b, int m) {
ini(m);
dft(a, 1); dft(b, 1);
for(int i=0; i<n; i++) a[i]=a[i]*b[i];
dft(a, -1);
}
}f;
int n1, n2, m, c[N];
ll a[N], b[N];
char s1[N], s2[N];
int main() {
freopen("in","r",stdin);
scanf("%s%s",s1,s2);
n1=strlen(s1); n2=strlen(s2);
for(int i=0; i<n1; i++) a[i] = s1[n1-i-1]-'0';
for(int i=0; i<n2; i++) b[i] = s2[n2-i-1]-'0';
m=n1+n2-1;
f.mul(a, b, m);
for(int i=0; i<m; i++) c[i]=a[i];
for(int i=0; i<m; i++) c[i+1]+=c[i]/10, c[i]%=10;
if(c[m]) m++;
for(int i=m-1; i>=0; i--) printf("%d",c[i]);
}