多项式算法合集

文章目录

• 多项式加减点乘点除
• 多项式乘法
• FFT：
• 前置
• 算法
• 代码实现
• NTT：
• 代码实现
• 多项式求逆：
• 求解过程
• 代码实现
• 多项式开方：
• 求解过程
• 代码实现
• 多项式除法和取模：
• 求解过程
• 代码实现
• 多项式求导与积分
• 代码实现
• 多项式Ln
• 求解过程
• 代码实现
• 多项式Exp
• 前置知识
• 泰勒展开
• 牛顿迭代
• 求解过程
• 代码实现
• 多项式多点求值
• 求解过程
• 代码实现
• 多项式快速插值
• 下降幂多项式乘法
• 其他技巧
• 多项式快速幂
• 分治NTT
• Bluestein
• 混合基NTT
• MTT
• 模板合集

头函数

#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_1x+a_2x^2+a_3x^3……a_nx^n$被称作该多项式的系数表示
我们可以通过带任意一个$x$都可以的到唯一的一个$f(x)$
但$oi$中一般$x$一般都只是一个不定元，不会带入特定值计算，比如用作表示下标之类的

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

---

虚数

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

复数更准确的定义是

$e^{ik}=cos(k)+isin(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,……n-1$
实际上对应的将平面单位圆周上的$n$个点
这些点将单位圆周均分成$n$块，且构成一个正$n$边形

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

单位根几个重要的性质：

1、消去引理

对于任何整数$ngeq0,kgeq 0,dgeq 0$
有$w_{dn}^{dk}=w_{n}^{k}$

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

2、折半引理

如果$ngeq 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、求和引理

对于任意整数$ngeq 1$和$n ot | k$
有$sum_{j=0}^{n-1}(w_{n}^k)^j=0$

考虑这其实是一个等比数列，$w_{n}^{k} ot= 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_1x+a_2x^2+a_3x^3……a_nx^n$
$B(x)=b_0+b_1x+b_2x^2+b_3x^3……b_mx^m$
我们要求一个$n+m$次多项式$C(x)=A(x)*B(x)$
更具体的$C(x)$满足$c_i=sum_{j=0}^{i}a_j*b_{i-j}$
也就是2个数列倒着乘的和，所谓的卷积

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

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

---

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

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

那么我们可以直接在$O(n)$的时间内求出
$C(x_1)=A(x_1)*B(x_1)，……，C(x_{n+1})=A(x_{n+1})*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}^na_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 n2 -1}a_{2i+1}(w_{n}^k)^{2i+1}$

$=sum_{i=0}^{frac n 2-1}a_{2i}(w_n^{2ki})+w_{n}^ksum_{i=0}^{frac n2 -1}a_{2i+1}(w_{n}^{2ki})$

$=sum_{i=0}^{frac n 2-1}a_{2i}(w_{frac n 2}^{k})+w_n^ksum_{i=0}^{frac n2 -1}a_{2i+1}(w_{frac n 2}^k)$

$=A_o(w_{frac n 2}^k)+{w_{n}^k}A_1(w_{frac n2 }^k)$

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

考虑对于$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 n2 -1}a_{2i+1}(w_{n}^{k+frac n2 })^{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 n2 -1}a_{2i+1}(-w_n^k)^{2i}$

$=sum_{i=0}^{frac n 2-1}a_{2i}(w_n^k)^{2i}-w_n^ksum_{i=0}^{frac n2 -1}a_{2i+1}(w_n^k)^{2i}$

$=sum_{i=0}^{frac n 2-1}a_{2i}(w_{frac n 2}^{k})-w_n^ksum_{i=0}^{frac n2 -1}a_{2i+1}(w_{frac n 2}^k)$

$=A_o(w_{frac n 2}^k)-{w_{n}^k}A_1(w_{frac n2 }^k)$

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

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

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

时间复杂度$T(n)=2*T(frac n2 )+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,……，g^{p-2}$互不相同的数
又由于费马定理，对一个素数$p$,有$a^{p-1}equiv 1(mod p)$
这个和单位根很相似

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

但这也就限制了$ntt$的模数必须是$k*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)'equiv 1(mod x^{lceil frac n 2 ceil})$

$A(x)B(x)'-A(x)B(x)equiv 0(mod x^{lceil frac n 2 ceil})$

