[学习笔记]多项式求逆

题意：如其名……

思路：

多项式的关键在于：用模的次数降次。
它的复杂度跟模数的次幂有关。
所以可以考虑对模数分治。参考
若多项式 $F$ 只有一项，直接求常数项的逆元（这也是判别多想式是否存在逆元的条件）。
设已知 $F(x)H(x)\equiv1\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}$
则 $F(x)G(x)\equiv1\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}$
所以 $G(x)-H(x)\equiv0\pmod{x^{\left\lceil\frac{n}{2}\right\rceil}}$
两边平方得： $G(x)^2+H(x)^2-2G(x)H(x)\equiv0\pmod{x^{n}}$
因为把平方看成两个多项式相乘 $\sum\limits_{j=0}^i(G(x)-H(x))^j(G(x)-H(x))^{i-j}$
而 $G(x)-H(x)$ 模 $x^{\left\lceil\frac{n}{2}\right\rceil}$ 为 $0$ ，即它的最小系数非负次数为 $\left\lceil\frac{n}{2}\right\rceil$ 。
当然分成的 $j$ 与 $i-j$ 必有一个小于 $\left\lceil\frac{n}{2}\right\rceil$ 。
因此两边同乘 $F(x)$ ，根据 $F(x)G(x)\equiv1\pmod{x^n}$ 得：
$G(x)\equiv2H(x)-F(x)H(x)^2\pmod{x^n}$

code:

ouo

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e6+5;
const ll mod=998244353;
ll ksm(ll a,ll b) {ll res=1;for(;b;b>>=1,a=a*a%mod)if(b&1)res=res*a%mod;return res;}
//ll ksm(ll x,ll y) {ll res=1;for(;y;y>>=1,x=x*x%mod)if(y&1)res=res*x%mod;return res;}
ll f[N],a[N],c[N],inv3,gen[2][N];
int up,l,rev[N],n;
void NTT(ll *a,int op) {
	for(int i=0;i<up;i++) {
		if(rev[i]>i)swap(a[i],a[rev[i]]);
	}
	for(int mid=1;mid<up;mid<<=1) {
		int len=mid<<1;ll w1=gen[op][len];
//		printf("!(%d,%d) %lld (%lld,%lld)\n",op,len,(w1+mod)%mod,ksm(3,(mod-1)/len),ksm(inv3,(mod-1)/len));
 		for(int l=0;l<up;l+=len) {
 			ll W=1;
			for(int i=0;i<mid;i++,W=W*w1%mod) {
				int p=l+i,q=p+mid;
				ll x=a[p],y=W*a[q];
				a[p]=(x+y)%mod;a[q]=(x-y)%mod;
			}
		}
	}
}

void gt_up(int len) {
	up=1,l=0;
	while(up<len) {up<<=1,l++;}
	for(int i=1;i<up;i++) {rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));}
}
ll ans[N];
void solve(int dep,ll *b) {	//b为该层(%x^{dep})得到的结果 
	if(dep==1) {b[0]=ksm(f[0],mod-2);return;}
	solve((dep+1)>>1,b);
	gt_up(dep<<1);
//	printf("dep=%d up=%d~~\n",dep,up);
	for(int i=0;i<dep;i++)a[i]=f[i];for(int i=dep;i<up;i++)a[i]=0;
	NTT(b,0);NTT(a,0);
//	printf("A:\n");for(int i=0;i<up;i++)printf("%d ",(a[i]%mod+mod)%mod);puts("");
//	printf("B:\n");for(int i=0;i<up;i++)printf("%d ",(b[i]%mod+mod)%mod);puts("");
	for(int i=0;i<up;i++) {c[i]=(2-a[i]*b[i]%mod)*b[i]%mod;}
	NTT(c,1);
	ll i_up=ksm(up,mod-2);
	for(int i=0;i<dep;i++) {b[i]=c[i]*i_up%mod;} for(int i=dep;i<up;i++)b[i]=0;
//	printf("S:\n");for(int i=0;i<up;i++)printf("%d ",(b[i]%mod+mod)%mod);puts("");
}
void gt_gen() {
	inv3=ksm(3,mod-2);
	gt_up(n<<1);
//	printf("!%d\n",up);
	gen[0][up]=ksm(3,(mod-1)/up);gen[1][up]=ksm(inv3,(mod-1)/up);
	for(int i=up;i;i>>=1) {
//		printf("i=%d %lld(%lld),%lld(%lld)\n",i,gen[0][i],ksm(3,(mod-1)/i),gen[1][i],ksm(inv3,(mod-1)/i));
		gen[0][i>>1]=gen[0][i]*gen[0][i]%mod;
		gen[1][i>>1]=gen[1][i]*gen[1][i]%mod;
	}
}
int main() {
	scanf("%d",&n);
	gt_gen();
	for(int i=0;i<n;i++) {scanf("%lld",&f[i]);}
	solve(n,ans);
	for(int i=0;i<n;i++) {printf("%lld ",(ans[i]+mod)%mod);}
	return 0;
}

ps.还有柯爱的
code:板套板，不会mtt，注意拆的系数最好 $\sqrt{p}$ ，我随便用的 $114514$ 它精度就G了。

