多项式的表示
系数表示法:即 (f(x)=sum limits_{i=0}^n a_i imes x^i) 的形式。
点值表示法:即给定 (n+1) 个互不相同的点的坐标来确定一个多项式。
复数
定义一个数 (i=sqrt{-1}) ,则所有形如 (a+bi)((a,b in mathbf{R}))的数称为复数,以复数构成的集合即为复数集 (mathbf{C}) 。
对于一个这样的数,称 (a) 为它的实部, (b) 为它的虚部。
复平面为一个平面直角坐标系,横轴为实轴,纵轴为虚轴。
每个复数都在复平面上对应着一个坐标为 ((a,b)) 的点,也对应着一个从原点到 ((a,b)) 的向量。
这个向量与实轴正方向的夹角即称为该复数的辐角。(若向量在实轴下方则辐角大于180度)
复数 (x=a+bi) 的模长定义为它在复平面上到原点的距离,记作 (|x|) ,即 (|x|=sqrt{a^2+b^2}) 。
定义一个复数在复平面上所对应的点关于实轴对称后,所得到的点所对应的复数与它互为共轭复数。
对于一个复数 (x=a+bi) ,它的共轭复数即为 (ar x=a-bi) 。
复数的运算法则
加法:((a+bi)+(c+di)=(a+c)+(b+d)i) 。
减法:((a+bi)-(c+di)=(a-c)+(b-d)i) 。
乘法:
((a+bi) imes (c+di)=ac+bci+adi+bdi^2=(ac-bd)+(bc+ad)i) 。
在复平面上,乘积的辐角等于两个乘数的辐角之和,模长等于两个乘数的模长之积。
除法:(dfrac{a+bi}{c+di}=dfrac{(a+bi)(c-di)}{(c+di)(c-di)}=dfrac{ac+bd}{c^2+d^2}+dfrac{bc-ad}{c^2+d^2}i) 。
欧拉公式:(e^{i heta}=cos heta +i sin heta) 。
单位根
在复数意义下,满足 (x^n=1) 的 (x) 称为 (n) 次单位根。
定义 (n) 次单位根 (omega_n=e^{frac{2 pi i}{n}}) ,则上面方程中 (x) 的解集就可以表示为 ({omega_n^k|k=0,1, cdots n-1}) 。
性质:(omega_n^k=omega_{2n}^{2k}) , (omega_{2n}^{k+n}=-omega_{2n}^k) 。
单位根可以用上面的欧拉公式展开计算。
在几何意义上,单位根就是把单位圆 (n) 等分后第一个在圆周上的数。
FFT
下面以多项式乘法为例。
考虑用点值表示法求多项式乘法的过程,
设一个多项式为 (f(x)) , 另一个多项式为 (g(x)) ,
那么,要找到 (n+1) 个在它们乘积的多项式上的点。
先在第一个多项式上找到 (n+1) 个点: ((x_1,f(x_1)),(x_2,f(x_2)), cdots ,(x_{n+1},f(x_{n+1}))) ,
同样地,在第二个多项式上找到 (n+1) 个点:((x_1,g(x_1)),(x_2,g(x_2)), cdots ,(x_{n+1},g(x_{n+1}))) ,
则可以得到在它们乘积的多项式上的 (n+1) 个点:((x_1,f(x_1) imes g(x_1)),(x_2,f(x_2) imes g(x_2)), cdots ,(x_{n+1},f(x_{n+1}) imes g(x_{n+1}))) 。
但是,这样做的复杂度明显太高。
考虑选点时选 (1,-1) 等幂具有特殊性质的数。
这里,就可以选择单位根。
考虑如何利用单位根的性质来优化上述过程。
设多项式 (f(x)=a_0+a_1 imes x+a_2 imes x^2+a_3 imes x^3 + cdots + a_{n-1} imes x^{n-1}) ,
按奇偶分组,得到 (f(x)=(a_0+a_2 imes x^2+a_4 imes x^4+ cdots +a_{n-2} imes x^{n-2})+(a_1 imes x+a_3 imes x^3+a_5 imes x^5 + cdots +a_{n-1} imes x^{n-1})) ,
设 (g(x)=a_0+a_2 imes x +a_4 imes x^2 + cdots +a_{n-2} imes x^{frac{n}{2}-1}) , (h(x)=a_1+a_3 imes x+a_5 imes x^2+ cdots +a_{n-1} imes x^{frac{n}{2}-1}) ,
那么, (f(x)=g(x^2)+xh(x^2)) 。
用 (omega_n^k) ((k< frac{n}{2}))代入,根据单位根的性质,得 (f(omega_n^k)=g(omega_n^{2k})+omega_n^k h(omega_n^{2k})=g(omega_{frac{n}{2}}^k)+omega_n^kh(omega_{frac{n}{2}}^k)) 。
接下来,用 (omega_n^{k+frac{n}{2}}) 代入,得 (f(omega_n^{k+frac{n}{2}})=g(omega_n^{2k+n})+omega_n^{k+frac{n}{2}}h(omega_n^{2k+n})=g(omega_n^{2k})-omega_n^k h(omega_n^{2k})=g(omega_{frac{n}{2}}^k)-omega_n^kh(omega_{frac{n}{2}}^k)) 。
这两个式子形式相同,可以一起算,然后分别递归求解。
注意要先补齐长度到 (2^m) 项。
但是,递归的常数极大,考虑能不能把递归转为迭代。
在递归时,因为是按照奇偶递归的,所以如果我们能求出递归到最后时的系数数组,就可以迭代实现了。
找一下规律,以一个 8 项多项式为例。
原数列: (a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7) ,
新数列: (a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7) 。
可以发现,新数列的下标就是原数列下标用二进制表示后翻转。
考虑怎么求翻转后的数,这个是可以 dp 求的。
从小到大求,这样在求到每个数时,设它的二进制有 (n) 位,则在此之前前 (n-1) 位翻转后的数一定已经求出来了,那么,只需把最后一位的贡献加上即可。
但是,只是这样还不够。我们还需要一个快速把点值表示转化为系数表示的方法。
设多项式 (f(x)=sum_{i=0}^{n-1} x^{a_i}) 经过前面的变换后得到的点值为 ({b_i}) ,即 (b_k= sum limits_{i=0}^{n-1} a_i imes omega_n^{ik}) ,
现在,我们要根据 (b) 还原出 (a) 。
直接给出结论: (a_k= frac{1}{n} sum limits_{i=0}^{n-1} b_i imes omega_n^{-ik}) 。
所以,取单位根为原来的倒数,直接对 (b) 跑一遍前面的过程,最后再除以 (n) 即可。
而单位根的倒数用前面的欧拉公式展开后形式是和单位根类似的,所以两种变化可以合在一起写。
void FFT(cp *a,int tp,int n)
{
for(int i=0;i<n;i++)
if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int l=1;l<n;l*=2)
{
cp w=(cp){cos(PI/(db)(l)),(db)(tp)*sin(PI/(db)(l))};
//前面枚举的是长度的一半,所以这里不用*2
for(int i=0;i<n;i+=l*2)
{
cp W=(cp){1,0};
for(int k=0;k<l;k++,W=W*w)
{
cp x=a[i+k],y=W*a[i+k+l];
a[i+k]=x+y,a[i+k+l]=x-y;
}
}
}
}
NTT
FFT 中用到了很多实数运算,所以精度不高。
主要的问题就是单位根是个复数,考虑能不能用一个整数代替。
单位根能用于 FFT ,原因就是性质优秀,所以,我们要找的就是满足单位根的一些性质的数。
而在模 (p) 意义下,就可以找到这样的数。
设 (p) 的原根为 (g) ,则 (g^{frac{p-1}{n}} pmod p) 就可以用于替代单位根。
对着单位根的几个性质推一推,可以发现都成立。
void NTT(int *a,int tp,int n)
{
for(int i=0;i<n;i++)
if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int l=1;l<n;l*=2)
{
int w=POW(tp?g:invg,(p-1)/(l*2));
for(int i=0;i<n;i+=l*2)
{
int W=1;
for(int k=0;k<l;k++,W=W*w%p)
{
int x=a[i+k],y=W*a[i+k+l]%p;
a[i+k]=(x+y)%p,a[i+k+l]=(x-y+p)%p;
}
}
}
}