多项式算法合集

文章目录

头函数

#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);
	a.resize(lim),b.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);
	a.resize(lim),b.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：

前置

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

对于一个nn次多项式f(x)f(x)
f(x)=a0+a1x+a2x2+a3x3……anxnf(x)=a_0+a_1x+a_2x^2+a_3x^3……a_nx^n被称作该多项式的系数表示
我们可以通过带任意一个xx都可以的到唯一的一个f(x)f(x)
但oioi中一般xx一般都只是一个不定元，不会带入特定值计算，比如用作表示下标之类的

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

虚数

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

复数更准确的定义是

eik=cos⁡(k)+isin⁡(k)e^{ik}=cos(k)+isin(k)

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

由于C++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);
    }
};

单位根

fftfft的根本
nn次单位根指的是满足wn=1w^n=1的复数
nn次单位根有nn个,分别表示为wnk，k=0,1,2,……n−1w_{n}^{k}，k=0,1,2,……n-1
实际上对应的将平面单位圆周上的nn个点
这些点将单位圆周均分成nn块，且构成一个正nn边形

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

单位根几个重要的性质：

1、消去引理

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

证明:wdndk=e2πidk/dn=e2πik/n=wnkw_{dn}^{dk}=e^{2pi idk/dn}=e^{2pi ik/n}=w_{n}^{k}

2、折半引理

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

3、求和引理

对于任意整数n≥1ngeq 1和n∤kn ot | k
有∑j=0n−1(wnk)j=0sum_{j=0}^{n-1}(w_{n}^k)^j=0

考虑这其实是一个等比数列，wnk≠1w_{n}^{k} ot= 1时
原式=(wnk)n−1wnk−1=(wnn)k−1wnk−1=0=frac{(w_{n}^{k})^{n}-1}{w_{n}^{k}-1}=frac{(w_{n}^n)^k-1}{w_{n}^{k}-1}=0

而当n∣kn|k时，原式显然为nn

这个性质在后面很重要

算法

考虑对于两个多项式
A(x)=a0+a1x+a2x2+a3x3……anxnA(x)=a_0+a_1x+a_2x^2+a_3x^3……a_nx^n
B(x)=b0+b1x+b2x2+b3x3……bmxmB(x)=b_0+b_1x+b_2x^2+b_3x^3……b_mx^m
我们要求一个n+mn+m次多项式C(x)=A(x)∗B(x)C(x)=A(x)*B(x)
更具体的C(x)C(x)满足ci=∑j=0iaj∗bi−jc_i=sum_{j=0}^{i}a_j*b_{i-j}
也就是2个数列倒着乘的和，所谓的卷积

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

而FFTFFT可以做到O(nlogn)O(nlogn)求出CC

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

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

那么我们可以直接在O(n)O(n)的时间内求出
C(x1)=A(x1)∗B(x1)，……，C(xn+1)=A(xn+1)∗B(xn+1)C(x_1)=A(x_1)*B(x_1)，……，C(x_{n+1})=A(x_{n+1})*B(x_{n+1})

现在我们考虑这样

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

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

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

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

DFT：

考虑我们带入nn次单位根wnw_n
A(wnk)=∑i=0nai(wnk)iA(w_n^k)=sum_{i=0}^na_i(w_n^k)^i
考虑将下标按照奇偶分类
A(wnk)=∑i=0n2−1a2i(wnk)2i+∑i=0n2−1a2i+1(wnk)2i+1A(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}

=∑i=0n2−1a2i(wn2ki)+wnk∑i=0n2−1a2i+1(wn2ki)=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})

=∑i=0n2−1a2i(wn2k)+wnk∑i=0n2−1a2i+1(wn2k)=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)

=Ao(wn2k)+wnkA1(wn2k)=A_o(w_{frac n 2}^k)+{w_{n}^k}A_1(w_{frac n2 }^k)

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

考虑对于A(wnk+n2)A(w_n^{k+frac n 2})
A(wnk+n2)=∑i=0n2−1a2i(wnk+n2)2i+wnk+n2∑i=0n2−1a2i+1(wnk+n2)2iA(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}

考虑单位根的消去引理

wnk+n2=−wnkw_{n}^{k+frac n 2}=-w_n^k

则A(wnk+n2)=∑i=0n2−1a2i(−wnk)2i+(−wnk)∑i=0n2−1a2i+1(−wnk)2iA(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}

=∑i=0n2−1a2i(wnk)2i−wnk∑i=0n2−1a2i+1(wnk)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}

=∑i=0n2−1a2i(wn2k)−wnk∑i=0n2−1a2i+1(wn2k)=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)

=Ao(wn2k)−wnkA1(wn2k)=A_o(w_{frac n 2}^k)-{w_{n}^k}A_1(w_{frac n2 }^k)

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

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

而A0A_0,A1A_1系数都只有n2frac n 2个，所以规模只有原来的一般
递归求解即可

时间复杂度T(n)=2∗T(n2)+O(n)=O(nlogn)T(n)=2*T(frac n2 )+O(n)=O(nlogn)

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

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

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

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

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：