$B(x)'-B(x)equiv 0(mod x^{lceil frac n 2 ceil})$

平方：

$B(x)'^2+B(x)^2-2B(x)'B(x)equiv 0(mod x^{ n})$
$A(x)(B(x)'^2+B(x)^2-2B(x)'B(x))equiv 0(mod x^{ n})$
$B(x)equiv 2B(x)'-A(x)B(x)'^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)^2equiv A(x)(mod x^n)$

满足$A[0]=1$

求解过程

倍增

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

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

$A(x)+B(x)'^2equiv 2B(x)B(x)'(mod x^n)$
$B(x)equiv frac{A(x)}{2B(x)'}+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^nB(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^nA(x)= x^mB(x)x^{n-m}C(x)+x^{n-m+1}*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'(x)=f'(x)g(x)$
$f'(x)=frac{g'(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+……+frac{f^n(x_0)}{n!}(x-x_0)^n+……$

考虑我们要构造一个函数$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'(x))equiv 0(mod x^{lceil frac n2 ceil})$
要求$f(g(x))equiv0(mod x^n)$

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

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

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

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

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

这就大功告成了

例：

比如多项式开根
就是解$B(x)^2-A(x)equiv 0(mod x^n)$
假设已知$B'(x)^2-A(x)equiv 0(mod x^{lceil frac n 2 ceil})$
这时候$f(B(x))=B(x)^2-A(x),$则$f^1(B(x))=2B(x)$
带入$B(x)equiv B'(x)-frac{B'(x)^2-A(x)}{2B'(x)}equiv frac{B'(x)^2+A(x)}{2B'(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'(x))-g(x)equiv 0(mod x^{lceil frac n 2 ceil})$

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

则$f(x)equiv f'(x)-frac{Ln(f'(x))-g(x)}{frac{1}{f'(x)}}equiv f'(x)-f'(x)Ln(f'(x))+f'(x)g(x)$

$equiv f'(x)(1-Ln(f'(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)……f(a_m)$

求解过程

考虑构造函数$P(x)=prod_{i=1}^{m}(x-a_i)$
显然$forall iin[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 iin[l,mid] ， f(x)\%P_1(x)=f(a_i)$
取模之后$f(x)$的次数减少了一半，继续递归求解即可

$P(x)$可以先分治$ntt$预处理出来
复杂度$T(n)=2*T(frac n 2)+O(nlogn)=O(nlog^2n)$

代码实现

第一次写$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_iprod_{j= ot i}frac{(x-x_j)}{x_i-x_j}$

化一下

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

考虑如何求出$G=prod_{j= ot i}(x_i-x_j)$
设$g(x)=prod_j(x-x_j)$
$j= ot i$就相当于除以了$x-x_i$

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

根据洛必达法则：
若
$lim_{x ightarrow a}f(x)=0,lim_{x ightarrow a}g(x)=0$

有
$lim_{x ightarrow a}frac{f(x)}{g(x)}=lim_{x ightarrow a}frac{f'(x)}{g'(x)}$

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

$f_{l,r}=sum_{i=l}^{r}frac{y_i}{g'(x_i)}prod_{j= ot i}(x-x_j)$

$=prod_{k=mid+1}^{r}(x-x_k)sum_{i=l}^{mid}frac{y_i}{g'(x_i)}prod_{j= ot i}^{[l,mid]}(x-x_j)+prod_{k=l}^{mid}(x-x_k)sum_{i=mid+1}^{r}frac{y_i}{g'(x_i)}prod_{j= ot 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'(x_i)$求出来再分治一波就完了

复杂度$O(nlog^2n)$

代码

---

下降幂多项式乘法

首先考虑对于$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^ne^x$

考虑对于$f(x)=sum_{i=0}^{infty}a_ix^{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_ji^{underline j}$

$=sum_{j=0}^{infty}a_jsum_{i=0}^{infty}frac{i^{underline j}}{i!}x^i=sum_{j=0}^{infty}a_jx^je^x$

$=e^xsum_{i=0}^{infty}a_ix^i$

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

代码

---

其他技巧

多项式快速幂

直接快速幂要多个$log$而且常数大（虽然$ln$和$exp$常数一样大死个仙人）
当$A[0]=1$时（$Ln$要求保证这个），$A(x)^k=Exp(ln(A(x))*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