快速傅里叶变换 & 快速数论变换

快速傅里叶变换 & 快速数论变换


[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},\)

\[A(\omega_n^k)=A_0(\omega_\frac{n}{2}^k)+\omega_n^kA_1(\omega_\frac{n}{2}^k)\\ A(\omega_n^{k+\frac{n}{2}})=A_0(\omega_\frac{n}{2}^k)-\omega_n^kA_1(\omega_\frac{n}{2}^k) \]

傅里叶逆变换

在单位复数根处插值
矩阵证明略
\(\omega_n^{-1}\)代替\(\omega_n\),计算结果每个元素除以\(n\)即可


实现

\(\omega\)可以预处理也可以递推,预处理精度更高
递归结束时每个元素所在的位置就是“二进制翻转”的位置,可以非递归的实现fft
加倍次数界,两个次数界为n的多项式相乘,次数界为2n-1,加倍到第一个大于等于的2的幂


注意:

  1. 我传入的参数是次数界n,最高次数n-1,数组中用0到n-1表示
  2. 取整用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 #include #include #include #include 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; }

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]);
}

posted @ 2017-03-30 21:46  Candy?  阅读(7465)  评论(0编辑  收藏  举报