多项式算法学习笔记

定义

形如 $\sum a_ix^{i}$ 的式子被称为多项式。定义多项式中每个单项式叫做这个多项式的项，定义多项式的次数为项的最高次数。显然 $(n+1)$ 个点可以唯一确定一个次数为 $n$ 的多项式。
接下来我们首先来看一下多项式的基本运算怎么实现（要不然谁用啊）

多项式加减

直接对齐相同次数然后系数直接加减就可以了吧。。。时间复杂度 $O(n)$ ，难度：普及组 T1（虽然这个级别已经不存在了233）

多项式乘法

首先我们考虑一下我们平时如何手算多项式算法。显然我们都是 $O(n^2)$ 暴力算的，所以就有了一个枚举所有可能的二元组然后统计答案的暴力做法。
当然在常规题目中我们会使用一种叫做 FFT（快速傅里叶变换）的 $O(nlogn)$ 的优秀算法来计算，详情请查看这篇文章。
哈哈你一定学完 FFT 然后回来了！我们来考虑一下如果多项式乘法在模意义下怎么做。
有一种叫做 NTT （快速数论变换）的东西，它把每次我们做快速傅里叶变换时用到的 $\omega$ 重新定义了一下。大概是令 $g$ 为模数 $p$ 的原根，我们就可以取 $w_n = g^{\frac{p-1}{n}}$ 。然后 IDFT 的时候直接乘逆元就可以了。
至于这个原根怎么找。。。一般你碰到模数为 $998244353$ 的时候就想下可不可以 NTT，因为 $998244353$ 是一个常用的 NTT 模数，一个原根是 $3$ 。
中心代码大致如下：

inline void NTT(int *A,int limit,int opt){
    FOR(i,0,limit-1) if(i < r[i]) std::swap(A[i],A[r[i]]);
    for(int mid = 1;mid < limit;mid <<= 1){
        int W = qpow(3,(ha-1)/(mid<<1));
        if(opt == -1) W = qpow(W,ha-2);
        for(int j = 0;j < limit;j += (mid<<1)){
            for(int k = 0,w = 1;k < mid;k++,w = w*W%ha){
                int x = A[j+k],y = w*A[j+mid+k]%ha;
                A[j+k] = (x+y)%ha;A[j+mid+k] = (x-y+ha)%ha;
            }
        }
    }
    if(opt == -1){ // IDFT
        int inv = qpow(limit,ha-2);
        FOR(i,0,limit-1) A[i] = A[i]*inv%ha;
    }
}

好了接下来才是关于多项式的一些基本的稍微硬核一点的操作。

确定多项式

上面提到了 $(n+1)$ 个点可以唯一确定一个次数为 $n$ 的多项式。
现在给定 $n$ 个点，请你确定这个多项式并且将 $k$ 代入求值。

做法一：高斯消元

我们暴力建方程然后暴力高斯消元！ $O(n^3)$ 。
什么你连高斯消元都不会？题目链接
其实就是暴力解方程啦~奉上高斯消元模板题的代码：

#include <algorithm>
#include <iostream>
#include <cstring>
#include <climits>
#include <cstdio>
#include <vector>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <queue>
#include <stack>
#include <map>
#include <set>

#define fi first
#define lc (x<<1)
#define se second
#define U unsigned
#define rc (x<<1|1)
#define Re register
#define LL long long
#define MP std::make_pair
#define CLR(i,a) memset(i,a,sizeof(i))
#define FOR(i,a,b) for(Re int i = a;i <= b;++i)
#define ROF(i,a,b) for(Re int i = a;i >= b;--i)
#define SFOR(i,a,b,c) for(Re int i = a;i <= b;i+=c)
#define SROF(i,a,b,c) for(Re int i = a;i >= b;i-=c)
#define DEBUG(x) std::cerr << #x << '=' << x << std::endl

const double EPS = 1e-7;
const int MAXN = 200+5;
double d[MAXN][MAXN],ans[MAXN];
int n;

