多项式初步

Author： RainAir
发布时间：February 28, 2019
2866views
4 comments
858 words
Categories：学习笔记

多项式算法学习笔记

定义

形如 $\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$。首先我们把这个式子进行一些化简。

$$ \begin{aligned} F(x)G(x)-1 &\equiv 0(\mod x^n) \\ F^2(x)G^2(x)+1-2F(x)G(x) &\equiv 0(\mod x^{2n}) \\ 2F(x)G(x)-F^2(x)G^2(x) &\equiv 1(\mod x^{2n}) \\ 2F(x)G(x)-F^2(x)G^2(x) &\equiv 1(\mod x^{2n}) \\ F(x)(2G(x)-F(x)G^2(x)) &\equiv 1(\mod x^{2n}) \\ \end{aligned} $$

不难发现最后一个式子的范围扩大到原来的 $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)$ 完成这个操作。
现在我们尝试对我们要构造的东西做一些改变。

$$ \begin{aligned} F(x) &= Q(x)G(x)+R(x) \\ F(\frac{1}{x}) &= Q(\frac{1}{x})G(\frac{1}{x})+R(\frac{1}{x}) \\ x^nF(\frac{1}{x}) &= x^{n-m}Q(\frac{1}{x})x^mG(\frac{1}{x})+x^{n-m+1}x^{m-1}R(\frac{1}{x}) \\ F_R(\frac{1}{x}) &= Q_R(\frac{1}{x})G_R(\frac{1}{x})+x^{n-m+1}R_R(\frac{1}{x}) \\ \end{aligned} $$

我们考虑两边对 $x^{n-m+1}$ 取模。

$$ \begin{aligned} F_R(\frac{1}{x}) &\equiv Q_R(\frac{1}{x})G_R(\frac{1}{x}) (\mod x^{n-m+1}) \\ Q_R(\frac{1}{x}) &\equiv F_R(\frac{1}{x})G_R^{-1}(\frac{1}{x}) \\ \end{aligned} $$

我们求一遍 $G_R$ 的逆，然后用 NTT 乘起来，就可以得到 $Q_R$ 了，最后
$$R(x) = F(x)-G(x)Q(x)$$
直接计算出来余数就可以了。
时间复杂度瓶颈在于多项式求逆，所以是 $O(nlogn)$。
代码就不贴了。就是上面的一些简单运用。

多项式求导/积分

都是非常简单的东西了。
$$f(x) = \sum_{i=0}^n a_ix^i$$
根据求导的加法法则和幂函数可以轻松得到：

$$ \begin{aligned} f'(x) &= \sum_{i=1}^n i\times a_ix^{i-1} \\ \int_l^r f(x)d(x) &= \sum_{i=1}^{n+1} \frac{a_{i-1}r^i}{i} - \sum_{i=1}^{n+1} \frac{a_{i-1}l^i}{i} \end{aligned} $$

补充一些同样适用于多项式的求导基本法则：
$(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)}$

多项式牛顿迭代

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

前置芝士1：Taylor 展开

Taylor 展开（泰勒展开）是一种用函数在某点的信息描述其附近取值的公式，也就是拟合某点及其周围的取值
具体来说我们从头开始构造这个拟合函数，设我们要使 $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$ 阶求导。
如果描述听不懂的话建议大家去看一下上面的视频

回到多项式牛顿迭代上

我们考虑使用和多项式求逆一样的 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})$，所以可以推一波式子：

$$ \begin{aligned} G(F_0(x))+G'(F_0(x))F(x)-G'(F_0(x))F_0(x) &= (\mod x^{2n}) \\ F(x) &= F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))} \end{aligned} $$

按照多项式求逆那样的套路迭代计算就可以了，我们来看一下一些经典的题目。

多项式开方

题目链接
给定一个多项式 $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)$，使用牛顿迭代法：