以下摘抄自miskcoomiskcoo神犇（markdownmarkdown太难写了QAQ）

在这里插入图片描述
I是对角矩阵，即对角线上都是1，其他都是0

代码实现

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：

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

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

首先引入原根的概念

对于一个素数pp，其原根gg定义为满足g0,g1,……，gp−2g^0,g^1,……，g^{p-2}互不相同的数
又由于费马定理，对一个素数pp,有ap−1≡1(mod p)a^{p-1}equiv 1(mod p)
这个和单位根很相似

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

但这也就限制了nttntt的模数必须是k∗2n+1k*2^n+1的形式
否则就要做任意模数nttntt
做法就是选出几个模数分别做之后CrtCrt合并答案（并不会）
代码和fftfft类似
由于wn−k=wnn−kw_n^{-k}=w_n^{n-k}
所以我们可以先做的时候不管正逆变换
最后把11~n−1n-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);
}

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

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

预处理原根

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%25\%到45%45\%不等

其实nttntt本身常数不算很大，运算常数大概也只有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−1n-1次多项式A(x)A(x)，求多项式B(x)B(x)满足：

A(x)B(x)≡1(mod xn)A(x)B(x)equiv 1(mod x^n)

求解过程

倍增：

若已知A(x)B(x)′≡1(mod x⌈n2⌉)A(x)B(x)'equiv 1(mod x^{lceil frac n 2 ceil})

A(x)B(x)′−A(x)B(x)≡0(mod x⌈n2⌉)A(x)B(x)'-A(x)B(x)equiv 0(mod x^{lceil frac n 2 ceil})

B(x)′−B(x)≡0(mod x⌈n2⌉)B(x)'-B(x)equiv 0(mod x^{lceil frac n 2 ceil})

平方：

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

复杂度T(n)=T(n2)+O(nlogn)=O(nlogn)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−1n-1次多项式A(x)A(x),求一个mod xnmod x^n意义下的多项式B(x)B(x)满足

B(x)2≡A(x)(mod xn)B(x)^2equiv A(x)(mod x^n)

满足A[0]=1A[0]=1

求解过程

倍增

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

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

A(x)+B(x)′2≡2B(x)B(x)′(mod xn)A(x)+B(x)'^2equiv 2B(x)B(x)'(mod x^n)
B(x)≡A(x)2B(x)′+B(x)2(mod xn)B(x)equiv frac{A(x)}{2B(x)'}+frac{B(x)}{2}(mod x^n)

复杂度T(n)=T(n2)+O(nlogn)=O(nlogn)T(n)=T(frac n 2)+O(nlogn)=O(nlogn)

代码实现

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

多项式除法和取模：

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

求解过程

考虑对一个最高次数为nn的多项式B(x)B(x)操作BR(x)=xnB(1x)B^R(x)=x^nB(frac 1 x)
会发现BR(x)B^R(x)只是B(x)B(x)的系数反转的柿子

考虑A(x)=B(x)C(x)+D(x)A(x)= B(x)C(x)+D(x)
A(x)A(x)最高项为nn,B(x)B(x)最高项为mm,则C(x)C(x)最高项为n−mn-m,D(x)D(x)为m−1m-1
两边同时乘一个xnx^n,并带入1xfrac 1 x
xnA(x)=xmB(x)xn−mC(x)+xn−m+1∗xm−1D(x)x^nA(x)= x^mB(x)x^{n-m}C(x)+x^{n-m+1}*x^{m-1}D(x)
AR(x)=BR(x)CR(x)+xn−m+1DR(x)A^R(x)=B^R(x)C^R(x)+x^{n-m+1}D^R(x)
考虑在mod xn−m+1mod x^{n-m+1}意义下，A,BA,B已知，CC最高为n−mn-m不影响，而DD被消去
则AR(x)≡BR(x)CR(x)(mod xn−m+1)A^R(x)equiv B^R(x)C^R(x)(mod x^{n-m+1})

多项式求逆就可以了

复杂度O(nlogn)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)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−1n-1次多项式g(x)g(x),求一个mod xnmod x^n意义下的多项式f(x)f(x)，满足：
f(x)≡Ln(g(x))(mod xn)f(x)equiv Ln(g(x))(mod x^n)

求解过程

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

复杂度O(nlogn)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(x0)+f1(x0)1!(x−x0)+f2(x0)2!(x−x0)2+……+fn(x0)n!(x−x0)n+……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)使g(x)g(x)完全拟合f(x)f(x)
那首先初始点x0x_0的值要和f(x0)f(x_0)一样
在此基础上只需要保证一阶导数，二阶导数……都完全相同即可
即∀i,fi(x)=gi(x)forall i,f^i(x)=g^i(x)
由于求gg第ii阶导数时为i!ai=fi(x)i!a_i=f^i(x)
即ai=fi(x)i!a_i=frac{f^i(x)}{i!}
所以得证

牛顿迭代

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

假设有函数ff和一个多项式g(x)mod xng(x)mod x^n
满足f(g(x))≡0(mod xn)f(g(x))equiv 0 ( mod x^n)
已知ff，求gg

