多项式

  1. Q: 什么是 (FFT)?
    A: (Fast Fourier Transform) 快速傅里叶变换。
  2. Q: 什么是 (DFT)?
    A: (Discrete Fourier Transform) 离散傅里叶变换。
  3. Q: 什么是 (IDFT)?
    A: (Inverse Discrete Fourier Transform) 逆·离散傅里叶变换。
  4. Q: 什么是 (NTT)?
    A: (number theory Transform) 数论变换。
  5. Q: 什么是 (MTT)?
    A: 任意模数数论变换。

(FFT)

或许有很多人会在这里强调好些引理,但我觉得在推导到那步时能够继续推下去就可以了。

没有必要过分强调。。。

这里将多项式和函数直接对应
其中 (f_{x_i}=f(x_i))
该方法用于在(O(nlog n))求离散卷积。
什么是离散卷积?
给出定义,令(y=f imes g)

[y_n=sum _{i=0}^{n}f_i imes g_{n-i} ]

或者这么表示。

[y_n=sum _{i=- infty}^{infty}f_i imes g_{n-i} ]

此时 所构成的多项式({y}) (=) (f imes g)
所以譬如我们都知道的十进制竖式乘法就是一种卷积。


引入概念:点值表达式

一个多项式可以表达成 (f=sumlimits_{i=0}^{n}a_i imes x^i),我们称为其为系数表达式
而此时的多项式我们可以理解成一个 (n) 次的函数
而在初中二年级我们就知道,任意 (n+1) 个点都可以确定一个 (n) 次的函数。
这就是为什么小学一些让你找规律的题其实并没有准确答案。
记点集为(egin{Bmatrix}(x_0,y_0),(x_1,y_1)cdots(x_n,y_n)end{Bmatrix})
(forall iin[0,n] f(x_i)=y_i)
由此我们将系数表达 (egin{Bmatrix} a_0,a_1cdots a_nend{Bmatrix}) 转换到点值表达 (egin{Bmatrix}(x_0,y_0),(x_1,y_1)cdots(x_n,y_n)end{Bmatrix}) 的操作称为一次傅里叶变换
而从点值重新转回系数表达的插值操作称为一次逆傅里叶变换


点值表达式的好处就在于,当两个点值表达式相乘时:

[f=egin{Bmatrix}(x_0,y_0),(x_1,y_1)cdots(x_n,y_n)end{Bmatrix} ]

[g=egin{Bmatrix}(x_0,z_0),(x_1,z_1)cdots(x_n,z_n)end{Bmatrix} ]

[h=f imes g ]

[h=egin{Bmatrix}(x_0,y_0),(x_1,y_1)cdots(x_n,y_n)end{Bmatrix} imes egin{Bmatrix}(x_0,z_0),(x_1,z_1)cdots(x_n,z_n)end{Bmatrix} ]

[forall iin[0,n], h(x_i)=f(x_i) imes g(x_i) ]

( h(x_i)=y_i imes z_i),(h) 的点值表达则为 (egin{Bmatrix}(x_0,y_0 imes z_0),(x_1,y_1 imes z_1)cdots(x_n,y_n imes z_n)end{Bmatrix})
我们实现了 (O(n)) 的离散卷积乘法。


但是我们还不知道如何将系数表达式转换为点值表达式
一般的我们可以,随便找 (n+1) 个数然后代入计算。
(O(2n)) 的复杂度我们并不允许。
我们突然意识到,随便取 (n+1) 个数是不可能的,这辈子都不可能。
这里我们引入概念,(n) 次单位根的概念:
满足 (x^n=1) 的所有的复数 (x)

引入复数概念

没啥特别的,高中基础吧。

[z=a+bi , (ain Re, bin Re, i^2=-1 ) ]

其中 记 (Re(z) = a) 表示复数 (z)的实部。
(Im(z)=b)表示复数 (z) 的虚部。
(mid z mid=sqrt{a^2+b^2}) 表示复数 (z) 的模长。
复数加法法则:实部相加,虚部相加。

