[Notes][多项式]杂记 · 多项式算法—多项式求逆多项式取模多项式开根…

最新推荐文章于 2021-02-05 11:11:33 发布

LowestJN

最新推荐文章于 2021-02-05 11:11:33 发布

阅读量4.8k

点赞数 3

CC 4.0 BY-SA版权

分类专栏：多项式 Notes 文章标签：多项式

本文链接：https://blog.youkuaiyun.com/Coldef/article/details/76020530

多项式同时被 2 个专栏收录

7 篇文章

订阅专栏

Notes

7 篇文章

订阅专栏

本文详细介绍了多项式算法，包括多项式加减法、乘法、求逆、开根、取模等操作，并探讨了它们在生成函数和多项式取模解线性递推数列中的应用。文章提供了算法的基本原理、复杂度分析以及相关资源链接。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

多项式

由若干个单项式相加组成的代数式叫做多项式
形如： $f(x)=\sum_{i=0}^{n}a_ix^i$ ，
$\deg~f(x)$ 称为 $f$ 的度，是 $f(x)$ 最高次项的次数。

生成函数

形如 $\sum_{i=0}^\infty a_ix^i$
生成函数又称母函数，往往和多项式算法联系起来起到优化转移的作用

多项式算法

加减法

多项式加减法比较简单
若 $f(x)=\sum_{i=0}^na_ix^i$ ， $g(x)=\sum_{i=0}^nb_ix^i$

则 $(f\pm g)(x)=\sum_{i=0}^{n}(a_i\pm b_i)x^i$

代码实现也比较简单

乘法

多项式乘法是所有多项式算法的基础，也是生成函数的卷积运算

若 $f(x)=\sum_{i=0}^na_ix^i$ ， $g(x)=\sum_{i=0}^nb_ix^i$

则 $(f\times g)(x)=\sum_{i=0}^{2n}\sum_{j=0}^{i}a_j\times b_{i-j} x^i$

朴素的多项式乘法是 $O(n^2)$ 的，但FFT（快速傅里叶变换）运用复数根的特性可以在 $O(n\log n)$ 的复杂度内实现多项式的系数表达和点值表达的转化，运用点值表达实现 $O(n)$ 的相乘，总复杂度就是 $O(n\log n)$ ，但是因为常数过大，往往小范围会使用暴力。

inline void FFT(E *a,int r){
  for(int i=0;i<n;i++) if(rev[i]>i) swap(a[rev[i]],a[i]);
  for(int i=1;i<n;i<<=1){
    E wn(cos(M_PI/i),r*sin(M_PI/i));
    for(int j=0;j<n;j+=(i<<1)){
      E w(1,0);
      for(int k=0;k<i;k++,w=w*wn){
    E x=a[j+k],y=w*a[j+k+i];
    a[j+k]=x+y; a[j+k+i]=x-y;
      }
    }
  }
  if(r==-1) for(int i=0;i<n;i++) a[i].real/=n;
}

运用生成函数的计数问题往往会对一个质数取模，如果这个质数可以表示成 $k\times 2^n+1$ ，其中 $k$ 为奇数 $n$ 大于多项式的度数，那么就可以用这个质数的原根代替复数根，实现多项式的系数对一个质数取模，这个算法就是NTT（快速数论变换），这种模数称为NTT模数。

inline void Pre(int n){
  num=n;
  int g=Pow(3,(P-1)/num);
  w[0][0]=w[1][0]=1; for(int i=1;i<num;i++) w[0][i]=1LL*w[0][i-1]*g%P;
  for(int i=1;i<num;i++) w[1][i]=w[0][num-i];
}

inline void NTT(int *a,int n,int r){
  for(int i=1;i<n;i++) if(rev[i]>i) swap(a[i],a[rev[i]]);
  for(int i=1;i<n;i<<=1)
    for(int j=0;j<n;j+=i<<1)
      for(int k=0;k<i;k++){
        int x=a[j+k],y=1LL*a[j+k+i]*w[r][num/(i<<1)*k]%P;
        a[j+k]=(x+y)%P; a[j+k+i]=(x+P-y)%P;
      }
  if(!r) for(int i=0,inv=Pow(n,P-2);i<n;i++) a[i]=1LL*a[i]*inv%P;
}

