多项式

“虚妄与实冢的交错，逆塑快速变换的坤仑。”

复数

形如 $a+bi,i=\sqrt{-1}$ 的数，其中 $a$ 为该复数的实根， $b$ 为该复数的虚根。一个复数 $x=a+bi$ 可以表示为一个向量 $(a,b)$ 。这样的数称为复数（$\text{complex}$）。

复数意义上的单位根

指的是坐标轴上的一个单位圆，被等分成 $n$ 份，有函数 $\omega_{k}^{n}$ 表示把这个单位圆等分成 $n$ 份取 $k$ 份，每一份大小为 $\frac{2\pi}{n}k$ ，被称为 $n$ 次单位根。

有性质：

$\omega^{i}_n\ne\omega^j_n,\forall i\ne j$ ；
$\omega^k_n=\cos \frac{2k\pi}{n}+i\sin\frac{2k\pi}{n}$ ；
$\omega_n^0=\omega_n^n=1$ ；
$\omega^{2k}_{2n}=\omega^{k}_n$ ；
$\omega_n^{k+\frac{\pi}{2}}=-\omega_n^k$ ；
$\omega_n^i\times\omega_n^j=\omega_n^{i+j}$ ；
$\omega_n^k=-\omega_n^{-k}$ 。

定义

形如：

$A(x)=a_0+a_1x+a_2x^2+...a_nx^n \nonumber$

的函数式。

性质

Ⅰ

用该多项式函数上任意 $n+1$ 个不同点，均可唯一一个 $n$ 个多项式。

这条性质也被称为是多项式的点表示法。

宏定义

幂级数

可列项相加的求和式称为级数，而对于和式 $\sum_{n=0}^{\infty}a_nx^n$ 中，每项均为非负整数次幂函数乘常数系数，这样的形式级数称为幂级数。

多项式环

在群论中有一个定义为环，而对于一般环 $R$ ，定义 $R$ 上的多项式环（$\text{polynomial ring}$）为 $R[x]$ 。

每个元素 $f$ 称为 $R$ 上的多项式，表示为：

$f=\langle f_0,f_1,f_2\dots f_n\rangle,(f_0,f_1,f_2,\dots,f_n\in R)$

直接将多项式定义为系数序列，表示为：

$f(x)=f_0+f_1x+f_2x^2+\dots+f_nx^n$

这也是我们熟知的多项式定义，这里设 $x$ 是一个形式符号，一个对系数位置的标识符。

当然，这里的多项式也包括了无穷项的存在：

$f(x)=f_0+f_1x+f_2x^2+\dots$

则得到形式幂级数环（$\text{formal power series ring}$）$R[[x]]$ ，其中每一个元素 $f$ 被称为形式幂级数（$\text{formal power series}$），简称为幂级数。

度

对于一个多项式 $f(x)$ ，其最高次项的次数被称为该多项式的度（$\text{Degree}$），也就是所说的次数，记作 $\deg f$ 。

运算

加减

$f(x)\pm g(x)=\sum_{i=0}a_ix\pm b_i x$ 。

乘法

首先定义多项式 $f(x),g(x)$ ：

$f(x)=a_0+a_1x+…+a_nx^n$

$g(x)=b_0+b_0x+…b_mx^m$

则有：

$Q(x)=f(x)\cdot g(x)=\sum\limits_{i=0}^n\sum\limits_{j=0}^{m}a_ib_jx^{i+j}=c_0+c_1x+…+c_{n+m}x^{n+m}$

多项式或是幂级数的乘法满足结合律，其加法也满足分配律。而如果 $R$ 为交换环或幺环，乘法则具有交换律和单位元。

而所谓的 $\text{FFT}$ 能够在 $\mathcal O(n2^n)$ 的时间内计算两个 $2^n$ 次多项式的乘积。

拉格朗日插值

这种东西可以在 $\mathcal O(n)$ 的时间内进行插值并得到一个 $n-1$ 次的多项式。

给出 $n$ 个点 $(x_i,y_i)$ ，求出这 $n$ 个点所在的最多 $n-1$ 次的多项式 $f(x)$ 中 $f(k)$ 的值。

差分法和待定系数法

详见 $\text{OI-Wiki}$ 。

插值法

给出表达式（具体推导我也不会）：

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

如果求多项式表达式的话，就需要 $\mathcal O(n^2)$ 来枚举，但如果只是求值的话，可以选择代入。

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

然后还需要求解逆元（用费马小定理），所以时间还是在 $\mathcal O(n^2)$ 。

拉格朗日插值之一

AC Code