点击查看代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef long double db;
const db pi=acos(-1.0);
const int N=1e6+5;
const ll mod=1e9+7;

ll ksm(ll a,ll b) {ll res=1;for(;b;b>>=1,a=a*a%mod)if(b&1)res=res*a%mod;return res;}
ll f[N],a[N],c[N],bas=sqrt(mod);
int up,l,rev[N],n;

struct cop {
	db img,rl;
	friend inline cop operator*(cop u,cop v) {return (cop){u.rl*v.img+u.img*v.rl,u.rl*v.rl-u.img*v.img};} 
	friend inline cop operator+(cop u,cop v) {return (cop){u.img+v.img,u.rl+v.rl};}
	friend inline cop operator-(cop u,cop v) {return (cop){u.img-v.img,u.rl-v.rl};}
}A[N],B[N],C[N],D[N],AC[N],BD[N],ABC[N];

void FFT(cop *a,int op) {
	for(int i=0;i<up;i++)if(rev[i]>i)swap(a[i],a[rev[i]]);
	for(int mid=1;mid<up;mid<<=1) {
		cop w1=(cop){op*sin(pi/mid),cos(pi/mid)};
		for(int len=mid<<1,l=0;l<up;l+=len) {
			cop W=(cop){0,1};
			for(int i=0;i<mid;i++,W=W*w1) {
				int p=l+i,q=p+mid;
				cop x=a[p],y=W*a[q];
				a[p]=x+y;a[q]=x-y;
			}
		}
	}
}

void gt_up(int len) {
	up=1,l=0;
	while(up<len) {up<<=1,l++;}
	for(int i=1;i<up;i++) {rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));}
}
ll ans[N];

void Mul(ll *f,ll *g) {
	for(int i=0;i<up;i++) {
		A[i]=(cop){0,f[i]/bas};B[i]=(cop){0,f[i]%bas};
		C[i]=(cop){0,g[i]/bas};D[i]=(cop){0,g[i]%bas};
	}
	FFT(A,-1),FFT(B,-1),FFT(C,-1),FFT(D,-1);
	for(int i=0;i<up;i++) {
		AC[i]=A[i]*C[i],ABC[i]=A[i]*D[i]+B[i]*C[i],BD[i]=B[i]*D[i];
	}
	FFT(AC,1),FFT(ABC,1),FFT(BD,1);
	for(int i=0;i<up;i++) {
		ll ac=(ll)round(AC[i].rl/up)%mod,abc=(ll)round(ABC[i].rl/up)%mod,bd=(ll)round(BD[i].rl/up)%mod;
		f[i]=(ac*bas%mod*bas+abc*bas+bd)%mod;
	}
}
void pw_Mul(ll *f) {
	for(int i=0;i<up;i++) {
		A[i]=(cop){0,f[i]/bas};B[i]=(cop){0,f[i]%bas};
	}
	FFT(A,-1),FFT(B,-1);
	for(int i=0;i<up;i++) {
		AC[i]=A[i]*A[i];ABC[i]=A[i]*B[i];BD[i]=B[i]*B[i];
	}
	FFT(AC,1),FFT(ABC,1),FFT(BD,1);
	for(int i=0;i<up;i++) {
		ll ac=(ll)round(AC[i].rl/up)%mod,abc=(ll)round(ABC[i].rl/up)%mod,bd=(ll)round(BD[i].rl/up)%mod;
		f[i]=(ac*bas%mod*bas+2*abc*bas%mod+bd)%mod;
//		printf("%lld ",(f[i]%mod+mod)%mod);
	}
//	puts(""); 
}

void solve(int dep,ll *b) {	 //b为该层(%x^{dep})得到的结果 
	if(dep==1) {b[0]=ksm(f[0],mod-2);return;}
	solve((dep+1)>>1,b);
	gt_up(dep<<1);
	for(int i=0;i<up;i++)c[i]=b[i],a[i]=(i<dep)?f[i]:0;
	pw_Mul(c);
	for(int i=dep;i<up;i++)c[i]=0;
	Mul(a,c);
	for(int i=0;i<dep;i++) {
		b[i]=(2*b[i]-a[i])%mod;
		if(b[i]<mod)b[i]+=mod;
	}
//	printf("S:\n");
//	for(int i=0;i<up;i++)printf("%lld ",(b[i]+mod)%mod);puts("");
}

int main() {
//	freopen("P4239_16.in","r",stdin);
//	freopen("ans.out","w",stdout);
	scanf("%d",&n);
	for(int i=0;i<n;i++) {scanf("%lld",&f[i]);}
	solve(n,ans);
	for(int i=0;i<n;i++) {printf("%lld ",(ans[i]+mod)%mod);}
	return 0;
}

注意
%.0lf以及浮点型转整形，是保留整数;除法不是下取整是向〇取整;四舍五入直接用round，负数就不能+0.5要-0.5转整型。

发表于 2022-06-28 20:15 Crazy!!! 阅读(145) 评论(0) 编辑收藏举报

搜索

我的标签

随笔分类

随笔档案

评论排行榜

题意：如其名……

思路：

code: