多项式汇总

一些模板与式子。

1.多项式乘法

梦开始的地方。

1.1. FFT 快速傅里叶变换

模板

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
const int N=4e6+10;
const double PI=acos(-1);
template<class T>
inline void read(T &x)
{
    x=0;int f=0; 
    char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=1;ch=getchar();}
    while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
    if(f)x=~x+1;
}
template<class T,class ...T1>
inline void read(T &x,T1 &...x1)
{
    read(x),read(x1...);
}
struct Complex
{
    double x,y;
    const Complex operator+(const Complex &t)const
    {
        return {x+t.x,y+t.y};
    }
    const Complex operator-(const Complex &t)const
    {
        return {x-t.x,y-t.y};
    }
    const Complex operator*(const Complex &t)const
    {
        return {x*t.x-y*t.y,x*t.y+y*t.x};
    }
}a[N],b[N];
int tot,bit,rev[N];
int n,m;
inline void FFT(Complex *a,int inv)
{
    for(int i=0;i<tot;i++)
        if(rev[i]<i)swap(a[rev[i]],a[i]);
    for(int mid=1;mid<tot;mid<<=1)
    {
        Complex w1={cos(PI/mid),sin(inv*PI/mid)};
        for(int i=0;i<tot;i+=mid*2)
        {
            Complex cur={1,0};
            for(int j=0;j<mid;j++,cur=cur*w1)
            {
                Complex x=a[i+j],y=cur*a[i+j+mid];
                a[i+j]=x+y,a[i+j+mid]=x-y;
            }
        }
    }
}
int main()
{
    read(n,m);
    for(int i=0;i<=n;i++)
    {
        int x;read(x);
        a[i].x=x;
    }
    for(int i=0;i<=m;i++)
    {
        int x;read(x);
        b[i].x=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. NTT 快速数论变换

模板

#include<iostream>
#include<cstdio>
#define ll long long
using namespace std;
const int N=4e6+10,mod=998244353,G=3,Gi=332748118;
template <class T>
inline void read(T &x)
{
    x=0;bool f=0;
    char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=1;ch=getchar();}
    while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
    if(f)x=-x;
}
template <class T,class ...T1>
inline void read(T &x,T1 &...x1)
{
    read(x),read(x1...);
}
int n,m;
int a[N],b[N];
int bit,tot,rev[N];
inline int adj(int x){return (x>=mod)?x-mod:x;}
inline int qpow(int x,int k)
{
	int res=1;
	while(k)
	{
		if(k&1)res=1ll*res*x%mod;
		x=1ll*x*x%mod;
		k>>=1;
	}
	return res;
}
inline void NTT(int *a,int inv)
{
	for(int i=0;i<tot;i++)
		if(rev[i]<i)swap(a[rev[i]],a[i]);
	for(int mid=1;mid<tot;mid<<=1)
	{
		int g=qpow(inv==1?G:Gi,(mod-1)/(mid<<1));
		for(int i=0;i<tot;i+=(mid<<1))
			for(int j=0,cur=1;j<mid;j++,cur=1ll*cur*g%mod)
			{
				int x=a[i+j],y=1ll*cur*a[i+j+mid]%mod;
				a[i+j]=adj(x+y),a[i+j+mid]=adj(x-y+mod);
			}
	}
	
}
int main()
{
	read(n,m);
	for(int i=0;i<=n;i++)
		read(a[i]);
	for(int i=0;i<=m;i++)
		read(b[i]);
	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);
	NTT(a,1),NTT(b,1);
	for(int i=0;i<tot;i++)
		a[i]=1ll*a[i]*b[i]%mod;
	NTT(a,-1);
	int inv=qpow(tot,mod-2);
	for(int i=0;i<n+m+1;i++)
		printf("%d ",1ll*a[i]*inv%mod);
	return 0;
}

2.多项式乘法逆

给定 $f(x)$，求 $f^{-1}(x)$。

倍增法

首先知道：

$$\begin{align} [x_0]f^{-1}(x)&=([x_0]f(x))^{-1}\notag \end{align}$$

考虑倍增。假设已经求出 $f(x)$ 在 $\pmod {x^{\lceil \frac{n}{2}\rceil}}$ 意义下逆元 $f_0^{-1}(x)$。推一下：

$$\begin{align} f(x)f_0^{-1}(x)&\equiv 1 \pmod {x^{\lceil \frac{n}{2}\rceil}}\notag \\ (f(x)f_0^{-1}(x)-1)^2&\equiv 0 \pmod {x^n}\notag \\ f(x)^2f_0^{-2}(x)-2f(x)f_0^{-1}(x)+1&\equiv 0 \pmod {x^n}\notag \\ \end{align}$$

同时除以 $f(x)$ 得到：

$$\begin{align} f(x)f_0^{-2}(x)-2f_0^{-1}(x)+f^{-1}(x)&\equiv 0 \pmod {x^n}\notag \\ \end{align}$$

整理一下，最终得到：

$$\begin{align} f^{-1}(x)&\equiv f_0^{-1}(x)(2-f(x)f_0^{-1}(x)) \pmod {x^n}\notag \\ \end{align}$$

递归计算即可，时间复杂度 $\mathcal O(n\log n)$。

模板

