Fork me on GitHub

多项式运算

Posted on | In

多项式运算

本文所讲所有的多项式运算全部基于vector的使用

1
typedef vector<int> Poly;

0.思想

自认为多项式最重要的思想就是倍增了

其次就是分治

建议先看看fft残疾人手册和ntt残疾人手册

1.多项式加法

$C(n)=A(n)+B(n)$

就是系数相加，非常简单的实现

1
2
3
4
5
6
Poly add(Poly a, Poly b) {
  a.resize(max(a.size(), b.size()));
  for (size_t i = 0; i < b.size(); i++)
    a[i] = add(a[i], b[i]);
  return a;
}

2.多项式减法

$C(n)=A(n)-B(n)$

系数相减，依旧简单的实现

1
2
3
4
5
6
Poly sub(Poly a, Poly b) {
  a.resize(max(a.size(), b.size()));
  for (size_t i = 0; i < b.size(); i++)
    a[i] = sub(a[i], b[i]);
  return a;
}

3.多项式乘法

$C(n)=\sum_{i=0}^nA(i)B(n-i)$

就是ntt

需要预处理原根的次幂

交换位置的时候动态处理就好，因为在倍增过程中ntt的长度不一样

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
int w[2][N];

void init() {
  for (int i = 1; i < (1 << 18); i <<= 1) {
    w[1][i] = w[0][i] = 1;
    int wn;
    wn = fast_pow(G, (Mod - 1) / (i << 1));
    for (int j = 1; j < i; j++) 
      w[1][i + j] = mul(wn, w[1][i + j - 1]);
    wn = fast_pow(G, Mod - 1 - (Mod - 1) / (i << 1));
    for (int j = 1; j < i; j++)
      w[0][i + j] = mul(wn, w[0][i + j - 1]);
  }
}

void transform(int *t, int len, int typ) {
  for (int i = 0, j = 0, k; j < len; j++) {
    if (j > i) swap(t[i], t[j]);
    for (k = (len >> 1); k & i; k >>= 1) i ^= k;
    i ^= k;
  }
  for (int i = 1; i < len; i <<= 1) {
    for (int j = 0; j < len; j += (i << 1)) {
      for (int k = 0; k < i; k++) {
        int x = t[j + k], y = mul(t[j + k + i], w[typ][i + k]);
        t[j + k] = add(x, y);
        t[j + k + i] = sub(x, y);
      }
    }
  }
  if (typ) return;
  int inv = fast_pow(len, Mod - 2);
  for (int i = 0; i < len; i++)
    t[i] = mul(t[i], inv);
}

Poly mul(const Poly &a, const Poly &b) {
  int len = a.size() + b.size() + 1;
  len = 1 << (int) ceil(log2(len));
  static Poly prea, preb;
  prea = a;
  preb = b;
  prea.resize(len);
  preb.resize(len); 
  transform(&prea[0], len, 1);
  transform(&preb[0], len, 1); 
  for (int i = 0; i < len; i++)
    prea[i] = mul(prea[i], preb[i]);
  transform(&prea[0], len, 0);
  clean(prea);
  return prea;
}

4.多项式求逆

对于多项式$A$，我们要求出$B$

满足$A\ast B=1(\mod x^{n+1})$

这里的$\mod x^n$是保留前$n+1$个系数的意思，不然$B(n)$会有无穷项

我们假设已经求出了$B’$满足$A\ast B’=1(\mod x^{\lceil\frac{n}{2}\rceil})$

那么因为存在$A\ast B=1(\mod x^{\lceil\frac{n}{2}\rceil})$

所以$A\ast (B-B’)=0(\mod x^{\lceil\frac{n}{2}\rceil})$

又因为A不是0，所以$(B-B’)=0(\mod x^{\lceil\frac{n}{2}\rceil})$

两边同时平方（模数也要平方）

$B^2-2\ast B\ast B’+B’^2=0(\mod x^n)$

因为有个$B^2$不好计算，所以把等式两边同时乘上$A$

得到$B-2\ast B’+A\ast B’^2=0(\mod x^n)$

即$B=2\ast B’-A\ast B’^2(\mod x^n)$

这样就可以递归求解了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
Poly inv(Poly a, int n) {
  if (n == 1) return Poly(1, fast_pow(a[0], Mod - 2));
  int len = 1 << ((int) ceil(log2(n)) + 1);
  Poly x = inv(a, (n + 1) >> 1), y;
  x.resize(len);
  y.resize(len);
  for (int i = 0; i < n; i++)
    y[i] = a[i];
  transform(&x[0], len, 1);
  transform(&y[0], len, 1);
  for (int i = 0; i < len; i++)
    x[i] = mul(x[i], sub(2, mul(x[i], y[i])));
  transform(&x[0], len, 0);
  x.resize(n);
  return x;
}

Poly inv(Poly a) {
  return inv(a, a.size());
}

5.多项式除法

