萌新刚学多项式,有一些地方很不理解

P5205 【模板】多项式开根

```cpp #include <iostream> #include <cstdio> #include <cstring> using namespace std; const long long MS=5000005; const int mod[3]={998244353,167772161,104857601}; inline int qpow(int x,int y,int p) { int res=1; for(;y>0;y>>=1) res=((y&1)?1ll*res*x%p:res),x=1ll*x*x%p; return res; } inline int getlen(int n) { int res=1; while(res<n) res<<=1; return res; } inline void NTT(int n,int a[],int tpe,int p) { for(int i=0,j=0;i<n;i++) { a[i]%=p; if(j<i) swap(a[i],a[j]); for(int l=n>>1;(j^=l)<l;l>>=1); } int g=tpe==1?3:qpow(3,p-2,p); for(int mid=1;mid<n;mid<<=1) { int len=mid<<1; int Wn=qpow(g,(p-1)/len,p); for(int l=0;l<n-len+1;l+=len) { for(int k=0,Wk=1;k<mid;k++,Wk=1ll*Wk*Wn%p) { int x=a[l+k],y=1ll*Wk*a[l+mid+k]%p; a[l+k]=(x+y)%p,a[l+mid+k]=(x-y+p)%p; } } } } inline void DFT(int n,int a[],int p){NTT(n,a,1,p);} inline void IDFT(int n,int a[],int p) { NTT(n,a,-1,p); int inv=qpow(n,p-2,p); for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%p; } int p_tmp[MS],p_tmp2[MS],p_tmp3[MS],p_tmp4[MS]; inline void PINV(int n,int a[],int res[],int p) // used tmp1~2 { int len; p_tmp[0]=res[0]=qpow(a[0],p-2,p); for(len=1;len<n<<1;len<<=1) // 为什么偏要求两倍? { int nxt=len<<1; for(int i=0;i<len;i++) p_tmp[i]=res[i],p_tmp2[i]=a[i]; DFT(nxt,p_tmp,p),DFT(nxt,p_tmp2,p); for(int i=0;i<nxt;i++) res[i]=((2ll-1ll*p_tmp[i]*p_tmp2[i]%p)*p_tmp[i]%p+p)%p; IDFT(nxt,res,p); for(int i=len;i<nxt;i++) res[i]=0; // 这个循环体不是根据已经推出的前 len 位,推前 nxt 位吗?那为啥要清空数组呢? } for(int i=0;i<len;i++) p_tmp[i]=p_tmp2[i]=0; for(int i=n;i<len;i++) res[i]=0; } inline void PSQRT(int n,int a[],int res[],int p) // used tmp1~4 { p_tmp3[0]=res[0]=1; int len,inv2=qpow(2,p-2,p); for(len=1;len<n<<1;len<<=1) // 为什么偏要求两倍? { int nxt=len<<1; for(int i=0;i<len;i++) p_tmp3[i]=a[i]; PINV(len,res,p_tmp4,p); DFT(nxt,p_tmp3,p),DFT(nxt,p_tmp4,p); for(int i=0;i<nxt;i++) p_tmp3[i]=1ll*p_tmp3[i]*p_tmp4[i]%p; IDFT(nxt,p_tmp3,p); for(int i=0;i<nxt;i++) res[i]=1ll*(p_tmp3[i]+res[i])*inv2%p; for(int i=len;i<nxt;i++) res[i]=0; // 这个循环体不是根据已经推出的前 len 位,推前 nxt 位吗?那为啥要清空数组呢?而且为什么这里清不清都没事,上面哪一处不清就会 WA 呢? } for(int i=0;i<len;i++) p_tmp3[i]=p_tmp4[i]=0; for(int i=n;i<len;i++) res[i]=0; } int n,a[MS],b[MS]; int main() { scanf("%d",&n); for(int i=0;i<n;i++) { scanf("%d",&a[i]); } PSQRT(n,a,b,mod[0]); for(int i=0;i<n;i++) { printf("%d ",b[i]); } printf("\n"); return 0; } ```
by lovely_ckj @ 2022-03-20 10:25:34