#include<bits/stdc++.h>
#define re register
typedef long long ll;
template<class T>
inline void read(T &x)
{
    x=0;
    char ch=getchar(),t=0;
    while(ch<'0'||ch>'9') t|=ch=='-',ch=getchar();
    while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    if(t) x=-x;
}
template<class T,class ...T1>
inline void read(T &x,T1 &...x1)
{
    read(x),read(x1...);
}
template<class T>
inline void write(T x)
{
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10+'0');
}
template<>
inline void write(char c)
{
    putchar(c);
}
template<>
inline void write(char *s)
{
    while(*s) putchar(*s++);
}
template<class T,class ...T1>
inline void write(T x,T1 ...x1)
{
    write(x),write(x1...);
}
template<class T>
inline bool checkMax(T &x,T &y)
{
    return x<y?x=y,1:0;
}
template<class T>
inline bool checkMin(T &x,T &y)
{
    return x>y?x=y,1:0;
}
const int MAXN=2e3+10;
const ll Mod=998244353;
ll N,K,Dx[MAXN],Dy[MAXN],ans,s1,s2;
inline ll qPow(ll a,ll b,ll p)
{
    ll res=1ll;
    while(b)
    {
        if(b&1) res=res*a%p;
        a=a*a%p;b>>=1;
    }
    return res;
}
inline ll inv(ll x)
{
    return qPow(x,Mod-2,Mod);
}
int main()
{
    // freopen("poly.in","r",stdin);
    // freopen("poly.out","w",stdout);
    read(N,K);
    for(int i=1;i<=N;++i) read(Dx[i],Dy[i]);
    for(int i=1;i<=N;++i)
    {
        s1=Dy[i]%Mod,s2=1;
        for(int j=1;j<=N;++j)
            if(i!=j) s1=s1*(K-Dx[j])%Mod,s2=s2*(Dx[i]-Dx[j])%Mod;
        ans+=s1*inv(s2)%Mod;
    }
    write((ans%Mod+Mod)%Mod);
    return 0;
}
/*

*/

线性插值

在一些特殊情况下，拉格朗日插值可以做到线性时间 $\mathcal O(n)$ 。需要满足已知点的横坐标为整数。

即当我们知道一个 $n$ 次多项式为 $f(x)$ ，且已知 $f(i),1\leq i\leq n,i\in\mathbb N$ ，就可以优化插值。

首先化简式子得到：

$\displaystyle f(x)=\sum_{i=1}^{n+1}y_i\prod_{i\neq j}\frac{x-x_j}{x_i-x_j}=\sum_{i=1}^{n+1}y_i\prod_{i\neq j}\frac{x-j}{i-j}$

然后对整个分式进行分别考虑，分子为：

$\displaystyle\frac{\prod_{j=1}^{n+1}(x-j)}{x-i}$

分母为：

$(-1)^{n+1-i}\times(i-1)!\times(n+1-i)!$

然后结合插值公式得到：

$\displaystyle f(x)=\sum_{i=1}^{n+1}y_i\times\frac{\prod_{j=1}^{n+1}(x-j)}{(x-i)\times(-1)^{n+1-i}\times(i-1)!\times(n+1-i)!}$

处理阶乘前缀和和阶乘逆，前后缀积，枚举一遍 $y_i$ ，求得结果，时间复杂度 $\mathcal O(n)$ 。

例题

CF622F

有一个结论，就是 $\displaystyle\sum_{i=1}^{n}i^k$ 是一个 $k+1$ 次的多项式，我不知道怎么证的，反正就是有这个结论。

然后捏，就线性筛出 $o^i,1\leq o\leq k+2$ 的值，得到 $k+2$ 个点，进行拉格朗日插值就可以在 $\mathcal O(n)$ 时间内计算出答案了。

AC Code