#include<iostream>
#include<cstdio>
using namespace std;
const int  N=4e6+10,mod=998244353,G=3,Gi=332748118;
template <class T>
inline void read(T &x)
{
    x=0;bool f=0;
    char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=1;ch=getchar();}
    while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
    if(f)x=-x;
}
template <class T,class ...T1>
inline void read(T &x,T1 &...x1)
{
    read(x),read(x1...);
}
int n;
int a[N],b[N],c[N];
int bit,tot,rev[N];
inline void init(int len)
{
	bit=0;
	while((1<<bit)<(len<<1))bit++;
	tot=1<<bit;
	for(int i=0;i<tot;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<bit-1);
}
inline int adj(int x){return (x>=mod)?x-mod:x;}
inline int qpow(int x,int k)
{
	int res=1;
	while(k)
	{
		if(k&1)res=1ll*res*x%mod;
		x=1ll*x*x%mod;
		k>>=1;
	}
	return res;
}
inline void NTT(int *a,int inv)
{
	for(int i=0;i<tot;i++)
		if(rev[i]<i)swap(a[rev[i]],a[i]);
	for(int mid=1;mid<tot;mid<<=1)
	{
		int g=qpow(inv==1?G:Gi,(mod-1)/(mid<<1));
		for(int i=0;i<tot;i+=(mid<<1))
			for(int j=0,cur=1;j<mid;j++,cur=1ll*cur*g%mod)
			{
				int x=a[i+j],y=1ll*cur*a[i+j+mid]%mod;
				a[i+j]=adj(x+y),a[i+j+mid]=adj(x-y+mod);
			}
	}
}
void poly_inv(int *a,int *b,int len)
{
	if(len==1)
	{
		b[0]=qpow(a[0],mod-2);
		return ;
	}
	poly_inv(a,b,(len+1)>>1);
	init(len);
	for(int i=(len+1)>>1;i<tot;i++)b[i]=0;
	for(int i=0;i<tot;i++)
		c[i]=(i<len?a[i]:0);
	NTT(b,1),NTT(c,1);
	for(int i=0;i<tot;i++)
		b[i]=1ll*adj(2-1ll*b[i]*c[i]%mod+mod)*b[i]%mod;
	NTT(b,-1);
	int inv=qpow(tot,mod-2);
	for(int i=0;i<tot;i++)
		b[i]=adj(1ll*b[i]*inv%mod+mod);
	for(int i=len;i<tot;i++)b[i]=0;
}
int main()
{
	read(n);
	for(int i=0;i<n;i++)
		read(a[i]);
	poly_inv(a,b,n);
	for(int i=0;i<n;i++)
		printf("%d ",b[i]);
	return 0;
}

3.多项式开根

给定 $g(x)$，求 $f(x)$ 满足：

$$\begin{align} f(x)^2\equiv g(x)\pmod x^n \end{align}$$

倍增法

首先知道：

$$\begin{align} [x_0]f(x)&=\sqrt{[x_0]g(x)}\notag \end{align}$$

考虑倍增。假设已经求出 $g(x)$ 在 $\pmod {x^{\lceil \frac{n}{2}\rceil}}$ 意义的平方根 $f_0(x)$。推一下：

$$\begin{align} f_0(x)^2&\equiv g(x) &\pmod {x^{\lceil \frac{n}{2}\rceil}} \notag\\ f_0(x)^2-g(x)&\equiv 0 &\pmod {x^{\lceil \frac{n}{2}\rceil}} \notag\\ (f_0(x)^2-g(x))^2&\equiv 0 &\pmod {x^n} \notag\\ (f_0(x)^2+g(x))^2&\equiv 4 f_0(x)^2g(x) &\pmod {x^n} \notag\\ \frac{(f_0(x)^2+g(x))^2}{4 f_0(x)^2} &\equiv f(x)^2 &\pmod {x^n} \notag\\ \frac{f_0(x)^2+g(x)}{2 f_0(x)} &\equiv f(x) &\pmod {x^n} \notag\\ \end{align}$$

最终得到：

$$\begin{align} 2^{-1}f_0(x)+2^{-1}f_0^{-1}(x)g(x) &\equiv f(x) &\pmod {x^n} \notag\\ \end{align}$$

需要一次多项式求逆，时间复杂度 $\mathcal O(n\log n)$。

多项式快速幂

4.多项式除法

给定一个 $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)$

做法

我们发现假如没有 $R(x)$，那么我们只需要对 $G(x)$ 做一次多项式求逆，再与 $F(x)$ 相乘即可。于是尝试去掉 $R(x)$ 的影响。

考虑一种变换，我们将 $F(x)$ 视作 $n$ 次多项式（实际次数小于等于 $n$），$F_R(x)=x^n F(\frac{1}{x})$，即 $[x_i]F_R(x)=[x_{n-i}]F(x)$。

开始推导：

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

注意 $Q_R(x)$ 的次数为 $n-m$，我们考虑对整个式子对 $x^{n-m+1}$ 取模并不会影响 $Q_R(x)$，并且可以消去 $R_R(x)$。

$$\begin{align} F_R(x) &\equiv Q_R(x) * G_R(x) &\pmod {x^{n-m+1}}\notag\\ Q_R(x) &\equiv F_R(x) * G_R^{-1}(x) &\pmod {x^{n-m+1}}\notag\\ \end{align}$$

于是可以得到 $Q(x)$，继而可以得到 $R(x)$。