多项式 学习笔记

记录多项式以及一些相关数论前置芝士qwq

拉格朗日插值

由 $n$ 个点 $(x_i,y_i)$ 可以唯一确定一个 $n-1$ 次多项式 $f(x)$，要求出 $f(k)$ 的值可以 $O(n^3)$ 高斯消元，而拉格朗日插值可以做到 $O(n^2)$ 求解。

由 $f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n$，$f(k)=a_0+a_{1}k+a_2{k}^2+\cdots +a_n{k}^n$ 得 $f(x)-f(k)=a_1(x-k)+a_2(x^2-k^2)+\cdots +a_n(x^n-k^n)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}x-k)$，则 $f(x)\equiv f(k)(\mathrm{m}\mathrm{o}\mathrm{d}x-k)$。

把 $n$ 个点都带入方程，可得同余方程组 $f(k)\equiv y_i(\mathrm{m}\mathrm{o}\mathrm{d}k-x_i))(1\le i\le n)$

套用 CRT 的求解，令 $M=\prod _{i=1}^n(k-x_i),m_i=\frac{M}{k-x_i}$，则有：
$f(k)=\sum _{i=1}^n{y}_i\times m_i\times m_i^{-1}$
$=\sum _{i=1}^n{y}_i\times \prod _{j\ne i}(k-x_j)\times m_i^{-1}$
$=\sum _{i=1}^n{y}_i\times \prod _{j\ne i}(k-x_i+x_i-x_j)\times \prod _{j\ne i}\frac{1}{x_i-x_j}$
$=\sum _{i=1}^n{y}_i\times \prod _{j\ne i}\frac{k-x_j}{x_i-x_j}$

中间在 $k-x_i$ 剩余系中转化了一下逆元，最后求出的就是拉格朗日插值的式子，在模意义和非模意义下通用。

$x_i$ 是连续整数的插值

对于 ${\forall }_{i=1}^n{x}_i=i$，可以将拉格朗日插值的式子进一步简化：

