多项式Ⅰ

多项式Ⅰ

update20181229

这篇屑文章咕咕了，我更改了一些代码写法（非递归->递归），从牛顿迭代的思想理解了一下，这里的一些思想可能比较繁琐，可以看看多项式Ⅱ

前言：本人只会瞎背板子，根本没怎么理解，而且只学了一两个最简单的（学一点更新一点

前置定义：多项式的度表示多项式的最高次数。

多项式求逆

• 定义

  [A(x)A^{-1}(x)equiv 1pmod {x^n} ]

  (n)是(A(x))的度，这里是两个多项式进行卷积之后对一个多项式((x^n))取模

  • 为什么要取模？

    如果不取模，除了只有常数项的多项式以外，其他多项式的逆元都是无穷级数。

  • 一个取模的感性理解

    对(x^n)取模相当于把次数大于(n)的项都扔掉

• 倍增

  若求出了(A(x)A'(x)equiv 1pmod {x^{lceilfrac{n}{2} ceil}})的(A'(x))

  [A(x)A^{-1}(x)equiv 1pmod {x^n} ]
  [A(x)A^{-1}(x)equiv 1 pmod {x^{lceilfrac{n}{2} ceil}} ]
  [A(x)(A^{-1}(x)-A'(x))equiv 0 pmod {x^{lceilfrac{n}{2} ceil}} ]
  [A^{-1}(x)-A'(x)equiv 0 pmod {x^{lceilfrac{n}{2} ceil}} ]
  [(A^{-1}(x)-A'(x))^2equiv 0 pmod {x^n} ]
  [(A^{-1}(x))^2-2A^{-1}(x)A'(x)+A'^2(x)equiv 0 pmod {x^n} ]
  [A^{-1}(x)-2A'(x)+A(x)A'^2(x)equiv 0 pmod {x^n} ]
  [A^{-1}(x)equiv 2A'(x)-A(x)A'^2(x)pmod {x^n} ]

  上面过程以已意会为主，谢谢。正经的觉得疑惑还是从定义去理解+讨论。

• 细节

  注意每次乘法时取长度取多不取少，比如那个三个多项式乘开开四倍长度。

• Code:

  #include <cstdio>
  #include <algorithm>
  #define ll long long
  const int N=(1<<18)+10;
  const int mod=998244353,G=3,Gi=332748118;
  #define add(x,y) ((x+y)%mod)
  #define Mul(x,y) (1ll*(x)*(y)%mod)
  int qp(int d,int k){int f=1;while(k){if(k&1)f=Mul(f,d);d=Mul(d,d);k>>=1;}return f;}
  int A[N],B[N],f[N],b[2][N],turn[N];
  void NTT(int *a,int len,int typ)
  {
      for(int i=0;i<len;i++)
          if(i<turn[i])
              std::swap(a[i],a[turn[i]]);
      for(int le=1;le<len;le<<=1)
      {
          int wn=qp(typ?G:Gi,(mod-1)/(le<<1));
          for(int p=0;p<len;p+=le<<1)
          {
              int w=1;
              for(int i=p;i<p+le;i++,w=Mul(w,wn))
              {
                  int tx=a[i],ty=Mul(w,a[i+le]);
                  a[i]=add(tx,ty);
                  a[i+le]=add(tx,mod-ty);
              }
          }
      }
      if(!typ)
      {
          int inv=qp(len,mod-2);
          for(int i=0;i<len;i++)
              a[i]=Mul(a[i],inv);
      }
  }
  void mul(int *a,int *b,int len)
  {
      int L=-1;for(int i=1;i<len;i<<=1,++L);
      for(int i=0;i<len;i++) turn[i]=turn[i>>1]>>1|(i&1)<<L,A[i]=B[i]=0;
      for(int i=0;i<len>>1;i++) A[i]=a[i],B[i]=b[i];
      NTT(A,len,1),NTT(B,len,1);
      for(int i=0;i<len;i++) A[i]=Mul(A[i],B[i]);
      NTT(A,len,0);
      for(int i=0;i<len;i++) a[i]=A[i];
  }
  void inv(int *a,int n)
  {
      int cur=0,len=2;
      b[cur][0]=qp(f[0],mod-2);
      while(len<=(n<<2))
      {
          cur^=1;
          for(int i=0;i<len>>1;i++) b[cur][i]=add(b[cur^1][i],b[cur^1][i]);
          mul(b[cur^1],b[cur^1],len);
          mul(b[cur^1],f,len);
          for(int i=0;i<len;i++) b[cur][i]=add(b[cur][i],mod-b[cur^1][i]);
          len<<=1;
      }
      for(int i=0;i<n;i++)
          a[i]=b[cur][i];
  }
  int main()
  {
      int n;scanf("%d",&n);
      for(int i=0;i<n;i++) scanf("%d",f+i);
      inv(f,n);
      for(int i=0;i<n;i++) printf("%d ",f[i]);
      return 0;
  }
  

多项式除法

[A(x)=B(x)D(x)+R(x) ]

注：(A(x))被除多项式，度为(n)，(B(x))除他的多项式，度为(m)，(D(x))商，度为(n-m)，余数为(R(x))，度不大于(m-1)

• 一种转换

  (A^r(x)=x^nA(frac{1}{x}))

  感性理解：把它的系数颠倒了，可以(O(n))处理

  • 有什么用？

    通过取模把余数扔掉。

• [A(x)=B(x)D(x)+R(x) ]
  [x^nA(x)=x^nB(x)D(x)+x^nR(x) ]
  [A^r(x)=B^r(x)D^r(x)+x^{n-m+1}R^r(x) ]
  [A^r(x)equiv B^r(x)D^r(x) pmod {x^{n-m+1}} ]

  然后求个逆把(D^r(x))弄出来，带回去就行了

• Code:

  #include <cstdio>
  #include <algorithm>
  const int N=(1<<20)+10;
  const int mod=998244353,Ge=3,Gi=332748118;
  #define add(x,y) ((x+y)%mod)
  #define Mul(x,y) (1ll*(x)*(y)%mod)
  int qp(int d,int k){int f=1;while(k){if(k&1)f=Mul(f,d);d=Mul(d,d);k>>=1;}return f;}
  int F[N],G[N],b[2][N],X[N],Y[N],A[N],B[N],turn[N];
  void NTT(int *a,int len,int typ)
  {
      for(int i=0;i<len;i++)
          if(i<turn[i])
              std::swap(a[i],a[turn[i]]);
      for(int le=1;le<len;le<<=1)
      {
          int wn=qp(typ?Ge:Gi,(mod-1)/(le<<1));
          for(int p=0;p<len;p+=le<<1)
          {
              int w=1;
              for(int i=p;i<p+le;i++,w=Mul(w,wn))
              {
                  int tx=a[i],ty=Mul(w,a[i+le]);
                  a[i]=add(tx,ty);
                  a[i+le]=add(tx,mod-ty);
              }
          }
      }
      if(!typ)
      {
          int inv=qp(len,mod-2);
          for(int i=0;i<len;i++)
              a[i]=Mul(a[i],inv);
      }
  }
  void polymul(int *a,int *b,int len)
  {
      int L=-1;for(int i=1;i<len;i<<=1,++L);
      for(int i=0;i<len;i++) turn[i]=turn[i>>1]>>1|(i&1)<<L,A[i]=B[i]=0;
      for(int i=0;i<len>>1;i++) A[i]=a[i],B[i]=b[i];
      NTT(A,len,1),NTT(B,len,1);
      for(int i=0;i<len;i++) A[i]=Mul(A[i],B[i]);
      NTT(A,len,0);
      for(int i=0;i<len;i++) a[i]=A[i];
  }
  void polyinv(int *a,int n)
  {
      int cur=0,len=2;
      b[cur][0]=qp(a[0],mod-2);
      while(len<=(n<<2))
      {
          cur^=1;
          for(int i=0;i<len>>1;i++) b[cur][i]=add(b[cur^1][i],b[cur^1][i]);
          polymul(b[cur^1],b[cur^1],len);
          polymul(b[cur^1],a,len);
          for(int i=0;i<len;i++) b[cur][i]=add(b[cur][i],mod-b[cur^1][i]);
          len<<=1;
      }
      for(int i=0;i<n;i++)
          a[i]=b[cur][i];
  }
  void Reverse(int *a,int len)
  {
      for(int i=0;i<<1<len;i++) std::swap(a[i],a[len-i-1]);
  }
  void polydiv(int *a,int *b,int n,int m)
  {
      for(int i=0;i<n;i++) X[i]=a[i];
      for(int i=0;i<m;i++) Y[i]=b[i];
      Reverse(a,n),Reverse(b,m);
      polyinv(b,n-m+1);
      int len=1;while(len<=n<<2) len<<=1;
      polymul(a,b,len);
      Reverse(a,n-m+1);
      for(int i=n-m+1;i<len;i++) a[i]=0;
      polymul(Y,a,len);
      for(int i=0;i<m-1;i++) b[i]=add(X[i],mod-Y[i]);
  }
  int main()
  {
      int n,m;scanf("%d%d",&n,&m);++n,++m;
      for(int i=0;i<n;i++) scanf("%d",F+i);
      for(int i=0;i<m;i++) scanf("%d",G+i);
      polydiv(F,G,n,m);
      for(int i=0;i<n-m+1;i++) printf("%d ",F[i]);puts("");
      for(int i=0;i<m-1;i++) printf("%d ",G[i]);
      return 0;
  }
  

多项式开根

• 像求逆那样进行倍增

  [T^2equiv F pmod {2^n} ]
  [T'^2equiv Fpmod {2^{lceilfrac{n}{2} ceil}} ]
  [(T+T')(T-T')equiv 0 pmod {2^{lceilfrac{n}{2} ceil}} ]
  [(T-T')equiv 0pmod {2^{lceilfrac{n}{2} ceil}} ]
  [T^2-2TT'+T'^2equiv 0pmod {2^n} ]
  [F-2TT'+T'^2equiv 0pmod {2^n} ]
  [Tequiv frac{F+T'^2}{2T'}pmod {2^n} ]

  或者再化简一步

  [Tequiv frac{F}{2T'}+frac{T'}{2}pmod {2^n} ]

• Code:

  void polysqrt(int *a,int n)
  {
      int len=2,cur=0;
      sq[cur][0]=1;
      while(len<=(n<<2))
      {
          cur^=1;
          for(int i=0;i<len>>1;i++) sq[cur][i]=mul(sq[cur^1][i],i2);
          for(int i=0;i<len>>1;i++) sq[cur^1][i]=add(sq[cur^1][i],sq[cur^1][i]);
          polyinv(sq[cur^1],len);
          polymul(sq[cur^1],a,len);
          for(int i=0;i<len;i++) sq[cur][i]=add(sq[cur][i],sq[cur^1][i]);
          len<<=1;
      }
      for(int i=0;i<n;i++) a[i]=sq[cur][i];
  }
  

一定还会更新的...