丑陋的多项式板子们

前言

这篇乘法写的贼烂,真想学懂建议看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))^2\equiv 0 \]

\[H^2(x) \equiv 2*(H\ast 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);
}
posted @ 2021-03-07 17:19  yzhx  阅读(77)  评论(0编辑  收藏  举报