Luogu P5205 【模板】多项式开根

前言

自己做出一道多项式题，虽调试艰辛，但是AC后自豪感爆棚。

思路：

题目

$B^2(x)=A(x)(mod{x}^n),$已知$A(x)$,求$B(x)$

首先，给出一个牛顿迭代公式$x=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$($x$为函数零点，$x_0$为已知点)

设$f(B(x))=B^2(x)-A(x),$

则$f^{\prime }(B(x))=2\ast B(x)$

由于上面的定义，所以我们现在求的就是$f(x)$的零点。

设已知$f(B^{\ast }(x))=(B^{\ast }{)}^2(x)-A(x)$,$(B^{\ast })$为$mod{x}^{n/2}$时的多项式解
那么根据上面给出的牛迭公式可得:$B(x)=(B^{\ast })(x)-\frac{(B^{\ast }{)}^2(x)-A(x)}{2\ast (B^{\ast })(x)}=\frac{(B^{\ast })(x)}{2}+\frac{A(x)}{2(B^{\ast })(x)}$

综上：$NTT$求逆元$\ast \log (n)$ 次迭代$B(x)$即可，复杂度:$O(n\ast \log n)$

代码：

#include<cstdio
#include<cctype>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=1<<18|10,mod=998244353;
const int g=3,G=332748118,inv2=499122177;//原根3,3、2的逆元 
template<class o>void qr(o&x) {
	char c=getchar();int f=1;x=0;
	while(!isdigit(c)){if(c=='-')f=-1;c=getchar();}
	while(isdigit(c))x=x*10+c-'0',c=getchar();
	x*=f;
}
template<class o>void qw(o x) {
	if(x/10)qw(x/10);
	putchar(x%10+'0');
}
template<class o> void pr1(o x) {
	if(x<0)x=-x,putchar('-');
	qw(x); putchar(' ');
}

ll power(ll a,ll b=mod-2) {
	ll c=1;
	while(b) {
		if(b&1) c=c*a%mod;
		a=a*a%mod; b=b>>1;
	}
	return c;
}
int m,tr[N];
ll a[N],f[N],A[N],B[N],tmp[N];
//B(x)=(B*)(x)/2 + A(x)/(2*(B*(x))) 
//a为原题中的A,f用于存2(B*),B用于解逆元,tmp为临时数组 

void NTT(ll *f,int e,int n) {
	for(int i=1;i<n;i++)
		if(i<tr[i]) swap(f[i],f[tr[i]]);
	for(int p=2,len;p<=n;p<<=1) {
		len=p>>1;
		ll tg=power(e,(mod-1)/p),buf,t;
		for(int k=0;k<n;k+=p) {
			buf=1;
			for(int l=k;l<k+len;l++) {
				t=f[l+len]*buf%mod;
				f[l+len]=(f[l]-t+mod)%mod;
				f[l]=(f[l]+t)%mod;
				buf=buf*tg%mod;
			}
		}
	}
}

void mult(ll *a,ll *b,int n) {//两个多项式大小均为n——2的次幂
	for(int i=0;i<n;i++) tmp[i]=b[i];
	for(int i=n;i<2*n;i++) tmp[i]=0;
	n<<=1;
	for(int i=1;i<n;i++) tr[i]=(tr[i>>1]>>1)|(i&1?n>>1:0);
	NTT(a,g,n); NTT(tmp,g,n); 	 for(int i=0;i<n;i++) a[i]=a[i]*tmp[i]%mod;
	NTT(a,G,n); ll Inv=power(n); for(int i=0;i<n;i++) a[i]=a[i]*Inv %mod;
}

void inv(int n) {//求f的逆元——最后存到B中 
	B[0]=inv2;
	for(int i=1;i<n;i++) B[i]=0;
	for(int p=2,len;p<=n;p<<=1) {
		len=p>>1;
		for(int i=0;i<len;i++)
			A[i]=2*B[i]%mod;
		mult(B,B,len); mult(B,f,p);
		for(int i=0;i<p;i++) B[i]=(A[i]-B[i]+mod)%mod;
		for(int i=p;i<2*p;i++) B[i]=0;
	}
	for(int i=0;i<n;i++) A[i]=0;
}

void Sqrt() {
	f[0]=2;//(B*)(x)=1
	int n; for(n=1;n<=m;n<<=1);
	for(int p=2;p<=n;p<<=1) {
		inv(p);
		mult(B,a,p);
		for(int i=0;i<p;i++)
			f[i]=(f[i]*inv2+2LL*B[i])%mod;
	}
	for(int i=0;i<=m;i++) pr1(f[i]*inv2%mod);
	puts("");
}

int main() {
	qr(m); m--;
	for(int i=0,x;i<=m;i++) 
		qr(x),a[i]=x;
	Sqrt(); return 0;
}