基础操作
FFT & NTT
基础知识。
考虑如何快速计算多项式乘法,发现常规的做法是没有前途的,考虑换种方法计算。
有一种计算方法是这样的(设要计算 (A imes B) ):
- 将 (A) 转化为点值表示。
- 将 (B) 转化为点值表示。
- (O(n)) 地计算点值表示的乘积 (C)。
- 将 (C) 转化为数值表示。
复杂度瓶颈为点值表示与数值表示的转化,朴素做是 (O(n^2))。
为了加速这个过程,FFT 引入了单位复根 (omega),NTT 引入了原根 (g)。
单位复根 (omega) 是 (x^n=1) 的解集,其满足如下性质:
- (omega^n_n=1);
- (omega^k_n=omega^{2k}_{2n});
- (omega^{k+n}_{2n}=-omega^k_{2n})。
考虑分治,对一个多项式 (f(x)) 奇偶性分组处理,比如一个多项式:
考虑奇偶性分组后的状态:
建立新函数:
那么待求解的 (f(x)) 可改为:
带入 (omega) 可得:
然后就可以分治做了。
诶等等,递归常数特别大呢!
为了不递归,我们考虑一开始就预处理出来每一项最后到了什么地方,这个做法是位逆序置换(又名蝴蝶置换)。
考虑找规律,发现将下标二进制表示后翻转即可,设 (R_i) 为 (i) 变到的地方,那么可以考虑这样求:
- 右移一位。
- 翻转。
- 再右移一位。
- 处理最高位。
那么 (R_i) 的递推式可以轻易的写成:
以上就是如何将数值表示转化为点值表示了,所以怎样转回来呢?
直接反转后除以 (len),或是取单位根的倒数跑即可,具体推导过程见 OI-Wiki。
模板:
void FFT(int n,cp* x,int op) {
for(int i=0;i<n;i++)
if(i<rk[i]) swap(x[i],x[rk[i]]);
for(int mid=1;mid<n;mid*=2) {
cp W(cos(Pi/mid),op*sin(Pi/mid));
for(int R=mid*2,j=0;j<n;j+=R) {
cp w(1,0);
for(int k=0;k<mid;k++,w=w*W) {
cp z1=x[j+k],z2=w*x[j+mid+k];
x[j+k]=z1+z2,x[j+mid+k]=z1-z2;
}
}
}
}
NTT 呢?
首先解释 (g),(g) 是原根,在 (mod p) 的情况下 (forall i ot =j,g^i otequiv g^jpmod p)。
我们用 (g^{frac{varphi(p)}{n}}) 代替 (omega) 进行计算,由 (varphi(p)=p-1) 且 (n=2^q),所以强制要求 (p=k imes 2^q+1),其中 (2^qge n)。
较 FFT 精度高,但却对模数有着致命的要求。
时间复杂度 (O(nlog n))。
分治 FFT & NTT
问题:求 (f_n),其中:
考虑 CDQ 分治,对于一个区间 ([l,r]),将其划分成 ([l,mid],[mid+1,r]) 两段。
先处理 ([l,mid]),再计算 ([l,mid]) 对 ([mid+1,r]) 的贡献,再处理 ([mid+1,r])。
中间的计算使用 NTT/FFT 解决,时间复杂度 (O(nlog^2n))
拉格朗日插值
问题,给定 (n) 个点 ((x_i,y_i)),确定一个多项式 (y=f(x)),求出 (f(k)) 的值。
显然的事实:
证明:(f(x)-f(y)=(a_0-a_0)+color{red}a_1(x-y)color{black}+a_2(x^2-y^2)+a_3(x^3-y^3)+cdots)。
那么我们就可以列一个线性同余方程组:
直接 CRT,所有模数的积为:
那么 (m_i=prod _{j ot=i} (x-y_i),m_i^{-1}=prod _{j ot=i} frac {1}{x_i-x_j})。
所以:
朴素实现复杂度是 (O(n^2)) 的。
在所取点值连续的时候,即 (x_i=i) 的时候,可以直接用前缀积优化到 (O(n))。
在需要动态加点的时候,可以使用重心拉格朗日插值法:
回顾之前的式子:
设 (g=prod^n_{i=1} k-x_i),那么:
设 (t_i=frac {y_i}{prod_{j ot=i} x_i-x_j}),那么:
加入一个点只需计算 (t_i) 即可。
求自然数幂和可以直接拉格朗日插值,时间复杂度 (O(n))。
多项式求导与积分
由于求导和积分满足可加性,所以拆成单项式分别求导或积分。
如果你学过微积分相关理论,会知道:
直接做即可,为数不多的 (O(n)) 时间复杂度。
多项式求逆
给定 (G(x)),求一个 (F(x)) 使得 (F(x)G(x)equiv0pmod {x^n})。
首先假设我们已经求出了一个 (f(x)) 使得 (f(x)G(x)equiv 0pmod{x^{n/2}}),考虑如何求出 (F(x))。
因为 (F(x)G(x)equiv 0pmod{x^n}),所以必然 (F(x)G(x)equiv 0 pmod{x^{n/2}}),两式相减,得 (F(x)-f(x)equiv0pmod{x^{n/2}}),此处 (G(x)) 被约掉了。
两边同时平方,得:
注意标红的部分,想一想,为什么?
考虑卷积的定义,在 ([0,n/2]) 的部分每一项都是 (0),而右边 ([n/2,n]) 的部分由 (a_i=sum^i_{j=1} a_ja_{i-j}) 得来,容易发现,(j) 与 (i-j) 中必有一个小于零。
两边同时乘以 (G(x)) 并移项,得:
直接使用 FFT 计算即可,时间复杂度 (O(nlog n))。
多项式开根
给定 (g(x)),求 (f(x)),使得:(f^2(x)equiv g(x)pmod{x^n})。
类似求逆,假设已经求出了 (f_0^2(x)equiv g(x)pmod{x^{n/2}})。
移项得:
两边平方得:
使用简单的数学公式 ((x-y)^2=(x+y)^2-4xy) 得:
将 (4f_0^2(x)) 移到左边,得:
左右两边同时开根,得:
直接计算即可,时间复杂度 (O(nlog n))。
多项式带余除法
问题:给定 (n) 次多项式 (F(x)) 和 (m) 次多项式 (G(x)),求 (Q(x)) 与 (R(x)),满足:
- (deg Q(x)=n-m),(deg R(x)<m);
- (F(x)=Q(x) imes G(x)+R(x))。
对于一个多项式 (f(x)),定义 (f_R(x)=x^nf(frac 1 x)),等价于将多项式 (f) 的系数翻转一次。
那么:
直接求出 (Q(x)),那么回代可得 (R(x)),时间复杂度 (O(nlog n))。
多项式对数函数
求 (B(x)equiv ln A(x)pmod{x^n})。
设 (F(x)=ln x),那么:
直接模拟即可,时间复杂度 (O(nlog n))。
多项式指数函数
求 (F(x)equiv e^{A(x) pmod{x^n}})。
考虑多项式牛顿迭代。
多项式牛顿迭代解决的问题是,给定函数 (G(x)),求 (G(F(x))equiv 0pmod {x^n})。
既然是迭代,那么当 (n=1) 的时候是要单独求出的,这里不讲。
类似求逆,假设已经求出 (G(F_0(x))equiv 0 pmod{x^{lceil frac n x ceil}}),考虑如何扩展。
把 (G(F(x))) 在 (F_0(x)) 处泰勒展开:
因为 (F(x)) 和 (F_0(x)) 的最后几项相同,所以 ((F(x)-F_0(x))^2) 的最低非 (0) 次数大于 (n),所以可以得到:
那么,将 (F(x)=e^{A(x)}) 变形得 (ln F(x)-A(x)=0)。
设 (G(F(x))=ln F(x)-A(x)),那么 (G'(F(x))=frac 1 {F(x)}),带入牛顿迭代公式得:
直接做,时间复杂度 (O(nlog n))。
参考代码
不想太长,所以放一份在剪贴板。