快速傅里叶变换 & 快速数论变换
[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{2pi i}{n}k})
主n次单位根 (omega_n = e^{frac{2pi i}{n}})
消去引理，折半引理，求和引理
(n)个(n)次单位复数根在乘法意义下形成一个群，与((Z_n,+))有相同的结构，因为(w(n,0)=w(n,n)=1 ightarrow w(n,j)*w(n,k)=w(n,(j+k) mod n))
FFT
离散傅里叶变换DFT
对于多项式(A(x)=sumlimits_{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)}: kge 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, 0le 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}}
eq 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^{2pi i})的等价，迭代到长度(l)时，(g_l=g^{frac{p-1}{l}})
或者(w_n=g_l=g_n^{frac{n}{l}}=g^{frac{p-1}{l}})
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]);
}