多项式算法合集

头函数

#define poly vector<int>
#define bg begin
#define pb push_back
const int mod=998244353,g=3;
inline int add(int a,int b){return a+b>=mod?a+b-mod:a+b;}
inline void Add(int &a,int b){a=add(a,b);}
inline int dec(int a,int b){return a>=b?a-b:a-b+mod;}
inline void Dec(int &a,int b){a=dec(a,b);}
inline int mul(int a,int b){return 1ll*a*b>=mod?1ll*a*b%mod:a*b;}
inline void Mul(int &a,int b){a=mul(a,b);}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));return res;}

多项式加减点乘点除

幼儿园小朋友应该都会了吧

inline poly operator +(poly a,poly b){
	poly c;int lim=max(a.size(),b.size());c.resize(lim);
	for(int i=0;i<lim;i++)c[i]=add(a[i],b[i]);return c;
}
inline poly operator -(poly a,poly b){
	poly c;int lim=max(a.size(),b.size());c.resize(lim);
	for(int i=0;i<lim;i++)c[i]=dec(a[i],b[i]);return c;
}
inline poly operator *(poly a,int b){
	for(int i=0;i<a.size();i++)Mul(a[i],b);return a;
}
inline poly operator /(poly a,int b){
	for(int i=0,inv=ksm(b,mod-2);i<a.size();i++)Mul(a[i],inv);
	return a;
}

---

多项式乘法

FFT：

前置

多项式的点值和系数表示法：

对于一个$n$次多项式$f(x)$
$f(x)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3\dots \dots a_n{x}^n$被称作该多项式的系数表示
我们可以通过带任意一个$x$都可以的到唯一的一个$f(x)$
但$oi$中一般$x$一般都只是一个不定元，不会带入特定值计算，比如用作表示下标之类的

而如果我们把不同的$n+1$带入进去得到的$n+1$个点值叫做点值表示法
由$n+1$个点也可以还原出唯一一个$n$次多项式

---

虚数

即$i，\sqrt{-1}$
考虑在平面直角坐标系内
将$y$轴用$i$来表示

复数更准确的定义是

$e^{ik}=\cos (k)+i\sin (k)$

这样平面上一个点就是$(a,bi)$的形式
而2个虚数相乘，就对应平面直角坐标系上2个向量模长相乘，极角相加

由于$C++$自带的复数很慢
所以我们一般手写复数结构体

const double pi=acos(-1);
struct complex{
    double x,y;
    complex(double _x=0,double _y=0):x(_x),y(_y){}
    friend inline complex operator +(const complex &a,const complex &b){
        return complex(a.x+b.x,a.y+b.y);
    }
    friend inline complex operator -(const complex &a,const complex &b){
        return complex(a.x-b.x,a.y-b.y);
    }
    friend inline complex operator /(const complex &a,const double &b){
        return complex(a.x/b,a.y/b);
    }
    friend inline complex operator *(const complex &a,const complex &b){
        return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
    }
};

---

单位根

$fft$的根本
$n$次单位根指的是满足$w^n=1$的复数
$n$次单位根有$n$个,分别表示为$w_n^k，k=0,1,2,\dots \dots n-1$
实际上对应的将平面单位圆周上的$n$个点
这些点将单位圆周均分成$n$块，且构成一个正$n$边形

更准确的表示为$w_n^k=e^{2\pi ik/n}$
即$w_n^k=\cos (2\pi k/n)+i\sin (2\pi k/n)$表示
即单位圆上的$n$点

单位根几个重要的性质：

1、消去引理

对于任何整数$n\ge 0,k\ge 0,d\ge 0$
有$w_{dn}^{dk}=w_n^k$

证明:$w_{dn}^{dk}=e^{2\pi idk/dn}=e^{2\pi ik/n}=w_n^k$

2、折半引理

如果$n\ge 0$且$n$为偶数，那么
对所有$n$次单位根平方，得到的集合是$n/2$个$n/2$次单位根的集合
说白了就是$w_n^{k+\frac{n}{2}}=-w_n^k$或者就是$w_n^{\frac{n}{2}}=-1$
证明：画个单位圆，$\frac{n}{2}$就是旋转了$180°$,自然取反

3、求和引理

对于任意整数$n\ge 1$和$n̸|k$
有${\sum }_{j=0}^{n-1}(w_n^k{)}^j=0$

考虑这其实是一个等比数列，$w_n^k̸=1$时
原式$=\frac{(w_n^k{)}^n-1}{w_n^k-1}=\frac{(w_n^n{)}^k-1}{w_n^k-1}=0$

而当$n|k$时，原式显然为$n$

这个性质在后面很重要

算法

考虑对于两个多项式
$A(x)=a_0+a_{1}x+a_2{x}^2+a_3{x}^3\dots \dots a_n{x}^n$
$B(x)=b_0+b_{1}x+b_2{x}^2+b_3{x}^3\dots \dots b_m{x}^m$
我们要求一个$n+m$次多项式$C(x)=A(x)\ast B(x)$
更具体的$C(x)$满足$c_i={\sum }_{j=0}^i{a}_j\ast b_{i-j}$
也就是2个数列倒着乘的和，所谓的卷积

