多项式全家桶【长期更新】

多项式

定义(表达式)

定义一个 n 次的多项式为:

F(x)=i=0nfixi

注意 :

  1. i=0 的那一项是常数项。
  2. 易得F(0)=f0
  3. 在OI中,多项式一般在同余条件下讨论mod p,在下文中会省略 "" 符号。
  4. 求多项式F(x),实际上就是把每一个fi求出来。

暴力全家桶

加法

给定两个多项式 F(x)G(x),求它们的和 H(x)
即: $ H(x) = F(x) + G(x) h_i = f_i + g_i O(n) $

乘法

即:H(x)=F(x)G(x)

hi=j=0ifjgij

就是合并同类项

复杂度 O(n2)

余数除法

已知F和H,求G和R

即:H(x)=F(x)G(x)+R(x)

就是小学奥数里的大除法。

O(n2)

求导和积分

若:G(x)=F(x)
则:gi=(i+1)×fi+1
若:H(x)=F(x)dx
则:hi=fi1i

求逆

F已知求G。

即:F(x)G(x)=1(mod p)

g0=f01(x=0,f01)j=0ifjgij=[i=0]f0gi=[i=0]j=1ifjgij 把第i项单独拿出来gi=[i=0]j=1ifjgijf0

递推求即可。

O(n2)

开根

G2(x)=F(x),已知F求G

0jigjgij=fi2gig0=fi0jigjgij就是把第一项和最后一项提出来,类比求逆gi=fi0jigjgij2g0

还是用递推,O(n2)

特别的,g02=f0 作为边界情况。

  • f00,则 G 的解数,取决于 f0 的平方根数量。(具体参见二次剩余)

  • f0=0,令 F(x)=xkH(x),满足 h(0)0,那么若 k 为奇数,则 G 无解;若 k 为偶数,记 P2(x)=H(x),那么 G(x)=xk2P(x)

这是因为开跟运算本身不一定有解。

求对数

给定多项式 F(x),求 G(x)lnF(x)(modxn),保证 f0=1

解:对原等式两边求导,得:

G(x)F(x)F1(x)(modxn)H(x)F1(x)(modxn)G(x)F(x)H(x)(modxn)(i+1)gi+1=j=0i(j+1)fj+1hij,g0=ln(f0)=0gi=j=0i1(j+1)fj+1hij1i

O(n2)

求指数

给定多项式 F(x),求 G(x)eF(x)(modxn),保证 f0=0

对原等式两边求导,得:

g0=ef0=e0=1

G(x)F(x)eF(x)(复合函数求导)F(x)G(x)(modxn)

(i+1)gi+1=j=0i(j+1)fj+1gij

此式子也可看成求对数的第三步和第四步,把F换成G,把G换成F

gi+1 时已求出 g0gi, 故可以递推求出 G(x)

O(n2)

求三角函数