$\{\begin{array}{l}\prod _{j\ne i}(i-j)=(-1{)}^{n-i}(i-1)!(n-i)!\\ \prod _{j\ne i}(k-j)=\frac{\prod _{j=1}^n(k-j)}{k-i}\end{array}$
$⟹f(k)=\sum _{i=1}^n{y}_i\times (-1{)}^{n-i}\frac{\prod _{j=1}^n(k-j)}{(i-1)!(n-i)!(k-i)}$

分子的和分母都可以预处理，于是可以做到 $O(n)$。

永远记不住的线性求逆元柿子：$inv[i]=(p-p/i)\times inv[p\%i]$

重心拉格朗日插值

当需要新加入一个点的时候，用朴素的拉插式子算是 $O(n^2)$ 的，但可以推柿子做到 $O(n)$：

$f(k)=\sum _{i=1}^n{y}_i\times \prod _{j\ne i}\frac{k-x_j}{x_i-x_j}$
$=\sum _{i=1}^n{y}_i\times \frac{\prod _{j=1}^{n}k-x_j}{(k-x_i)\prod _{j\ne i}(x_i-x_j)}$
$=\prod _{i=1}^n(k-x_i)\sum _{i=1}^n\frac{y_i}{(k-x_i)\prod _{j\ne i}(x_i-x_j)}$

设 $g(k)=\prod _{i=1}^n(k-x_i)$，$t_i=\prod _{j\ne i}\frac{1}{x_i-x_j}$，则 $f(k)=g(k)\sum _{i=1}^n\frac{y_i{t}_i}{k-x_i}$，每次新增一个点 $n$ 可以 $O(1)$ 更新 $g(k)$，$O(n)$ 求 $t_n$，$O(n)$ 更新 $t_1$ 到 $t_{n-1}$，再 $O(n)$ 更新 $f(k)$，这个点就插入成功了。

快速傅里叶变换

前置芝士：

• 多项式的系数表示法：即常见的函数形式，形如 $f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n$；
• 多项式的点值表示法：对于一个 $n-1$ 次多项式，$n$ 个不同的 $x$ 代入能得到 $n$ 个不同的点 $(x_i,y_i)$，这 $n$ 个点唯一确定了该多项式；
• 多项式的系数表示和点值表示可以相互转化。

快速傅里叶变换（FFT）就是一种 $O(n\log n)$ 将系数表示转化为点值表示，进而求出两个多项式的乘积的方法。

离散傅里叶变换

傅里叶把多项式先填到最高次项（系数可以为 $0$）为 $n=2^k$ 的形式，然后规定点值表示中的点为 $n$ 个模长为 $1$ 的复数，使得这些复数在平面直角坐标系中把单位圆 $n$ 等分。

为了方便，将这些点以 $(1,0)$ 为起点逆时针编号为 $0$ 到 $n-1$，第 $k$ 个点的坐标为 $(\cos \frac{k}{n}2\pi ,\sin \frac{k}{n}2\pi )$，表示的复数为 $w_n^k$，其中 $w_n^1$ 叫作 $n$ 次单位根。

这种点值选取的方式就是离散傅里叶变换（DFT）。

根据复数和三角函数的一些知识，得到（设 $a$ 为常数）：

• $w_n^k=w_{an}^{ak}$
• $w_n^k=-w_n^{k+\frac{n}{2}}$
• $w_n^k=w_n^{k+an}$

快速傅里叶变换

在 DFT 的基础上，我们推一些柿子：将 $f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^{n-1}$ 按 $x$ 的幂次的奇偶拆开并化简：

$f(x)=(a_0+a_2{x}^2+\cdots +a_{n-2}{x}^{n-2})+(a_{1}x+a_3{x}^3+\cdots +a_{n-1}{x}^{n-1})$
$=(a_0+a_2(x^1{)}^2+\cdots +a_{n-2}(x^{\frac{n-2}{2}}{)}^2)+x(a_1+a_3(x^1{)}^2+\cdots a_{n-2}(x^{\frac{n-2}{2}}{)}^2)$

令 $g(x)=a_0+a_{2}x+\cdots +a_{n-2}{x}^{\frac{n-2}{2}}$，$h(x)=a_1+a_{3}x+\cdots +a_{n-1}{x}^{\frac{n-2}{2}}$，则 $f(x)=g(x^2)+xh(x^2)$。

把 DFT 中选择的点代入：

• $k<\frac{n}{2}$，$f(w_n^k)=g(w_n^{2k})+w_n^{k}h(w_n^{2k})=g(w_{\frac{n}{2}}^k)+w_n^{k}h(w_{\frac{n}{2}}^k)$
• $k\ge \frac{n}{2}$，$f(w_n^k)=g(w_n^{2k+n})+w_n^{k+\frac{n}{2}}h(w_n^{2k+n})=g(w_{\frac{n}{2}}^k)-w_n^{k}h(w_{\frac{n}{2}}^k)$

那么这就是分治了，但是由于大常数选手的递归版跑得过慢等原因，代码就不放了/kk

递归本身常数就是比较大的，考虑能不能优化。模拟一下每次选奇偶点分前后两半的过程，发现如果把 $i$ 的二进制表示前后翻转所得的数为 $j$，那么递归后原本在位置 $i$ 的数会跑到位置 $j$，这就是位逆序变换。因此可以预处理出每个数的最终位置，一点点向上合并即可。

实现的时候，复数可以用 STL 的 complex，复数建议手写因为 complex 常数很大，预处理 $i$ 翻转二进制位所得的数 $po{s}_i$，就可以写出 FFT 了：

inline void fft(com *a,com *omg){
	ff(i,0,n-1) if(pos[i]<i) swap(a[pos[i]],a[i]);
	for(int len=2;len<=n;len<<=1){
		int m=len>>1;
		for(int i=0;i<n;i+=len) ff(j,0,m-1){
			com tmp=a[i+j+m]*omg[n/len*j];
			a[i+j+m]=a[i+j]-tmp,a[i+j]+=tmp;
		}	
	}
} 

离散傅里叶逆变换

我们用 FFT 成功加速了系数表示转点值表示的过程，但目的是想加速多项式乘法。FFT 是 $O(n\log n)$ 的，两个点值表示的多项式相乘是 $O(n)$ 的，想求乘法我们还需要快速将点值表示转化为系数表示，就需要离散傅里叶逆变换（IDFT）。

假设我们对 $f(x)$ 进行 DFT之后得到了一系列纵坐标 $(y_0,y_1,\cdots ,y_{n-1})$，构造函数 $g(x)=y_0+y_{1}x+\cdots +y_{n-1}{x}^{n-1}$，把 $w_n^0,w_n^{-1},\cdots ,w_n^{-(n-1)}$ 代入得到 $(z_0,z_1,\cdots ,z_{n-1})$：

$z_k=\sum _{i=0}^{n-1}{y}_i\times (w_n^{-k}{)}^i$
$=\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j\times (w_n^i{)}^j\times (w_n^{-k}{)}^i$
$=\sum _{j=0}^{n-1}{a}_j\sum _{i=0}^{n-1}(w_n^{j-k}{)}^i$
而 $\sum _{i=0}^{n-1}(w_n^{j-k}{)}^i=\{\begin{array}{l}j=kn\\ j\ne k\frac{(w_n^{j-k}{)}^n-1}{w_n^{j-k}-1}=0\end{array}$

故当且仅当 $j=k$，$z_k=a_j\times n$。

所以对 $g(x)$ 代入 $w_n^0,w_n^{-1},\cdots ,w_n^{-(n-1)}$ 之后求 DFT 的结果除以 $n$ 就可以得到 $f(x)$ 的各项系数，而形如 $w_n^{-1}$ 的就是共轭复数。

学了一大顿之后就又可以写高精乘法了：（亲测不预处理 $w_n^k$ 会更快）

#include<bits/stdc++.h>
#define ll long long
#define rep(i,s,e) for(int i=(s);i<=(e);++i)
#define Rep(i,s,e) for(int i=(e);i>=(s);--i)
using namespace std;
inline int read(){
    int x=0,f=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){if(ch=='-') f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    return x*f;
}
const int N=1<<21;
const double pai=acos(-1);
int n1,n2,n=1,m,pos[N],ans[N];
char s1[N],s2[N];
struct qwq{
    double x,y;
}a[N],b[N];
inline qwq operator + (const qwq &a,const qwq &b){
    return {a.x+b.x,a.y+b.y};
}
inline qwq operator - (const qwq &a,const qwq &b){
    return {a.x-b.x,a.y-b.y};
}
inline qwq operator * (const qwq &a,const qwq &b){
    return {a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y};
}
inline void fft(qwq a[],int flag){
    rep(i,0,n-1) if(pos[i]<i) swap(a[pos[i]],a[i]);
    for(int len=2;len<=n;len<<=1){
        int m=len>>1;
        qwq x={cos(2*pai/len),flag*sin(2*pai/len)};
        for(int i=0;i<n;i+=len){
            qwq now={1,0};
            rep(j,0,m-1){
                qwq tmp=a[i+j+m]*now;
                a[i+j+m]=a[i+j]-tmp,a[i+j]=a[i+j]+tmp;
                now=now*x;
            }
        }
    }
}
signed main(){
    scanf("%s",s1),n1=strlen(s1);
    scanf("%s",s2),n2=strlen(s2);
    rep(i,0,n1-1) a[i].x=s1[n1-i-1]-'0';
    rep(i,0,n2-1) b[i].x=s2[n2-i-1]-'0';
    while(n<n1+n2) n<<=1,++m;
    rep(i,0,n-1){
        int tmp=0;
        rep(j,0,m-1) if(i&(1<<j)) tmp|=(1<<m-j-1);
        pos[i]=tmp;
    }
    fft(a,1),fft(b,1);
    rep(i,0,n-1) a[i]=a[i]*b[i];
    fft(a,-1);
    rep(i,0,n-1){
        ans[i]+=floor(a[i].x/n+0.5);
        if(ans[i]>=10) ans[i+1]+=ans[i]/10,ans[i]%=10;
    }
    while(n>1&&ans[n-1]==0) --n;
    Rep(i,0,n-1) putchar(ans[i]+'0');
    putchar('\n');
}

快速数论变换

感觉 OI Wiki 上原根之外的前置知识都可以直接背下来？

原根

阶：使得 $a^n\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}m)$ 成立的最小正整数 $n$ 叫做 $a$ 模 $m$ 的阶，符号 ${\delta }_m(a)$。