#include<bits/stdc++.h>
#define re register
typedef long long ll;
template<class T>
inline void read(T &x)
{
    x=0;
    char ch=getchar(),t=0;
    while(ch<'0'||ch>'9') t|=ch=='-',ch=getchar();
    while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    if(t) x=-x;
}
template<class T,class ...T1>
inline void read(T &x,T1 &...x1)
{
    read(x),read(x1...);
}
template<class T>
inline void write(T x)
{
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10+'0');
}
template<>
inline void write(char c)
{
    putchar(c);
}
template<>
inline void write(char *s)
{
    while(*s) putchar(*s++);
}
template<class T,class ...T1>
inline void write(T x,T1 ...x1)
{
    write(x),write(x1...);
}
template<class T>
inline bool checkMax(T &x,T &y)
{
    return x<y?x=y,1:0;
}
template<class T>
inline bool checkMin(T &x,T &y)
{
    return x>y?x=y,1:0;
}
const int MAXN=1e6+10;
const int Mod=1e9+7;
int N,K,ans;
int Pre[MAXN],Suf[MAXN],Inv[MAXN],Fac[MAXN];
int pri[MAXN],Tot,f[MAXN];
bool is[MAXN];
inline int qPow(int a,int b,int p)
{
    int res=1;
    while(b)
    {
        if(b&1) res=1ll*res*a%p;
        a=1ll*a*a%p;b>>=1;
    }
    return res;
}
inline void sieve(int limit)
{
    f[1]=1;
    for(int i=2;i<=limit;++i)
    {
        if(!is[i]) pri[++Tot]=i,f[i]=qPow(i,K,Mod);
        for(int j=1;j<=Tot&&1ll*i*pri[j]<=limit;++j)
        {
            is[i*pri[j]]=1;
            f[i*pri[j]]=1ll*f[i]*f[pri[j]]%Mod;
            if(i%pri[j]==0) break;
        }
    }
    for(int i=2;i<=limit;++i) f[i]=(f[i-1]+f[i])%Mod;
}
int main()
{
    // freopen("poly.in","r",stdin);
    // freopen("poly.out","w",stdout);
    read(N,K);
    sieve(K+2);
    if(N<=K+2) return printf("%d",f[N])&0;
    Pre[0]=Suf[K+3]=1;
    for(int i=1;i<=K+2;++i) Pre[i]=1ll*Pre[i-1]*(N-i)%Mod;
    for(int i=K+2;i>=1;--i) Suf[i]=1ll*Suf[i+1]*(N-i)%Mod;
    Fac[0]=Inv[0]=Fac[1]=Inv[1]=1;
    for(int i=2;i<=K+2;++i)
    {
        Fac[i]=1ll*Fac[i-1]*i%Mod;
        Inv[i]=1ll*(Mod-Mod/i)*Inv[Mod%i]%Mod;
    }
    for(int i=2;i<=K+2;++i) Inv[i]=1ll*Inv[i-1]*Inv[i]%Mod;
    for(int i=1;i<=K+2;++i)
    {
        int P=1ll*Pre[i-1]*Suf[i+1]%Mod;
        int Q=1ll*Inv[i-1]*Inv[K+2-i]%Mod;
        int Mul=((K+2-i)&1)?-1:1;
        ans=(ans+1ll*(Q*Mul+Mod)%Mod*P%Mod*f[i]%Mod)%Mod;
    }
    write(ans);
    return 0;
}
/*
4 3
*/

快速傅里叶变换

俗称为 $\text{FFT}$ ，用于处理多项式乘法。

在满足性质：

用该多项式函数上任意 $n+1$ 个不同点，均可唯一一个 $n$ 个多项式。

是使用 $\text{FFT}$ 的前提与条件。

例题

求解 $C(x)=A(x)B(x)$ ，其中 $A(x)$ 为 $n$ 次式，$B(x)$ 为 $m$ 次式。

首先，这个问题在现实中的计算时间应该是 $\mathcal O(nm)$ 的，这个不必赘述。而 $\text{FFT}$ 能将复杂度降低到 $\mathcal O(n\log n)$ 。

正变换：系数表示法 -> 点表示法

利用点表示法，我们需要选择 $n+m+1$ 个点，表示为 $x_1,x_2…x_{n+m+1}$ ，因为有性质 $C(x_i)=A(x_i)B(x_i),\forall i\in A,B$ ，则我们求解 $n+m+1$ 个特殊点，然后通过点表示法转化回系数表示法。

有性质决定，我们选择的特殊点是复数的单位根。

有 $A(x)=a_0+a_1x+a_2x^2…+a_{n-1}x^{n-1}$ ，我们需要快速求出 $A(\omega_{k}^{n}),k\in[0,n-1]$ ，则是这道题的关键。

我们把 $A(x)$ 分开奇偶，表示为：

$A(x)=(a_0+a_2x^2+…a_{n-2}x^{n-2})+(a_1x+a_3x^{3}+…a_{n-1}x^{n-1})$

并将前者用 $A_1(x)$ 表示，后者用 $A_2(x)$ 来表示。

对于 $A_1(x)$ ，我们用 $x^2$ 换元 $x$ ，得到：

$A_1(x)=a_0+a_2x+a_4x^2+…+a_{n-2}x^{\frac{n}{2}-1}$

同样有：

$A_2(x)=a_1+a_3x+a_5x^2+…+a_{n-1}x^{\frac{n}{2}-1}$

带回原式则有：