前置知识:[欧拉公式(eix)](https://www.cnblogs.com/water-flower/p/18651768#- 定义五:在复数中的定义(欧拉公式)-) ,泰勒展开

由欧拉公式:

eix=cosx+isinx

eix=cosxisinx

两式相加得

cosx=eix+eix2

两式相减得

sinx=eixeix2i

现在我们把三角函数换成了指数函数,用指数的方法来推导三角函数。

O(n2)

快速傅里叶变换

终于进入正题了。

前置内容

  1. 多项式插值:n + 1 个点值可确定一个 n 次多项式。

    很好理解。用待定系数法高斯消元。

    系数矩阵满秩,所以不会是无数解。

    有如果无解,说明有两行的系数,上面的系数都是下面的系数的k倍。但是对于不同的x,{xi} 的集合不会是这种比例关系。

  2. 复数运算:参见数学书

  3. 单位根
    P.S. 以下为了方便书写,均将 ω 简写为 w,且将 ωni 表示为 wi

    定理:任何复系数一元n次多项式方程在复数域上至少有一个根。

    大多数数学家都认为这是对的。

    推论:任何复系数一元 n 次多项式方程在复数域上恰好有 n 个根。

    设现有一根 xi

    若将原 n 次多项式因式分解,必定存在一项形如 (xxi)。把这一项去掉,得到一个 n-1 次的多项式。而这个新的多项式必有一根。再把这个根去掉。重复这个操作。因为可以降幂 n 次,所以有 n 个根。

    单位根定义xn1=0n 个根,记作 w0,w1...wn1

    单位根几何意义:建立复数域的坐标系。做单位圆。运算用向量的运算。这里以 n=6 为例。

    用三角函数来表达这些单位根可得: wi=cos2iπn+isin2iπn

    首先显然有:wkn+i=wi,这相当与多转 k 圈。

    通过化简,(wi)2=(cos2iπn+isin2iπn)2=cos22iπn+isin22iπn=w2i,这个可以看为把和x轴的夹角翻倍。

    同理(大概是用一些二项式定理和欧拉定理展开后乱搞),(wi)n=win=w0

    所以一个单位根的 n 次方都是 1+0i

    下面给出一些关于单位根的运算性质 :

    • wi+n=wi

    • wiwj=wi+j

    • (wi)j=wij 次方转下标

    • wi=wi

    • wknki=wni

      单位根反演

      (应该只在IDFT的证明中有用到)

      这里用 wi 代表 win

      i=0n1wik=n[nk]

      证明:

      i=0n1wik=i=0n1wki

      该式为等比数列,所以:

      i=0n1wki={n if wk=11wkn1wk if wk!=1

      观察到,wkn=wnk=1,并且 wk=1n|k

      所以上式化简为:i=0n1wik=n[nk]

FFT

这玩意可以做到 O(nlogn) 的多项式乘法。老牛逼了。

核心思路

根据多项式插值,我们根据 F(x)G(x) 的 n + 1 个点值,然后用插值法,可以插出这个多项式。

通过优秀的选点(单位根),可以做到 O(nlogn)

现在我们有两个要做的事情:

  1. DFT:输入一个 n1 次多项式的系数列,快速得到其在 n 次单位根处的点值列,也就是求 F(wi)
  2. IDFT:在优秀的时间复杂度内插值出多项式,也就是用 H(xi) 去求 hi(点值求系数)。

具体做法

先来解决DFT:

(1)F(wi)=i=0n1fiwki(2)=i=0n21f2iwk2i+i=0n21f2i+1wk2i+1把i按照奇偶分组(3)=i=0n21f2iwk2i+wk1i=0n21f2i+1wk2i

DFT(F)=DFT({f0,f2fn2})k+wkDFT({f1,f3fn1})kF(x)=F1(x)+xF2(x)

带入 wk(k<n2)F(wk)=F1(wk)+wkF2(wk)

带入 wk+n2(k<n2)F(wk)=i=0n21f2iwk2i+wk+n21i=0n21f2i+1wk2i

wk+n22i=w2k+ni=w2ki=wk2iwk+n2=wk 这里 wk+n2n2 为底, wk 以 n 为底。

所以带入 wk+n2(k<n2)F(wk)=F1(wk)wkF2(wk)

于是可以递归求解。

复杂度是 T(n)=2T(n2)+O(n) ,是 O(nlogn)

可以发现,DFT能做到高效的原理实际上是利用的单位根的运算性质。

再来看IDTF:

根据单位根反演可以推到出这个公式:

nhk=i=0n1wkiH(xi)

证明:

(4)i=0n1wkiH(xi)=i=0n1wkij=0n1wijhj(5)=j=0n1hji=0n1wkiwij(6)=j=0n1hji=0n1wijk(7)=j=0n1hjn[n|jk](8)=nhk(只有 j=k 的那一项有值)

这个公式与DTF相比,左边多一个 n ,右边的 k 变成了 -k,其余一致。

常数优化:非递归FFT

注意到,第 k 个点和第 k + len / 2 个点 由 第 k 个点和第 k + len/2 共同转移来,并且每次len除2。

若执行到了第 k 次,那么把从右往左数第 k 位是 0 的数往前提。

所以对于最后一排,最后一位是 0 的排在前面,如果最后一位都是 0,那么倒数第二排是 0 的排在前面...以此类推。

容易发现,这相当于倒着比较二进制数的大小。

考虑最后一行的初始化:第 i 位就应该为 fi表示二进制取反。

写成代码是这样的:

  for (int i = 0; i < len; ++i) {
    rev[i] = rev[i >> 1] >> 1;
    if (i & 1) rev[i] |= len >> 1;\\变换最高位从左往右第一位?神秘
  }
  for (int i = 0; i < len; ++i) if (i < rev[i]) swap(y[i], y[rev[i]]); \\这样才能刚好交换一次

之后就可以从低位向高位递推了。

附上代码:$$

typedef complex<double > com;
const int N=(1e6+6) * 4;
const double PI = acos(-1);
int n,m,rev[N];
com f[N],g[N];

void change(int n) { 
	int L = log2(n);
	for(int i=0;i<n;++i)rev[i] = (rev[i>>1]>>1) | ((i&1) << (L-1));
}
void FFT(com *f, int n,int op) {
    change(n)
	for(int i=0;i<n;++i) if(rev[i] < i) swap(f[i], f[rev[i] ]);
	for(int mid = 1;mid < n;mid <<= 1){
		int j = mid << 1; 
		com nxt(cos(PI / mid), sin(PI / mid) * op);
		for(int st = 0;st < n; st += j){
			com w(1.0, 0.0);
			for(int i = st;i < st + (j>>1); ++i, w *= nxt){
				com tmp1 = f[i],tmp2 = w * f[i + mid];
				f[i] = tmp1 + tmp2;
				f[i + mid] = tmp1 - tmp2;
			}
		}
	}
}

NTT

除法意味着精度误差,所以很多人不喜欢除法,而使用在模数的逆元来替代除法。

FFT用的是单位根的点值,而NTT则运用原根,并且可以避免除法运算,被广泛使用。

原根定义:对于模数 m,满足 gcd(g,m)=1gk1(modm) 最小正整数 k 等于 φ(m) 的一个数 g。注意,一个模数 m 可能有多个原根。

原根有和单位根一样优良的性质。设 g 表示 p 的一个原根,则令 wi=giφ(p)n

  • wi+n=wi
  • wiwj=wi+j
  • wij=wij
  • 1gi=gi
  • giφ(p)n = (gφ(p)n)i

但是这不代表我们可以直接将他套用进 FFT 中。因为这里必须要求 nφ(p),所以这对 n 产生了很大的限制。

幸运的是我们发现质数 998244353 的欧拉函数值 998244352=119×223,因此当 n223 时,都能够很好的进行这个算法,223O(nlogn) 算法能通过的数据范围略大一点,所以大多数时候可以直接使用这个NTT算法,这也是 998244353 成为常用模数的原因。(尽管有些时候题目与多项式无关)

998244353 的原根有3114514这使得这个数字充满了美好的寓意

代码(摘自 oi-wiki.org)

void ntt(int *x, int lim, int opt) {
  int i, j, k, m, gn, g, tmp;
  for (i = 0; i < lim; ++i)
    if (r[i] < i) swap(x[i], x[r[i]]);
  for (m = 2; m <= lim; m <<= 1) {
    k = m >> 1;
    gn = qpow(3, (P - 1) / m);
    for (i = 0; i < lim; i += m) {
      g = 1;
      for (j = 0; j < k; ++j, g = 1ll * g * gn % P) {
        tmp = 1ll * x[i + j + k] * g % P;
        x[i + j + k] = (x[i + j] - tmp + P) % P;
        x[i + j] = (x[i + j] + tmp) % P;
      }
    }
  }
  if (opt == -1) {
    reverse(x + 1, x + lim);
    int inv = qpow(lim, P - 2);
    for (i = 0; i < lim; ++i) x[i] = 1ll * x[i] * inv % P;
  }
}

快速多项式全家桶

加法和乘法和求导积分略过。

前置知识:牛顿迭代

对于一个方程 F(x)=0 求近似解,我们先瞎猜一个解 x ,算出 F(x)(x,F(x)) 处的切线,算切线与 x 轴的交点,把交点 x1 当成新的解并重复此操作。

当我们把 x 换成一个函数,我们发现牛顿迭代又可以用于求解函数。

假如一个函数 F(x) 满足一个关系 G(F(x))=0,这里会把 F(x) 当作未知量,并且 G 是一个泛函数。

为了方便表示,设 FkF(modx2k) ,能得到 Fk 满足 G(Fk)%x2k=0,我们希望用 Fk 求出 Fk+1,直到达到精度要求。

那么如果我们先预估一个函数 F0(x),就能得到 F1(x)=F0(x)G(F0(x)))G(F0(x)) ,并不断继续重复以提高精度。

证明:

因为 G(Fk+1(x))%x2k+1=0,将其在 Fk 处泰勒展开得:

i=0+\infinG(i)(Fk(x))i!(F(x)k+1Fk(x))i0(modx2k+1)

为了求 Fk+1 我们将其对 x2k+1 次方取模

因为 Fk2k 次的,所以 (Fk+1(x)Fk(x))i 的最低次项是 2ki 次的,所以只有 i 等于 0,1 时才有值。化简上式:

i=01G(i)(Fk(x))i!(Fk+1(x)Fk(x))i=0G(Fk(x))+G(Fk(x))(Fk+1(x)Fk(x))=0Fk+1(x)=G(Fk(x))Fk(x)G(Fk(x))G(Fk(x))=Fk(x)G(Fk(x))G(Fk(x))

求逆

F已知求G。

即:F(x)G(x)=1(mod p)

我们要想办法把这玩意化成可以牛顿迭代的形式:

可以得到 H(G)=1GF=0。因为 F 已知,这里被当作常数,所以 H 就是关于 G 的函数。

迭代的初值 : g0=f01 (带入 x=0,f01 用普通的求逆元的方法即可)

使用牛顿迭代的公式得到 : Gk+1=GkH(Gk)H(Gk)=Gk1GkF1Gk2

我们对这个式子化简(分子分母同乘 Gk2 )可得: Gk+1=2GkFGk2.

时间复杂度 T(n)=T(n/2)+O(nlogn)=O(nlogn).

余数除法

已知F和H,求G和R

即:H(x)=F(x)G(x)+R(x)

思路:如果没有 R(x),这个题就用逆元做,所以要想办法避开 R.

用普通的方法很难求,但是观察到,R 的次数比 H 低,于是考虑把函数系数翻转后,发现R的后几项的系数为0,再用取模的方式把 R 干掉。

做法:设 n=degH,m=degF

rH(x)=xnH(1x),即把 H(x) 系数翻转之后得到的函数,同理有 rF(x)=xmF(1x)rG(x)=xnmG(1x)rR(x)=xmR(1x)

容易得知 rF(x)rQ(x)rG(x)(modxnm+1)

证明:

H(x)=F(x)G(x)+R(x)xnH(1x)=xmF(1x)xnmG(1x)+xnR(1x)rH=rFrG+xnmrRrHrFrG(modxnm+1) 哦,这真是一个关键点!

rG 的次数刚好是 nm 次的,所以可以用求逆的方法算出 rG ,于是可以算出 G,于是可以算出 R

复杂度同求逆 O(nlogn).

开根

G2(x)=F(x),已知F求G

转成可以牛顿迭代的形式:H(G)=G2F=0.

Gk+1=GkGk2F2GkGk(x)2+F(x)2Gk(x)(modx2k+1)

O(nlogn)

特别的,g02=f0 作为边界情况。

  • f00,则 G 的解数,取决于 f0 的平方根数量。(具体参见二次剩余)

  • f0=0,令 F(x)=xkH(x),满足 h(0)0,那么若 k 为奇数,则 G 无解;若 k 为偶数,记 P2(x)=H(x),那么 G(x)=xk2P(x)

这是因为开跟运算本身不一定有解。

ln

给定多项式 F(x),求 G(x)lnF(x)(modxn),保证 f0=1

解:对原等式两边求导,得:

G(x)F(x)F1(x)(modxn)

g0=0

O(nlogn)

exp

给定多项式 F(x),求 G(x)eF(x)(modxn),保证 f0=0

g0=ef0=e0=1

对两边取 ln

lnGF=0Gk+1=GklnGkF1Gk=Gk(1lnGk+F)

O(nlogn)

求三角函数

cosx=eix+eix2

sinx=eixeix2i

现在我们把三角函数换成了指数函数,用指数的方法来推导三角函数。

O(nlogn)

多项式与分治

给定长为 n 的序列 ai,求多项式 i=1n(xai)

n105

i=1n(xai)=i=1mid(xai)×i=mid+1n(xai)

看作两个多项式相乘,用 FFT

T(n)=2T(n/2)+O(nlogn)=O(nlogn2)

参考资料:上课的课件,大佬的博客

posted @   花子の水晶植轮daisuki  阅读(19)  评论(1编辑  收藏  举报
相关博文:
阅读排行:
· 在鹅厂做java开发是什么体验
· 百万级群聊的设计实践
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
· 永远不要相信用户的输入:从 SQL 注入攻防看输入验证的重要性
· 浏览器原生「磁吸」效果!Anchor Positioning 锚点定位神器解析
https://blog-static.cnblogs.com/files/zouwangblog/mouse-click.js
点击右上角即可分享
微信分享提示