一些性质：

• $\forall a^n\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}m),{\delta }_m(a)\mid n⟹{\delta }_m(a)\mid \varphi (m)$
• ${\forall }_{i,j\in [1,{\delta }_m(a)],i\ne j}{a}^i\not\equiv a^j(\mathrm{m}\mathrm{o}\mathrm{d}m)$
• $\gcd (a,m)=1,{\delta }_m(a^k)=\frac{{\delta }_m(a)}{\gcd (k,{\delta }_m(a))}$

原根：若 $\gcd (a,m)=1,{\delta }_m(a)=\varphi (m)$，则 $a$ 是 $m$ 的原根。

• 判定定理：${\forall }_{p\mid \varphi (m)}{a}^{\frac{\varphi (m)}{p}}\not\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}m)⟺a$ 是 $m$ 的原根；
• 存在定理：只有 $2,4,p^a,2p^a$ 才存在原根，其中 $p$ 为奇素数；
• 原根个数：若 $m$ 有原根，则其原根个数为 $\varphi (\varphi (m))$；
• $m$ 的最小原根 $g$ 不超过 $m^{\frac{1}{4}}$，所有其它原根均为 $g^k(\gcd (k,\varphi (m)=1))$。

快速数论变换

FFT 的缺点在于浮点数运算会产生精度问题，实际上在模数是一些特定质数 $p=k\times 2^n+1$ 时可以将单位根替换为原根进行计算。这些质数被称为 NTT 模数，常见的是 $998244353$，原根为 $g=3$。

