```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