推柿子跟求逆一样,分治(倍增)的思想:不想写了
推出\((F-G)^2 \equiv0\pmod{x^n}\)
所以\(G=\dfrac{F^2+A}{2F}\)
边界处要用二次剩余的Cipolla算法。
因此只要会多项式求逆、乘法,二次剩余即可。
#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.这个取掉了一个模调了半天。