考虑如果我们直接拆开暴力做是$O(n^2)$的
当然也有一种分治乘法可以做到$n^{lo{g}_{2}3}$(大概就是大常数$n\sqrt{n}$)

而$FFT$可以做到$O(nlogn)$求出$C$

---

下面为了方便假设$n=m$，$n̸=m$也没有区别

考虑如果我们有$n+1$个值$x_1,\dots \dots ，x_{n+1}$
如果已经求出$A(x_1)，\dots \dots A(x_{n+1})$
和$B(x_1),\dots \dots B(x_{n+1})$
也就是分别求出$x_1,\dots \dots ，x_n+1$的点值

那么我们可以直接在$O(n)$的时间内求出
$C(x_1)=A(x_1)\ast B(x_1)，\dots \dots ，C(x_{n+1})=A(x_{n+1})\ast B(x_{n+1})$

现在我们考虑这样

先对$A(x)$,$B(x)$求出$n+1$个点的点值，乘起来得到$C(x)$的点值
又对于一个$n$次多项式，如果我们知道其$n+1$个不同点值
就一定可以还原出一个唯一的多项式
于是最后再由点值表示还原出原来的$C(x)$

注意由于实际乘出来$C(x)$最高系数是$2n-1$
所以我们需要带入$2n$个点求值
于是我们会将$A(x),B(x)$补0到$2n$次项
实际上由于$fft$的特殊性
我们会将项数补充到2的整数次幂次

当然非二的整数次幂项的多项式乘法也是可以做的（见下面补充的混合基$fft$和$Bluestein$）

第一步将系数转成点值是正变换，称作$DFT$,单次复杂度$O(nlogn)$
由点值还原系数为逆变换，称作$IDFT$,单次复杂度$O(nlogn)$
于是总复杂度就是$O(nlogn)$
整个操作被称为快速傅里叶变换$(FFT)$

---

DFT：

考虑我们带入$n$次单位根$w_n$
$A(w_n^k)={\sum }_{i=0}^n{a}_i(w_n^k{)}^i$
考虑将下标按照奇偶分类
$A(w_n^k)={\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i}(w_n^k{)}^{2i}+{\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i+1}(w_n^k{)}^{2i+1}$

$={\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i}(w_n^{2ki})+w_n^k{\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i+1}(w_n^{2ki})$

$={\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i}(w_{\frac{n}{2}}^k)+w_n^k{\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i+1}(w_{\frac{n}{2}}^k)$

$=A_o(w_{\frac{n}{2}}^k)+w_n^k{A}_1(w_{\frac{n}{2}}^k)$

这样$A_0$和$A_1$都只有$\frac{n}{2}$项了
可以继续递归去做
尽管现在复杂度并没有变化

考虑对于$A(w_n^{k+\frac{n}{2}})$
$A(w_n^{k+\frac{n}{2}})={\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i}(w_n^{k+\frac{n}{2}}{)}^{2i}+w_n^{k+\frac{n}{2}}{\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i+1}(w_n^{k+\frac{n}{2}}{)}^{2i}$

考虑单位根的消去引理

$w_n^{k+\frac{n}{2}}=-w_n^k$

则$A(w_n^{k+\frac{n}{2}})={\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i}(-w_n^k{)}^{2i}+(-w_n^k){\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i+1}(-w_n^k{)}^{2i}$

$={\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i}(w_n^k{)}^{2i}-w_n^k{\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i+1}(w_n^k{)}^{2i}$

$={\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i}(w_{\frac{n}{2}}^k)-w_n^k{\sum }_{i=0}^{\frac{n}{2}-1}{a}_{2i+1}(w_{\frac{n}{2}}^k)$

$=A_o(w_{\frac{n}{2}}^k)-w_n^k{A}_1(w_{\frac{n}{2}}^k)$

我们发现唯一不同的就是第二项的符号
也就是说如果我们知道$A_o(w_{\frac{n}{2}}^k)$和$A_1(w_{\frac{n}{2}}^k)$
我们就可以同时知道$A(w_n^{k+\frac{n}{2}})$和$A(w_n^k)$

考虑对于$k\in [0,n-1],$我们求$A(w_n^k)$
就只需要知道$k\in [0,\frac{n}{2}-1],A_0(w_{\frac{n}{2}}^k)$和$A_1(w_{\frac{n}{2}}^k)$
就可以在$O(n)$的时间内得到$A$

而$A_0$,$A_1$系数都只有$\frac{n}{2}$个，所以规模只有原来的一般
递归求解即可

时间复杂度$T(n)=2\ast T(\frac{n}{2})+O(n)=O(nlogn)$

这里也是之所以要把项数补充到$2^k$
因为每一次都把$n$项分成$\frac{n}{2}$项
如果$n$是奇数，那就没法分开了

由于递归常数比较大，而一般$fft$有关的题时间瓶颈就在$DFT$上面
$DFT$不知道为什么 常数也很大
写的差的$fft$甚至可以跑$1e6$的数据$T$掉

所以我们考虑迭代做
由于每个数最终在的位置和原来不一样
所以我们要预处理出最终的位置上