[Eg:(1+3i)+(2-i)=(3+2i) ]

而加法的实质意义是平行四边形法则。
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-1000.jpg %}
复数乘法法则:相当于拆括号。

[Eg:(1+3i) imes(2+2i)=1 imes 2+3i imes 2+1 imes 2i+3i imes 2i=2-6+6i+2i=-4+8i ]

而乘法的实质意义是旋转相似。
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-1001.jpg %}
每一个复数 (z)((Re(z),Im(z))) 在平面上表示出。
发现对于任意的一个复数都可以在平面上表示。
于是乎,我们将这个平面称为复平面。
复数相乘几何意义:模长相乘,副角相加。
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-1002.jpg %}

引入单位根概念

2次单位根:(1),(-1)
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-0800.jpg %}
3次单位根:(1),(omega),(omega ^{2})
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-0801.jpg %}
4次单位根:(1),(i),(-1),(,-i)
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-0802.jpg %}
我们发现 (2^k)这种单位根有很好的轴对称以及中心对称性,为了简化问题,下文所说的 (n) 均默认 (n=2^k,kin Z)
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-0803.jpg %}
我们沿逆时针一圈把单位根标序号,对于第 (i) 个我们记录为 (omega_{n}^{i-1}).

于是我们生成了

(omega_{n}^{0}),(omega_{n}^{1}cdots omega_{n}^{n-1})

(omega_{n}^{k}) 性质何在?
首先我们这么理解为 (omega_{n})(k) 次幂
由欧拉定理:

[e^{ivarTheta}=cos {varTheta} + isin{varTheta} ]

其中 (varTheta) 为复数的辐角。
根据其定义,有以下性质。

[omega^{k}_{n}=e^{ ifrac{2pi k}{n}}=cos {frac{2pi k}{n}}+isin{frac{2pi k}{n}} ]

  1. [omega_{2n}^{2k}=omega_{n}^{k} ]

  2. [omega_{n}^{k}=-omega_{n}^{k+frac{n}{2}} ]

我们已经确定了选那些数了,剩下的只剩带入求出该点的函数值!

快速傅里叶变换

我们如何快速求出他们 (omega_{n}^{0}) , (omega_{n}^{1}cdots omega_{n}^{n-1}) 的函数值?

[f(x)=a_0+a_1x+a_2x^2+cdots +a_{n-1}x^{n-1}, (x=omega_{n}^{k}) ]

按次数奇偶划分。

[f(x)=(a_0+a_2x^2+a_4x^4+cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+a_5x^5+cdots+a_{n-1}x^{n-1}), (x=omega_{n}^{k}) ]

划分成两个多项式。

[f(x)=f_1(x^2)+xf_2(x^2), (x=omega_{n}^{k}) ]

发现:

[f(omega_{n}^{k})=f_1(omega_{n}^{2k})+omega_{n}^{k}f_2(omega_{n}^{2k})\ 而f(omega_{n}^{k+frac{n}{2}})=f_1((-omega_{n}^{k})^2)-omega_{n}^{k}f_2((-omega_{n}^{k})^2)\Leftrightarrow\ f(omega_{n}^{k+frac{n}{2}})=f_1(omega_{n}^{2k})-omega_{n}^{k}f_2(omega_{n}^{2k}) ]

发现什么了吗? 在 (omega_{n}^{k})(omega_{n}^{k+frac{n}{2}}) 的函数值,只差个负号。

[所以我们可以只求出 [0,frac{n}{2})的函数值就可以了。 ]

而我们每次这么操作时所需求解的区间就会折半,分治的思想。
所以快速傅里叶复杂度保证在 (O(nlog n)) 了。

逆·快速傅里叶变换

我们引进范德蒙德矩阵