根据原根的性质可知 $g^{kn}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 且 $i\in [0,p)$ 范围内 $g_i\mathrm{m}\mathrm{o}\mathrm{d}p$ 两两不同，于是单位根满足的性质原根全部满足，就可以类似 FFT 来求解，IDFT 的求共轭复数变为求逆元即可。

模板代码：

#include<bits/stdc++.h>
#define ll long long
#define rep(i,s,e) for(int i=(s);i<=(e);++i)
#define Rep(i,s,e) for(int i=(e);i>=(s);--i)
using namespace std;
inline int read(){
	int x=0,f=1;
	char ch=getchar();
	while(ch>'9'||ch<'0'){if(ch=='-') f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
	return x*f;
}
const int N=3e6+5,mod=998244353;
int n1,n2,m,n=1,p[N];
ll a[N],b[N];
inline ll ksm(ll x,int y){
    ll res=1;
    for(;y;y>>=1){
        if(y&1) res=res*x%mod;
        x=x*x%mod;
    }
    return res;
}
inline void NTT(ll a[],int flag){
    rep(i,0,n-1) if(p[i]<i) swap(a[i],a[p[i]]);
    for(int len=2;len<=n;len<<=1){
        int m=len>>1;ll x=ksm(3,(mod-1)/len);
        if(flag==-1) x=ksm(x,mod-2);
        for(int i=0;i<n;i+=len){
            ll now=1;
            rep(j,0,m-1){
                ll u=a[i+j],v=a[i+j+m]*now%mod;
                a[i+j]=(u+v)%mod,a[i+j+m]=(u-v+mod)%mod;
                now=now*x%mod;
            }
        }
    }
    if(flag==-1){
        int inv=ksm(n,mod-2);
        rep(i,0,n-1) a[i]=a[i]*inv%mod;
    }
}
signed main(){
    n1=read(),n2=read();
    rep(i,0,n1) a[i]=read();
    rep(i,0,n2) b[i]=read();
    while(n<(n1+n2+2)) n<<=1,++m;
    rep(i,1,n){
        int now=0;
        rep(j,0,m-1) if(i&(1<<j)) now|=(1<<m-j-1);
        p[i]=now;
    }
    NTT(a,1),NTT(b,1);
    rep(i,0,n-1) a[i]=a[i]*b[i]%mod;
    NTT(a,-1);
    rep(i,0,n1+n2) printf("%lld ",a[i]);
    putchar('\n');
}

快速莫比乌斯变换 & 快速沃尔什变换

感谢可爱的猫告诉我基于分治的 and 和 or 卷积叫 FWT，基于 dp 卷积的叫 FMT。好强大的樱雪喵。

与卷积、或卷积和异或卷积的本质都是仿照 FFT 将多项式 $A,B$ 转化为点值表示 $FWT[A],FWT[B]$，根据 $FWT[C]=FWT[A]\times FWT[B]$ 得到 $FWT[C]$，再逆运算得到多项式 $C$（以下称逆运算为 $UFWT[A]$）

OR 卷积

要求 $c_k=\sum _{i|j=k}{a}_i{b}_j$，考虑根据 $i|k=k,j|k=k⟺(i|j)|k=k$，设 $FWT[A{]}_k=\sum _{i|k=k}{a}_i$。

考虑分治，把长为 $2^n$ 的数组写成分成前后两段，前一段二进制最高位全为 $0$，后一段全为 $1$，更低位前后两段的情况是一样的。把这两段分别记为 $A0$ 和 $A1$，则 $A1$ 的子集也包括最高位为 $0$ 的子集，即：
$$\begin{gathered} FWT[A]=merge(FWT[A0],FWT[A0]+FWT[A1]) \\ UFWT[A]=merge(UFWT[A0],UFWT[A1]-UFWT[A0]) \end{gathered}$$

AND 卷积

要求 $c_k=\sum _{i\&j=k}{a}_i{b}_j$，考虑根据 $i\&k=k,j\&k=k⟺(i\&j)\&k=k$，设 $FWT[A]=\sum _{i\&k=k}{a}_i$。

同理 OR 卷积可知
$$\begin{gathered} FWT[A]=merge(FWT[A0]+FWT[A1],FWT[A0]) \\ UFWT[A]=merge(UFWT[A0]-UFWT[A1],UFWT[A1]) \end{gathered}$$

XOR 卷积

要求 $c_k=\sum _{i\mathrm{xor}j=k}{a}_i{b}_j$，设 $d(x)$ 为 $x$ 二进制下 $1$ 个数的奇偶性，根据 $d(i\&k)\mathrm{xor}d(j\&k)=d((i\mathrm{xor}j)\&k)$，设 $FWT[A]=\sum _i(-1{)}^{d(i\&k)}{a}_i$。

由于只有 $1\&1=1$，会使得 $d()$ 的奇偶性改变，所以可得：
$$\begin{gathered} FWT[A]=merge(FWT[A0]+FWT[A1],FWT[A0]-FWT[A1]) \\ UFWT[A]=merge(\frac{UFWT[A0]+UFWT[A1]}{2},\frac{UFWT[A0]-UFWT[A1]}{2}) \end{gathered}$$

于是就可以在 $O(n\log n)$ 的复杂度内分别求出三个卷积。

#include<bits/stdc++.h>
#define ll long long
#define rep(i,s,e) for(int i=(s);i<=(e);++i)
#define Rep(i,s,e) for(int i=(e);i>=(s);--i)
using namespace std;
inline int read(){
    int x=0,f=1;
    char ch=getchar();
    while(ch>'9'||ch<'0'){if(ch=='-') f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    return x*f;
}
const int N=(1<<17)+5,mod=998244353;
int n;
ll a[N],b[N],c[N],x[N],y[N],inv;
inline ll ksm(ll x,int y){
    ll res=1;
    for(;y;y>>=1){
        if(y&1) res=res*x%mod;
        x=x*x%mod;
    }
    return res;
}
inline void OR(ll a[],int flag){
    for(int len=2;len<=n;len<<=1){
        int m=len>>1;
        for(int i=0;i<n;i+=len) rep(j,0,m-1){
            int u=a[i+j],v=a[i+j+m];
            a[i+j+m]=flag==1?(u+v)%mod:(v-u+mod)%mod;
        }
    }
}
inline void AND(ll a[],int flag){
    for(int len=2;len<=n;len<<=1){
        int m=len>>1;
        for(int i=0;i<n;i+=len) rep(j,0,m-1){
            int u=a[i+j],v=a[i+j+m];
            a[i+j]=flag==1?(u+v)%mod:(u-v+mod)%mod;
        }
    }
}
inline void XOR(ll a[],int flag){
    for(int len=2;len<=n;len<<=1){
        int m=len>>1;
        for(int i=0;i<n;i+=len) rep(j,0,m-1){
            int u=a[i+j],v=a[i+j+m];
            a[i+j]=flag==1?(u+v)%mod:(u+v)%mod*inv%mod;
            a[i+j+m]=flag==1?(u-v+mod)%mod:(u-v+mod)%mod*inv%mod;
        }
    }
}
signed main(){
    n=1<<read(),inv=ksm(2,mod-2);
    rep(i,0,n-1) a[i]=x[i]=read();
    rep(i,0,n-1) b[i]=y[i]=read();
    OR(a,1),OR(b,1);
    rep(i,0,n-1) c[i]=a[i]*b[i]%mod,a[i]=x[i],b[i]=y[i];
    OR(c,-1);
    rep(i,0,n-1) printf("%lld ",c[i]);putchar('\n');
    AND(a,1),AND(b,1);
    rep(i,0,n-1) c[i]=a[i]*b[i]%mod,a[i]=x[i],b[i]=y[i];
    AND(c,-1);
    rep(i,0,n-1) printf("%lld ",c[i]);putchar('\n');
    XOR(a,1),XOR(b,1);
    rep(i,0,n-1) c[i]=a[i]*b[i]%mod;
    XOR(c,-1);
    rep(i,0,n-1) printf("%lld ",c[i]);putchar('\n');
}

坑填的差不多，大概就这样？