丑陋的多项式板子们

前言

这篇乘法写的贼烂，真想学懂建议看oi-wiki。
y1s1，但是感觉后面的多项式其他的基础处理还是讲的可以的。
(代码轻压，都是跟 WernerYin 学的)
特别适合蒟蒻yzhx复习

乘法

FFT

主要思路是将“单位根”带入参与运算的两个多项式 (DFT)，然后进行点值直接相乘，再反带回来 (IDFT) 得到最后的答案多项式

所谓“单位根”，是指一些复数平面的单位圆上的n等分点，即：n次单位根
记为：(omega_n)

比如：

而它们有什么性质呢？

CODE

//头文件略
const int _=5e6+5;
struct Comp{
   db x,y;
   Comp(db _x=0, db _y=0): x(_x), y(_y){};
   friend Comp operator + (Comp a, Comp b) {return Comp(a.x+b.x, a.y+b.y); }
   friend Comp operator - (Comp a, Comp b) {return Comp(a.x-b.x, a.y-b.y); }
   friend Comp operator * (Comp a, Comp b) {return Comp(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x); }
}F[_],G[_];
int rev[_];
void change(Comp *f,int l)
{
   for(re int i=1;i<=l;i++) {rev[i]=rev[i>>1]>>1; if(i&1) rev[i]|=l>>1; }
   for(re int i=1;i<=l;i++) if(i<rev[i]) swap(f[i],f[rev[i]]);
}
const db pi=acos(-1);
void FFT(Comp *f,int l,int on)
{
   change(f,l);
   for(re int h=2;h<=l;h<<=1) {
      Comp wn(cos(2*pi/h),sin(2*pi*on/h));
      for(re int j=0;j<l;j+=h) {
         Comp w(1,0);
         for(re int k=j;k<j+h/2;k++) {
            Comp u=f[k], v=f[k+h/2]*w;
            f[k]=u+v; f[k+h/2]=u-v; w=w*wn;
         }
      }
   }
   if(on==-1) for(re int i=0;i<l;i++) f[i].x/=l;
}
int n,m,ans[_];
int main()
{
   n=read(), m=read();
   for(re int i=0;i<=n;i++) F[i].x=read(), F[i].y=0;
   for(re int i=0;i<=m;i++) G[i].x=read(), G[i].y=0;
   n++, m++; int len=1;
   while(len<n*2 || len<m*2) len<<=1;
   FFT(F,len,1), FFT(G,len,1);
   for(re int i=0;i<len;i++) F[i]=F[i]*G[i];
   FFT(F,len,-1);
   for(re int i=0;i<len;i++) ans[i]=int(F[i].x+0.5);
   len=n+m-1;
   for(re int i=0;i<len;i++) cout<<ans[i]<<' ';
}

NTT

把单位根替换为指定模数下的原根即可（不懂原根戳这里）
其他操作一样，直接贴代码了

CODE

in void NTT(ll *f,int l,int o)
{
    change(f,l);
    for(re int h=2;h<=l;h<<=1) {
        ll wn=qpow((o==1 ? g: ig), (mod-1)/h);
        for(re int i=0;i<l;i+=h) {
            ll w=1;
            for(re int k=i;k<i+h/2;k++) {
                ll u=f[k], v=f[k+h/2]*w%mod; w=w*wn%mod;
                f[k]=(u+v)%mod, f[k+h/2]=(u-v+mod)%mod;
            }
        }
    }
    if(o==-1) {
        ll inv=qpow(l,mod-2);
        for(re int i=0;i<l;i++) f[i]=f[i]*inv%mod;
    }
}

求逆

多项式求逆，其实和数字求逆元一样，主要是充当模意义下的除法

（以下n表多项式的项数，不是次数！！！）
即：求出(G(x)),满足：(F(x)ast G(x) equiv1 pmod {x^{n/2}})

令 : (F(x)ast G(x) equiv1 pmod {x^{n/2}})
(F(x)ast H(x) equiv1 pmod x^{n})
显然 (H(x)) 的长度相较于 (G(x)) 更长

所以主要思路是：找出H(x)关于G(x)的表达式，并不断更新

上面两式相减得（后面的mod n不写了）：

[F(x)ast (G(x)-H(x))equiv 0 ]

因为F(x)已知非零：

[(G(x)-H(x))equiv 0 ]

[(G(x)-H(x))^2equiv 0 ]

[H^2(x) equiv 2*(Hast G)(x)-G^2(x) ]

再把 F(x) 乘上去，因为 (F(x)ast H(x)equiv1 pmod x^{n})
(G(x)在%n意义下不为1)
有：

[H(x) equiv 2*G（x)-G^2(x)ast F(x) ]

然后递归处理即可，递归的边界是 n=1时 G(x)为F(0)的逆元

CODE

void SolveInv(ll *f,ll *h,int m)
{
   if(m==1) { h[0]=qpow(f[0],mod-2); return ;}
   int l=1, mid=m+1>>1; while(l<2*m) l<<=1;
   SolveInv(f,h,mid);
   for(re int i=0;i<m;i++) a[i]=f[i]; memset(a+m,0,(l-m)<<3);
   for(re int i=1;i<l;++i){rev[i]=rev[i>>1]>>1; if(i&1) rev[i]|=l>>1;} 
   //以下三次NTT位数相同，可提前预处理rev
   NTT(a,l,1), NTT(h,l,1);
   for(re int i=0;i<l;++i) h[i]=h[i]*(2-a[i]*h[i]%mod+mod)%mod;
   NTT(h,l,0); memset(h+m,0,(l-m)<<3);
}

开根

求G(x),满足：(G^2(x)equiv F(x) pmod{n/2})

思路：

主要也是利用和求逆一样的思想

令 : $$G^2(x)equiv F(x) pmod{(n/2)/2}$$
$$H^2(x)equiv F(x) pmod{n/2}$$
两式相减化简后再将(H^2(x)equiv F(x)) 带入，可得：

[H(x)equiv (1/2)ast(G(x)+F(x)/G(x)) ]

注意到中间有个加号，所以只需要NTT时处理F(x)与G(x)的逆即可

CODE

#define mem(x) memset(x+m,0,(l-m)<<3);
in void Sqrt(ll *f,ll *h,int m)
{
   if(m==1) {h[0]=1; return ;}
   int l=1; while(l<(m<<1)) l<<=1;; Sqrt(f,h,m+1>>1);
   memset(b,0,l<<3);Inv(h,b,m);
   for(re int i=1;i<l;++i) rev[i]=rev[i>>1]>>1|(i&1?l>>1:0);
   for(re int i=0;i<m;i++) c[i]=f[i]; mem(c);
   NTT(c,l), NTT(b,l); for(re int i=0;i<l;++i) b[i]=b[i]*c[i]%mod;
   NTT(b,l,0); for(re int i=0;i<l;++i) h[i]=(h[i]+b[i])%mod*inv2%mod;
   mem(h);
}