$$ \begin{aligned} B(x) & \equiv B_0(x)-\frac{B_0^2(x)-A(x)}{2B_0(x)} \\ & \equiv \frac{F_0^2(x)+A(x)}{2F_0(x)} \end{aligned} $$

然后就可以直接迭代计算了。时间复杂度 $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)}$
然后代入牛顿迭代公式，得到

$$ \begin{aligned} B(x) & \equiv B_0(x) - \frac{G(B_0(x))}{G'(B_0(x))} \\ B(x) & \equiv B_0(x) - B_0(x)(lnB_0(x)-A(x)) \\ B(x) & \equiv B_0(x)(1-lnB_0(x)+A(x)) \end{aligned} $$

因为 $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;
}

Last modification：October 3rd, 2020 at 08:21 am

如果觉得我的文章对你有用，请随意赞赏

4 comments

gxy
July 23rd, 2023 at 06:31 pm

Orz

Reply
菜鸡
January 14th, 2020 at 03:42 pm

orz

Reply
蔡鸡
October 4th, 2019 at 04:05 pm

神rainair

Reply
memset0
March 3rd, 2019 at 05:26 pm

rainair 太强啦！！！

Reply

多项式初步

RainAir • 2019 年 02 月 28 日

多项式算法学习笔记

定义

多项式加减

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

多项式乘法

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;
}

做法二：拉格朗日插值

#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;
}

扩展阅读：拉格朗日插值

多项式求逆

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.......

我们考虑两边对 $x^{n-m+1}$ 取模。

$$ \begin{aligned} F_R(\frac{1}{x}) &\equiv Q_R(\frac{1}{x})G_R(\frac{1}{x}) (\mod x^{n-m+1}) \\ Q_R(\frac{1}{x}) &\equiv F_R(\frac{1}{x})G_R^{-1}(\frac{1}{x}) \\ \end{aligned} $$

多项式求导/积分

都是非常简单的东西了。
$$f(x) = \sum_{i=0}^n a_ix^i$$
根据求导的加法法则和幂函数可以轻松得到：

$$ \begin{aligned} f'(x) &= \sum_{i=1}^n i\times a_ix^{i-1} \\ \int_l^r f(x)d(x) &= \sum_{i=1}^{n+1} \frac{a_{i-1}r^i}{i} - \sum_{i=1}^{n+1} \frac{a_{i-1}l^i}{i} \end{aligned} $$

多项式牛顿迭代

前置芝士1：Taylor 展开

回到多项式牛顿迭代上

$$ \begin{aligned} G(F_0(x))+G'(F_0(x))F(x)-G'(F_0(x))F_0(x) &= (\mod x^{2n}) \\ F(x) &= F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))} \end{aligned} $$

按照多项式求逆那样的套路迭代计算就可以了，我们来看一下一些经典的题目。

多项式开方

$$ \begin{aligned} B(x) & \equiv B_0(x)-\frac{B_0^2(x)-A(x)}{2B_0(x)} \\ & \equiv \frac{F_0^2(x)+A(x)}{2F_0(x)} \end{aligned} $$

然后就可以直接迭代计算了。时间复杂度 $O(nlogn)$。

多项式求 ln

多项式求 exp

$$ \begin{aligned} B(x) & \equiv B_0(x) - \frac{G(B_0(x))}{G'(B_0(x))} \\ B(x) & \equiv B_0(x) - B_0(x)(lnB_0(x)-A(x)) \\ B(x) & \equiv B_0(x)(1-lnB_0(x)+A(x)) \end{aligned} $$

因为 $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;
}

多项式算法学习笔记

定义

多项式加减

多项式乘法

确定多项式

做法一：高斯消元

做法二：拉格朗日插值

多项式求逆

多项式除法/取模

多项式求导/积分

多项式牛顿迭代

前置芝士1：Taylor 展开

回到多项式牛顿迭代上

多项式开方

多项式求 ln

多项式求 exp

4 comments