int main(){
    scanf("%d",&n);
    FOR(i,1,n) FOR(j,1,n+1) scanf("%lf",&d[i][j]);
    FOR(i,1,n){
        int r = i;
        FOR(j,i+1,n) if(std::fabs(d[r][i]) < std::fabs(d[j][i])) r = j;
        if(std::fabs(d[r][i]) < EPS){
            puts("No Solution");
            return 0;
        }
        if(i != r) std::swap(d[r],d[i]);
        double div = d[i][i];
        FOR(j,i,n+1) d[i][j] /= div;
        FOR(j,i+1,n){
            div = d[j][i];
            FOR(k,i,n+1) d[j][k] -= d[i][k]*div;
        }
    }
    ans[n] = d[n][n+1];
    ROF(i,n-1,1){
        ans[i] = d[i][n+1];
        FOR(j,i+1,n) ans[i] -= (d[i][j]*ans[j]);
    }
    FOR(i,1,n) printf("%.2f\n",ans[i]);
    return 0;
}

做法二：拉格朗日插值

背下式子：
设该多项式为 $f(x)$ ，第 $i$ 个点的坐标为 $(x_i,y_I)$ ，需要求在 $k$ 点的取值。

$f(k) = \sum_{i=0}^n y_i \prod_{i \neq j}\frac{k-x_j}{x_i-x_j}$

证明显然：因为这个式子满足了如果代入 $x_i$ 就肯定会返回 $y_i$ 。根据 $n+1$ 个点唯一确定一个多项式，就可以知道这个多项式等价于原多项式了。

#include <algorithm>
#include <iostream>
#include <cstring>
#include <climits>
#include <cstdio>
#include <vector>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <queue>
#include <stack>
#include <map>
#include <set>

#define fi first
#define lc (x<<1)
#define se second
#define U unsigned
#define rc (x<<1|1)
#define Re register
#define LL long long
#define MP std::make_pair
#define CLR(i,a) memset(i,a,sizeof(i))
#define FOR(i,a,b) for(Re int i = a;i <= b;++i)
#define ROF(i,a,b) for(Re int i = a;i >= b;--i)
#define SFOR(i,a,b,c) for(Re int i = a;i <= b;i+=c)
#define SROF(i,a,b,c) for(Re int i = a;i >= b;i-=c)
#define DEBUG(x) std::cerr << #x << '=' << x << std::endl

const int MAXN = 2000+5;

const int ha = 998244353;

inline int qpow(int a,int n=ha-2){
    int res = 1;
    while(n){
        if(n & 1) res = 1ll*res*a%ha;
        a = 1ll*a*a%ha;
        n >>= 1;
    }
    return res;
}

int n,k,x[MAXN],y[MAXN],ans;

int main(){
    scanf("%d%d",&n,&k);
    FOR(i,1,n) scanf("%d%d",x+i,y+i);
    FOR(i,1,n){
        int s1 = 1,s2 = 1;
        FOR(j,1,n){
            if(j == i) continue;
            s1 = 1ll*s1*((1ll*k-x[j]+ha)%ha)%ha;s2 = 1ll*s2*((1ll*x[i]-x[j]+ha)%ha)%ha;
        }
        s2 = qpow(s2);s1 = 1ll*s1*s2%ha*y[i]%ha;
        ans = (ans+s1)%ha;
    }
    printf("%d\n",ans);
    return 0;
}

扩展阅读：维基百科-拉格朗日插值

多项式求逆

题目链接
给定一个多项式 $F(x)$ ，让你求出一个多项式 $G(x)$ ，满足 $F(x) * G(x) \equiv 1(\mod x^n)$
我们可以考虑倍增地构造多项式 $G$ 。首先我们把这个式子进行一些化简。
$图片说明$
不难发现最后一个式子的范围扩大到原来的 $2$ 倍，并且又凑成了一开始的形式。
我们考虑迭代处理，首先当规模为 $1$ 时显然答案系数就是对应系数的逆元。
我们设对于规模为 $n$ 的问题得到的多项式为 $G_0$ ，我们想要推出扩大一倍规模的答案 $G_1$ ，根据上面的推导，可以有