不过有些丧心病狂的题的模数不是NTT模数，这个时候需要CRT合并。

多项式求逆

比如BZOJ4555

最后会得到 $f(x)=f(x)g(x)+1$ ，这是一个非齐次线性递推，因为求的第n项不大，只要直接求出 $f(x)$ 就行了

很容易可以得到 $f(x)={1\over 1-g(x)}$ ，如果令 $h(x)=1-g(x)$ ，那么 $f(x)=(h(x))^{-1}$ ，所以就要求 $h(x)$ 的逆元
求逆元的话网上讲解也有很多，鏼爷集训队论文里也有，大概是用了倍增的做法，摘一下
求 $A(x)$ 的逆元，即是求 $B(x)$ ，使 $A(x)B(x)\equiv 1\pmod {x^{2^k}}$

如果我们已经求出 $G(x)$ ， $A(x)G(x)\equiv 1\pmod {x^n}$
那么

A (x) G (x) \equiv 1 (mod x n)

$A(x)G(x)\equiv 1\pmod {x^n}$

A (x) G (x) - 1 \equiv 0 (mod x n)

$A(x)G(x)-1\equiv 0\pmod {x^n}$

(A (x) G (x) - 1) 2 \equiv 0 (mod x 2 n)

$(A(x)G(x)-1)^2\equiv 0\pmod {x^{2n}}$

A 2 (x) G 2 (x) - 2 A (x) G (x) + 1 \equiv 0 (mod x 2 n)

$A^2(x)G^2(x)-2A(x)G(x)+1\equiv 0\pmod {x^{2n}}$

A (x) (2 G (x) - A (x) G 2 (x)) \equiv 1 (mod x 2 n)

$A(x)(2G(x)-A(x)G^2(x))\equiv 1\pmod {x^{2n}}$

n=0 n = 0 $n=0$ 时，求一下

A(x) A ( x ) $A(x)$ 常数项的乘法逆元就可以了

那么用 $O(n\log n)$ 的复杂度把 $B(x)$ 赋成 $2G(x)-A(x)G^2(x)$ 进入下一层运算就可以了
总复杂度 $T(n)=T(n/2)+\mathcal{O}(n\log n)=\mathcal{O}(n\log n)$