$A(x)=A_1(x^2)+xA_2(x^2)$

带回 $\omega_k^n$ ，其中 $k\in [0,\frac{n}{2}-1]$ 。则有：

$A(\omega^k_n)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})=A_1(\omega_{\frac{n}{2}}^k)+\omega_n^kA_2(\omega_{\frac{n}{2}}^k)$

同理得到：（$k+n$ 相当于是 $k$ ，因为 $n$ 就表示为一圈）

$A(\omega_n^{k+\frac{n}{2}})=A_1(\omega^{k}_{\frac{n}{2}})-\omega_n^kA_2(\omega_{\frac{n}{2}}^k)$

那么，对于一个 $n$ 次式的 $A(x)$ ，我们把其拆分为 $2$ 个 $\frac{n}{2}$ 的子多项式，从而求出 $A(x)$，那么最后直到拆分成 $0$ 次式，最后递归操作，时间复杂度为 $\mathcal O(n\log n)$ 。

反变换：点表示法 -> 系数表示法

设现在我们已经求出了 $n+1$ 个点 $(\omega_n^k,A(w_n^k))$ ，我们设原多项式为：

$A(x)=c_0+c_1x+…c_{n-1}x^{n-1}$

存在定理：$c_k=\sum\limits_{i=0}^{n-1}y_i(\omega_n^{-k})^i$

转化未知数得到：$B(x)=y_0+y_1x+…y_{n-1}x^{n-1}$ ，则有 $c_k=B(\omega_n^{-k})$ 。（这个 $B(x)$ 是我们相乘之后得到的答案的点表示法）所以我们需要的是通过反变换求出 $B(\omega_n^0)\sim B(\omega_n^{-(n-1)})$ 即是答案。然后根据单位根的性质七，我们实际上就是需要求解 $B(\omega_n^0)\sim B(\omega_n^{n-1})$ 。

现在我们再次组织 $c_k$ ：

$\begin{align} c_k \nonumber &= \sum\limits_{i=0}^{n-1}y_i(\omega_n^{-k})^i \nonumber \\&= \sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega^i_n)^j)(\omega_n^{-k})^i \nonumber \\&=\ \sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i) \nonumber \\&= na_k \nonumber \end{align}$

最后得到 $c_k=na_k$ ，从而得到 $a_k=\frac{c_k}{n}$ ，得到答案中的系数。

迭代法的变迁

实际上，我们不需要使用递归来实现多项式的正变换，而是可以使用递推的方式来解决。

有一个性质：对于系数 $i$ ，它在正变换的末尾（递归的结尾部分）的标号会将其二进制位取反：

例如：对于 $8$ 次式 $a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7,a_8$ ；

变换之后的系数顺序是 $a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7$ 。

用二进制表示为：

$000,001,010,011,100,101,110,111$

和

$000,100,010,110,001,101,011,111$

这样就浅显易懂了。从而，在实现时，我们可以预处理出当前系数会跳到的位置，然后从底端递推到原式以减少迭代消耗的时间。

P3803 【模板】多项式乘法（FFT）

AC Code

typedef long long ll;
const int MAXN=3e6+1;
const double Pi=acos(-1);
template<class T>
inline void read(T &x)
{
    x=0;
    char ch=gh(),t=0;
    while(ch<'0'||ch>'9') t|=ch=='-',ch=gh();
    while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=gh();
    if(t) x=-x;
}
template<class T,class ...T1>
inline void read(T &x,T1 &...x1)
{
    read(x),read(x1...);
}
template<class T>
inline void write(T x)
{
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10+'0');
}
int N,M;
struct Complex
{
    double x,y;
    Complex operator+(const Complex& t) const
    {
        return {x+t.x,y+t.y};
    }
    Complex operator-(const Complex& t) const
    {
        return {x-t.x,y-t.y};
    }
    Complex operator*(const Complex& t) const
    {
        return {x*t.x-y*t.y,x*t.y+y*t.x};
    }
}a[MAXN],b[MAXN];
int Rev[MAXN],Bit,Tot;
inline void Fft(Complex a[],int inv)
{
    for(int i=0;i<Tot;++i)
        if(i<Rev[i]) std::swap(a[i],a[Rev[i]]);
    for(int mid=1;mid<Tot;mid<<=1)
    {
        auto w_1=Complex({std::cos(Pi/mid),inv*std::sin(Pi/mid)});
        for(int i=0;i<Tot;i+=mid*2)
        {
            auto w_k=Complex({1,0});
            for(int j=0;j<mid;++j,w_k=w_k*w_1)
            {
                auto x=a[i+j],y=w_k*a[i+j+mid];
                a[i+j]=x+y,a[i+j+mid]=x-y;
            }
        }
    }
}
int main()
{
    // freopen("fft.in","r",stdin);
    // freopen("fft.out","w",stdout);
    read(N,M);
    for(int i=0;i<=N;++i) scanf("%lf",&a[i].x);
    for(int i=0;i<=M;++i) scanf("%lf",&b[i].x);
    while((1<<Bit)<N+M+1) ++Bit;
    Tot=1<<Bit;
    for(int i=0;i<Tot;++i) Rev[i]=(Rev[i>>1]>>1)|((i&1)<<(Bit-1));
    Fft(a,1),Fft(b,1);
    for(int i=0;i<Tot;++i) a[i]=a[i]*b[i];
    Fft(a,-1);
    for(int i=0;i<=N+M;++i) printf("%d ",(int)(a[i].x/Tot+0.5));
    return 0;
}
/*
1 2
1 2
1 2 1
*/