$G_1 = 2G(x)-F(x)G^2(x) = G(x)(2-F(x)G(x))$

所以我们直接两遍 NTT 就能扩大一倍的问题规模了。
时间复杂度： $T(n) = 2T(n/2) + n log n \Rightarrow O(nlogn)$ 。
这里就只给迭代写法啦 qwq 貌似迭代确实比递归快很多呢。

inline void Inv(int *A,int *B,int n){
    B[0] = qpow(A[0],ha-2);
    int len;
    for(len = 1;len < (n<<1);len <<= 1){
        int d = len<<1;
        FOR(i,0,len-1) a[i] = A[i],b[i] = B[i];
        FOR(i,0,d-1) pos[i] = (pos[i>>1]>>1)|((i&1)?len:0);
        NTT(a,d,1);NTT(b,d,1);
        FOR(i,0,d-1) B[i] = ((2-a[i]*b[i]%ha+ha)%ha*b[i]%ha)%ha;
        NTT(B,d,-1);
        FOR(i,len,d-1) B[i] = 0;
    }
    FOR(i,0,len-1) a[i] = b[i] = 0;
    FOR(i,n,len-1) B[i] = 0;
}

多项式除法

题目链接
给定一个 $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)$ 。
显然如果能整除的话我们直接求个逆就可以了。那现在有余数怎么办啊 Orz.......
别急，我们首先定义一种操作 R(reverse)，满足： $A_R(x) = x^nA(\frac{1}{x})$ 。
不难看出 $R$ 操作的本质是多项式翻转。所以我们就可以 $O(n)$ 完成这个操作。
现在我们尝试对我们要构造的东西做一些改变。
$图片说明$
我们考虑两边对 $x^{n-m+1}$ 取模。
$图片说明$
我们求一遍 $G_R$ 的逆，然后用 NTT 乘起来，就可以得到 $Q_R$ 了，最后 $R(x) = F(x)-G(x)Q(x)$ 直接计算出来余数就可以了。
时间复杂度瓶颈在于多项式求逆，所以是 $O(nlogn)$ 。
代码就不贴了。就是上面的一些简单运用。
多项式求导/积分
都是非常简单的东西了。

$f(x) = \sum_{i=0}^n a_ix^i$

根据求导的加法法则和幂函数可以轻松得到：
$图片说明$
补充一些同样适用于多项式的求导基本法则：

$(cA(x))' = cA'(x)$
$(A(x)+B(x))' = A'(x)+B'(x)$
$(A(x)B(x))' = A'(x)B(x)+A(x)B'(x)$
$(\frac{A(x)}{B(x)})' = \frac{A'(x)B(x)-A(x)B'(x)}{B^2(x)}$
写这些主要是为了介绍后面的牛顿迭代 $qwq$

多项式牛顿迭代

处理多项式复合的高级神器 Orz
问题大概是这样的：给出 $G(x)$ ，求 $F(x)$ 满足 $G(F(x)) \equiv 0 (\mod x^n)$ 。
首先为了让像我这样萌萌哒的弱校初中选手也能看懂，我决定放一些前置芝士。

前置芝士1：Taylor 展开

这个大家都会 qwq。

Taylor 展开（泰勒展开）是一种用函数在某点的信息描述其附近取值的公式，也就是拟合某点及其周围的取值啦qwq。
具体来说我们从头开始构造这个拟合函数，设我们要使 $P$ 及其周围一段区间拟合这个原函数。首先我们可以使得我们构造出的函数和待拟合的函数在 $P$ 的取值相等。但是发现这只有一个值？一点也不像怎么办？我们可以使得他们的二阶倒数，三阶导数.....一直相等下去，然后我们就会发现这两个函数在那附近神奇的重合了！我们看个 GIF 来详细的了解这个过程。

所以我们就可以得到 Taylor 展开来拟合多项式的式子了：

$f(x) \approx g(x) = g(0) + \sum_{i=1} \frac{f^{(i)}(0)}{i!}x^i$

