[学习笔记]多项式开根

思路：

推柿子跟求逆一样，分治（倍增）的思想：不想写了
推出\((F-G)^2 \equiv0\pmod{x^n}\)
所以\(G=\dfrac{F^2+A}{2F}\)
边界处要用二次剩余的Cipolla算法。
因此只要会多项式求逆、乘法，二次剩余即可。

code

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e6+5;
const ll mod=998244353;
const int inf=1e9;
ll I_2,inv2;
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;}
struct cop {
	ll rl,img;
	inline friend cop operator *(cop x,cop y) {
		cop res;
		res.rl=(x.rl*y.rl+I_2*x.img%mod*y.img)%mod;
		res.img=(x.img*y.rl+x.rl*y.img)%mod;
		return res;
	}
};

ll cop_ksm(cop a,ll b) {
	cop res=(cop){1,0};
	for(;b;b>>=1,a=a*a) if(b&1)res=res*a;
	return (res.rl+mod)%mod;
}

//二次剩余 
ll _SQ(ll a) {
	if(!a){return 0;}
	ll mod_2=(mod-1)>>1;
	ll b;
	while(1) {
		b=rand();
		I_2=(b*b-a)%mod;
		if(ksm(I_2,mod_2)!=1)break;
	}
	ll x1=cop_ksm((cop){b,1},(mod+1)>>1);
	ll x2=mod-x1;
	return min(x1,x2);
}

ll a[N],inv3,gen[2][N],f[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];
 		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));}
}

void poly_inv(int dep,ll *f,ll *b) {	// b = f mod(x^dep) 
	if(dep==1) {
		b[0]=ksm(f[0],mod-2);
		return;
	}
	poly_inv((dep+1)>>1,f,b);
	gt_up(dep<<1);
	for(int i=0;i<dep;i++)a[i]=f[i];for(int i=dep;i<up;i++)a[i]=0;
//	printf("a:\n");for(int i=0;i<up;i++)printf("%lld ",(a[i]%mod+mod)%mod);puts("");
//	printf("b:\n");for(int i=0;i<up;i++)printf("%lld ",(b[i]%mod+mod)%mod);puts("");
	NTT(b,0);NTT(a,0);
	for(int i=0;i<up;i++) {a[i]=(2-a[i]*b[i]%mod)*b[i]%mod;}
	NTT(a,1);
	ll i_up=ksm(up,mod-2);
//	printf("!%d %lld\n",dep,i_up);
	for(int i=0;i<dep;i++) {b[i]=a[i]*i_up%mod;} for(int i=dep;i<up;i++)b[i]=0;
//	printf("aa:\n");for(int i=0;i<up;i++)printf("%lld ",(a[i]%mod+mod)%mod);puts("");
//	printf("EB:\n");for(int i=0;i<up;i++)printf("%lld ",(b[i]%mod+mod)%mod);puts("");
}

ll c[N];
void poly_sqrt(int dep,ll *f,ll *b) {
	if(dep==1) {b[0]=_SQ(f[0]);return;}
	poly_sqrt((dep+1)>>1,f,b);
	gt_up(dep<<1);int tup=up;
	for(int i=0;i<up;i++) {c[i]=0;}
	poly_inv(dep,b,c);
	up=tup;
	for(int i=0;i<up;i++) {a[i]=(i<dep)?f[i]:0;}
	NTT(c,0);NTT(a,0);
	for(int i=0;i<up;i++) a[i]=a[i]*c[i]%mod;
	NTT(a,1);
	for(int i=0,iup=ksm(up,mod-2);i<dep;i++)b[i]=(b[i]+iup*a[i])%mod*inv2%mod;
	for(int i=dep;i<up;i++)b[i]=0;
}

void gt_gen(int len) {
	inv3=ksm(3,mod-2);
	gt_up(len);
	gen[0][up]=ksm(3,(mod-1)/up);gen[1][up]=ksm(inv3,(mod-1)/up);
	for(int i=up;i;i>>=1) {
		gen[0][i>>1]=gen[0][i]*gen[0][i]%mod;
		gen[1][i>>1]=gen[1][i]*gen[1][i]%mod;
	}
}
ll ans[N];
int main() {
	srand(time(0));
	inv2=ksm(2,mod-2);
	scanf("%d",&n);
	gt_gen(n<<1);
	for(int i=0;i<n;i++) {scanf("%lld",&f[i]);}
	poly_sqrt(n,f,ans);
	for(int i=0;i<n;i++)printf("%lld ",(ans[i]+mod)%mod);
	return 0;
}

emmm.这个取掉了一个模调了半天。

标签：cop,多项式,ll,笔记,res,rl,开根,gen,mod
来源： https://www.cnblogs.com/bestime/p/16425065.html