多项式求逆

多项式求逆

定义

(displaystyle f(x) =sum^{n-1}_{k=0}a_kx^k)(g(x) =sum^{n-1}_{k=0}b_kx^k),使得
(displaystyle f(x)g(x)equiv 1 (mod x^n))

(displaystyle f(x)g(x)) 的前(n)项中只有常数项为(1),其余项均为(0),称为

(f(x))(displaystyle x^n)的逆元,或简称求(f(x))的逆元。

多项式可逆当且仅当其常数项可逆,证明见“方法一”(displaystyle b_0=frac{1}{a_0})

方法1:分治FFT

从常数项出发,比较系数得(displaystyle b_0=frac{1}{a_0})(displaystyle sum^k_{i=0}b_ia_{k-i}=0)那么有

[b_k=sum^{k-1}_{i=0}b_i(-frac{a_{k-i}}{a_0}) ]

从而可以通过分治FFT 计算,时间复杂度(displaystyle Theta(n log^2 n))

方法2:倍增法

考虑倍增,从局部的逆元(低次)出发,设(displaystyle g_0(x) = g(x) mod x^k),那么有

[f(x)g_0(x)-1equiv 0 (mod x^k) ]

从而有

[f(x)g_0(x)-1)^2equiv 0 (mod x^{2k}) ]

可以发现模的次数增加了,此时展开左边可得

[g(x)equiv g_0(x)(2-f(x)g_0(x))(mod x^{2k}) ]

从而可以算出(g(x))的前(2k)项。

通过摊还分析得时间复杂度为(displaystyle Theta(n log n))

实现

多项式求逆的实现不难,两种方法都可以,但第二种更快,也更容易实现,故此处以方法二为例。
我们每次递归算出模(displaystyle x^{frac{n}{2}})的逆元,记为(g_0),再套用(displaystyle g(x)equiv g_0(x)(2-f(x)g_0(x))(mod x^{2k})),用FFT计算即可

void solve(LL *a,LL *b,int p){/*a:seq,b:inv*/
	if(p==1){
		b[0]=fpm(a[0],MOD-2);
		return;
	}
	solve(a,b,(p+1)>>1);
	lim=1; L=0;
	while(lim<(p<<1)){
		lim<<=1;
		L++;
	}
	for(int i=0;i<lim;i++)
		rev[i]=rev[i>>1]>>1|(1&i)<<(L-1);
	for(int i=0;i<p;i++)
		tmp[i]=a[i];
	for(int i=p;i<lim;i++)
		tmp[i]=0;
	NTT(tmp,1); NTT(b,1);
	for(int i=0;i<lim;i++)
		b[i]=(2-b[i]*tmp[i]%MOD+MOD)%MOD*b[i]%MOD;
	NTT(b,-1);
	for(int i=p;i<lim;i++)
		b[i]=0;
}

例题

[P4238 【模板】多项式求逆](%3Ca href="https://www.luogu.org/problem/P4238"%3Ehttps://www.luogu.org/problem/P4238%3C/a%3E)

这是纯粹的模板题,直接给代码

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int INF=1e9+7,MAXN=3e6+10/*Min:2^20+10*/;
const LL G=3,MOD=998244353,INV=332748118;
inline LL fpm(LL base,LL p){
	LL ret=1;
	while(p){
		if(p&1)
			ret=ret*base%MOD;
		base=base*base%MOD;
		p>>=1;
	}
	return ret%MOD;
}
int N,M,lim=1,L,rev[MAXN];
inline void NTT(LL *a,int type){
	for(int i=0;i<lim;i++)
		if(i<rev[i])
			swap(a[i],a[rev[i]]);
	for(int mid=1;mid<lim;mid<<=1){
		int len=mid<<1/*n*/;
		LL Wn=fpm(G,(MOD-1)/len);
		for(int j=0;j<lim;j+=len){
			LL w=1,t1,t2;
			for(int k=0;k<mid;k++){
				t1=a[j+k],t2=w*a[j+k+mid]%MOD;
				a[j+k]=(t1+t2)%MOD;
				a[j+k+mid]=(t1-t2+MOD)%MOD;
				w=w*Wn%MOD;
			}
		}
	}
	if(type==-1){
		LL lim_inv=fpm(lim,MOD-2);
		reverse(a+1,a+lim);
		for(int i=0;i<lim;i++)
			a[i]=a[i]*lim_inv%MOD;
	}
}
LL tmp[MAXN];
void solve(LL *a,LL *b,int p){/*a:seq,b:inv*/
	if(p==1){
		b[0]=fpm(a[0],MOD-2);
		return;
	}
	solve(a,b,(p+1)>>1);
	lim=1; L=0;
	while(lim<(p<<1)){
		lim<<=1;
		L++;
	}
	for(int i=0;i<lim;i++)
		rev[i]=rev[i>>1]>>1|(1&i)<<(L-1);
	for(int i=0;i<p;i++)
		tmp[i]=a[i];
	for(int i=p;i<lim;i++)
		tmp[i]=0;
	NTT(tmp,1); NTT(b,1);
	for(int i=0;i<lim;i++)
		b[i]=(2-b[i]*tmp[i]%MOD+MOD)%MOD*b[i]%MOD;
	NTT(b,-1);
	for(int i=p;i<lim;i++)
		b[i]=0;
}
LL a[MAXN],b[MAXN];
int main(){
	scanf("%d",&N);
	for(int i=0;i<N;i++){
		scanf("%lld",a+i);
		a[i]=(a[i]+MOD)%MOD;
	}
	solve(a,b,N);
	for(int i=0;i<N;i++)
		printf("%lld ",b[i]);
	return 0;
}
[P4239 【模板】多项式求逆(加强版)](%3Ca href="https://www.luogu.org/problem/P4239"%3Ehttps://www.luogu.org/problem/P4239%3C/a%3E)

这道题是在任意模数下的多项式求逆,需要将上一个板子的NTT改成MTT

原文地址:https://www.cnblogs.com/guoshaoyang/p/11296025.html