其中 $f^{(n)}(x)$ 表示对原函数图像 $x$ 这个点进行 $n$ 阶求导。
如果描述听不懂的话建议大家去看一下上面的视频的...毕竟视频里有图还有可爱的 $\pi$ ~

回到多项式牛顿迭代上

我们考虑使用和多项式求逆一样的 trick：倍增来构造答案多项式。
首先当 $n=1$ 时， $G(F(x)) \equiv 0 (\mod x)$ ，这个就是要单独求出来的。
现在我们假设已经求出了

$G(F_0(x)) \equiv 0(\mod x^n)$

考虑如何扩展到 $x^{2n}$ 下，我们把 $G(F(x))$ 放到 $F_0(x)$ 上做 Taylor 展开

$G(F(x)) = G(F_0(x)) + \sum_{i=1} \frac{G^{(i)}}{i!}(F(x)-F_0(x))^i$

因为倍增构造，我们可以发现 $F(x)$ 和 $F_0(x)$ 后 $n$ 项相同，所以 $(F(x)-F_0(x))^2$ 的最低非 $0$ 项次数大于 $2n$ 。所以可以得到

$G(F(x)) \equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x)) (\mod x^{2n})$

然后因为我们要求的 $G(F(x))$ 需要满足 $G(F(x)) \equiv 0(mod x^{2n})$ ，所以可以推一波式子：
$图片说明$
按照多项式求逆那样的套路迭代计算就可以了，我们来看一下一些经典的题目。

多项式开方

题目链接
给定一个多项式 $A(x)$ ，要求求出一个多项式 $B(x)$ ，满足 $B^2(x) \equiv A(x)(\mod x^n)$
简单的变换： $B^2(x)-A(x) \equiv 0(\mod x^n)$ 。
设 $G(B(x)) = B^2(x)-A(x)$ ，则 $G'(B(x)) = 2B(x)$ ，使用牛顿迭代法：
$图片说明$
然后就可以直接迭代计算了。时间复杂度 $O(nlogn)$ 。

多项式求 ln

题目链接
给定一个多项式 $A(n)$ ，求一个多项式 $B(n)$ ，满足 $B(n) = ln A(n)$ 。
我们设 $G(x) = F(A(x))$ ，定义 $F(x) = ln(x)$ 。
那么显然 $G'(x) = F'(A(x))A'(x) = \frac{A'(x)}{A(x)}$ 。
那么求逆就可以算出 $G'$ 。然后因为求导和不定积分互逆，所以对 $G'$ 求一个不定积分就可以了。
代码不贴了（主要是后面有大杂烩）

多项式求 exp

题目链接
多项式大杂烩，练习之前学的怎么样了qwq
给定一个多项式 $A(x)$ ，要你求出一个 $B(x)$ ，满足 $B(x) \equiv e^{A(x)} (\mod x^n)$ 。
我们考虑计算出

$B(x) \equiv e^{A(x)} (\mod x^n)$

变形一下可得

$lnB(x)-A(x) \equiv 0 (\mod x^n)$

设 $G(B(x)) = lnB(x)-A(x)$ ，然后直接套用牛顿迭代方法：
首先求导：得到 $G'(B(x)) = \frac{1}{B(x)}$
然后代入牛顿迭代公式，得到
$图片说明$
因为 $A(0)=0$ ，所以 $B(x)$ 的常数项是 $1$ 。（有些牛顿迭代题目要用二次剩余算系数但是我还没接触过QAQ）
然后就是板子大全了qwq
放一下代码吧：

#include <algorithm>
#include <iostream>
#include <cstring>
#include <climits>
#include <cstdio>
#include <vector>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <queue>
#include <stack>
#include <map>
#include <set>

#define Re register
#define LL long long
#define U unsigned
#define FOR(i,a,b) for(Re int i = a;i <= b;++i)
#define ROF(i,a,b) for(Re int i = a;i >= b;--i)
#define SFOR(i,a,b,c) for(Re int i = a;i <= b;i+=c)
#define SROF(i,a,b,c) for(Re int i = a;i >= b;i-=c)
#define CLR(i,a) memset(i,a,sizeof(i))
#define BR printf("--------------------\n")
#define DEBUG(x) std::cerr << #x << '=' << x << std::endl
#define int LL
const int MAXN = 5e5 + 5;
const int ha = 998244353;

inline int qpow(int a,int n){
    int res = 1;
    while(n){
        if(n & 1) res = res*a%ha;
        a = a*a%ha;
        n >>= 1;
    }
    return res;
}

int a[MAXN],b[MAXN],c[MAXN],d[MAXN],r[MAXN],g[MAXN],f[MAXN],x[MAXN],n;

inline void NTT(int *A,int limit,int opt){
    FOR(i,0,limit-1) if(i < r[i]) std::swap(A[i],A[r[i]]);
    for(int mid = 1;mid < limit;mid <<= 1){
        int W = qpow(3,(ha-1)/(mid<<1));
        if(opt == -1) W = qpow(W,ha-2);
        for(int j = 0;j < limit;j += (mid<<1)){
            for(int k = 0,w = 1;k < mid;k++,w = w*W%ha){
                int x = A[j+k],y = w*A[j+mid+k]%ha;
                A[j+k] = (x+y)%ha;A[j+mid+k] = (x-y+ha)%ha;
            }
        }
    }
    if(opt == -1){
        int inv = qpow(limit,ha-2);
        FOR(i,0,limit-1) A[i] = A[i]*inv%ha;
    }
}

inline void Inv(int *A,int *B,int n){
    B[0] = qpow(A[0],ha-2);
    int len;
    for(len = 1;len < (n<<1);len <<= 1){
        int d = len<<1;
        FOR(i,0,len-1) a[i] = A[i],b[i] = B[i];
        FOR(i,0,d-1) r[i] = (r[i>>1]>>1)|((i&1)?len:0);
        NTT(a,d,1);NTT(b,d,1);
        FOR(i,0,d-1) B[i] = (2-a[i]*b[i]%ha+ha)%ha*b[i]%ha;
        NTT(B,d,-1);FOR(i,len,d-1) B[i] = 0;
    }
    FOR(i,0,len-1) a[i] = b[i] = 0;
    FOR(i,n,len-1) B[i] = 0;
}

inline void Direv(int *A,int *B,int n){
    FOR(i,1,n-1) B[i-1] = i*A[i]%ha;B[n-1] = 0;
}

inline void Inter(int *A,int *B,int n){
    FOR(i,1,n-1) B[i] = A[i-1]*qpow(i,ha-2)%ha;B[0] = 0;
}

inline void getln(int *A,int *B,int n){
    int *C = c,*D = d;
    Direv(A,C,n);Inv(A,D,n);
    int len = 0,limit = 1;
    while(limit <= n) limit <<= 1,len++;
    FOR(i,1,limit-1) r[i] = (r[i>>1]>>1)|((i&1)<<(len-1));
    NTT(C,limit,1);NTT(D,limit,1);
    FOR(i,0,limit-1) C[i] = C[i]*D[i]%ha;
    NTT(C,limit,-1);Inter(C,B,n);
}

inline void Exp(int *A,int *B,int n){
    B[0] = 1;int len;
    for(len = 1;len < (n<<1);len <<= 1){
        int d = len<<1;
        getln(B,f,len);
        FOR(i,0,len-1) f[i] = (A[i]+(i==0)-f[i]+ha)%ha;
        FOR(i,0,d-1) r[i] = (r[i>>1]>>1)|((i&1)?len:0);
        NTT(B,d,1);NTT(f,d,1);
        FOR(i,0,d-1) B[i] = B[i]*f[i]%ha;
        NTT(B,d,-1);
        FOR(i,len,d-1) B[i] = 0;
    }
    FOR(i,0,len-1) f[i] = 0;FOR(i,n,len-1) B[i] = 0;
}

signed main(){
    scanf("%lld",&n);
    FOR(i,0,n-1) scanf("%lld",x+i);
    int len;for(len=1;len<=n;len<<=1);Exp(x,g,len);
    FOR(i,0,n-1) printf("%lld ",g[i]);
    return 0;
}

多项式相关简单入门