FFT 优化高精度乘法

我们将一个高精度表示的数组转化为一个多项式：

$A=\overline{a_{n-1}a_{n-2}\dots a_0}=a_{n-1}\times 10^{n-1}+a_{n-2}\times 10^{n-2}+\dots a_0\times 10^0$ 。

那么，我们把 $10$ 看作是未知数 $x$ ，则 $a$ 就是需要求的系数。

这样做同样是有根据的，因为有 $C(x)=A(x)B(x)$ ，所以才会有 $C(10)=A(10)B(10)=AB$ 。

然后就可做了。

AC Code

using namespace std;
const int MAXN=3e6+1;
const double Pi=acos(-1);
template<class T>
inline void read(T &x)
{
    x=0;
    char ch=gh(),t=0;
    while(ch<'0'||ch>'9') t|=ch=='-',ch=gh();
    while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=gh();
    if(t) x=-x;
}
template<class T,class ...T1>
inline void read(T &x,T1 &...x1)
{
    read(x),read(x1...);
}
template<class T>
inline void write(T x)
{
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10+'0');
}
template<class T>
inline bool checkMax(T &x,T &y)
{
    return x<y?x=y,1:0;
}
template<class T>
inline bool checkMin(T &x,T &y)
{
    return x>y?x=y,1:0;
}
struct Cl
{
    double x,y;
    Cl operator+(const Cl &a) const
    {
        return {x+a.x,y+a.y};
    }
    Cl operator-(const Cl &a) const
    {
        return {x-a.x,y-a.y};
    }
    Cl operator*(const Cl &a) const
    {
        return {x*a.x-y*a.y,x*a.y+y*a.x};
    }
}a[MAXN],b[MAXN];
char str1[MAXN],str2[MAXN];
int res[MAXN],N,M;
int Rev[MAXN],Bit,Tot;
inline void Fft(Cl a[],int inv)
{
    for(int i=0;i<Tot;++i)
        if(i<Rev[i]) std::swap(a[i],a[Rev[i]]);
    for(int mid=1;mid<Tot;mid<<=1)
    {
        auto w1=Cl({cos(Pi/mid),inv*sin(Pi/mid)});
        for(int i=0;i<Tot;i+=mid*2)
        {
            auto wk=Cl({1,0});
            for(int j=0;j<mid;++j,wk=wk*w1)
            {
                auto x=a[i+j],y=wk*a[i+j+mid];
                a[i+j]=x+y,a[i+j+mid]=x-y;
            }
        }
    }
}
int main()
{
    // freopen(".in","r",stdin);
    // freopen(".out","w",stdout);
    scanf("%s%s",str1,str2);
    N=strlen(str1)-1,M=strlen(str2)-1;
    for(int i=0;i<=N;++i) a[i].x=str1[N-i]-'0';
    for(int i=0;i<=M;++i) b[i].x=str2[M-i]-'0';
    while((1<<Bit)<N+M+1) ++Bit;
    Tot=1<<Bit;
    for(int i=0;i<Tot;++i) Rev[i]=((Rev[i>>1]>>1))|((i&1)<<(Bit-1));
    Fft(a,1),Fft(b,1);
    for(int i=0;i<Tot;++i) a[i]=a[i]*b[i];
    Fft(a,-1);
    int k=0;
    for(int i=0,t=0;i<Tot||t;++i)
    {
        t+=a[i].x/Tot+0.5;
        res[k++]=t%10;
        t/=10;
    }
    while(k>1&&!res[k-1]) --k;
    for(int i=k-1;i>=0;--i) write(res[i]);
    return 0;
}
/*
114514 
1919810
*/