多项式
FFT
需要用到复数的内容，然而我并不会，所以写一点。
然后现在就会了对吧。
思考
考虑我们如果要算两个多项式 \(f,g\) 的乘法，可以怎么做。
首先可以思考暴力的想法，\(\mathcal O(n\times m)\) 是显然的。但是似乎没有太大的优化前途。
于是我们不难想到从答案进行入手，既然答案是一个 \(N=n+m\) 次的多项式，那我直接把这个多项式给插值出来不久行了吗？
看起来很对，但实际上发现变成了 \(\mathcal O(N^2)\)，变慢了？？？
关键性质
但是根据拉格朗日插值，我们可以任意选 \(N\) 个值带进去！！！选 \(1\sim N\)，好像没什么用。
这时候，就需要惊人的想象力，我们可以选择 \(N\) 次单位根的若干次方！
做法
多项式求值
本质就是求一个向量和矩阵的乘法：
设 \(\alpha\) 为 \(N\) 次单位根，那么不难发现矩阵里面的元素就只有 \(N\) 种了，但是似乎还是很难算。
于是我们又需要一些性质：
当 \(N\) 为偶数时，\(\alpha^k=-\alpha^{k+\frac{N}{2}}\)。
这似乎意味着我们可以将 \(f(\alpha^k),f(\alpha^{k+\frac{N}{2}})\) 一起算，只需要算一遍，那每次不就只需要算 \(\frac{N}{2}\) 个值了吗？
但是如果我们又直接暴力计算 \(f(\alpha^k)\)，那就没什么区别了。
但如果将 \(f\) 的奇次项和偶次项分开看：
尝试将 $\alpha^k $ 带入：
同理
注意到 \(N\) 次单位根下的 \(\alpha_{N}^{2k}\) 其实就是 \(\frac{N}{2}\) 次单位根下的 \(\alpha_{\frac{N}{2}}^k\)，那这就变成了一个子问题！
注意到 \(f'(x)\) 和 \(f''(x)\) 的形式与 \(f(x)\) 一样，只是次数减半，那么就可以递归进行处理。
又注意到
那么两个东西就可以一起算出来了！
还原多项式
但是现在我们只会算单位根的值，插值求多项式系数似乎还是 \(\mathcal O(N^2)\) 的。注意到求值本值是乘了一个矩阵，那还原就用求的值乘上逆矩阵就行了！！！
而
的逆矩阵一看就很特殊，我们惊人的发现他的逆矩阵就是：
每一项 \(\alpha^k\) 都变成了它的共轭 \(\alpha^{N-k}\)！
那就再带回去求一遍值就行了。
复杂度 \(\mathcal O(n\times \log n)\)。
实现
因为要用到浮点数，常数较大，所以递归版特别慢。而思考它的本质后，就可以研究出非递归版的写法。
首先，FFT 先递归将每一层的奇次项和偶次项分开，那么观察递归到出口时，实际上就是将当前二进制进行翻转后的值赋给他。因此可以直接一层循环实现。
然后，再是合并，我们可以使用“蝶形运算优化”，即将 \(\alpha^k\) 和 \(\alpha^{k+\frac{N}{2}}\) 的点值原地进行计算，比较好理解。
代码如下：
#include<bits/stdc++.h>
using namespace std;
const int maxn=5e6+5;
const double pi2=2.0*acos(-1);
int n,m;
struct cp{
double a,b;
cp operator + (const cp &x){
return {a+x.a,b+x.b};
}
cp operator - (const cp &x){
return {a-x.a,b-x.b};
}
cp operator * (const cp &x){
return {a*x.a-b*x.b,a*x.b+b*x.a};
}
}w[maxn],f[maxn],g[maxn];
void fft(cp a[],int n,int ty){
for(int i=0,j=0;i<n;i++){
if(i<j){
swap(a[i],a[j]);
}
for(int k=n/2;j^=k,j<k;k>>=1);
}
w[0]={1,0};
for(int p=1;p<n;p<<=1){
double c=1.0/(p*2);
cp now={cos(pi2*c),ty*sin(pi2*c)};
for(int j=p*2-2;j>=0;j-=2){
w[j]=w[j>>1];
w[j+1]=w[j]*now;
}
for(int i=0;i<n;i+=(p<<1)){
for(int j=i;j<i+p;j++){
cp x=a[j],y=a[j+p]*w[j-i];
a[j]=x+y,a[j+p]=x-y;
}
}
}
if(ty==-1){
double c=1.0/n;
for(int i=0;i<n;i++){
a[i].a*=c,a[i].b*=c;
}
}
}
int main(){
scanf("%d%d",&n,&m);
for(int i=0;i<=n;i++){int x;scanf("%d",&x);f[i].a=x;}
for(int i=0;i<=m;i++){int x;scanf("%d",&x);g[i].a=x;}
int len=n+m+1;
if((1<<__lg(len))!=len){
len=(1<<(__lg(len)+1));
}
fft(f,len,1);
fft(g,len,1);
for(int i=0;i<len;i++){
f[i]=f[i]*g[i];
}
fft(f,len,-1);
for(int i=0;i<=n+m;i++)printf("%d ",int(f[i].a+0.5));
return 0;
}
NTT
运用
对多项式乘法后的每个系数取模。
一般来说对于质数取模。
思路
注意到原根与单位根很类似。我们可以将模 \(m\) 意义下的 \(N\) 次单位根 \(\omega^N\) 表示为原根 \(g^{\frac{m-1}{N}}\)，然后性质这些也类似！
注意到 \(\varphi(m)=m-1\) 每一次被除 \(2\)，所以 NTT 模数一般会有很多 \(2\) 作为 \(m-1\) 的因子，比如 \(998244353=119\times 2^{23}+1\)。
代码
代码和 FFT 差不多，不过卡一下常数总是好的。
int n,a[N],ff[N],f[N],g[N];
ll jc[N],ny[N],w[N];
int poww(int x,int y,int p){
int res=1;
while(y>0){
if(y&1)res=(ll)res*x%p;
x=(ll)x*x%p;y>>=1;
}
return res;
}
int __mod(int x){
return x>=mod?x-mod:x;
}
void ntt(int a[],int n,int ty){
for(int i=0,j=0;i<n;i++){
if(i<j)swap(a[i],a[j]);
for(int k=n>>1;(j^=k)<k;k>>=1);
}
w[0]=1;
for(int p=1;p<n;p<<=1){
int t=p<<1;
int now=poww(ty==1?3:332748118,(mod-1)/t,mod);
for(int i=t-2;i>=0;i-=2){
w[i]=w[i>>1];w[i+1]=w[i]*now%mod;
}
for(int i=0;i<n;i+=t){
for(int j=i;j<i+p;j++){
int x=a[j],y=w[j-i]*a[j+p]%mod;
a[j]=__mod(x+y),a[j+p]=__mod(x-y+mod);
}
}
}
if(ty==-1){
int c=poww(n,mod-2,mod);
for(int i=0;i<n;i++){
a[i]=(ll)a[i]*c%mod;
}
}
}
任意模数多项式乘法
一般来说值域范围比较大，FFT 会爆精度，模数不是 NTT 模数。这时候就需要用到一些技巧。
三模NTT
可以用 9 遍 NTT，即任取三个 \(10^9\) 级别的 NTT 模数，然后分别算出来答案，再用 CRT 还原回真实值后再对任意模数取模。
优点是比较好想，缺点是常数比较大。
拆系数FFT
8 FFT
我们将 \(f,g\) 的每一项系数 \(a_i\) 拆成 \(x\times 2^{15}+y\)，然后再对应进行相乘。
即将原多项式拆成两部分 \(f_1,f_2\)，其中 \(f=f_1\times 2^{15}+f_2\)，\(g\) 也同理。
那么 \(f*g=f_2*g_2+(f_1*g_2+f_2*g_1)\times 2^{15}+f_1*g_1\times 4^{15}\)。
而这个做法需要 8 遍 FFT，肯定比不过 9 遍 NTT。
5 FFT
注意到 FFT 引入了复数，那么我们不妨将 \(f_1,f_2\) 的结果统一放在一个多项式的实部和虚部，即引入 \(f'=f_1+f_2I\) 那么我们惊奇的发现 \(f'*g_1\) 和 \(f'*g_2\) 的实部和虚部就是我们想要的所有东西。这个做法需要 5 FFT，已经很快了。
还有 4 FFT 的做法，暂时还没搞，到时候再补。
FWT
xor/or/and 卷积
定义 \(\otimes\in{\operatorname{or},\operatorname{xor},\operatorname{and}}\)，求
根据 FFT，依然考虑转换成点乘。即构造：
最后再逆变换回去。
几种操作的变换方式分别为：
or 卷积
for(int p=1;p<n;p<<=1){
int c=p<<1;
for(int i=0;i<n;i+=c){
for(int j=i;j<i+p;j++){
a[j+p]=ty==1?__mod(a[j+p]+a[j]):__mod(a[j+p]-a[j]+mod);
}
}
}
and 卷积
for(int p=1;p<n;p<<=1){
int c=p<<1;
for(int i=0;i<n;i+=c){
for(int j=i;j<i+p;j++){
a[j]=ty==1?__mod(a[j]+a[j+p]):__mod(a[j]-a[j+p]+mod);
}
}
}
xor 卷积
for(int p=1;p<n;p<<=1){
int c=p<<1;
for(int i=0;i<n;i+=c){
for(int j=i;j<i+p;j++){
int x=a[j],y=a[j+p];
a[j]=__mod(x+y),a[j+p]=__mod(x-y+mod);
}
}
}
if(ty==-1){
ll c=poww(n,mod-2);
for(int i=0;i<n;i++){
a[i]=c*a[i]%mod;
}
}
子集卷积
实际上可以转换为
这样看起来就很像 or 卷积了。但是还要满足 \(i\cap j=\emptyset\)，再进行变形即可发现：
那么将 or 卷积扩域，再记录一维 \(i\) 表示 \(|T|\)。每次单独对 \(|T|\) 相同的数跑 or 卷积，最后再合并。
for(int i=0;i<siz;i++){
cin>>a[__builtin_popcount(i)][i];
}
for(int i=0;i<siz;i++){
cin>>b[__builtin_popcount(i)][i];
}
for(int i=0;i<=n;i++){
fwt(a[i],siz,1);fwt(b[i],siz,1);
}
for(int j=0;j<=n;j++){
for(int k=0;k<=j;k++){
for(int i=0;i<siz;i++){
c[j][i]=__mod(c[j][i]+(long long)a[k][i]*b[j-k][i]%mod);
}
}
}
for(int i=0;i<=n;i++){
fwt(c[i],siz,-1);
}
全家桶
这里的多项式乘法大多都是 NTT。
多项式求逆
一般是对于多项式 \(F(x)\)，求一个多项式 \(G(x)\)，满足 \(F(x)*G(x) \equiv 1 \mod x^n\)，系数对 NTT 模数 \(p\) 取模。
暴力
发现可以暴力地从小到大递推来做，但是复杂度是 \(\mathcal O(nm)\) 的。但是可以从暴力做法得出 \(F(x)\) 有逆的充要条件为 \(F(x)\) 的常数项在 \(\mod p\) 意义下有逆。
倍增
可以使用倍增来解决这个问题。
首先，让 \(n=2^k\)，然后考虑倍增。
假设现在算出了 \(F(x)*G_m(x)\equiv 1 \mod x^m\)，考虑如何计算 \(F(x)*G_{2m}(x)\equiv 1 \mod x^{2m}\)。
需要进行一些推导。
首先，需要将 \(G_m(x),G_{2m}(x)\) 放在一起，有：
想让模数变为 \(x^{2m}\)，不难想到将同余式两边平方。
\(F(x)*G_{2m}(x)\equiv 1 \mod x^{2m}\) 这个条件还没用，将 \(F(x)\) 乘进去，神奇的事情发生了！
右边都是已知的多项式，直接 NTT 就行了，复杂度单次倍增是 \(\mathcal O(m\log m)\) 的，总时间复杂度即为 \(\mathcal O(n\log n)\)。
当然，还可以优化一点常数。
这里要注意一个细节，就是 \(G_m(x)^2*F(x)\) 最高次是 \(m-1+m-1+2m-1=4m-3\) 次，需要带入 \(4m\) 个点值。
优化
以上做法每次需要 3 遍长度为 \(4m\) 的 NTT。但是可以进行优化。
观察一下原式：
最关键就是求 \(G_m(x)^2*F(x)\)，将其分为两步，先计算 \(C(x)=F(x)*G_m(x)\mod x^{2m}\)，再计算 \(D(x)=C(x)*G_m(x)\mod x^{2m}\)。
如何优化？发现 \(C(x)\) 的 \(C[0]=1\)，\(C[1],C[2]...C[m-1]=0\)。这意味着只需要算对 \(C[m]\sim C[2m-1]\) 就行了。
那么可以直接跑长度为 \(2m\) 的 NTT，乘出来的 \(C[2m]\sim C[3m-1]\) 会累加到 \(C[0]\sim C[m-1]\)，但是没关系，因为我们会直接将其赋值。
对于 \(D(x)\) 同理，因为 \(G_{2m}(x)\) 的前 \(m\) 项与 \(G_m(x)\) 相同，所以 \(D(x)\) 的前 \(m\) 位同样可以不管。
总体只需要 5 次长度为 \(2m\) 的 NTT。
多项式除法
给定一个 \(n\) 次多项式 \(F(x)\) 和一个 \(m\) 次多项式 \(G(x)\) ，求出多项式 \(Q(x)\), \(R(x)\)，满足以下条件：
- \(Q(x)\) 次数为 \(n-m\)，\(R(x)\) 次数小于 \(m\)
- \(F(x) = Q(x) * G(x) + R(x)\)
在模一个 NTT 模数 \(p\) 的情况下计算。通俗来讲，就是计算类似于带余除法的商和余数（大除法）。
首先，可以直接大除法 \(\mathcal O((n-m)m)\) 解决。
如何优化？需要求解 \(Q(x),R(x)\) 两个东西，不好同时求解，考虑消掉一个进行求解。
直觉消去 \(R(x)\) 方便一些。
方法
接下来的操作就有些妙了。
定义
即
就是将多项式进行了翻转。接下来的推导比较神奇：
而现在就可以将 \(R_R(x)\) 消掉了，要求的 \(Q_R(x)\) 的最高次项为 \(x^{n-m}\)，所以可以将两边同时模一个 \(x^{n-m+1}\)。
求出 \(Q_R(x)\) 后，将其系数翻转就可得到 \(Q(x)\)，将其带回原式，即可得到：
需要求逆和多项式乘法，时间复杂度 \(\mathcal O(n\log n)\)。
常系数齐次线性递推
求一个递推数列 \(f_i\) 的第 \(n\) 项，满足：
矩阵快速幂
不难想到可以直接使用矩阵快速幂，复杂度 \(\mathcal O(k^3\log n)\)。
考虑如何优化，需要挖掘一些性质。
快速幂+多项式除法
对于任意的 \(f_j\)，我们都可以将其表示为
即将 \(f_j\) 展开后算出系数之和。
而展开的过程不难发现与多项式的降幂很类似。
具体的，将 \(f_i\) 看做 \(x^i\) 的系数。那么展开即可表示为：
而对于原式，将其减去上式即达到了降幂的效果，不难发现，最终结果即原式对于 \(x^k-\sum_{i=1}^k a_ix^{k-i}\) 这个多项式取模。而这个多项式被称为特征多项式。
那么要计算 \(x^n \mod x^k-\sum_{i=1}^k a_ix^{k-i}\)，自然就可以写过快速幂，需要用到多项式除法（取模）。
如果模数 \(p\) 为 NTT 模数，自然可以做到 \(\mathcal O(k\log k\log n)\)。如果不是 NTT 模数，那么取模可以直接使用大除法，多项式乘法也可以暴力，复杂度 \(\mathcal O(k^2\log n)\)。
多项式开根
求 \(G(x)^2\equiv F(x)\mod x^n\)。
暴力
根求逆一样，开根显然也可以递推，首先 \(G[0]^2\equiv F[0]\mod p\)，即 \(F[0]\) 在模 \(p\) 意义下需要有二次剩余。然后就可以递推了。
倍增
发现暴力递推的过程与多项式求逆十分神似，于是同样考虑倍增。
假设当前求得：
要推出 \(G_{2m}\)：
而根据递推的过程 \(G_{2m}(x)\) 的低位与 \(G_m(x)\) 相同，就有：
接下来与推导多项式求逆就一样了：
需要求逆，复杂度与多项式求逆一样，\(\mathcal O(n\log n)\)
多项式多点求值
对于 \(n\) 次多项式 \(F(x)\)，求 \(F(a_1),F(a_2),F(a_3)...F(a_m)\) 的点值。
分治
暴力带点求值肯定不行，那考察一下性质，
注意到 \(F(a_i)\equiv F(x)\mod x-a_i\)，但是一个一个求甚至比暴力还慢，所以需要加速这个过程。
考虑分治，对于区间 \([l,r]\)，当前算好了 \(F'(x)\)，将 \(F'(x)\mod \prod_{i=l}^r (x-a_i)\)，然后再递归左右区间。最后 \(l=r\) 时即为 \(F(a_l)\) 的答案。
\(\prod _{i=l}^r (x-a_i)\) 可以先用分治 NTT 预处理，复杂度 \(\mathcal O(n\log^2 n)\)。多项式取模会取 \(2n\) 次，复杂度 \(\mathcal O(n\log n)\)。
所以总时间复杂度 \(\mathcal O(n\log^2 n)\)。
多项式 ln
恶补极限，导数，积分中......
给定 \(F(x)\)，求 \(G(x)\equiv \ln F(x)\mod x^n\)。
需要满足 \(F[0]=1\) 才可 ln。
设 \(f(x)=\ln x\)，则要求：
分别求导，可得：
然后就可以求出 \(G'(x)\)，最后再积分回去。
多项式求导即：
积分即：
所以：
多项式 exp
给定 \(A(x)\)，求 \(B(x)\equiv e^{A(x)}\mod x^n\)。同样 \(A[0]=0\) 为 \(B(x)\) 存在的充要条件。
发现十分神奇，\(e\) 不是超越数吗？怎么 \(A(x)\) 次方之后变成多项式了？还有多项式次方是什么？
但是这些都不重要。
首先，带有 \(e\) 肯定不好办，先取一个对数。
牛顿迭代
然后就不会了。但是先别急，需要学习牛顿迭代。
并不需要知到牛顿迭代如何证明，为什么可以用，直接用公式就行了。
先确定 \(G(x)=\ln x-A(x)\)。其中 \(x\) 表示一个多项式，\(A(x)\) 看做常数。
那么套用公式：
多项式快速幂
给定 \(A(x),k\)，求 \(B(x)\equiv A(x)^k\mod x^n\)。\(k\) 是一个𡘙的数。
二进制快速幂
显然可以直接使用普通的快速幂，时间复杂度 \(O(n\log n\log k)\)，但如果 \(k\) 特别大，仍然不够优秀。
使用 ln 和 exp
注意到
因此，只需要使用多项式 ln 和 exp 即可。
但是有一个问题。\(\ln A(x)\mod x^n\) 存在的条件是 \(A[0]=1\)，但是 \(A(x)\) 的常数项显然不一定为 \(1\)。
所以要进行一些调整。具体而言，找到一个 \(C(x)\)，使得 \(C[0]=1\)，且：
容易发现 \(q\) 即为 \(A(x)\) 最低的一位，使得 \(A[q]\neq 0\)。\(p=A[q]\)。那么：
这样就行了。注意 \(p^k\) 中的 \(k\) 要模 \(\varphi(p)\)。而 \(qk\ge n\) 时，\(B(x)\) 即为 \(0\)。
快速插值
分治
即快速拉格朗日插值。可做到 \(\mathcal O(n\log^2 n)\)。
首先回顾拉格朗日插值的式子：
考虑如何快速计算。
首先，注意到 \(\prod_{j \neq i}x-x_j\) 实际上比较好算，虽然有下标不相等的限制，但同样可以用分治 NTT 求出。而分母就不是很好看，所以提到外面当做系数：
设
那么
\(G_i(x_i)\) 看起来很像多项式多点求值，但是每个 \(x_i\) 带入的多项式长得不一样。需要想办法让他们一样。
注意到对于 \(x_i\)，\(\forall j\neq i,G_j(x_i)=0\)，所以可以将他们加起来，\(G_i(x_i)\) 的值也不会发生变化。
而 \(G(x)\) 同样是可以分治 NTT 求解的，还需要多项式多点求值出 \(G(x_1),G(x_2)...G(x_{n})\)。
原式变为：
注意到 \(F(x)\) 同样可以分治 NTT 求解，所以总复杂度为 \(\mathcal O(n\log^2 n)\)。
优化
上面做法带有 \(4\) 倍多的常数，比较慢，本质上在算 \(G(x)\) 时会做两次乘法合并。但是根据求导的乘法法则，容易发现：
可以优化 \(\frac{1}{4}\) 左右的常数。
代码
先贴一个求逆，除法，开根的代码，接下来还要写一堆东
西，准备好接受审判了。
取模（多项式除法）有时候除数多项式是一样的，所以可以提前算好 \(G_R(x)^{-1}\)，多算几位，可以减少常数。当然，如果次数很少或者每次不同，先算 \(G_R(x)^{-1}\) 需要的位数肯定是快一些的。
因为多项式多点求值需要存 \(O(n)\) 个多项式，所以把代码改成了 vector，记得 resize，然后就和数组一样用。
浙公网安备 33010602011771号