多项式从入门到入坟

参考

MCAdam的FFT和NTT学习笔记

NTT与多项式全家桶 - command_block

基本概念

多项式一般有三种表达方式：

系数表达式：即用一个系数数组 $A$ 表示这个多项式。
即 $F(x)=\sum_{i=0}^n A_ix^i$
点值表达式：因为 $n+1$ 个点确定一个多项式，所以可以用一个 $n+1$ 个点表示一个多项式，显然点值表达式不唯一。
下降幂表达式，可以通过斯特林数与系数表达式相互转化。

卷积

当两个多项式（设系数数组为 $F,G$ 分别为 $n$ 次和 $m$ 次），当二者相乘的时候，得到的多项式的第 $i$ 项系数为： $\sum_{j+k=i}F_j \times G_k$

该形式就是卷积。

正常两个多项式相乘需要 $O(n^2)$ 的复杂度，而利用多项式算法可以优化成 $O(n\log_n)$ 。

FFT

接下来看如何优化卷积。

复数

基础概念请自行翻阅必修 2。

一个复数可以表示成 $a+bi$ 的形式。

加法： $(a+bi)+(c+di)=(a+c)+(b+d)i$
减法： $(a+bi)-(c+di)=(a-c)+(b-d)i$
乘法： $(a+bi)\times(c+di)=(ac-bd)+(ad \times bc)i$

两个复数相乘，模长相乘，幅角相加

struct cmplx{
	double x,y;
	cmplx(){}
	cmplx(double a,double b){x=a,y=b;}
	cmplx operator + (const cmplx &A)const{return cmplx(x+A.x,y+A.y);}
	cmplx operator - (const cmplx &A)const{return cmplx(x-A.x,y-A.y);}
	cmplx operator * (const cmplx &A)const{return cmplx(x*A.x-y*A.y,x*A.y+y*A.x);}
};

单位根

显然， $x^n=1$ 在复数域内有且只有 $n$ 个解，我们称这些解为 $n$ 次单位根，其幅角为 $\frac{a}{n}\times 2\pi(0 \le a <n,a \in Z)$

记第 $k$ 个单位根为 $\omega_n^k$ 。

单位根有以下性质(可以用几何意义来理解)：

$\omega_n^i \times \omega_n^j=\omega_n^{i+j}$

$\omega_{pn}^{pk}=\omega_n^k$

$\omega_n^{k+n/2}=-\omega_n^k$

$\omega_n^k=\omega_n^{k\bmod n}$

而因为

$\omega_n^1= cos\frac{2\pi}{n}+sin\frac{2\pi}{n}i$

$\omega_n^k=(\omega_n^1)^k$

就可以得到所有单位根。

用于快速把系数表达式转成点值表达式。

以下的 $n$ 都满足 $2^k=n,k \in N$ 。

对于一个多项式 $A(x)=\sum_{i=0}^{n-1}x^iF(i)$ ，
$A(x)=F(0)+xF(1)+x^2F(2)+...+x^{n-1}F(n-1)$ 。

我们把它奇偶分开。

$A1(x)=F(0)+xF(2)+...+x^{\frac{n}{2}-1}F(n-2)$

$A2(x)=F(1)+xF(3)+...+x^{\frac{n}{2}-1}F(n-1)$

显然 $A(x)=A1(x^2)+xA2(x^2)$

设 $k<\frac{n}{2}$

把 $\omega_n^k$ 带入

$A(\omega_n^k)=A1((\omega_n^k)^2)+\omega_n^kA2((\omega_n^k)^2)$

$=A1(\omega_n^{2k})+\omega_n^kA2(\omega_n^{2k})$

$=A1(\omega_{n/2}^k)+\omega_n^kA2(\omega_{n/2}^k)$

把 $\omega_n^{k+n/2}$ 带入

$A(\omega_n^{k+n/2})=A1((\omega_n^{k+n/2})^2)+\omega_n^{k+n/2}A2((\omega_n^{k+n/2})^2)$

$=A1(\omega_n^{2k+n})+\omega_n^{k+n/2}A2(\omega_n^{2k+n})$

$\because \omega_n^k=\omega_n^{k\ mod\ n}$

$\therefore A(\omega_n^{k+n/2})=A1(\omega_n^{2k})+\omega_n^{k+n/2}A2(\omega_n^{2k})$

$=A1(\omega_{n/2}^k)+\omega_n^{k+n/2}A2(\omega_{n/2}^k)$

$=A1(\omega_{n/2}^k)-\omega_n^kA2(\omega_{n/2}^k)$

综上所述：

$A(\omega_n^k)=A1(\omega_{n/2}^k)+\omega_n^kA2(\omega_{n/2}^k)$

$A(\omega_n^{k+n/2})=A1(\omega_{n/2}^k)-\omega_n^kA2(\omega_{n/2}^k)$

就可以用分治在 $O(n\log_n)$ 的时间内求出 $n$ 个点值，就得到 $A$ 的点值表达式。

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const int N=1e6+10;
const double pi=acos(-1);
struct cmplx
{
    double x,y;
    cmplx (double a=0,double b=0){ x=a,y=b; }
    cmplx operator + (cmplx const &a) const{ return cmplx(x+a.x,y+a.y); }
    cmplx operator - (cmplx const &a) const{ return cmplx(x-a.x,y-a.y); }
    cmplx operator * (cmplx const &a) const{ return cmplx(x*a.x-y*a.y,x*a.y+y*a.x); }
}v[N],sav[N];
inline void DFT(cmplx *f,int len)
{
    if(len==1) return;
    cmplx *f1=f,*f2=f+len/2;//用指针，改f1/f2，f也会跟着改；改f，f1/f2也会改（它们用同样的地址） 
    for(int i=0;i<len;i++)//所以要开临时的数组 
        sav[i]=f[i];
    for(int i=0;i<len/2;i++)
        f1[i]=sav[i<<1],f2[i]=sav[i<<1|1];//奇偶分开 
    DFT(f1,len/2),DFT(f2,len/2);//分治 
    cmplx bas=cmplx(1,0),buf=cmplx(cos(2*pi/len),sin(2*pi/len));//算单位根 
    for(int i=0;i<len/2;i++,bas=bas*buf)
        sav[i]=f1[i]+bas*f2[i],sav[i+len/2]=f1[i]-bas*f2[i];//套那两个式子 
    for(int i=0;i<len;i++)
        f[i]=sav[i];
}
int main()
{
    int n,m=1;
    scanf("%d",&n);
    for(int i=0;i<n;i++)
        scanf("%lf",&v[i].x);//输入系数 
    while(m<n) m<<=1;//补全 
    DFT(v,m);
    return 0;
}

注意： $n=2^k$ ， $n$ 为项数！

IDFT

其实是DFT的逆运算。

设 $F(x)$ DFT 之后得到 $G(x)$ 。

单位根反演乱搞可得：

F (x) [p] = \frac{1}{n} \sum_{i = 0}^{n - 1} (ω_{n}^{- p})^{i} G (x) [i] = \frac{1}{n} G (ω_{n}^{- p})

$F(x)[p]=\frac{1}{n}\sum_{i=0}^{n-1}(\omega_n^{-p})^iG(x)[i]=\frac{1}{n}G(\omega_n^{-p})$

长得很像 DFT，所以代码也差不多（只差一个符号）。

因为DFT和IDFT长得太像，所以把他们合并成FFT。

而且三角函数可以预处理。

通过找规律，初始序列最后变换到的位置是可以算出来的,这就叫蝴蝶变换：

for(int i=0;i<len;i++) cir[i]=(cir[i>>1]>>1)|(i&1?len>>1:0);

所以我们就得到了FFT,下面是迭代版的。

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const int N=4e6+10;
int cir[N];
double tsin[N],tcos[N];
const double pi=acos(-1);
struct cmplx
{
    double x,y;
    cmplx (double a=0,double b=0){ x=a,y=b; }
    cmplx operator + (cmplx const &a) const{ return cmplx(x+a.x,y+a.y); }
    cmplx operator - (cmplx const &a) const{ return cmplx(x-a.x,y-a.y); }
    cmplx operator * (cmplx const &a) const{ return cmplx(x*a.x-y*a.y,x*a.y+y*a.x); }
}F[N],G[N];
inline void FFT(cmplx *f,int len,double tag)
{
    for(register int i=0;i<len;i++) if(i<cir[i]) swap(f[i],f[cir[i]]);
    for(register int l=2;l<=len;l<<=1)
    {
        int size=l>>1;
        cmplx buf=cmplx(tcos[l],tag*tsin[l]);
        for(register int i=0;i<len;i+=l)
        {
            cmplx bas=cmplx(1,0); 
            for(register int j=i;j<i+size;j++,bas=bas*buf)
            {
                cmplx tmp=bas*f[j+size];
                f[j+size]=f[j]-tmp;
                f[j]=f[j]+tmp;
            }
        }
    }
}
int main()
{
    int n,m,len=1;
    scanf("%d%d",&n,&m);
    for(register int i=0;i<=n;i++)
        scanf("%lf",&F[i].x);
    for(register int i=0;i<=m;i++)
        scanf("%lf",&G[i].x);
    while(len<=n+m) len<<=1;
    for(register int i=len;i>=1;i>>=1) tsin[i]=sin(2*pi/i),tcos[i]=cos(2*pi/i);
    for(register int i=0;i<len;i++) cir[i]=(cir[i>>1]>>1)|((i&1)?len>>1:0);
    FFT(F,len,1),FFT(G,len,1);
    for(register int i=0;i<len;i++) F[i]=F[i]*G[i];
    FFT(F,len,-1);
    for(register int i=0;i<=n+m;i++) printf("%d ",(int)(F[i].x/len+0.49));
    return 0;
}

在 OI 中我们讨厌实数的运算，于是有了 NTT。

我们首先要找到一个东西替代单位根，于是找到了原根。