说白了就是用来解f(g(x))≡0(mod xn)f(g(x))equiv 0(mod x^n)之类的方程

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

假设已经知道f(g′(x))≡0(mod x⌈n2⌉)f(g'(x))equiv 0(mod x^{lceil frac n2 ceil})
要求f(g(x))≡0(mod xn)f(g(x))equiv0(mod x^n)

考虑将f(g(x))f(g(x))在g′(x)g'(x)处泰勒展开
f(g(x))=f(g′(x))+f1(g′(x))1!(g(x)−g′(x))+f2(g′(x))2!(g(x)−g(x)′)2+……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)≡g′(x)(mod x⌈n2⌉)g(x)equiv g'(x)(mod x^{lceil frac n 2 ceil})
所以g(x)−g′(x)g(x)-g'(x)最低项次数一定大于⌈n2⌉lceil frac n 2 ceil

则在mod xnmod x^n意义下，整个式子从(g(x)−g(x)′)2(g(x)-g(x)')^2开始都为00了
则f(g(x))≡f(g′(x))+f1(g′(x))(g(x)−g′(x))(mod xn)f(g(x))equiv f(g'(x))+{f^1(g'(x))}(g(x)-g'(x))(mod x^n)

又f(g(x))≡0(mod xn)f(g(x))equiv 0 (mod x^n)

g(x)≡g′(x)−f(g′(x))f1(g′(x))(mod xn)g(x)equiv g'(x)-frac{f(g'(x))}{f^1(g'(x))}(mod x^n)

这就大功告成了

例：

比如多项式开根
就是解B(x)2−A(x)≡0(mod xn)B(x)^2-A(x)equiv 0(mod x^n)
假设已知B′(x)2−A(x)≡0(mod x⌈n2⌉)B'(x)^2-A(x)equiv 0(mod x^{lceil frac n 2 ceil})
这时候f(B(x))=B(x)2−A(x),f(B(x))=B(x)^2-A(x),则f1(B(x))=2B(x)f^1(B(x))=2B(x)
带入B(x)≡B′(x)−B′(x)2−A(x)2B′(x)≡B′(x)2+A(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−1n-1次多项式g(x)g(x),求一个mod xnmod x^n意义下的多项式f(x)f(x)，满足：
f(x)≡eg(x)(mod xn)f(x)equiv e^{g(x)}(mod x^n)

也就是Ln(f(x))≡g(x)(mod xn)Ln(f(x))equiv g(x)(mod x^n)

求解过程

倍增：

考虑原问题，即求解方程Ln(f(x))−g(x)≡0(mod xn)Ln(f(x))-g(x)equiv 0(mod x^n)

同样假设已经知道Ln(f′(x))−g(x)≡0(mod x⌈n2⌉)Ln(f'(x))-g(x)equiv 0(mod x^{lceil frac n 2 ceil})

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

则f(x)≡f′(x)−Ln(f′(x))−g(x)1f′(x)≡f′(x)−f′(x)Ln(f′(x))+f′(x)g(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)

≡f′(x)(1−Ln(f′(x))+g(x))(mod xn)equiv f'(x)(1-Ln(f'(x))+g(x)) (mod x^n)

复杂度T(n)=T(n2)+O(nlogn)=O(nlogn)T(n)=T(frac n 2)+O(nlogn)=O(nlogn)

代码实现

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

多项式多点求值

已知一个nn次多项式f(x)f(x)，求f(a1)……f(am)f(a_1)……f(a_m)

求解过程

考虑构造函数P(x)=∏i=1m(x−ai)P(x)=prod_{i=1}^{m}(x-a_i)
显然∀i∈[1,m],P(ai)=0forall iin[1,m],P(a_i)=0

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

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

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

考虑分治优化
P1(x)=∏i=lmid(x−ai),P2(x)=∏i=mid+1r(x−ai)P_1(x)=prod_{i=l}^{mid}(x-a_i),P_2(x)=prod_{i=mid+1}^{r}(x-a_i)
∀i∈[l,mid]，f(x)%P1(x)=f(ai)forall iin[l,mid] ， f(x)\%P_1(x)=f(a_i)
取模之后f(x)f(x)的次数减少了一半，继续递归求解即可

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

代码实现

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

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

多项式快速插值

写在这里了

下降幂多项式乘法

写在这里了

其他技巧

多项式快速幂

直接快速幂要多个loglog而且常数大（虽然lnln和expexp常数一样大死个仙人）
当A[0]=1A[0]=1时（LnLn要求保证这个），A(x)k=Exp(ln(A(x))∗k)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;
}

对于A[0]≠1A[0] ot=1也可以很简单的做
写在这里了

MTTMTT

写在这里了

模板合集

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);
	a.resize(lim),b.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);
	a.resize(lim),b.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;int inv2=ksm(2,mod-2);
	for(int lim=4;lim<(deg<<2);lim<<=1){
		c=a,c.resize(lim>>1);
		init_rev(lim),d=Inv(b,lim>>1),
		c=c*d,b.resize(lim);
		for(int i=0;i<lim;i++)b[i]=mul(inv2,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
【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/stargazer-cyk/p/12328281.html