void Inv(int *a,int *b,int n){
    if(n==1){
        b[0]=Pow(a[0],P-2); return ;
    }
    Inv(a,b,n>>1);
    for(int i=0;i<n;i++) tmp[i]=a[i],tmp[i+n]=0;
    int L=0; while(!(n>>L&1)) L++;
    for(int i=1;i<(n<<1);i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
    FFT(tmp,n<<1,1); FFT(b,n<<1,1);
    for(int i=0;i<(n<<1);i++)
        tmp[i]=(1LL*b[i]*2+P-1LL*tmp[i]*b[i]%P*b[i]%P)%P;
    FFT(tmp,n<<1,0);
    for(int i=0;i<n;i++) b[i]=tmp[i],b[i+n]=0;
}

应用的话，就是一些解一些现行递推式，或者求伯努利数
预处理伯努利数：http://blog.youkuaiyun.com/coldef/article/details/72876762

学习自：http://blog.miskcoo.com/2015/05/polynomial-inverse

多项式开根

就是对 $A(x)$ ，求 $B(x)$ ，使 $B^2(x)\equiv A(x)\pmod {x^n}$
同样是用倍增的方法。
已求得 $G(x)$ ，使 $G^2(x)\equiv A(x)\pmod {x^n}$
那么

B 2 (x) \equiv G 2 (x) (mod x n)

$B^2(x)\equiv G^2(x)\pmod {x^n}$

[B 2 (x) - G 2 (x)] 2 \equiv 0 (mod x 2 n)

$[B^2(x)-G^2(x)]^2\equiv 0\pmod {x^{2n}}$

B 4 (x) - 2 B 2 (x) G 2 (x) + G 4 (x) \equiv 0 (mod x 2 n)

$B^4(x)-2B^2(x)G^2(x)+G^4(x)\equiv 0\pmod {x^{2n}}$

B 4 (x) + 2 B 2 (x) G 2 (x) + G 4 (x) \equiv 4 B 2 (x) G 2 (x) (mod x 2 n)

$B^4(x)+2B^2(x)G^2(x)+G^4(x)\equiv 4B^2(x)G^2(x)\pmod {x^{2n}}$

[B 2 (x) + G 2 (x)] 2 \equiv [2 B (x) G (x)] 2 (mod x 2 n)

$[B^2(x)+G^2(x)]^2\equiv [2B(x)G(x)]^2\pmod {x^{2n}}$

B 2 (x) + G 2 (x) \equiv 2 B (x) G (x) (mod x 2 n)

$B^2(x)+G^2(x)\equiv 2B(x)G(x)\pmod {x^{2n}}$

A (x) + G 2 (x) \equiv 2 B (x) G (x) (mod 2 2 n)

$A(x)+G^2(x)\equiv 2B(x)G(x)\pmod {2^{2n}}$

B (x) \equiv A ( x ) + G 2 ( x ) 2 G ( x ) (mod 2 2 n)

$B(x)\equiv {A(x)+G^2(x)\over 2G(x)}\pmod {2^{2n}}$

n=0 n = 0 $n=0$ 的时候就是

B(x) B ( x ) $B(x)$ 的常数项就是1（我也不知道为什么，Manchery说的）
时间复杂度

T(n)=T(n/2)+多项式求逆的复杂度+O(nlogn)=O(nlogn) T ( n ) = T ( n / 2 ) + 多项式求逆的复杂度 + O ( n log ⁡ n ) = O ( n log ⁡ n ) $T(n)=T(n/2)+多项式求逆的复杂度+\mathcal{O}(n\log n)=\mathcal{O}(n\log n)$

void Sqrt(int *a,int *b,int n){
    if(n==1) return void(b[0]=1);
    Sqrt(a,b,n>>1);
    memset(invb,0,sizeof(int)*n); Inv(b,invb,n);
    int L=0; while(!(n>>L&1)) L++;
    for(int i=1;i<(n<<1);i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
    for(int i=0;i<n;i++) tmp[i]=a[i],tmp[i+n]=0;
    FFT(invb,n<<1,1); FFT(tmp,n<<1,1); 
    for(int i=0,inv2=P+1>>1;i<(n<<1);i++) tmp[i]=1LL*tmp[i]*inv2%P*invb[i]%P;
    FFT(tmp,n<<1,0);
    for(int i=0,inv2=P+1>>1;i<n;i++) b[i]=(1LL*b[i]*inv2+tmp[i])%P;
}

学习自： http://www.cnblogs.com/acha/p/6472798.html

多项式取模

此取模不是指系数取模，而是多项式模多项式

UPD:
暴力的多项式取模比较简单，
设模多项式为 $M(x)=\sum_{i=0}^n a_i x^i$ ，化一下可以得到 $x^n=\sum_{i=0}^{n-1}{a_i\over a_n}x^i$ ，将所有大于 $n$ 的项带入就行

对于 $A(x)$ ， $B(x)$ ，设 $\deg A(x)=n , \deg B(x)=m(m \le n)$
求 $D(x)$ ， $R(x)$ ，满足 $A(x)=D(x)B(x)+R(x)$
只要求出 $D(x)$ ，通过多项式加减乘就可以得到 $R(x)$
把 ${1\over x}$ 带入多项式

A (1 x) = D (1 x) B (1 x) + R (1 x)

$A({1\over x})=D({1\over x})B({1\over x})+R({1\over x})$

两边乘 $x^n$

x n A (1 x) = x n - m D (1 x) x m B (1 x) + x n R (1 x)

$x^nA({1\over x})=x^{n-m}D({1\over x})x^mB({1\over x})+x^nR({1\over x})$

设 $f^R(x)=x^{\deg f(x)}f({1\over x})$

则

A R (x) = D R (x) B R (x) + x n - deg R (x) R R (x)

$A^R(x)=D^R(x)B^R(x)+x^{n-\deg R(x)}R^R(x)$

发现 $\deg R(x)< m$

那么

A R (x) \equiv D R (x) B R (x) (mod x n - m + 1)

$A^R(x)\equiv D^R(x)B^R(x)(\mod x^{n-m+1})$

A R ( x ) B R ( x ) \equiv D R (x) (mod x n - m + 1)

${A^R(x)\over B^R(x)}\equiv D^R(x)(\mod x^{n-m+1})$

这样多项式求个逆就可以求出 $D(x)$ 啦
带回去加加减减就可以得到 $R(x)$

inline void Div(int *a,int n,int *b,int m){
  static int tmp[N],A[N],B[N];
  for(int i=0;i<m;i++) t[i]=b[m-i-1];
  for(int i=0;i<n;i++) A[i]=a[n-i-1];
  int nn=1,d=n-m+1; for(;nn<d<<1;nn<<=1);
  for(int i=n;i<nn;i++) A[i]=0;
  for(int i=d;i<nn;i++) t[i]=0;
  for(int i=0;i<nn;i++) B[i]=0;
  Inv(t,B,nn);
  for(int i=d;i<nn;i++) B[i]=0;
  Mul(A,B,max(n,d));
  for(int i=d;i<=n<<1;i++) A[i]=0;
  for(int i=0;i<d;i++) if(i>d-i-1) swap(A[i],A[d-i-1]);
  for(int i=0;i<m;i++) t[i]=b[i];
  Mul(t,A,max(d,m));
  for(int i=0;i<n;i++) a[i]=(a[i]+P-t[i])%P;
}

学习自：http://blog.miskcoo.com/2015/05/polynomial-division

~~大致就是这些算法~~

~~还有多项式求exp什么的…还不会…先挖个坑~~

多项式exp

可以参考金策的2015年国家集训队论文

B (x) = e A (x)

$B(x)=e^{A(x)}$

两边同求导

B' (x) = e A (x) A' (x)

$B'(x)=e^{A(x)}A'(x)$

B' (x) = B (x) A' (x)

$B'(x)=B(x)A'(x)$

对于每一项的系数

(i + 1) b i + 1 = \sum j + k = i b j a k + 1 (k + 1)

$(i+1)b_{i+1}=\sum_{j+k=i}b_ja_{k+1}(k+1)$

也就是说

b i = 1 i \sum j = 1 i j a j b i - j

$b_i={1\over i}\sum_{j=1}^ija_jb_{i-j}$

分治FFT

牛顿迭代不会~~

void exp(int *a,int *b,int l,int r){
  if(l==r){
    add(b[l],1LL*inv[l]*a[l]%P); return ;
  }
  int mid=l+r>>1;
  exp(a,b,l,mid); 
  static int tmpa[N],tmpb[N];
  for(int i=l;i<=mid;i++) tmpa[i-l]=b[i];
  for(int i=1;i<=r-l;i++) tmpb[i]=a[i];
  int m,L=0;
  for(m=1;m<r-l+1;m<<=1,L++); m<<=1;
  for(int i=1;i<m;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
  NTT(tmpa,m,1); NTT(tmpb,m,1);
  for(int i=0;i<m;i++) tmpa[i]=1LL*tmpa[i]*tmpb[i]%P;
  NTT(tmpa,m,0);
  for(int i=mid+1;i<=r;i++) add(b[i],1LL*inv[i]*tmpa[i-l]%P);
  for(int i=0;i<m;i++) tmpa[i]=tmpb[i]=0;
  exp(a,b,mid+1,r);
}