其实我们就是要找到一个东西像单位根一样转 $n$ 次就回来的东西，这不就是原根吗，转 $p-1$ 回来，那我们设置转一次是 $\frac{p-1}{n}$ ，那不就刚好转 $n$ 次就回来了吗。

这需要 $n\mid p-1$ ，而 $n=2^k$ ，所以 $p-1$ 中需要很多的 $2$ 。

一般我们都是用 $998244353$ 为模数：

$998244353$ 大小合适，且是质数
$998244353-1=2^{23}\times 119$

#include<iostream>
#include<stdio.h>
#define p 998244353
#define g 3
#define invg 332748118
#define N 4000005
#define int long long
using namespace std;
inline int read(){
	int x=0,f=1;
	char ch=getchar();
	while(ch<'0' || ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	return x*f;
}
int n,m,F[N],G[N],cir[N];
int power(int x,int b){
	int res=1;
	while(b){
		if(b&1) res=res*x%p;
		b>>=1;
		x=x*x%p;
	}
	return res;
}
void NTT(int *f,int lenth,int tag){
	for(int i=0;i<lenth;++i) if(i<cir[i]) swap(f[i],f[cir[i]]);
	for(int len=2;len<=lenth;len<<=1){
		int sze=len>>1,buf=power((tag?g:invg),(p-1)/len);
		for(int i=0;i<lenth;i+=len){
			int bas=1;
			for(int j=i;j<i+sze;j++,bas=bas*buf%p){
				int tmp=bas*f[j+sze]%p;
				f[j+sze]=(f[j]-tmp+p)%p;
				f[j]=(f[j]+tmp)%p;
			}
		}
	}
}
signed main(){
	n=read(),m=read();
	for(int i=0;i<=n;++i) F[i]=read();
	for(int i=0;i<=m;++i) G[i]=read();
	int len=1;
	while(len<=n+m) len*=2;
	for(int i=0;i<len;++i) cir[i]=(cir[i>>1]>>1)|((i&1)?len>>1:0);
	NTT(F,len,1),NTT(G,len,1);
	for(int i=0;i<len;++i) F[i]=F[i]*G[i]%p;
	NTT(F,len,0);
	for(int i=0;i<n+m+1;++i) printf("%lld ",F[i]*power(len,p-2)%p); 
	return 0;
}

封装后的版本

#define mo 998244353
#define g 3
#define invg 332748118
int cir[N],rt;
int power(int x,int b){
    int res=1;
    while(b){
        if(b&1) res=res*x%mo;
        b>>=1;
        x=x*x%mo;
    }
    return res;
}
void clear(int *A,int x,int y){
    for(int i=x;i<y;++i) A[i]=0;
}
void copy(int *A,int *B,int x,int y){
    for(int i=x;i<y;++i) A[i]=B[i];
}
void mul(int *A,int *B,int x,int y){
    for(int i=x;i<y;++i) A[i]=A[i]*B[i]%mo;
}
void NTT(int *f,int lenth,int tag){
    if(lenth!=rt)for(int i=0;i<lenth;++i)cir[i]=(cir[i>>1]>>1)|((i&1)?lenth>>1:0);
    rt=lenth;
    for(int i=0;i<lenth;++i) if(i<cir[i]) swap(f[i],f[cir[i]]);
    for(int len=2;len<=lenth;len<<=1){
        int sze=len>>1,buf=power((tag?g:invg),(mo-1)/len);
        for(int i=0;i<lenth;i+=len){
            int bas=1;
            for(int j=i;j<i+sze;j++,bas=bas*buf%mo){
                int tmp=bas*f[j+sze]%mo;
                f[j+sze]=(f[j]-tmp+mo)%mo;
                f[j]=(f[j]+tmp)%mo;
            }
        }
    }
    if(tag==1) return;
    int invn=power(lenth,mo-2);
    for(int i=0;i<lenth;++i) f[i]=f[i]*invn%mo;
}
int tmul[N];
void Polymul(int *A,int *B,int len1,int len2){
    #define sav tmul
    int len=1;
    while(len<(len1<<1)) len<<=1;
    copy(sav,B,0,len);
    NTT(A,len,1),NTT(sav,len,1);
    for(int i=0;i<len;++i) A[i]=A[i]*sav[i]%mo;
    NTT(A,len,0);
    clear(A,len2,len),clear(sav,0,len);
    #undef sav
}

多项式求逆

传送门

设 $G(x)$ 为答案，当 $F(x)$ 只有1项时，答案为 $F(x)$ 的逆元，其他情况考虑倍增。

设 $mod\ x^{\lceil \frac{n}{2}\rceil}$ 时的答案为 $H(x)$

∵ F (x) G (x) \equiv 1, F (x) H (x) \equiv 1 (m o d x^{⌈ \frac{n}{2} ⌉})

$\because F(x)G(x)\equiv1,F(x)H(x)\equiv1\ (mod\ x^{\lceil \frac{n}{2}\rceil})$

∴ G (x) - H (x) \equiv 0 (m o d x^{⌈ \frac{n}{2} ⌉})

$\therefore G(x)-H(x)\equiv0\ (mod\ x^{\lceil \frac{n}{2}\rceil})$

(G (x) - H (x))^{2} \equiv 0 (m o d x^{n})

$(G(x)-H(x))^2\equiv0\ (mod\ x^n)$

G (x)^{2} - 2 G (x) H (x) + H (x)^{2} \equiv 0 (m o d x^{n})

$G(x)^2-2G(x)H(x)+H(x)^2\equiv0\ (mod\ x^n)$

两边同乘 $F(x)$

G (x) - 2 H (x) + H (x)^{2} F (x) \equiv 0 (m o d x^{n})

$G(x)-2H(x)+H(x)^2F(x)\equiv0\ (mod\ x^n)$

G (x) = 2 H (x) - H (x)^{2} F (x)

$G(x)=2H(x)-H(x)^2F(x)$

int inv1[N],inv2[N],inv3[N];
void Polyinv(int *A,int len){
	#define sav1 inv1
	#define sav2 inv2
	#define sav3 inv3
    int lim=1;while(lim<len) lim<<=1;
    sav1[0]=power(A[0],mo-2);
    for(int l=2;l<=lim;l<<=1){
        int r=l<<1;copy(sav3,A,0,l);
        for(int i=0;i<(l>>1);++i) sav2[i]=(sav1[i]<<1)%mo;
        NTT(sav1,r,1),mul(sav1,sav1,0,r);
        NTT(sav3,r,1),mul(sav1,sav3,0,r);
        NTT(sav1,r,0),clear(sav1,l,r);
        for(int i=0;i<l;++i) sav1[i]=(sav2[i]-sav1[i]+mo)%mo;
    }
    copy(A,sav1,0,len);
    clear(sav1,0,lim<<1),clear(sav2,0,lim<<1),clear(sav3,0,lim<<1);
    #undef sav1
    #undef sav2
    #undef sav3
}

分治FFT

传送门

这题看似是 $F=F*G$ 。
但仔细想想却不是，因为如果按照这样 $G[0]$ 应该是 $0$ ，但题目又说 $F[0]=1$ ，所以并不能直接卷。

所以我们把 $F 和 G$ 卷一下。

F (x) * G (x) = \sum_{k = 0}^{\infty} \sum_{i + j = k} F [i] G [j] x^{k}

$F(x)*G(x)=\sum_{k=0}^{\infty}\sum_{i+j=k}F[i]G[j]x^k$

当 $k=0$ ，那一项为 $0（G[0]=0）$
当 $k>0$ ，为 $\sum_{i=0}^kF[k-i]G[i]=\sum_{i=1}^kF[k-i]G[i]=F[k]$ 。

所以 $F*G$ 和 $F$ 只有在 $k=0$ 时差了个1，所以 $F=F*G+1$ 。

移个项得： $F=\dfrac{1}{1-G}$ ，写个求逆就行了。

但是分治 FFT 的板子怎么能写求逆呢？而且求逆不是能完全代替分治 FFT 的。

分治 FFT 的本质是解方程，是一种半在线卷积。

考虑 $F=F*G+1$ ，模仿 CDQ 分治的过程，先算左边，再算左边对右边的贡献，最后再算右边。

具体地说，当我们在处理 $[l,r]$ ，先 solve(l,mid)，然后把 $F[l,mid]$ 和 $G[0,r-l]$ 卷起来贡献到 $F[mid+1,r]$ ，再 solve(mid+1,r)。

但是有些题会出现自己卷自己的情况，例如 #50. 【UR #3】链式反应。

这题如果不管一些边界条件和细节的话，其实就是 $F=\frac{1}{2}AF^2$ 。

首先这不能求逆解决了，考虑分治 FFT，分两种情况考虑：

$l\not=0$ ，这时候 $[0,r-l]$ 肯定已经算出来了，直接贡献即可。
$l=0$ ，我们可以先只把 $[l,mid]$ 卷 $[l,mid]$ 的贡献算上，剩下的留给情况 1。所以情况 1 时不仅要考虑 $[l,mid]$ 对 $[mid+1,r]$ 的贡献，还要考虑 $[0,r-l]$ 对 $[mid+1,r]$ 的贡献，但因为 $[l,mid]$ 卷 $[0,r-l]$ 和 $[-,r-l]$ 卷 $[l,mid]$ 是一样的，所以我们需要做的只是在情况 1 时的贡献乘一个 2 而已。

void solve(int l,int r){
	if(l==r) return;
	int mid=l+r>>1;
	solve(l,mid);
	if(l!=0){
		for(int i=l;i<=mid;++i) tmp1[i-l]=F[i];
		for(int i=0;i<=r-l;++i) tmp2[i]=F[i];
		polyMul(tmp1,tmp2,mid-l+1,r-l+1,r-l+1);
		polyMul(tmp1,P,r-l+1,r-l+1,r-l+1);
		for(int i=mid+1;i<=r;++i) (F[i]+=tmp1[i-l-1]*fac[i-1]%mo*ifac[i]%mo)%=mo;
		clear(tmp1,0,r-l+1),clear(tmp2,0,r-l+1);
	}
	else{
		for(int i=l;i<=mid;++i) tmp1[i-l]=F[i];
		for(int i=l;i<=mid;++i) tmp2[i-l]=F[i];
		polyMul(tmp1,tmp2,mid-l+1,mid-l+1,r-l+1);
		polyMul(tmp1,P,r-l+1,r-l+1,r-l+1);
		for(int i=mid+1;i<=r;++i) (F[i]+=tmp1[i-l-1]*fac[i-1]%mo*ifac[i]%mo*inv2%mo)%=mo;
		clear(tmp1,0,r-l+1),clear(tmp2,0,r-l+1);
	}
	solve(mid+1,r);
}

其实形如 $F=A*F$ ，其中 $A$ 是已知的，我们可以认为这是一次的类型。形如 $F=A*F^2$ ，这种是二次的类型，以此类推，其实二次就是自己卷自己的类型，所以上文讲了一次和二次的解法。

但方程不仅有几次，还有几元，比如 #428. 【集训队作业2018】普通的计数题。这题本质是一个二元二次： $F=A*G,G=F*G$ 。多元和一元的解法没多大区别，但是要注意，这里 $l\not=0$ 时的贡献不能直接乘 2 了，因为原来我们直接乘 2 是因为 $F[l,mid]*F[0,r-l]=F[0,r-l]*F[l,mid]$ ，而这里 $F[l,mid]*G[0,r-l]=G[l,mid]*F[0,r-l]$ 并不成立！所以我们只能重新算一遍，而不能直接乘以 2！

至于三次的解法其实很简单，例如 $F=F^3$ ，设 $F_2=F^2,F=F_2*F$ 当二元二次解就好了。但是例如 $F=F^8$ ，其实不需要分别记 $F_2,F_3\ldots F_8$ ，其实只需要记 $F_2,F_4,F_8$ 就好了，利用二进制拆分的原理可以优化一下常数。

多项式除法

传送门

给定 $n$ 次多项式 $F(x)$ 和 $m$ 次的 $G(x)$ ,求 $n-m$ 次的 $Q(x)$ 和小于 $m$ 次的 $R(x)$ 。

满足 $F(x)\equiv G(x)Q(x)+R(x)\ (mod\ x^n)$

思路是想办法让 $R(x)$ 被模掉，但只有高次的会被模掉，而 $R(x)$ 是最低次的。但是我们可以把所有多项式都反转，这样 $R(x)$ 就是最高次的了！

设 $F_0(x)$ 是 $F(x)$ 反转之后的多项式。

显然 $F_0(x)=\sum_{i=0}^{n}F[n-i]x^i$ 。

但我们需要更好的表达方式，发现 $F_0(x)=x^nF(\dfrac{1}{x})$ ，展开一下即可证明。

把上面的式子用 $\dfrac{1}{x}$ 代替，并两边乘 $x^n$ 。

x^{n} F (\frac{1}{x}) = x^{n} G (\frac{1}{x}) Q (\frac{1}{x}) + x^{n} R (\frac{1}{x})

$x^nF(\dfrac{1}{x})=x^nG(\dfrac{1}{x})Q(\dfrac{1}{x})+x^nR(\dfrac{1}{x})$

左边就是 $F_0(x)$ ，而 $G(x)$ 是m次， $Q(x)$ 是n-m次，加起来刚好 n 次， $R(x)$ 可以当成 m-1 次的，所以：

F_{0} (x) = G_{0} (x) Q_{0} (x) + x^{n - m + 1} R_{0} (x)

$F_0(x)=G_0(x)Q_0(x)+x^{n-m+1}R_0(x)$

这下 $x^{n-m+1}R_0(x)$ 是最高次了，所以：

F_{0} (x) \equiv G_{0} (x) Q_{0} (x) (m o d x^{n - m + 1})

$F_0(x)\equiv G_0(x)Q_0(x)\ (mod\ x^{n-m+1})$

求出 $Q(x)$ 后代回去得出 $R(x)$ 。

int div1[N],div2[N],rr[N];
void rev(int *A,int len){
	#define sav rr
	copy(sav,A,0,len);
	for(int i=0;i<len;++i) A[i]=sav[len-i-1];
	#undef sav
}
void Polydiv(int *A,int *B,int len1,int len2){
	#define sav1 div1
	#define sav2 div2
	int len0=len1-len2+1;
	rev(A,len1),copy(sav1,A,0,len0),rev(A,len1);
	rev(B,len2),copy(sav2,B,0,len0),rev(B,len2);
	Polyinv(sav2,len0),Polymul(sav2,sav1,len0,len0),rev(sav2,len0);
	Polymul(B,sav2,len1,len1);
	for(int i=0;i<len2-1;++i) B[i]=(A[i]-B[i]+mo)%mo;
	copy(A,sav2,0,len0),clear(A,len0,len1),B[len2-1]=0;
	clear(sav1,0,len1),clear(sav2,0,len1);
	#undef sav1
	#undef sav2
}

多项式求导与积分

不多说，大家都会。

int invcnt=1,inv[N];
inline void get_inv(int len){
    inv[1]=1; if(invcnt>=len) return;
    for(int i=invcnt+1;i<=len;++i) inv[i]=(mo-mo/i)*inv[mo%i]%mo;
}
inline void Polyder(int *A,int len){
    for(int i=1;i<=len;i++) A[i-1]=A[i]*i%mo; A[len]=0;
}
inline void Polyint(int *A,int len){
    get_inv(len);
    for(int i=len;i>=1;i--) A[i]=A[i-1]*inv[i]%mo; A[0]=0;
}

多项式对数函数

传送门

给定 $F(x)$ ，求 $G(x)$ 使得 $G(x) \equiv \ln F(x)\ \pmod {x^n}$

因为 $\ln$ 求导后很好看，所以两边求导。

G^{'} (x) \equiv \frac{F^{'} (x)}{F (x)}