据说找规律得到了$O(n)$预处理的方法
代码如下：
没看懂，选择全文背诵
当然也可以$nlogn$模拟最终位置

int rev[N<<2];
inline void init(int lim){
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)*(lim>>1));
}

我们先把每个数放到最终应该在的地方然后一步步迭代回去就是了

inline void DFT(complex f[],int lim){
    for(int i=0;i<lim;i++)if(rev[i]>i)swap(f[i],f[rev[i]]);
    for(int mid=1;mid<lim;mid<<=1){
        complex now=plx(cos(pi/mid),sin(pi/mid));
        for(int i=0;i<lim;i+=(mid<<1)){
            complex w=plx(1,0);
            for(int j=0;j<mid;j++,w=w*now){
                plx p0=f[i+j],p1=w*f[i+j+mid];
                f[i+j]=p0+p1,f[i+j+mid]=p0-p1;
            }
        }
    }
}

---

IDFT：

以下摘抄自$miskcoo$神犇（$markdown$太难写了QAQ）

---

代码实现

const double pi=acos(-1);
struct plx{
    double x,y;
    plx(double _x=0,double _y=0):x(_x),y(_y){}
    friend inline plx operator +(const plx &a,const plx &b){
        return plx(a.x+b.x,a.y+b.y);
    }
    friend inline plx operator -(const plx &a,const plx &b){
        return plx(a.x-b.x,a.y-b.y);
    }
    friend inline plx operator /(const plx &a,const double &b){
        return plx(a.x/b,a.y/b);
    }
    friend inline plx operator *(const plx &a,const plx &b){
        return plx(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
    }
};
inline void fft(plx f[],int lim,int kd){//kd表示在做正变换还是逆变换
    for(int i=0;i<lim;i++)if(rev[i]>i)swap(f[i],f[rev[i]]);
    for(int mid=1;mid<lim;mid<<=1){
        plx now=plx(cos(pi/mid),kd*sin(pi/mid));
        for(int i=0;i<lim;i+=(mid<<1)){
            plx w=plx(1,0);
            for(int j=0;j<mid;j++,w=w*now){
                plx p0=f[i+j],p1=w*f[i+j+mid];
                f[i+j]=p0+p1,f[i+j+mid]=p0-p1;
            }
        }
    }
    if(kd==-1)for(int i=0;i<lim;i++)f[i]=f[i]/lim;
}

---

NTT：

由于$fft$涉及复数和实数运算，实际会出现精度误差
在整数运算时难免会出锅

所以我们考虑一个能在模意义下的变换
这就是$ntt$

首先引入原根的概念

对于一个素数$p$，其原根$g$定义为满足$g^0,g^1,\dots \dots ，g^{p-2}$互不相同的数
又由于费马定理，对一个素数$p$,有$a^{p-1}\equiv 1(modp)$
这个和单位根很相似

我们考虑单位根之所以能够做$fft$，是因为$w$的三个性质
考虑如果我们对于一个素数$p=2^n\ast k+1$,我们令$g_n=g^k$
这样我们就可以满足$g_n^0,g_n^1,g_n^2\dots \dots g_n^{n-1}$互不相同且满足这些性质（不想证了）

但这也就限制了$ntt$的模数必须是$k\ast 2^n+1$的形式
否则就要做任意模数$ntt$
做法就是选出几个模数分别做之后$Crt$合并答案（并不会）
代码和$fft$类似
由于$w_n^{-k}=w_n^{n-k}$
所以我们可以先做的时候不管正逆变换
最后把$1$~$n-1$反转一下就可以了

代码实现

inline void ntt(poly &f,int lim,int kd){
    for(int i=0;i<lim;i++)if(i<rev[i])swap(f[i],f[rev[i]]);
    for(int mid=1;mid<lim;mid<<=1){
        int now=ksm(g,(mod-1)/(mid<<1));
        for(int i=0;i<lim;i+=(mid<<1)){
            int w=1;
            for(int j=0;j<mid;j++,w=mul(w,now)){
                int a0=f[i+j],a1=mul(w,f[i+j+mid]);
                f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
            }
        }
    }
    if(kd==-1&&(reverse(f.bg()+1,f.bg()+lim),1))for(int i=0,inv=ksm(lim,mod-2);i<lim;i++)f[i]=mul(f[i],inv);
}

由于$ntt$中每次乘起来很耗费常数
我也不知道为什么，但就是很耗时间

于是我们可以预处理出原根优化常数

预处理原根

const int N=1000005,C=21;
int *w[22];
inline void init_w(){
	for(int i=1;i<=C;i++)
	w[i]=new int[1<<(i-1)];
	int wn=ksm(g,(mod-1)/(1<<C));
	w[C][0]=1;
	for(int i=1;i<(1<<(C-1));i++)w[C][i]=mul(w[C][i-1],wn);
	for(int i=C-1;i;i--)
	for(int j=0;j<(1<<(i-1));j++)
	w[i][j]=w[i+1][j<<1];
}

速度比不预处理快了差不多$25\%$到$45\%$不等

其实$ntt$本身常数不算很大，运算常数大概也只有5、6左右
不过下标不连续可能会导致慢一些

inline void ntt(poly &f,int lim,int kd){
	for(int i=0;i<lim;i++)if(i>rev[i])swap(f[i],f[rev[i]]);
	for(int mid=1,l=1;mid<lim;mid<<=1,l++)
	for(int i=0;i<lim;i+=(mid<<1))
		for(int j=0,a0,a1;j<mid;j++)
			a0=f[i+j],a1=mul(f[i+j+mid],w[l][j]),
			f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
	if(kd==-1&&(reverse(f.begin()+1,f.begin()+lim),1))
		for(int inv=ksm(lim,mod-2),i=0;i<lim;i++)Mul(f[i],inv);
}

乘法

可以在比较小的时候暴力加循环展开 优化常数

inline poly operator *(poly a,poly b){
	int deg=a.size()+b.size()-1,lim=1;
	if(deg<=128){
		poly c(deg,0);
		for(int i=0;i<a.size();i++)
		for(int j=0;j<b.size();j++)
			Add(c[i+j],mul(a[i],b[j]));
		return c;
	}
	while(lim<deg)lim<<=1;init(lim);
	a.resize(lim),ntt(a,lim,1);
	b.resize(lim),ntt(b,lim,1);
	for(int i=0;i<lim;i++)Mul(a[i],b[i]);
	ntt(a,lim,-1),a.resize(deg);
	return a;
}

---

多项式求逆：

已知一个$n-1$次多项式$A(x)$，求多项式$B(x)$满足：

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

求解过程

倍增：

若已知$$A(x)B(x{)}^{\prime }\equiv 1(mod{x}^{\lceil \frac{n}{2}\rceil })$$

$$A(x)B(x{)}^{\prime }-A(x)B(x)\equiv 0(mod{x}^{\lceil \frac{n}{2}\rceil })$$

$$B(x{)}^{\prime }-B(x)\equiv 0(mod{x}^{\lceil \frac{n}{2}\rceil })$$

平方：

$$B(x{)}^{\prime 2}+B(x{)}^2-2B(x{)}^{\prime }B(x)\equiv 0(mod{x}^n)$$
$$A(x)(B(x{)}^{\prime 2}+B(x{)}^2-2B(x{)}^{\prime }B(x))\equiv 0(mod{x}^n)$$
$$B(x)\equiv 2B(x{)}^{\prime }-A(x)B(x{)}^{\prime 2}(mod{x}^n)$$

复杂度$T(n)=T(\frac{n}{2})+O(nlogn)=O(nlogn)$

注意次数，最高到3倍，开的4倍

代码实现

inline poly Inv(poly a,int deg){
	poly c,b(1,ksm(a[0],mod-2));
	for(int lim=4;lim<(deg<<2);lim<<=1){
		init(lim);
		c=a,c.resize(lim>>1);
		c.resize(lim),ntt(c,lim,1);
		b.resize(lim),ntt(b,lim,1);
		for(int i=0;i<lim;i++)Mul(b[i],dec(2,mul(b[i],c[i])));
		ntt(b,lim,-1),b.resize(lim>>1);
	}b.resize(deg);return b;
}

---

多项式开方：

已知一个$n-1$次多项式$A(x)$,求一个$mod{x}^n$意义下的多项式$B(x)$满足

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

满足$A[0]=1$

求解过程

倍增

首先当$n=1$时要满足$A[0]=1$（否则要二次剩余解，老子不会）

假设已知$$B(x{)}^{\prime 2}\equiv A(x)(mod{x}^{\lceil \frac{n}{2}\rceil })$$
由于$\lceil \frac{n}{2}\rceil$以上的项是不会有影响的
所以$$B(x{)}^{\prime }\equiv B(x)(mod{x}^{\lceil \frac{n}{2}\rceil })$$
移项平方得：
$$B(x{)}^2+B(x{)}^{\prime 2}-2B(x)B(x{)}^{\prime }\equiv 0(mod{x}^n)$$

$$A(x)+B(x{)}^{\prime 2}\equiv 2B(x)B(x{)}^{\prime }(mod{x}^n)$$
$$B(x)\equiv \frac{A(x)}{2B(x{)}^{\prime }}+\frac{B(x)}{2}(mod{x}^n)$$

复杂度$T(n)=T(\frac{n}{2})+O(nlogn)=O(nlogn)$

代码实现

inline poly Sqrt(poly a,int deg){
	poly b(1,1),c,d;
	for(int lim=4;lim<(deg<<2);lim<<=1){
		c=a,c.resize(lim>>1);
		init(lim),d=Inv(b,lim>>1),
		c.resize(lim),ntt(c,lim,1);
		d.resize(lim),ntt(d,lim,1);
		for(int i=0;i<lim;i++)Mul(c[i],d[i]);
		ntt(c,lim,-1),b.resize(lim>>1);
		for(int i=0;i<(lim>>1);i++)b[i]=mul(inv2,add(b[i],c[i]));
	}b.resize(deg);return b;
}

---

多项式除法和取模：

给定一个$n$次多项式$A(x)$和一个$m$次多项式$B(x)$
求一个$n-m$次多项式$C(x)$和$m-1$次多项式$D(x)$满足：
$$A(x)=B(x)C(x)+D(x)$$

求解过程

考虑对一个最高次数为$n$的多项式$B(x)$操作$B^R(x)=x^{n}B(\frac{1}{x})$
会发现$B^R(x)$只是$B(x)$的系数反转的柿子

考虑$$A(x)=B(x)C(x)+D(x)$$
$A(x)$最高项为$n$,$B(x)$最高项为$m$,则$C(x)$最高项为$n-m$,$D(x)$为$m-1$
两边同时乘一个$x^n$,并带入$\frac{1}{x}$
$$x^{n}A(x)=x^{m}B(x)x^{n-m}C(x)+x^{n-m+1}\ast x^{m-1}D(x)$$
$$A^R(x)=B^R(x)C^R(x)+x^{n-m+1}{D}^R(x)$$
考虑在$mod{x}^{n-m+1}$意义下，$A,B$已知，$C$最高为$n-m$不影响，而$D$被消去
则$$A^R(x)\equiv B^R(x)C^R(x)(mod{x}^{n-m+1})$$

多项式求逆就可以了

复杂度$O(nlogn)$

代码实现

inline poly operator /(poly a,poly b){
	int lim=1,deg=a.size()-b.size()+1;
	reverse(a.bg(),a.end());
	reverse(b.bg(),b.end());
	while(lim<=deg)lim<<=1;
	b=Inv(b,lim),b.resize(deg);
	a=a*b,a.resize(deg);
	reverse(a.bg(),a.end());
	return a;
}
inline poly operator %(poly a,poly b){
	poly c=a-(a/b)*b;
	c.resize(b.size()-1);
	return c;
}

---

多项式求导与积分

我怕是自学的是一个假的微积分
假装自己会的差不多了

复杂度$O(n)$

代码实现

inline poly deriv(poly a){
	for(int i=0;i<a.size()-1;i++)a[i]=mul(a[i+1],i+1);
	a.pob();return a;
}
inline poly integ(poly a){
	a.pb(0);
	for(int i=a.size()-1;i;i--)a[i]=mul(a[i-1],inv[i]);
	a[0]=0;
	return a;
}

---

多项式Ln

已知一个$n-1$次多项式$g(x)$,求一个$mod{x}^n$意义下的多项式$f(x)$，满足：
$$f(x)\equiv Ln(g(x))(mod{x}^n)$$

求解过程

由于有公式
若$$f(x)=ln(g(x))$$
则$$g^{\prime }(x)=f^{\prime }(x)g(x)$$
$$f^{\prime }(x)=\frac{g^{\prime }(x)}{g(x)}$$
且要满足$g[0]=1$否则老子不会
求导求逆最后再积分就可以了

复杂度$O(nlogn)$

代码实现

/*
if f(x)=Ln(g(x))
then g'(x)=f'(x)g(x)
*/
inline poly Ln(poly a,int lim){
	a=integ(deriv(a)*Inv(a,lim)),a.resize(lim);
	return a;
}

---

多项式Exp

前置知识

泰勒展开
牛顿迭代

泰勒展开

$$f(x)=f(x_0)+\frac{f^1(x_0)}{1!}(x-x_0)+\frac{f^2(x_0)}{2!}(x-x_0{)}^2+\dots \dots +\frac{f^n(x_0)}{n!}(x-x_0{)}^n+\dots \dots$$

考虑我们要构造一个函数$g(x)$使$g(x)$完全拟合$f(x)$
那首先初始点$x_0$的值要和$f(x_0)$一样
在此基础上只需要保证一阶导数，二阶导数……都完全相同即可
即$\forall i,f^i(x)=g^i(x)$
由于求$g$第$i$阶导数时为$i!a_i=f^i(x)$
即$a_i=\frac{f^i(x)}{i!}$
所以得证

牛顿迭代

在多项式中一般用来解这类问题：

假设有函数$f$和一个多项式$g(x)mod{x}^n$
满足$f(g(x))\equiv 0(mod{x}^n)$
已知$f$，求$g$

说白了就是用来解$f(g(x))\equiv 0(mod{x}^n)$之类的方程

首先在$n=1$的时候即常数时单独求

假设已经知道$f(g^{\prime }(x))\equiv 0(mod{x}^{\lceil \frac{n}{2}\rceil })$
要求$f(g(x))\equiv 0(mod{x}^n)$

考虑将$f(g(x))$在$g^{\prime }(x)$处泰勒展开
$$f(g(x))=f(g^{\prime }(x))+\frac{f^1(g^{\prime }(x))}{1!}(g(x)-g^{\prime }(x))+\frac{f^2(g^{\prime }(x))}{2!}(g(x)-g(x{)}^{\prime }{)}^2+\dots \dots$$

首先显然有$g(x)\equiv g^{\prime }(x)(mod{x}^{\lceil \frac{n}{2}\rceil })$
所以$g(x)-g^{\prime }(x)$最低项次数一定大于$\lceil \frac{n}{2}\rceil$

则在$mod{x}^n$意义下，整个式子从$(g(x)-g(x{)}^{\prime }{)}^2$开始都为$0$了
则$$f(g(x))\equiv f(g^{\prime }(x))+f^1(g^{\prime }(x))(g(x)-g^{\prime }(x))(mod{x}^n)$$

又$f(g(x))\equiv 0(mod{x}^n)$

$$g(x)\equiv g^{\prime }(x)-\frac{f(g^{\prime }(x))}{f^1(g^{\prime }(x))}(mod{x}^n)$$

这就大功告成了

例：

比如多项式开根
就是解$B(x{)}^2-A(x)\equiv 0(mod{x}^n)$
假设已知$B^{\prime }(x{)}^2-A(x)\equiv 0(mod{x}^{\lceil \frac{n}{2}\rceil })$
这时候$f(B(x))=B(x{)}^2-A(x),$则$f^1(B(x))=2B(x)$
带入$B(x)\equiv B^{\prime }(x)-\frac{B^{\prime }(x{)}^2-A(x)}{2B^{\prime }(x)}\equiv \frac{B^{\prime }(x{)}^2+A(x)}{2B^{\prime }(x)}$
就是我们推出来的式子

---

已知一个$n-1$次多项式$g(x)$,求一个$mod{x}^n$意义下的多项式$f(x)$，满足：
$$f(x)\equiv e^{g(x)}(mod{x}^n)$$

也就是$$Ln(f(x))\equiv g(x)(mod{x}^n)$$

求解过程

倍增：

考虑原问题，即求解方程$$Ln(f(x))-g(x)\equiv 0(mod{x}^n)$$

同样假设已经知道$Ln(f^{\prime }(x))-g(x)\equiv 0(mod{x}^{\lceil \frac{n}{2}\rceil })$

令$P(f(x))=Ln(f(x))-g(x)$,则$P^1(f(x))=\frac{1}{f(x)}$

则$$f(x)\equiv f^{\prime }(x)-\frac{Ln(f^{\prime }(x))-g(x)}{\frac{1}{f^{\prime }(x)}}\equiv f^{\prime }(x)-f^{\prime }(x)Ln(f^{\prime }(x))+f^{\prime }(x)g(x)$$

$$\equiv f^{\prime }(x)(1-Ln(f^{\prime }(x))+g(x))(mod{x}^n)$$

复杂度$T(n)=T(\frac{n}{2})+O(nlogn)=O(nlogn)$

代码实现

inline poly exp(poly a,int deg){
    poly b(1,1),c;a.resize(deg<<1);
    for(int lim=2;lim<(deg<<1);lim<<=1){
        c=Ln(b,lim);
        for(int i=0;i<lim;i++)c[i]=dec(a[i],c[i]);
        Add(c[0],1),b=b*c;
        b.resize(lim);
    }b.resize(deg);return b;
}

---

多项式多点求值

已知一个$n$次多项式$f(x)$，求$f(a_1)\dots \dots f(a_m)$

求解过程

考虑构造函数$P(x)={\prod }_{i=1}^m(x-a_i)$
显然$\forall i\in [1,m],P(a_i)=0$

假设$f(x)=P(x)g(x)+A(x)$
$A(x)=f(x)\%P(x)$

那显然$f(a_i)=A(a_i)$
但由于$P(x)$是$m$次的，没有起到优化的作用

而考虑对于$k，P(x)=(x-a_k)$
则$f(x)\%P(x)$后就只剩下一个常数项，即$f(a_i)$的值了
但是这样一次就$nlogn$了

考虑分治优化
$P_1(x)={\prod }_{i=l}^{mid}(x-a_i),P_2(x)={\prod }_{i=mid+1}^r(x-a_i)$
$\forall i\in [l,mid]，f(x)\%P_1(x)=f(a_i)$
取模之后$f(x)$的次数减少了一半，继续递归求解即可

$P(x)$可以先分治$ntt$预处理出来
复杂度$T(n)=2\ast T(\frac{n}{2})+O(nlogn)=O(nlo{g}^{2}n)$

代码实现

第一次写$T$掉了，预处理了一波单位根就跑过去了
也可以在比较小的时候暴力秦九韶展开，然并没写

poly f[N<<2];
int a[N];
int n,m;
#define lc (u<<1)
#define rc ((u<<1)|1)
#define mid ((l+r)>>1)
void build(int u,int l,int r,int *v){
	if(l==r){f[u].pb(dec(0,v[l])),f[u].pb(1);return;}
	build(lc,l,mid,v),build(rc,mid+1,r,v);
	f[u]=f[lc]*f[rc];
}
void calc(int u,int l,int r,poly now,int *v){
	if(l==r){v[l]=now[0];return;}
	calc(lc,l,mid,now%f[lc],v),calc(rc,mid+1,r,now%f[rc],v);
}
#undef lc
#undef rc
#undef mid

---

多项式快速插值

考虑传统的拉格朗日插值法插多项式是$O(n^2)$的
即构造函数$$f(x)=\sum _{i=1}^n{y}_i\prod _{j≠i}\frac{(x-x_j)}{x_i-x_j}$$

化一下

$$f(x)=\sum _{i=1}^n\frac{y_i}{{\prod }_{j≠i}(x_i-x_j)}\prod _{j≠i}(x-x_j)$$

考虑如何求出$G={\prod }_{j≠i}(x_i-x_j)$
设$g(x)={\prod }_j(x-x_j)$
$j≠i$就相当于除以了$x-x_i$

那就变成了$\frac{g(x_i)}{(x-x_i)}$
但是这个分子分母就都是0了，没法直接求

根据洛必达法则：
若
$$\underset{x\to a}{\lim }f(x)=0,\underset{x\to a}{\lim }g(x)=0$$

有
$$\underset{x\to a}{\lim }\frac{f(x)}{g(x)}=\underset{x\to a}{\lim }\frac{f^{\prime }(x)}{g^{\prime }(x)}$$

同时取导得到$G=g^{\prime }(x_i)$
接下来考虑对整个式子分治
设$f_{l,r}$表示分治$[l,r]$得到的答案

$$f_{l,r}=\sum _{i=l}^r\frac{y_i}{g^{\prime }(x_i)}\prod _{j≠i}(x-x_j)$$

$$=\prod _{k=mid+1}^r(x-x_k)\sum _{i=l}^{mid}\frac{y_i}{g^{\prime }(x_i)}\prod _{j≠i}^{[l,mid]}(x-x_j)+\prod _{k=l}^{mid}(x-x_k)\sum _{i=mid+1}^r\frac{y_i}{g^{\prime }(x_i)}\prod _{j≠i}^{[mid+1,r]}(x-x_j)$$

$$=\prod _{i=mid+1}^r(x-x_i)f_{l,mid}+\prod _{i=l}^{mid}(x-x_i)f_{mid+1，r}$$

先分治$ntt$求出$g$，多点求值把$g^{\prime }(x_i)$求出来再分治一波就完了

复杂度$O(nlo{g}^{2}n)$

代码

---

下降幂多项式乘法

首先考虑对于$x^{\underline{n}}$构建$EGF$
$x^{\underline{n}}={\sum }_{i=n}^{\infty }\frac{i^{\underline{n}}}{i!}{x}^i={\sum }_{i=n}^{\infty }\frac{1}{(i-n)!}{x}^i=x^n{e}^x$

考虑对于$f(x)={\sum }_{i=0}^{\infty }{a}_i{x}^{\underline{i}}$的点值构造$EGF$

$g(x)={\sum }_{i=0}^{\infty }\frac{f(i)}{i!}{x}^i={\sum }_{i=0}^{\infty }\frac{x^i}{i!}{\sum }_{j=0}^{\infty }{a}_j{i}^{\underline{j}}$

$={\sum }_{j=0}^{\infty }{a}_j{\sum }_{i=0}^{\infty }\frac{i^{\underline{j}}}{i!}{x}^i={\sum }_{j=0}^{\infty }{a}_j{x}^j{e}^x$

$=e^x{\sum }_{i=0}^{\infty }{a}_i{x}^i$

所以只需要用普通多项式的系数乘个$e^x$就得到了点值的$EGF$
点值还原原多项式只需要乘一个$e^{-x}$即可

代码

---

其他技巧

多项式快速幂

直接快速幂要多个$log$而且常数大（虽然$ln$和$exp$常数一样大死个仙人）
当$A[0]=1$时（$Ln$要求保证这个），$$A(x{)}^k=Exp(ln(A(x))\ast k)$$
洛谷板子读入时取模

inline poly ksm(poly a,int deg,int k){
	a=exp(Ln(a,deg)*k,deg),a.resize(deg);
	return a;
}

---

分治NTT

---

Bluestein

---

混合基NTT

---

MTT

---

模板合集

const int mod=998244353,g=3;
inline int add(int a,int b){return a+b>=mod?a+b-mod:a+b;}
inline void Add(int &a,int b){a=add(a,b);}
inline int dec(int a,int b){return a>=b?a-b:a-b+mod;}
inline void Dec(int &a,int b){a=dec(a,b);}
inline int mul(int a,int b){return 1ll*a*b>=mod?1ll*a*b%mod:a*b;}
inline void Mul(int &a,int b){a=mul(a,b);}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));return res;}
const int N=100005,C=17;
int *w[18];
int rev[N<<2];
#define poly vector<int> 
#define pb push_back
#define bg begin
inline void init_rev(int lim){
	for(int i=0;i<lim;i++)rev[i]=(rev[i>>1]>>1)|((i&1)*(lim>>1));
}
inline void init_w(){
	for(int i=1;i<=C;i++)
	w[i]=new int[1<<(i-1)];
	int wn=ksm(g,(mod-1)/(1<<C));
	w[C][0]=1;
	for(int i=1;i<(1<<(C-1));i++)w[C][i]=mul(w[C][i-1],wn);
	for(int i=C-1;i;i--)
	for(int j=0;j<(1<<(i-1));j++)
	w[i][j]=w[i+1][j<<1];
}
inline void ntt(poly &f,int lim,int kd){
	for(int i=0;i<lim;i++)if(i>rev[i])swap(f[i],f[rev[i]]);
	for(int mid=1,l=1;mid<lim;mid<<=1,l++)
	for(int i=0;i<lim;i+=(mid<<1))
		for(int j=0,a0,a1;j<mid;j++)
			a0=f[i+j],a1=mul(f[i+j+mid],w[l][j]),
			f[i+j]=add(a0,a1),f[i+j+mid]=dec(a0,a1);
	if(kd==-1&&(reverse(f.bg()+1,f.bg()+lim),1))
		for(int inv=ksm(lim,mod-2),i=0;i<lim;i++)Mul(f[i],inv);
}
inline poly operator +(poly a,poly b){
	poly c;int lim=max(a.size(),b.size());c.resize(lim);
	for(int i=0;i<lim;i++)c[i]=add(a[i],b[i]);return c;
}
inline poly operator -(poly a,poly b){
	poly c;int lim=max(a.size(),b.size());c.resize(lim);
	for(int i=0;i<lim;i++)c[i]=dec(a[i],b[i]);return c;
}
inline poly operator *(poly a,int b){
	for(int i=0;i<a.size();i++)Mul(a[i],b);return a;
}
inline poly operator /(poly a,int b){
	for(int i=0,inv=ksm(b,mod-2);i<a.size();i++)Mul(a[i],inv);
	return a;
}
inline poly operator *(poly a,poly b){
	int deg=a.size()+b.size()-1,lim=1;
	if(deg<=128){
		poly c(deg,0);
		for(int i=0;i<a.size();i++)
		for(int j=0;j<b.size();j++)
			Add(c[i+j],mul(a[i],b[j]));
		return c;
	}
	while(lim<deg)lim<<=1;init(lim);
	a.resize(lim),ntt(a,lim,1);
	b.resize(lim),ntt(b,lim,1);
	for(int i=0;i<lim;i++)Mul(a[i],b[i]);
	ntt(a,lim,-1),a.resize(deg);
	return a;
}
inline poly Inv(poly a,int deg){
	poly c,b(1,ksm(a[0],mod-2));
	for(int lim=4;lim<(deg<<2);lim<<=1){
		init(lim);
		c=a,c.resize(lim>>1);
		c.resize(lim),ntt(c,lim,1);
		b.resize(lim),ntt(b,lim,1);
		for(int i=0;i<lim;i++)Mul(b[i],dec(2,mul(b[i],c[i])));
		ntt(b,lim,-1),b.resize(lim>>1);
	}b.resize(deg);return b;
}
inline poly Sqrt(poly a,int deg){
	poly b(1,1),c,d;
	for(int lim=4;lim<(deg<<2);lim<<=1){
		c=a,c.resize(lim>>1);
		init(lim),d=Inv(b,lim>>1),
		c.resize(lim),ntt(c,lim,1);
		d.resize(lim),ntt(d,lim,1);
		for(int i=0;i<lim;i++)Mul(c[i],d[i]);
		ntt(c,lim,-1),b.resize(lim>>1);
		for(int i=0;i<(lim>>1);i++)b[i]=mul(inv[2],add(b[i],c[i]));
	}b.resize(deg);return b;
}
inline poly operator /(poly a,poly b){
	int lim=1,deg=a.size()-b.size()+1;
	reverse(a.bg(),a.end());
	reverse(b.bg(),b.end());
	while(lim<=deg)lim<<=1;
	b=Inv(b,lim),b.resize(deg);
	a=a*b,a.resize(deg);
	reverse(a.bg(),a.end());
	return a;
}
inline poly operator %(poly a,poly b){
	poly c=a-(a/b)*b;
	c.resize(b.size()-1);
	return c;
}
inline poly deriv(poly a){
	for(int i=0;i<a.size()-1;i++)a[i]=mul(a[i+1],i+1);
	a.pob();return a;
}
inline poly integ(poly a){
	a.pb(0);
	for(int i=a.size()-1;i;i--)a[i]=mul(a[i-1],inv[i]);
	a[0]=0;
	return a;
}
inline poly Ln(poly a,int lim){
	a=integ(deriv(a)*Inv(a,lim)),a.resize(lim);
	return a;
}
inline poly exp(poly a,int deg){
	poly b(1,1),c;int n=a.size();
	for(int lim=2;lim<(deg<<1);lim<<=1){
		c=Ln(b,lim);
		for(int i=0;i<lim;i++)c[i]=dec(i<n?a[i]:0,c[i]);
		Add(c[0],1),b=b*c;
		b.resize(lim);
	}b.resize(deg);return b;
}
inline poly ksm(poly a,int deg,int k){
	a=exp(Ln(a,deg)*k,deg),a.resize(deg);
	return a;
}
poly f[N<<2];
#define mid ((l+r)>>1)
#define lc (u<<1)
#define rc ((u<<1)|1)
inline void build(int u,int l,int r,int *a){
	if(l==r){
		f[u].clear();
		f[u].pb(dec(0,a[l]));
		f[u].pb(1);return;
	}build(lc,l,mid,a),build(rc,mid+1,r,a);
	f[u]=f[lc]*f[rc];
}
inline void calc(int u,int l,int r,poly g,int *a){
	if(l==r){a[l]=g[0];return;}
	calc(lc,l,mid,g%f[lc],a);
	calc(rc,mid+1,r,g%f[rc],a);
}
inline void getans(poly a,int *b,int num){
	build(1,1,num,b),calc(1,1,num,a,b);
}
#undef mid
#undef lc
#undef rc