在一些迷惑操作下,我理解了![](//图.tk/0) 首先,循环的边界有误。循环体是知道前 $\dfrac{len}{2}$ 位,要求前 $len$ 位。 然后求出前 $2n$ 位是因为有“虚假信息”。
by lovely_ckj @ 2022-03-20 10:42:29


```cpp #include <iostream> #include <cstdio> #include <cstring> using namespace std; const long long MS=5000005; const int mod[3]={998244353,167772161,104857601}; inline int qpow(int x,int y,int p) { int res=1; for(;y>0;y>>=1) res=((y&1)?1ll*res*x%p:res),x=1ll*x*x%p; return res; } inline int getlen(int n) { int res=1; while(res<n) res<<=1; return res; } inline void NTT(int n,int a[],int tpe,int p) { for(int i=0,j=0;i<n;i++) { a[i]%=p; if(j<i) swap(a[i],a[j]); for(int l=n>>1;(j^=l)<l;l>>=1); } int g=tpe==1?3:qpow(3,p-2,p); for(int mid=1;mid<n;mid<<=1) { int len=mid<<1; int Wn=qpow(g,(p-1)/len,p); for(int l=0;l<n-len+1;l+=len) { for(int k=0,Wk=1;k<mid;k++,Wk=1ll*Wk*Wn%p) { int x=a[l+k],y=1ll*Wk*a[l+mid+k]%p; a[l+k]=(x+y)%p,a[l+mid+k]=(x-y+p)%p; } } } } inline void DFT(int n,int a[],int p){NTT(n,a,1,p);} inline void IDFT(int n,int a[],int p) { NTT(n,a,-1,p); int inv=qpow(n,p-2,p); for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%p; } int p_tmp[MS],p_tmp2[MS],p_tmp3[MS],p_tmp4[MS]; inline void PINV(int n,int a[],int res[],int p) // used tmp1~2 { int len; p_tmp[0]=res[0]=qpow(a[0],p-2,p); for(len=2;len<n<<1;len<<=1) // 知道前 len/2 位,要求前 len 位 { int lim=len<<1; for(int i=0;i<len;i++) p_tmp[i]=res[i],p_tmp2[i]=a[i]; DFT(lim,p_tmp,p),DFT(lim,p_tmp2,p); for(int i=0;i<lim;i++) res[i]=((2ll-1ll*p_tmp[i]*p_tmp2[i]%p)*p_tmp[i]%p+p)%p; IDFT(lim,res,p); for(int i=len;i<lim;i++) res[i]=0; } for(int i=0;i<len;i++) p_tmp[i]=p_tmp2[i]=0; for(int i=n;i<len;i++) res[i]=0; } inline void PSQRT(int n,int a[],int res[],int p) // used tmp1~4 { p_tmp3[0]=res[0]=1; int len,inv2=qpow(2,p-2,p); for(len=2;len<n<<1;len<<=1) // 知道前 len/2 位,要求前 len 位 { int lim=len<<1; for(int i=0;i<len;i++) p_tmp3[i]=a[i]; PINV(len,res,p_tmp4,p); DFT(lim,p_tmp3,p),DFT(lim,p_tmp4,p); for(int i=0;i<lim;i++) p_tmp3[i]=1ll*p_tmp3[i]*p_tmp4[i]%p; IDFT(lim,p_tmp3,p); for(int i=0;i<len;i++) res[i]=1ll*(p_tmp3[i]+res[i])*inv2%p; } for(int i=0;i<len;i++) p_tmp3[i]=p_tmp4[i]=0; for(int i=n;i<len;i++) res[i]=0; } int n,a[MS],b[MS]; int main() { scanf("%d",&n); for(int i=0;i<n;i++) { scanf("%d",&a[i]); } PSQRT(n,a,b,mod[0]); for(int i=0;i<n;i++) { printf("%d ",b[i]); } printf("\n"); return 0; } ```
by lovely_ckj @ 2022-03-20 10:42:49


@[lovely_ckj](/user/251130) 解答一下 `inv` 的问题 两倍长:无论是递归还是迭代,最后一步需要求 $F(x)G_0(x)^2$,其中 $F(x)$ 是长为 $n$ 的,$G_0(x)$ 是长为 $n / 2$ 的,所以总长可以是 $2n$。 清空数组:你当前求的是模 $x^{len}$ 意义下的 $G_0(x)$,所以要清空。
by int1 @ 2022-03-20 10:47:55


`sqrt` 没写过,但是这些细节上的了解是类似的
by int1 @ 2022-03-20 10:48:59


inv 可以做到空间 $n$
by MoYuFang @ 2022-03-20 10:58:55


多项式乘法封装成一个函数是不是更好?
by MoYuFang @ 2022-03-20 11:00:04


若求 $F(x)=\frac{1}{G(x)}$,迭代公式就是 $F_*(x)=2F(x)-G(x)F^2(x)$。注意到 $G(x)F(x)\% x^{n/2}=1$,实际上我们只用求出 $G(x)F(x)$ 的后 $n/2$ 位,因为 $G(x)$ 有 $n$ 位 $F(x)$ 有 $n/2$ 位,可以用长度为 $n$ 的 $\text{NTT}$,这样循环卷积的溢出只会使前 $n/2$ 位错误,但后 $n/2$ 不受影响。取出 $G(x)F(x)$ 的后 $n/2$ 位后再跟 $F(x)$ 乘一次,$\text{NTT}$ 长度同样为 $n$,这样就能得到 $G(x)F^2(x)$ 的后 $n/2$ 位了。
by MoYuFang @ 2022-03-20 11:12:57


```cpp #define _for(i, a, b) for(int i = (a); i < (b); ++i) void ntt(int *a, int n); void intt(int *a, int n); int lb(int x){ int lbn = -1; while(x>0) ++lbn, x>>=1; return lbn; } void inv(int *a, int n){ //2*6ntt(n) = 12ntt(n) static int b[maxn], c[maxn]; re int nn = 1<<lb(n-1)+1;//将 n 扩展成 2 的幂次 _for(i, n, nn) a[i] = 0; b[0] = get_inv(a[0]); for(re int n = 2; n <= nn; n <<= 1){ _for(i, n/2, n) b[i] = 0; _for(i, 0, n) c[i] = a[i]; ntt(b, n); ntt(c, n); _for(i, 0, n) c[i] = (ll)b[i]*c[i]%mod; intt(c, n); _for(i, n/2, n) c[i-n/2] = c[i], c[i] = 0; ntt(c, n); _for(i, 0, n) c[i] = (ll)b[i]*c[i]%mod; intt(b, n); intt(c, n); _for(i, n/2, n) b[i] = (mod-c[i-n/2])%mod; } _for(i, 0, nn) a[i] = b[i]; } ```
by MoYuFang @ 2022-03-20 11:18:34


@[int1](/user/119009) @[MoYuFang](/user/474113) thx!
by lovely_ckj @ 2022-03-20 11:29:54


| 下一页