$G'(x) \equiv \dfrac{F'(x)}{F(x)}$

G (x) \equiv \int \frac{F^{'} (x)}{F (x)}

$G(x) \equiv \int\dfrac{F'(x)}{F(x)}$

int ln1[N];
inline void Polyln(int *A,int len){
	#define sav ln1
	copy(sav,A,0,len),Polyder(A,len),Polyinv(sav,len);
	Polymul(A,sav,len,len),Polyint(A,len-1),clear(sav,0,len);
	#undef sav
}

多项式牛顿迭代

求 $F(x)$ 满足 $G(F(x)) \equiv 0\ (mod\ x^n)$

考虑递归求解，假设我们已经求出 $\mod x^{\lceil\frac{n}{2}\rceil}$ 的解 $F_0(x)$ 。

在 $F_0(x)$ 处套泰勒展开：

G (F (x)) = G (F_{0} (x)) + (G (F_{0} (x)))^{'} (F (x) - F_{0} (x))

$G(F(x))=G(F_0(x))+(G(F_0(x)))'(F(x)-F_0(x))$

+ \frac{(G (F_{0} (x)))^{″}}{2} (F (x) - F_{0} (x))^{2} + . . .

$\qquad \qquad \quad+\dfrac{(G(F_0(x)))''}{2}(F(x)-F_0(x))^2+...$

已知 $F(x)$ 和 $F_0(x)$ 的后 $\lceil\dfrac{n}{2}\rceil$ 项是一样的，所以 $F(x)-F_0(x)$ 的大于等于 2 的次方不为零的项已经在 n 次以上了，在 $\mod x^n$ 情况下就是 0 ，所以：

G (F (x)) = G (F_{0} (x)) + (G (F_{0} (x)))^{'} (F (x) - F_{0} (x))

$G(F(x))=G(F_0(x))+(G(F_0(x)))'(F(x)-F_0(x))$

移个项得：

F (x) = F_{0} (x) - \frac{G (F_{0} (x))}{(G (F_{0} (x)))^{'}}

$F(x)=F_0(x)-\dfrac{G(F_0(x))}{(G(F_0(x)))'}$

多项式指数函数

传送门

给定 $A(x)$ ，求 $F(x)$ 使得 $F(x)\equiv e^{A(x)}\ (mod\ x^n)$

考虑使用牛顿迭代求解。

F (x) = e^{A (x)}

$F(x)=e^{A(x)}$

\ln (F (x)) - A (x) = 0

$\ln(F(x))-A(x)=0$

设 $G(F(x))=\ln(F(x))-A(x)$ 。
但是我们还要知道 $(G(F_0(x)))'$ 的值，所以把 $A(x)$ 当做常数，把 $F(x)$ 当成自变量，两边求导。

(G (F (x)))^{'} = \frac{1}{F (x)}

$(G(F(x)))'=\dfrac{1}{F(x)}$

现在可以套牛顿迭代了。

F (x) = F_{0} (x) - \frac{G (F_{0} (x))}{(G (F_{0} (x)))^{'}}

$F(x)=F_0(x)-\dfrac{G(F_0(x))}{(G(F_0(x)))'}$

F (x) = F_{0} (x) - \frac{\ln (F_{0} (x)) - A (x)}{\frac{1}{F_{0} (x)}}

$F(x)=F_0(x)-\dfrac{\ln(F_0(x))-A(x)}{\dfrac{1}{F_0(x)}}$

F (x) = F_{0} (x) (1 - \ln (F_{0} (x)) + A (x))

$F(x)=F_0(x)(1-\ln(F_0(x))+A(x))$

int exp_1[N],exp_2[N];
void Polyexp(int *A,int len){
	#define sav1 exp_1
	#define sav2 exp_2
	int lim=1;while(lim<len) lim<<=1;
	sav1[0]=1;
	for(int l=2;l<=lim;l<<=1){
		copy(sav2,sav1,0,l>>1),Polyln(sav2,l);
		for(int i=0;i<l;++i) sav2[i]=(A[i]-sav2[i]+mo)%mo;
		sav2[0]=(sav2[0]+1)%mo;
		Polymul(sav1,sav2,l,l);
	}
	copy(A,sav1,0,len),clear(sav1,0,lim),clear(sav2,0,lim);
	#undef sav1
	#undef sav2
}

多项式快速幂

传送门

直接套快速幂的话会 log 方的复杂度，但我们有了指数和对数，就可以将其优化。

A (x)^{k} = (e^{\ln A (x)})^{k} = e^{k \ln A (x)}

$A(x)^k=(e^{\ln A(x)})^k=e^{k\ln A(x)}$

void Polypow(int *A,int len,int k){
	Polyln(A,len);
	for(int i=0;i<len;++i) A[i]=A[i]*k%mo;
	Polyexp(A,len);
}

多项式开根

传送门

求 $F(x)$ ，使得 $F(x)\equiv \sqrt{A(x)}\ (mod\ x^n)$

同样考虑牛顿迭代。

F^{2} (x) = A (x)

$F^2(x)=A(x)$

G (F (x)) = F^{2} (x) - A (x) = 0

$G(F(x))=F^2(x)-A(x)=0$

F (x) = F_{0} (x) - \frac{F_{0}^{2} (x) - A (x)}{2 F_{0} (x)} = \frac{F_{0}^{2} (x) + A (x)}{2 F_{0} (x)}

$F(x)=F_0(x)-\dfrac{F_0^2(x)-A(x)}{2F_0(x)}=\dfrac{F_0^2(x)+A(x)}{2F_0(x)}$

int sqrt_1[N],sqrt_2[N];
void Polysqrt(int *A,int len){
	#define sav1 sqrt_1
	#define sav2 sqrt_2
	int lim=1;while(lim<len) lim<<=1;
	sav1[0]=1;
	for(int l=2;l<=lim;l<<=1){
		for(int i=0;i<(l<<1);++i) sav2[i]=(sav1[i]<<1)%mo;
		Polyinv(sav2,l);
		NTT(sav1,l,1),mul(sav1,sav1,0,l),NTT(sav1,l,0);
		for(int i=0;i<l;++i) sav1[i]=(sav1[i]+A[i])%mo;
		Polymul(sav1,sav2,l,l);
	}
	copy(A,sav1,0,len),clear(sav1,0,lim),clear(sav2,0,lim);
	#undef sav1
	#undef sav2
}

最后来个全的！

#include<iostream>
#include<stdio.h>
#define N 
#define int long long
using namespace std;
inline int read(){
	int x=0,f=1;
	char ch=getchar();
	while(ch<'0' || ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	return x*f;
}
#define mo 998244353
#define g 3
#define invg 332748118
int cir[N],rt;
int power(int x,int b){
    int res=1;
    while(b){
        if(b&1) res=res*x%mo;
        b>>=1;
        x=x*x%mo;
    }
    return res;
}
void clear(int *A,int x,int y){
    for(int i=x;i<y;++i) A[i]=0;
}
void copy(int *A,int *B,int x,int y){
    for(int i=x;i<y;++i) A[i]=B[i];
}
void mul(int *A,int *B,int x,int y){
    for(int i=x;i<y;++i) A[i]=A[i]*B[i]%mo;
}
void NTT(int *f,int lenth,int tag){
    if(lenth!=rt)for(int i=0;i<lenth;++i)cir[i]=(cir[i>>1]>>1)|((i&1)?lenth>>1:0);
    rt=lenth;
    for(int i=0;i<lenth;++i) if(i<cir[i]) swap(f[i],f[cir[i]]);
    for(int len=2;len<=lenth;len<<=1){
        int sze=len>>1,buf=power((tag?g:invg),(mo-1)/len);
        for(int i=0;i<lenth;i+=len){
            int bas=1;
            for(int j=i;j<i+sze;j++,bas=bas*buf%mo){
                int tmp=bas*f[j+sze]%mo;
                f[j+sze]=(f[j]-tmp+mo)%mo;
                f[j]=(f[j]+tmp)%mo;
            }
        }
    }
    if(tag==1) return;
    int invn=power(lenth,mo-2);
    for(int i=0;i<lenth;++i) f[i]=f[i]*invn%mo;
}
int tmul[N];
void Polymul(int *A,int *B,int len1,int len2){
    #define sav tmul
    int len=1;
    while(len<(len1<<1)) len<<=1;
    copy(sav,B,0,len);
    NTT(A,len,1),NTT(sav,len,1);
    for(int i=0;i<len;++i) A[i]=A[i]*sav[i]%mo;
    NTT(A,len,0);
    clear(A,len2,len),clear(sav,0,len);
    #undef sav
}
int inv1[N],inv2[N],inv3[N];
void Polyinv(int *A,int len){
    #define sav1 inv1
    #define sav2 inv2
    #define sav3 inv3
    int lim=1;while(lim<len) lim<<=1;
    sav1[0]=power(A[0],mo-2);
    for(int l=2;l<=lim;l<<=1){
        int r=l<<1;copy(sav3,A,0,l);
        for(int i=0;i<(l>>1);++i) sav2[i]=(sav1[i]<<1)%mo;
        NTT(sav1,r,1),mul(sav1,sav1,0,r);
        NTT(sav3,r,1),mul(sav1,sav3,0,r);
        NTT(sav1,r,0),clear(sav1,l,r);
        for(int i=0;i<l;++i) sav1[i]=(sav2[i]-sav1[i]+mo)%mo;
    }
    copy(A,sav1,0,len);
    clear(sav1,0,lim<<1),clear(sav2,0,lim<<1),clear(sav3,0,lim<<1);
    #undef sav1
    #undef sav2
    #undef sav3
}
int invcnt=1,inv[N];
inline void get_inv(int len){
    inv[1]=1; if(invcnt>=len) return;
    for(int i=invcnt+1;i<=len;++i) inv[i]=(mo-mo/i)*inv[mo%i]%mo;
}
inline void Polyder(int *A,int len){
    for(int i=1;i<len;i++) A[i-1]=A[i]*i%mo; A[len-1]=0;
}
inline void Polyint(int *A,int len){
    get_inv(len);
    for(int i=len;i>=1;i--) A[i]=A[i-1]*inv[i]%mo; A[0]=0;
}
int ln1[N];
inline void Polyln(int *A,int len){
    #define sav ln1
    copy(sav,A,0,len),Polyder(A,len),Polyinv(sav,len);
    Polymul(A,sav,len,len),Polyint(A,len-1),clear(sav,0,len);
    #undef sav
}
int exp_1[N],exp_2[N];
void Polyexp(int *A,int len){
    #define sav1 exp_1
    #define sav2 exp_2
    int lim=1;while(lim<len) lim<<=1;
    sav1[0]=1;
    for(int l=2;l<=lim;l<<=1){
        copy(sav2,sav1,0,l>>1),Polyln(sav2,l);
        for(int i=0;i<l;++i) sav2[i]=(A[i]-sav2[i]+mo)%mo;
        sav2[0]=(sav2[0]+1)%mo;
        Polymul(sav1,sav2,l,l);
    }
    copy(A,sav1,0,len),clear(sav1,0,lim),clear(sav2,0,lim);
    #undef sav1
    #undef sav2
}
void Polypow(int *A,int len,int k){
    Polyln(A,len);
    for(int i=0;i<len;++i) A[i]=A[i]*k%mo;
    Polyexp(A,len);
}
int sqrt_1[N],sqrt_2[N];
void Polysqrt(int *A,int len){
    #define sav1 sqrt_1
    #define sav2 sqrt_2
    int lim=1;while(lim<len) lim<<=1;
    sav1[0]=1;
    for(int l=2;l<=lim;l<<=1){
        for(int i=0;i<(l<<1);++i) sav2[i]=(sav1[i]<<1)%mo;
        Polyinv(sav2,l);
        NTT(sav1,l,1),mul(sav1,sav1,0,l),NTT(sav1,l,0);
        for(int i=0;i<l;++i) sav1[i]=(sav1[i]+A[i])%mo;
        Polymul(sav1,sav2,l,l);
    }
    copy(A,sav1,0,len),clear(sav1,0,lim),clear(sav2,0,lim);
    #undef sav1
    #undef sav2
}
signed main(){
	
	return 0;
}

posted @ 2022-06-04 18:12 CHiSwsz 阅读(54) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

相关博文：

· [WC2019] 数树

· 2024 syzx 冬季训练 5 - 数据结构

· 多项式小记

· 多项式全家桶

· 多项式运算

阅读排行：
· 无需6万激活码！GitHub神秘组织3小时极速复刻Manus，手把手教你使用OpenManus搭建本
· C#/.NET/.NET Core优秀项目和框架2025年2月简报
· Manus爆火，是硬核还是营销？
· 一文读懂知识蒸馏
· 终于写完轮子一部分：tcp代理了，记录一下

xzzduang的博客

多项式从入门到入坟

参考

基本概念

卷积

FFT

注意： $n=2^k$ ， $n$ 为项数！

多项式求逆

分治FFT

多项式除法

多项式求导与积分

多项式对数函数

多项式牛顿迭代

多项式指数函数

多项式快速幂

多项式开根

最后来个全的！

公告

搜索

常用链接

随笔档案

阅读排行榜

推荐排行榜

xzzduang的博客

多项式从入门到入坟

参考

基本概念

卷积

FFT

注意：n=2kn=2kn=2^k，nnn 为项数！

多项式求逆

分治FFT

多项式除法

多项式求导与积分

多项式对数函数

多项式牛顿迭代

多项式指数函数

多项式快速幂

多项式开根

最后来个全的！

公告

搜索

常用链接

随笔档案

阅读排行榜

推荐排行榜

注意： $n=2^k$ ， $n$ 为项数！