多项式学习笔记
\(\LaTeX\) 过多警告 ⚠
学多项式了,由于肯定会忘所以来写笔记了,同时边学边写算加深记忆吧。这边还是用上次网络流笔记的格式。
下文中角度均以弧度制表示。
一、前置芝士知识
\(~~~~\) 1.多项式: 称一个关于 \(x\) 的式子:
\(~~~~\) 为 \(\pmb{n}\) 次多项式,\(a_i\) 为 \(\pmb{i}\) 次项系数(是一个常数),\(n\) 为该多项式的次数。
\(~~~~\) 由上面的定义,我们也可以把 \(n\) 次多项式记为一个 \(n\) 次函数 \(y=f(n)\) 。
\(~~~~\) 而我们知道,给出 \((n+1)\) 个点的坐标,此时可以确定一个 \(n\) 次函数的解析式。(由于一个 \(n+1\) 元的非无解线性方程组有唯一解的必要条件是有 \(n\) 个本质不同的方程)
\(~~~~\) 综上,我们可以用两种方式来记一个多项式:
\(~~~~ ~~~~\) 1.1 系数表示法:写成上面那个式子的形式。
\(~~~~ ~~~~\) 1.2 点值表示法:给出 \((n+1)\) 个互不相同,且能形成非无解的方程组的点的坐标。
\(~~~~\) 2.多项式运算
\(~~~~\) 为了方便,我们假定下面进行运算的两个多项式 \(A(x)=\sum_{i=0}^n a_i\times x^i\),\(B(x)=\sum_{i=0}^nb_i\times x^i\) 的次数均为 \(n\) 。(不足的一个不足部分的系数就是 \(0\) )
\(~~~~ ~~~~\) 2.1 多项式加减法:
\(~~~~ ~~~~\) 乘法分配律,显然。
\(~~~~ ~~~~\) 2.2 多项式乘法:
\(~~~~ ~~~~\) 这个也很显然吧,还是乘法分配律,多用几遍而已。
\(~~~~ ~~~~\) 2.3 点值表示法下运算
\(~~~~ ~~~~\) 首先,取可以表示多项式 \(A\) 的 \((n+1)\) 个点,并取相同横坐标的,能表示多项式 \(B\) 的 \((n+1)\) 个点。
\(~~~~ ~~~~\) 之后,将两个横坐标相等的点的纵坐标直接 加/减/乘 就可以得到结果的多项式的点值表示。(由于多项式的运算也可以看作两个多项式值的运算。)
\(~~~~ ~~~~\) 需要注意的是,由于乘法产生的多项式最高次数为 \(2n\) ,所以初始应该取 \(2n+1\) 个点。
\(~~~~\) 3.复数
\(~~~~\) 定义常数 \(\pmb i\) 为虚根,满足:
\(~~~~\) 则所有形如:
\(~~~~\) 的数 \(z\) 被称作复数,同时所有复数组成的集合为复数集,记作 \(C\).
\(~~~~\) 特殊地,当 \(\pmb{b \not = 0}\) 时,这样的数称为虚数,这样的数不能在数轴上表示,因为找不到它的相对大小。
\(~~~~\) 称复数 \(z=a+bi\) 中,\(a\) 为 \(z\) 的实部,\(b\) 为 \(z\) 的虚部。
\(~~~~\) 4.复平面
\(~~~~\) 一个平面直角坐标系,由互相垂直的横轴(即实轴)与纵轴(即虚轴)组成。
\(~~~~\) 对于一个复数 \(z=a+bi\) ,它可以在复平面上表示为一个从 \((0,0)\) 指向 \((a,b)\) 的向量。
\(~~~~\) 而一个从 \((0,0)\) 指向 \((a,b)\) 的向量也一定唯一对应一个 \(z=a+bi\) 的向量。这个很显然吧。
\(~~~~\) 特殊地,当向量指向的点纵坐标为 \(0\) 时,它表示一个实数。
\(~~~~\) 记以横轴为始边,复数 \(\pmb z\) 对应的向量 \(\pmb Z\) 为终边的有向正角为复数 \(z\) 的幅角。
\(~~~~\) 注:正角为始边逆时针旋转使得其与终边重合的角度 \(\theta\),负角为始边顺时针旋转使得其与终边重合的角度 \(\theta\) 的相反数。
\(~~~~\) 举个例子:
\(~~~~\) 上图的 \(\theta\) 就是幅角。
\(~~~~\) 5.相关定义
\(~~~~\) 以下设复数 \(z=a+bi\) 。
\(~~~~ ~~~~\) 5.1 复数的模: 一个复数的模为它在复平面上对应向量的长度,\(z\) 的模记作 \(|z|\) 。距离公式得:
\(~~~~ ~~~~\) 5.2 共轭复数: 两个复数若实部相等,虚部互为相反数,则它们互为共轭复数。记 \(z\) 的共轭复数为 \(\overline z\) 。
\(~~~~ ~~~~\) 设 \(z\) 的幅角为 \(\theta_1\) ,\(\overline z\) 的幅角为 \(\theta_2\) 。则:
\(~~~~ ~~~~\) 5.3 复平面上共轭复数: 两共轭复数对应的向量在复平面上关于横轴对称。因此不难发现上面关于幅角和的性质:
\(~~~~\) 6 复数加减法
\(~~~~ ~~~~\) 6.1 代数加减 :设 \(z_1=a+bi,z_2=x+yi\) ,则 \(z_1\pm z_2=(a\pm x)+(b \pm y)i\) 。
\(~~~~ ~~~~\) 6.2 复平面上加减 :既然两个复数已经被表示为向量了,那向量加减直接使用平行四边形定则即可。什么,你说你不会? 物理必修一第二章欢迎您
\(~~~~\) 7 复数乘法
\(~~~~ ~~~~\) 7.1 代数乘法: 仍然设 \(z_1=a+bi\),\(z_2=x+yi\) ,那么:
\(~~~~ ~~~~\) 代入
\(~~~~ ~~~~\) 所以 \(z_0\) 的向量可以表示为 \((ax-by,ay+bx)\) 。
\(~~~~ ~~~~\) 如果两个复数互为共轭复数,则 \(z\times \overline z=(a+bi)\times (a-bi)=a^2+b^2\) ,而这肯定是一个实数。
\(~~~~ ~~~~\) 7.2 复平面上乘法: 设 \(z_0=z_1\times z_2\) ,\(z_0,z_1,z_2\) 的幅角分别是 \(\theta_0,\theta_1,\theta_2\) 。则 \(\theta_0=\theta_1+\theta_2\) 。
\(~~~~ ~~~~\) \(\mathcal{Proof.}\)
\(~~~~ ~~~~\) 首先:
\(~~~~ ~~~~\) 由反正切的加法:
\(~~~~ ~~~~\) 至于为什么其实类比一下 \(\tan(\alpha+\beta)\) 就好了。
\(~~~~ ~~~~\) 所以:
\(~~~~ ~~~~\) \(\mathcal{Q.E.D}\)
\(~~~~ ~~~~\) 而模长的关系为 \(|z_0|=|z_1|\times |z_2|\) 。这个仿照上面去证就好了。
\(~~~~\) 8.复数的指数幂
\(~~~~\) 由欧拉公式:
\(~~~~\) 所以当 \(\theta=\pi\) 时:
\(~~~~\) 9.单位根
\(~~~~\) 下文的 \(n\) 均为 \(2\) 的正整数幂,即满足 \(n=2^k,k\in \text N^*\) 。
\(~~~~ ~~~~\) 9.1 代数定义 在 \(C\) (复数域)下,满足 \(x^n=1\) 的 \(x\) 被称为 \(n\) 次单位根,而由代数基本定理,\(n\) 次单位根共 \(n\) 个。
\(~~~~ ~~~~\) 9.2 复平面上定义 :在复平面上,以原点为圆心,\(1\) 为半径作单位圆。以原点为起点,单位圆的 \(n\) 等分点为终点作 \(n\) 个向量且其中一个向量表示的数为 \(1\) 的所有向量表示的复数。
\(~~~~ ~~~~\) 9.3 本原单位根 在 \(n\) 次单位根中,设幅角为正且最小的一个为 \(\omega_n\) ,则 \(\omega_n\) 为 \(n\) 次本原单位根。
\(~~~~ ~~~~\) 9.4 单位圆上的向量: 若我们在单位圆上有一个幅角为 \(\theta\) 的向量,那么这个向量的终点坐标为 \((\cos \theta,\sin \theta)\) 。
\(~~~~ ~~~~\) 9.4 单位根的值: 由上可知,在 \(n\) 次单位根中,除本原单位根外其余 \(n-1\) 个单位根分别是 \(\omega_n^2,\omega_n^3,\dots,\omega_n^n\) (幅角叠加,模长始终为 \(1\) ),需要注意的是 \(\omega_n^n=1\) 。
\(~~~~\) 若想计算 \(n\) 次单位根的值,肯定要先算出本原单位根。显然 \(n\) 次本原单位根的幅角为 \(\dfrac{2\pi}{n}\) ,那么 \(n\) 次本原单位根的终点坐标为 \((\cos \dfrac{2\pi}{n},\sin \dfrac{2\pi}{n})\) ,它表示的值就是 \(\cos \dfrac{2\pi}{n}+i\times \sin \dfrac{2\pi}{n}\) 。而根据欧拉公式,它可以被记为 \(e^{i \frac{2\pi}{n}}\) ,故任意 \(n\) 次单位根的值 \(\omega_n^k=e^{ik\frac{2\pi}{n}}\) 。
\(~~~~ ~~~~\) 9.5 单位根的性质:
\(~~~~ ~~~~ ~~~~\) 9.5.1: \(\omega_{an}^{ak}=\omega_n^k\) 。证明:代入求值的公式,发现 \(a\) 可以约掉。或者从几何的角度来看,\(n\) 等分单位圆的第 \(k\) 个点就是 \(an\) 等分单位圆的第 \(ak\) 个点。
\(~~~~ ~~~~ ~~~~\) 9.5.2: \(\omega_n^{k+\frac{n}{2}}=-\omega_n^{\frac{n}{2}}\) 。证明:代入求值公式。或者从几何的角度来看,这就相当于单位圆的 \(n\) 等分点的第 \(k\) 个多绕半周单位圆,此时它们表示的值互为相反数。
二、快速傅里叶变换
\(~~~~\) 还记得上面的多项式乘法吗,我们来分析它们的时间复杂度。
\(~~~~\) 若用系数表示法,则两式依次相乘为 \(\mathcal{O(n^2)}\) 。
\(~~~~\) 若用点值表示法,则选点需要 \(\mathcal{O(n^2)}\) (称为求值) ,还原回系数表示法需要 \(\mathcal{O(n^3)}\) (称为插值) 。
\(~~~~\) 因此下面优将分别优化求值和插值。
快速傅里叶变换(FFT)
\(~~~~\) 考虑一个 \(n\) 项,\(n-1\) 次的多项式 \(A(x)=\sum_{i=0}^{n-1} a_i\times x^i\) 。
\(~~~~\) 将它转化为点值表示,一种方法是我们分别带 \(n\) 次单位根进去,但这样是 \(n^2\) 的。
\(~~~~\) 现在将它的各个单项式按次数分为奇次和偶次:
\(~~~~\) 令:\(A_0(x),A_1(x)\) 为两个 \((\frac{n}{2}-1)\) 次多项式,满足:
$~~~~ $则 \(A(x)=A_0(x^2)+x\times A_1(x^2)\) 。
\(~~~~\) 所以,我们的问题变成计算 \(A_0(x^2)\) 和 \(A_1(x^2)\) 的点值,这样可以达到 \(\mathcal O(n)\) 计算 \(A(x)\) 。
\(~~~~\) 此时我们尝试带一个 \(n\) 次单位根 \(\omega_n^k\) :
\(~~~~\) 那顺手再用一下第二个性质(一.9.5.2
),代入 \(\omega_n^{k+\frac{n}{2}}\)
\(~~~~\) 你似乎发现了,代入 \(\omega_n^k\) 与 \(\omega_n^{k+\frac{n}{2}}\) 只会造成 \(A_1\) 的系数符号不同。
\(~~~~\) 那这意味着我们在第一个式子依次取 \(\omega_n^k,k\in [0,\frac n 2-1]\) 时,第二个式子可以马上帮我们算出取 \(\omega_n^k,k\in [\frac n 2,n-1]\) 时的点值。
\(~~~~\) 也就是说我们可以每次缩小问题规模为原来的一半,直到只剩下常数项。
\(~~~~\) 那么时间复杂度由于问题规模每次缩小为原来一半,类比线段树,为 \(\mathcal O(n \log n)\) 。
离散傅里叶逆变换 (IDFT)
\(~~~~\) 上面的过程最后你得到的其实是点值表示法的多项式,所以你必须把它还原为系数表示法。
\(~~~~\) 设: 数列 \(\{b_i\}\) 表示一个 \((n-1)\) 次多项式进行 FFT 后得到的点值,即:
\(~~~~\) 这里直接给一个结论:
\(~~~~\) (你发现这东西跟上面的式子非常相似)
\(~~~~\) \(\mathcal{Proof.}\)
\(~~~~\) 设有另一数列 \(\{c_i\}\) ,满足:
\(~~~~\) 注意是 \(b_i\),因为我们是根据 \(b\) 来推 \(a\) ,也就是化成跟上面类似的形式。
\(~~~~\) 考虑怎么算后面一坨。
\(~~~~\) 设等比数列求和 \(S(x)=\sum_{i=0}^{n-1} x^i\) ,将 \(\omega_n^k\) 代入:
\(~~~~\) \(k<n\) ,分两种情况讨论:
\(~~~~ ~~~~\) 当 \(k\not=0\) 时,原式为 \(0\) ;
\(~~~~ ~~~~\) 当 \(k=0\) 时,直接代入 \(S(\omega_n^0)=S(1)=n\) 。
\(~~~~\) 好,返回上面的那个式子
\(~~~~ ~~~~\) 当 \(j\not = k\) 时,\(\sum\) 后面一坨是 \(0\) ;
\(~~~~ ~~~~\) 当 \(j=k\) 时,\(\sum\) 后面一坨是 \(n\) 。
\(~~~~\) 再代入 \(c_k\) :
\(~~~~\) \(\mathcal{Q.E.D}\).
快速傅里叶逆变换 (IFFT)
\(~~~~\) 上面的式子当然不代也没什么事,毕竟我们只需要 \(\sf{FFT}\) 求出 \(\{c_i\}\) 即可。
\(~~~~\) 具体来说,\(\omega_n^{-ik}=\omega_n^{-ik+n}\) (两次 一.9.5.2
)
\(~~~~\) 所以
\(~~~~\) 我们先求出 \(c_k\) 本身在 \(b_i\) 下的点值,然后把 \(c[1,n]\) (注意不是 \(c[0,n]\) ,因为 \(\omega_n^0=\omega_n^n=1\))倒置一下就好了。
\(~~~~\) 然后你就做到了 \(\mathcal O(n \log n)\) 计算多项式乘法。
FFT的优化 ——迭代实现
\(~~~~\) 你发现原序列为 \([0,7]\) 的上升序列,而之后的序列为原序列每个二进制数的逆序。
\(~~~~\) 所以提前逆序一下二进制数,可以避免递归。然后向上合并即可。
代码
查看代码
#include <cmath>
#include <cstdio>
#include <algorithm>
using namespace std;
int n,m,Lg=0;
const double PI=acos(-1);
struct Complex{
double r,i;
Complex(){}
Complex(double R,double I){r=R,i=I;}
Complex operator +(Complex &x){return Complex(r+x.r,i+x.i);}
Complex operator -(Complex &x){return Complex(r-x.r,i-x.i);}
Complex operator *(Complex &x){return Complex(r*x.r-i*x.i,r*x.i+i*x.r);}
void operator +=(Complex &x){r+=x.r,i+=x.i;}
void operator *=(Complex &x){double tmp=r;r=r*x.r-i*x.i,i=tmp*x.i+i*x.r;}
}F[10000005],G[10000005],H[10000005];
int To[10000005];
void Swap(Complex &a,Complex &b){Complex t=a;a=b;b=t;}
void FFT(Complex *A,int op)
{
for(int i=0;i<n;i++) if(i<To[i]) Swap(A[To[i]],A[i]);
for(int i=1;i<n;i<<=1)// i=每层子问题/2
{
Complex W=Complex(cos(PI/i),sin(PI/i)*op); // 本原单位根 \omega_i
for(int j=0;j<n;j+=i<<1)
{
Complex w=Complex(1,0);
for(int k=0;k<i;k++,w*=W)
{
Complex x=A[j+k],y=w*A[i+j+k];
A[j+k]=x+y;A[i+j+k]=x-y;
}
}
}
}
double Fabs(double x){return x>0?x:-x;}
int main() {
scanf("%d %d",&n,&m);
for(int i=0;i<=n;i++) scanf("%lf",&F[i].r);
for(int i=0;i<=m;i++) scanf("%lf",&G[i].r);
for(m+=n,n=1;n<=m;n<<=1,Lg++);
for(int i=0;i<n;i++) To[i]=(To[i>>1]>>1)|((i&1)<<(Lg-1));
FFT(F,1);FFT(G,1);
for(int i=0;i<n;i++) H[i]=F[i]*G[i];
FFT(H,-1);
for(int i=0;i<=m;i++) printf("%.0f ",Fabs(H[i].r)/n);
return 0;
}
参考资料
题解P3803 【模板】多项式乘法(FFT)————一扶苏一
三、快速数论变换
\(~~~~\) \(\sf{FFT}\) 是挺好,但它涉及到了实数运算,因此掉精度(残缺的字符串)和无法取模都是它的问题。换句话说如果你的多项式系数都很大,那你就没办法了。
\(~~~~\) 因此我们考虑用什么东西来替换单位根,那么就要知道替换的东西要满足哪些性质。
\(~~~~\) 来看单位根的性质:
-
所有 \(\omega_n^k(0\leq k\leq n-1)\) 互不相同:保证带进去时点的坐标不同。
-
\(\omega_{an}^{ak}=\omega_{n}^k\) :\(\sf FFT\) 使用的性质,可以参见上面。
-
\(\omega_n^{k+\frac{n}{2}}=-\omega_n^{\frac{n}{2}}\) :\(\sf FFT\) 使用的性质,可以参见上面。
-
\(\sum_{i=0}^{n-1} \omega_n^{ki}\) 在 \(k \not =0\) 时值为 \(0\) ,在 \(k=0\) 时值为 \(n\) :\(\sf IFFT\) 使用的性质,可以参见上面。
\(~~~~\) 综合上面的性质,我们可以使用原根来替代单位根 。
原根及其性质
\(~~~~\) 考虑一个质数 \(p=nq+1\) 。\(n\) 为 \(2\) 的幂。则 \(p\) 的原根 \(g\) 为使得 \(g^i(0\leq i\leq \varphi(p)=p-1)\) 互不相同的一个正整数。由定义可得,一个数的原根可能有很多个。
\(~~~~\) 由于定义,这些性质原根也是可以满足的。
求一个数的原根
\(~~~~\) 我们使用枚举法,即依次枚举一个数,判断其是否是 \(p\) 的原根。一般来说,质数的原根都很小。
\(~~~~\) 而检验时,我们使用如下的性质:
\(~~~~ ~~~~\) 对于一个数 \(g\) ,最小的满足 \(g^k \equiv 1 \pmod p\) 的正整数 \(k\) 一定是 \(p-1\) 的约数。
\(~~~~\) 故如果我们在 \(p-1\) 的所有约数中找不到一个 \(q\) 满足 \(g^q\equiv 1\pmod p\) ,就说明这是原根。
代码
查看代码
bool check(ll gg,ll x)
{
for(int i=2;i*i<=x-1;i++)
{
if((x-1)%i==0&&(qpow(gg,i,x)==1||qpow(gg,(x-1)/i,x)==1)) return false;
}
return true;
}
void GetG(ll x)
{
ll ret=0;
for(ll i=2;i<x;i++) if(check(i,x)){ret=i;break;}
}
NTT
\(~~~~\) 我们只需要把原来 \(\sf FFT\) 中所有的 \(\omega_n\) 改为 \(g^q\) 即可。但注意除法需要改为求逆元。
代码
查看代码
const int MOD=998244353,g=3,gi=332748118;
ll qpow(ll a,ll b)
{
ll ret=1;
while(b)
{
if(b&1) ret=ret*a%MOD;
a=a*a%MOD;b>>=1;
}
return ret;
}
void NTT(ll *S,int op)
{
for(int i=0;i<n;i++) if(i<To[i]) swap(S[i],S[To[i]]);
for(int i=1;i<n;i<<=1)
{
ll W=qpow(op==1?g:gi,(MOD-1)/(i<<1));
for(int j=0;j<n;j+=i<<1)
{
ll w=1;
for(int k=0;k<i;k++,w=w*W%MOD)
{
ll x=S[j+k],y=w*S[i+j+k]%MOD;
S[j+k]=(x+y)%MOD;S[i+j+k]=(x-y)%MOD;
}
}
}
if(op==-1)
{
ll Inv=qpow(n,MOD-2);
for(int i=0;i<N;i++) S[i]=S[i]*Inv%MOD;
}
}
参考资料
四、三模NTT/MTT
咕,咕咕咕 (拟声)
五、多项式求逆
\(~~~~\) 多项式做除法的话,就必须像普通数在取余的情况下做除法一样求得逆元。
\(~~~~\) 一般这东西朗读并背诵全文就行了。
本部分的一小点前置知识
\(~~~~\) 多项式的度数: 多项式最高次项的次数(或者说就是多项式的次数),多项式 \(A(x)\) 的度为 \(deg_A\) 。
\(~~~~\) 多项式的带余除法: 设有多项式 \(A(x),B(x)\) ,则存在唯一的 \(Q(x),R(x)\) 满足 \(A(x)=B(x)Q(x)+R(x)\) ,且 \(deg_R<deg_B\) 称 \(Q(x)\) 为商,\(R(x)\) 为余数。同时有:
\(~~~~\) 多项式的逆元: 若存在多项式 \(B(x)\) ,满足 \(A(x) \times B(x) \equiv 1 \pmod{x^n}\) ,则称 \(B(x)\) 为 \(A(x)\) 的逆元,记作 \(A^{-1}(x)\) 。后面的模数 \(x^n\) 其实相当于舍弃掉乘积在 \(x^{n+1}\) 及以后的所有系数。
如何求逆元
\(~~~~\) 假设我们现在求到了 \(A(x)\) 在 \(\bmod x^{\lceil \frac{n}{2} \rceil}\) 下的逆元 \(B'(x)\) ,现在想求 \(B(x)=A^{-1}(x)\) 。
\(~~~~\) 因此有:
\(~~~~\) 显然,一个在更广泛的模意义下的逆元乘起来它的系数从低到高依然是 \(1,0,0,0,0,\dots\) ,因此有:
\(~~~~\) 下 \(-\) 上:
\(~~~~\) 显然 \(A \bmod x^{\lceil \frac{n}2\rceil}\) 不为 \(0\) ,故 \(B-B' \bmod x^{\lceil \frac{n}2\rceil}=0\)
\(~~~~\) 因此:
\(~~~~\) 考虑如何上升到 \(x^n\) ,很难不难想到进行平方:
\(~~~~\) \(\mathcal{Proof.}\)
\(~~~~\) 设平方后的多项式为 \(F(x)\) ,系数序列为 \(b_i\)
\(~~~~\) 对于 \(i=0\) ,显然跟之前一样,即 \(b_i=1\) 。
\(~~~~\) 对于 \(1\leq i \leq \lceil \dfrac{n}{2} \rceil\) :由于原来其本身和次数小于它的所有非常数项的系数是 \(0\) ,卷积后肯定也是 \(0\) 。
\(~~~~\) 对于 \(\lceil \dfrac{n}{2} \rceil < i \leq n\) ,\(b_i=\sum_{j=0}^i a_j\times a_{i-j}\) ,而由于原先的 \(0\) 是过半的,因此不论 \(j\) 取什么值,都有,各加数都为 \(0\) 。
\(~~~~\) 综上,对于 \(i=0\) ,\(b_i=1\) ;否则:\(b_i=0\) 。
\(~~~~\) \(\mathcal{Q.E.D.}\)
\(~~~~\) 展开平方:
\(~~~~\) 然而如果直接把 \(B^2\) 甩过去,会出现多项式开根,而多项式开根又需要多项式求逆,所以这两个东西都不可做 ,所以我们两边同时卷上 \(A\) :
\(~~~~\) 而 \(AB\equiv 1\pmod{x^n}\):
\(~~~~\) 每层递归计算就行了,出口就是多项式只剩下常数项时,这时就是普通的求逆元了。
代码
查看代码
#include <cstdio>
#include <algorithm>
#define ll long long
using namespace std;
int To[1000005],N;
const ll MOD=998244353;
ll qpow(ll a,ll b) {......}//快速幂
const int g=3,gi=qpow(g,MOD-2);
void NTT(ll *S,int op)//NTT
{
......
if(op==-1)
{
ll Inv=qpow(N,MOD-2);
for(int i=0;i<N;i++) S[i]=S[i]*Inv%MOD;
}
}
int n;
ll F[1000005],G[1000005],C[1000005];
void GetInv(int deg,ll *A,ll *B)
{
if(deg==1){B[0]=qpow(A[0],MOD-2);return;}
GetInv((deg+1)>>1,A,B);
for(N=1;N<=(deg<<1);N<<=1);
for(int i=0;i<N;i++) To[i]=(To[i>>1]>>1)|((i&1)*(N>>1)),C[i]=A[i];
for(int i=deg;i<N;i++) C[i]=0;
NTT(C,1);NTT(B,1);
for(int i=0;i<N;i++) B[i]=B[i]*(((2ll-B[i]*C[i]%MOD)+MOD)%MOD)%MOD;
NTT(B,-1);
for(int i=deg;i<N;i++) B[i]=0;
}
int main() {
scanf("%d",&n);
for(int i=0;i<n;i++) scanf("%lld",&F[i]);
GetInv(n,F,G);
for(int i=0;i<n;i++) printf("%lld ",G[i]);
return 0;
}
六、多项式开根
\(~~~~\) 上面是不是提到了多项式开根?
\(~~~~\) 求多项式 \(g(x)\) 使得 \(g^2(x)\equiv f(x) \pmod {x^n}\) 即多项式开根。
\(~~~~\) 一般情况下,被开方的多项式常数项为完全平方数。(否则就没法直接 \(\sf NTT\) 做了)
求解过程
\(~~~~\) 跟求逆差不多,我们假设 \(g^2_0(x) \equiv f(x) \pmod{2^{\lceil \frac x 2 \rceil}}\) ,即上一层的解。
\(~~~~\) 那么:
\(~~~~\) 继续沿用上面的套路,两边平方,具体证明见上。
\(~~~~\) 进行一个变形:
\(~~~~\) 两边开根:
\(~~~~\) 然后我们就得到了 \(g\) 的式子。
\(~~~~\) 具体实现的时候仍然递归,但由于每层还要求逆,所以时间复杂度为 \(\mathcal{O(n \log^2 n)}\) 。
代码
不知道为什么,我的代码不开O2 过洛谷模板,9s多,开了 3s多。
所以仅供参考。具体实现方式还有开根时每层循环来迭代的,可能会快一些。
查看代码
#include <cstdio>
#include <algorithm>
#define ll long long
using namespace std;
const ll MOD=998244353;
ll To[5000005],N;
ll qpow(ll a,ll b){...}//快速幂
const ll g=3,gi=qpow(g,MOD-2);
const ll Inv2=qpow(2,MOD-2);
inline ll Add(ll a,ll b){return ((a+b)%MOD+MOD)%MOD;}
inline ll Mul(ll a,ll b){return a*b%MOD;}
void NTT(ll *S,ll op){...}//多项式乘法
ll C[5000005];
void GetInv(int deg,ll *A,ll *B){...}//求逆元
ll inv[5000005],D[5000005];//注意新开数组存
void GetSqrt(int deg,ll *A,ll *B)
{
if(deg==1){B[0]=1;return;}
GetSqrt((deg+1)>>1,A,B);
for(N=1;N<=(deg<<1);N<<=1);
for(int i=0;i<N;i++) To[i]=(To[i>>1]>>1)|((i&1)*(N>>1)),D[i]=inv[i]=0;
for(int i=0;i<deg;i++) D[i]=A[i];
GetInv(deg,B,inv);
NTT(inv,1);NTT(D,1);NTT(B,1);
for(int i=0;i<N;i++) B[i]=Mul(Add(B[i],Mul(D[i],inv[i])),Inv2);
NTT(B,-1);
for(int i=deg;i<N;i++) B[i]=0;
}
ll F[5000005],G[5000005];
int main() {
int n;read(n);
for(int i=0;i<n;i++) scanf("%lld",&F[i]);
GetSqrt(n,F,G);
for(int i=0;i<n;i++) printf("%lld ",G[i]);
return 0;
}
七、多项式除法
\(~~~~\) 一般说的多项式除法指带余除法,即给定多项式 \(F(x),G(x)\) ,求 \(Q(x),R(x)\) 满足 \(F(x)=G(x)*Q(x)+R(x)\) ,其中 \(deg_R<deg_G\) 。
\(~~~~\) 这个东西其实一般也是朗读并全文背诵。
\(~~~~\) 顺便提一下,一般人工做这个东西的时候都是用大除法。
求解过程
\(~~~~\) 由上面的定义,我们知道:
\(~~~~\) 显然更换代入的 \(x\) 上式仍然成立:
\(~~~~\) 设: \(deg_F=n,deg_G=m\),两边同时乘上 \(x^n\) :
\(~~~~\) 我们发现一件事:\(x^n\times \sum_{i=0}^{n} a_i\times \dfrac{1}{x^i}=\sum_{i=0}^n a_{n-i}x^i\) ,即得到了原来的多项式系数序列翻转后的多项式。
\(~~~~\) 把上面的式子都拆一下,不难确认:\(deg_Q=n-m,deg_R<m\) ,因此:
\(~~~~\) 记 \(x^n\times F(\dfrac{1}{x})\) 为 \(F_R(x)\) ,表示翻转后的 \(F(x)\) ,则:
\(~~~~\) 来思考一个问题,我们为什么要这么拆?想到原来为什么不可做的原因就是 \(R(x)\) 无法确定。
\(~~~~\) 而现在,\(R_R(x)\) 乘上了一个 \(x^{n-m+1}\) ,可以借此消掉它:
\(~~~~\) 所以:
\(~~~~\) 显然,\(n-m=deg_{G_R}< n-m+1\) 所以这是没有影响的,直接对 \(G_R(x)\) 求逆元即可。
\(~~~~\) 最后,\(R(x)=F(x)-G(x)*Q(x)\) ,直接求即可。
代码
这题细节挺多的,稍不注意就会打挂,这里代码放完整版,个人错过的点注释在下面了。
查看代码
#include <cstdio>
#include <algorithm>
#define ll long long
using namespace std;
int N,To[4000005];
const ll MOD=998244353;
inline ll Add(ll a,ll b){return (((a+b)%MOD)+MOD)%MOD;}
inline ll Mul(ll a,ll b){return a*b%MOD;}
ll qpow(ll a,ll b)
{
ll ret=1;
while(b)
{
if(b&1) ret=Mul(ret,a);
b>>=1;a=Mul(a,a);
}
return ret;
}
const ll gg=3,ggi=qpow(gg,MOD-2);
void NTT(ll *S,ll op)
{
for(int i=0;i<N;i++) if(i<To[i]) swap(S[i],S[To[i]]);
for(int i=1;i<N;i<<=1)
{
ll W=qpow(op==1?gg:ggi,(MOD-1)/(i<<1));
for(int j=0;j<N;j+=i<<1)
{
ll w=1;
for(int k=0;k<i;k++,w=Mul(w,W))
{
ll x=S[j+k],y=Mul(S[i+j+k],w);
S[j+k]=Add(x,y);S[i+j+k]=Add(x,-y);
}
}
}
if(op==-1)
{
ll Inv=qpow(N,MOD-2);
for(int i=0;i<N;i++) S[i]=Mul(S[i],Inv);
}
}
ll C[4000005],A1[4000005],A2[4000005];
void mul(ll *A,ll *B,ll n,ll m)
{
for(N=1;N<=n+m;N<<=1);
for(int i=0;i<N;i++) To[i]=(To[i>>1]>>1)|((i&1)*(N>>1));
for(ll i=0;i<N;i++) A1[i]=A[i],A2[i]=B[i];
NTT(A1,1);NTT(A2,1);
for(ll i=0;i<N;i++) A[i]=Mul(A1[i],A2[i]);
NTT(A,-1);
}
void GetInv(int deg,ll *A,ll *B)
{
if(deg==1){B[0]=qpow(A[0],MOD-2);return;}
GetInv((deg+1)>>1,A,B);
for(N=1;N<=(deg<<1);N<<=1);
for(int i=0;i<N;i++) To[i]=(To[i>>1]>>1)|((i&1)*(N>>1)),C[i]=A[i];
for(int i=deg;i<N;i++) C[i]=0;
NTT(C,1);NTT(B,1);
for(int i=0;i<N;i++) B[i]=Mul(Add(2ll,-Mul(C[i],B[i])),B[i]);
NTT(B,-1);
for(int i=deg;i<N;i++) B[i]=0;
}
ll f[4000005],g[4000005],q[4000005],inv[4000005];
void Div(int n,int m,ll *F,ll *G,ll *Q)
{
for(int i=0;i<n;i++) Q[i]=F[i];
for(int i=0;i<m;i++) g[i]=G[i];
reverse(Q,Q+n);reverse(g,g+m);
for(int i=n-m+1;i<m;i++) g[i]=0;
GetInv(n-m+1,g,inv);//注意即使g只有m次,这里求在mod x^{n-m+1} 意义下的逆元。
mul(Q,inv,n,n-m+1);//最好封装一下乘法,不然很容易挂
reverse(Q,Q+n-m+1);
for(int i=n-m+1;i<=N;i++) Q[i]=0;
for(int i=0;i<=N;i++) inv[i]=0;
}
void Rest(int n,int m,ll *F,ll *G,ll *Q,ll *R)
{
for(int i=0;i<n;i++) f[i]=F[i];
for(int i=0;i<m;i++) g[i]=G[i];
for(int i=0;i<n-m+1;i++) q[i]=Q[i];
mul(g,q,m,n-m+1);
for(int i=0;i<m-1;i++) R[i]=Add(f[i],-g[i]);//这里最后多项式加减直接对应系数加减就好了
}
ll F[4000005],G[4000005],Q[4000005],R[4000005];
int main() {
int n,m;scanf("%d %d",&n,&m);n++;m++;
for(int i=0;i<n;i++) scanf("%lld",&F[i]);
for(int i=0;i<m;i++) scanf("%lld",&G[i]);
Div(n,m,F,G,Q); for(int i=0;i<n-m+1;i++) printf("%lld ",Q[i]);puts("");
Rest(n,m,F,G,Q,R); for(int i=0;i<m-1;i++) printf("%lld ",R[i]);
return 0;
}
来自 2021/08/23 的吐槽:果然忘完了,回来一看不仅没填坑而且还一堆错,为让之前的各位读者或许在某些地方感到困惑而致歉。
2021/08/26 的更新
八、多项式对数函数(多项式 ln)
\(~~~~\) 求对数函数啊~
前置芝士
\(~~~~\) 多项式求导:\(f(x)\) 的导函数记为 \(f'(x)\)
\(~~~~\) 对于函数 \(f(x)=x^a\),有 \(f'(x)=ax^{a-1}\)。
\(~~~~\) 对于函数 \(f(x)=g(x)+h(x)\),有 \(f'(x)=g'(x)+h'(x)\)。
\(~~~~\) 对于复合函数, \(f(g(x))'=f'(g(x))g'(x)\)
\(~~~~\) 对数函数的导函数: \(\ln'(x)=\dfrac{1}{x}\)。
\(~~~~\) 多项式积分: \(f(x)\) 的积分记为 \(\int f(x) dx\)
\(~~~~\) 对于函数 \(f(x)=x^a\),有 \(\int f(x) dx=\dfrac{1}{a+1}x^{a+1}\)。
\(~~~~\) 对于函数 \(f(x)=g(x)+h(x)\) ,有 \(\int f(x) dx =\int g(x) dx+\int h(x) dx\)。
\(~~~~\) 由上述定义可以积分和求导互为逆运算。
求解过程
\(~~~~\) 记 \(G(x) \equiv \ln (F(x)) \pmod {x^n}\) ,则我们要求的即为 \(G(x)\)。
\(~~~~\) 对右边求导: \(\ln(F(x))'=ln'(F(x))F'(x)=\dfrac{F'(x)}{F(x)}\) ,故:
\(~~~~\) 那么求出 \(F'(x)\times F^{-1}(x)\) 后再积分一下就是 \(G(x)=\ln F(x)\) 了。
代码
查看代码
int A1[500005],A2[500005];
void mul(int *A,int *B,int n,int m)
{
for(N=1;N<=n+m;N<<=1);
for(int i=0;i<N;i++) To[i]=(To[i>>1]>>1)|((i&1)*(N>>1));
for(int i=0;i<N;i++) A1[i]=A[i],A2[i]=B[i];
NTT(A1,1);NTT(A2,1);
for(int i=0;i<N;i++) A[i]=Mul(A1[i],A2[i]),A1[i]=A2[i]=0;
NTT(A,-1);
}
int C[500005];
void GetInv(int *A,int *B,int deg)
{
if(deg==1){B[0]=qpow(A[0],MOD-2);return;}
GetInv(A,B,(deg+1)>>1);
for(N=1;N<=(deg<<1);N<<=1);
for(int i=0;i<N;i++) To[i]=(To[i>>1]>>1)|((i&1)*(N>>1)),C[i]=A[i];
for(int i=deg;i<N;i++) C[i]=0;
NTT(C,1);NTT(B,1);
for(int i=0;i<N;i++) B[i]=Mul(Add(2ll,-Mul(C[i],B[i])),B[i]);
NTT(B,-1);
for(int i=deg;i<N;i++) B[i]=0;
}
void GetDev(int *F,int *G,int deg)
{
for(int i=1;i<deg;i++) G[i-1]=1ll*F[i]*i%MOD;
G[deg-1]=0;
}
void GetInvDev(int *F,int *G,int deg)
{
for(int i=1;i<deg;i++) G[i]=1ll*F[i-1]*qpow(i,MOD-2)%MOD;
G[0]=0;
}
int a[500005],b[500005];
void GetLn(int *F,int *G,int deg)
{
GetDev(F,a,deg);GetInv(F,b,deg);
mul(a,b,deg,deg);
GetInvDev(a,G,deg);
}