[left [ egin{matrix} 1 & x_0 & x_0^2 & cdots & x_0^{n-1} \ 1 & x_1 & x_1^2 & cdots & x_1^{n-1} \ vdots & vdots & vdots & ddots & vdots\ 1 & x_{n-1} & x_{n-1}^2 & cdots & x^{n-1}_{n-1} end{matrix} ight] left [ egin{matrix} a_0 \ a_1 \ vdots\ a_{n-1} end{matrix} ight] = left [ egin{matrix} y_0 \ y_1 \ vdots\ y_{n-1} end{matrix} ight] ]

这里的 ((x_i,y_i)) 即为点值表达。
(a_i) 为系数表达。
我们记录

[ V=left [ egin{matrix} 1 & x_0 & x_0^2 & cdots & x_0^{n-1} \ 1 & x_1 & x_1^2 & cdots & x_1^{n-1} \ vdots & vdots & vdots & ddots & vdots\ 1 & x_{n-1} & x_{n-1}^2 & cdots & x^{n-1}_{n-1} end{matrix} ight]]

在得知 ((x_i,y_i)) 后我们只要算出 (V) 的逆矩阵 (V^{-1})

[left [ egin{matrix} a_0 \ a_1 \ vdots\ a_{n-1} end{matrix} ight] = V^{-1} left [ egin{matrix} y_0 \ y_1 \ vdots\ y_{n-1} end{matrix} ight]]

好了,现在大力求逆矩阵就完了,所以我们开始高斯消元。
ちょっとまって
这尼玛不还是 (O(n^3))?
哦,我们突然发现(真的是发现)好像(V^{-1}_{i , j}=dfrac{omega _{n}^{-ij}}{n})

你丫不是扯淡?

证明

[P=V imes V^{-1} , P_{ i, j}=sum _{k=0}^{n-1}frac{omega_{n}^{ki} imesomega _{n}^{-jk}}{n}=sum _{k=0}^{n-1}frac{omega_{n}^{k(i-j)}}{n} ]

发现(omega_{n}^{k(i-j)})为等比数列
(i eq j)

[原式=dfrac{1-w_{n}^{(i-j)n}}{1-omega_{n}^{i-j}}\ 根据定义,omega_{n}^{(i-j)n}=(omega_{n}^{n})^{(i-j)}=1\原式=0 ]

(i=j) 时,原式 (=dfrac{n imes omega_{n}^{0}}{n}=1)
综上:得出 (P)(n) 阶单位阵。
根据逆矩阵定义

[V^{-1}_{i , j}=dfrac{omega _{n}^{-ij}}{n} ]

这时我们只要在 (FFT) 的过程中记录一个 (flag) 在系数上(×-1) 即可

void fft(Complex* a,int len,int f)
{
    if(len==1) return;
    Complex a0[len/2],a1[len/2];
    for(int i=0;i<len;i++,i++)
    {
        a0[i/2]=a[i];
        a1[i/2]=a[i+1];
    }
    fft(a0,len/2,f); fft(a1,len/2,f);
    Complex wn(cos(f*2.0*pie/len),sin(f*2.0*pie/len));
    Complex w(1,0);
    for(int i=0;i<(len/2);i++)
    {
        a[i]=a0[i]+w*a1[i];
        a[i+len/2]=a0[i]-w*a1[i];
        w=w*wn;
    }
}

还没完!
(FFT)还没完,我们还要减小常数

1.小trick

w*a1[i] 被算了两遍,没意义,所以

for(int i=0;i<(len/2);i++)
{
    Complex t=w*a1[i];
    a[i]=a0[i]+t;
    a[i+len/2]=a0[i]-t;
    w=w*wn;
}

说实话,仅仅这么改屁用没有。

2.我不用递归啦

观察最后得到的数列:
{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-1003.jpg %}

{% fb_img https://cdn.jsdelivr.net/gh/protons-z/cdn@6.1.3/images/Sol/2020-04-1004.jpg %}
观察化成二进制后正好反转了。
最后的数列我们得到了。
得到递推式:

for(int i=0;i<=lim;i++) 
    rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));

然后从底层向上迭代。
具体如下:

const double p=acos(-1.0);
inline void fft(Complex *a,int f)
{
	for(int i=0;i<lim;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int mid=1;mid<lim;mid=mid<<1)
	{
		Complex wn(cos(1.0*p/mid),f*sin(1.0*p/mid));
		for(int r=mid<<1,j=0;j<lim;j+=r)
		{
			Complex w(1.0,0.0);
			for(int k=0;k<mid;k++,w=w*wn)
			{
				Complex x=a[k+j],y=w*a[k+j+mid];
				a[k+j]=x+y;
				a[k+j+mid]=x-y;
			}
		}
	}
}

(FFT) 到此结束了。

(NTT)

由于(FFT) 炸精还慢,在模数特别时,(NTT)诞生了。

(mod m) 意义下,当 ((m,a)=1) 时,使得 (a^requiv 1 pmod m) 的最小的 (r) ,叫做 (a) 关于 (m) 的阶。记为 (delta_{n}(a)=r)
性质:

  1. ((a,m)=1),且 (a^n=1pmod m),则 (delta_{m}(a)mid n)
  2. (delta_{m}(a)mid varphi(m))

原根

(delta_{m}(a)=varphi(m)) 则称 (a)(mod m)的一个原根。
性质:

  1. 原根存在 (Leftrightarrow) (m=2,4,p^e,2p^e)
  2. (m)(varphi(varphi(m))) 个原根。
  3. (g)(m) 的一个原根,那么所有的原根为(g,g^2,g^3, cdots ,g^{varphi(i)})

发现有些性质和单位根相同,所以我们就使用原根代替单位根,实现快速数论卷积。
实现与 (FFT) 类似。


(MTT)

我们的 (NTT) 十分依赖模数,这使 (NTT) 很鸡肋。

理论基础(雾(

一般认为,两个FFT跑得比三个NTT稍微快一点。

  1. 有一种方法是将非 (NTT) 模数拆成 (3) 个原根模数,最后 (Crt) 合并。但是这是九次 (NTT) 慢死你。。。
  2. 根据command_block大佬的博客,这里提供一种五次 (FFT) 的思路。
    首先是拆系数我们将 (f) 拆成 (f=a_1 imes d + b_1,d=2^{15})(g=a_2 imes d + b_2,d=2^{15})
    (f imes g=a_1 imes a_2 imes d^2+(a_1 imes b_2+a_2 imes b_1) imes d+b_1 imes b_2)
    emmm 四次 (DFT),三次 (IDFT),相当于 (7)(FFT) 慢死了。。
    ちょっとまって
    我们看:这种结构让我们想到了虚数乘法时的样子。
    我们设 (A_1) 代表 (a_1),(B_1) 代表 (b_1),(A_2) 代表 (a_2),(B_2) 代表 (b_2)
    设复多项式 (F=A_1+iA_2) , (G=B_1+iB_2)
    (P_1=F imes G=(A_1 imes B_1-A_2 imes B_2)+i(A_1 imes B_2+A_2 imes B_1))
    再设 (H=A_1-iA_2)
    (P_2=H imes G=(A_1 imes B_1+A_2 imes B_2)+i(A_1 imes B_2-A_2 imes B_1))
    那么 (P_1+P_2=2(A_1 imes B_1+iA_1 imes B_2))
    (A_1 imes B_1=Re(P_1+P_2)) , (A_1 imes B_2=Im(P_1+P_2))
    所以现在可以将(A_1 imes A_2,A_1 imes B_2,A_2 imes B_1,B_1 imes B_2)
    所以将 (F,G,H) 转为点值表达花费 (3)(FFT)
    (P_1,P_2) 转回系数表达花费 (2)(FFT)
    总记 (5)(FFT)
    (5FFT<9NTT-FFT) 优化胜利啦{% emoji_coda 2233/daxiao.png %}

(FWT)

待更。。。

原文地址:https://www.cnblogs.com/Phoenix41/p/14076165.html