对于n次多项式$A$和m次多项式$B$，构造出小于n-m次多项式$C$和小于m次的多项式$D$

满足$A=B*C+D$

我们考虑一下多项式$A$的系数反转形式记做$RevA$

所以$RevA(x)=x^nA(\frac{1}{x})$

我们把原来的式子左右同时乘上$x^n$

得到$A(\frac{1}{x})\ast x^n=B(\frac{1}{x})\ast C(\frac{1}{x})\ast x^n+D(\frac{1}{x})\ast x^n$

所以$RevA(x)=RevB(x)\ast RevC(x)+x^{n-m+1}RevD(x)$

在$\mod x^{n-m+1}$的意义下发现$RevD$就没了，然而并不影响$RevA,RevB,RevC$

所以我们只需要求出$RevB$在$\mod x^{n-m+1}$意义下的逆元，然后就可以算出$RevC$和$RevD$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Poly rev(Poly a) {
  clean(a);
  reverse(a.begin(), a.end());
  return a;
}

void div(Poly a, Poly b, Poly &c, Poly &d) {
  int n = a.size(), m = b.size();
  Poly ra = rev(a), rb = rev(b);
  ra.resize(n - m + 1);
  rb.resize(n - m + 1);
  c = mul(ra, inv(rb));
  c.resize(n - m + 1);
  c = rev(c);  
  d = sub(a, mul(b, c));
  clean(c), clean(d);
}

6.多项式求导

因为多项式的形式很简单，所以求导直接可以模拟出来

1
2
3
4
5
6
7
Poly deri(Poly a) {
  int n = a.size();
  for (int i = 1; i < n; i++)
    a[i - 1] = mul(a[i], i);
  a.resize(n - 1);
  return a;
}

7.多项式积分

积分也直接模拟吧。。

1
2
3
4
5
6
7
8
Poly inte(Poly a) {
  int n = a.size();
  a.resize(n + 1);
  for (int i = n; i >= 1; i--)
    a[i] = mul(a[i - 1], fast_pow(i, Mod - 2));
  a[0] = 0;
  return a;
}

8.多项式对数函数

对于n次多项式$A$我们需要求出n次多项式$B$满足$B=\ln(A)$

两边同时求导

$B’=\frac{A’}{A}$

所以把$A$的导和逆乘起来然后积分回去就可以了。。

1
2
3
4
5
6
Poly ln(Poly a) {
  int len = a.size();
  a = inte(mul(deri(a), inv(a)));
  a.resize(len);
  return a;
}

9.多项式指数函数

对于n次多项式$A$我们需要求出来n次多项式$B$满足$B=e^A$

所以$\ln(B)=A $

然后我们构造一个函数$F(B)=\ln(B)-A$

我们要求函数$F(A)=0$的时候的解

根据泰勒展开和牛顿迭代：

$f(x)=f(x_0)+f’(x_0)(x-x_0)$

我们可以得到迭代公式：$x=x_0-\frac{f(x_0)}{f’(x_0)}$

对应到求exp上面来：$B=B_0-\frac{F(B_0)}{F’(B_0)}$

展开得到：$B=B_0-\frac{\ln(B_0)-A}{\frac{1}{B_0}}=B_0(1-\ln(B_0)+A)$

然后就可以根据上面的这个式子进行倍增，因为每一次计算长度会变成原来的两倍

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Poly exp(Poly a, int n) {
  if (n == 1) return Poly(1, 1);
  Poly x = exp(a, (n + 1) >> 1), y;
  x.resize(n);
  y = ln(x);
  for (int i = 0; i < n; i++)
    y[i] = sub(a[i], y[i]);
  y[0]++;
  x = mul(x, y);
  x.resize(n);
  return x;
}

Poly exp(Poly a) {
  return exp(a, a.size());
}

10.多项式开根

给定n次多项式$A$，求$B$满足$B^2=A(\mod x^{n+1})$

和求exp的思路比较类似，同样用到了牛顿迭代法

构造函数$F(B)=B^2-A$，求$F(B)=0$的解

带入迭代公式：$x=x_0-\frac{f(x_0)}{f’(x_0)}$

得到$B=B_0-\frac{B_0^2-A}{2B_0}=\frac{B_0^2+A}{2B_0}=\frac{1}{2}\ast (B_0+\frac{A}{B_0})$

依旧是倍增进行求解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
Poly sqrt(Poly a, int n) {
  if (n == 1) return Poly(1, 1);
  Poly x = sqrt(a, (n + 1) >> 1), y;
  x.resize(n);
  y.resize(n);
  for (int i = 0; i < n; i++)
    y[i] = a[i];
  x = add(mul(inv(x), y), x);
  int inv2 = fast_pow(2, Mod - 2);
  for (int i = 0; i < n; i++)
    x[i] = mul(x[i], inv2);
  x.resize(n);
  return x;
}

Poly sqrt(Poly a) {
  return sqrt(a, a.size());
}