多项式全家桶（半）

最新推荐文章于 2021-04-07 13:42:45 发布
最新推荐文章于 2021-04-07 13:42:45 发布
阅读量911 收藏
点赞数
分类专栏： 多项式及生成函数
版权声明：本文为博主原创文章，遵循 CC 4.0 BY 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_45483201/article/details/113747432
版权

快速傅里叶变换(FFT)

多项式表示

系数表示法：

一个 n n n次多项式可以用 n + 1 n + 1 n+1个系数表示出来： f ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a n − 1 x n − 1 + a n x n f(x) = a_0 + a_1 x + a_2 x ^ 2 + \dots + a_{n - 1} x ^{n- 1} + a_n x ^n f(x)=a0​+a1​x+a2​x2+⋯+an−1​xn−1+an​xn。

点值表示法：

通过线性代数，高斯消元我们可以知道，一个 n n n次多项式可以通过 n + 1 n + 1 n+1个点联立方程组解得：

f ( x ) = { ( x 0 , f ( x 0 ) , ( x 1 , f ( x 1 ) ) , ( x 2 , f ( x 2 ) ) ) … ( x n − 1 , f ( x n − 1 ) ) , ( x n , f ( x n ) ) } f(x) = \{(x_0, f(x_0), (x_1, f(x_1)), (x_2, f(x_2))) \dots (x_{n - 1}, f(x_{n - 1})), (x_n, f(x_n)) \} f(x)={(x0​,f(x0​),(x1​,f(x1​)),(x2​,f(x2​)))…(xn−1​,f(xn−1​)),(xn​,f(xn​))}。

有离散傅里叶变换 D F T DFT DFT（把一个多项式从系数表示变成点值表示）， I D F T IDFT IDFT（把一个多项式从点值表示变成系数表示），

而 F F T FFT FFT，就是通过选取某些特殊 x x x点，来加速 D F T , I D F T DFT, IDFT DFT,IDFT的一种方法。

点值表示法的多项式相乘：

F ( x ) = f ( x ) g ( x ) F(x) = f(x) g(x) F(x)=f(x)g(x)

F ( x ) = { ( x 0 , f ( x 0 ) g ( x 0 ) ) , ( x 1 , f ( x 1 ) g ( x 1 ) ) , ( x 2 , f ( x 2 ) g ( x 2 ) ) … ( x n − 1 , f ( x n − 1 ) g ( x n − 1 ) ) , ( x n , f ( x n ) g ( x n ) ) } F(x) = \{(x_0, f(x_0)g(x_0)), (x_1, f(x_1)g(x_1)), (x_2, f(x_2)g(x_2)) \dots (x_{n - 1}, f(x_{n - 1})g(x_{n - 1})), (x_n, f(x_n)g(x_n)) \} F(x)={(x0​,f(x0​)g(x0​)),(x1​,f(x1​)g(x1​)),(x2​,f(x2​)g(x2​))…(xn−1​,f(xn−1​)g(xn−1​)),(xn​,f(xn​)g(xn​))}

由此我们想要得到两个多项式相乘的系数，只需要先对两个多项式进行 D F T DFT DFT，然后对应的点值相乘，再做一次 I D F T IDFT IDFT，即可求得系数。

引入复数

两个复数相乘的结果为，模长相乘，辐角相加，证明如下：

有两复数 A ( a cos ⁡ θ 1 , a sin ⁡ θ 1 i ) , B ( b cos ⁡ θ 2 , b sin ⁡ θ 2 i ) A(a \cos \theta_1, a \sin \theta_1 i), B(b \cos \theta_2, b \sin \theta_2i) A(acosθ1​,asinθ1​i),B(bcosθ2​,bsinθ2​i)，用极角 + 模长来表示。

两复数相乘有 A × B = ( a b ( cos ⁡ θ 1 cos ⁡ θ 2 − sin ⁡ θ 1 sin ⁡ θ 2 ) , a b ( sin ⁡ θ 1 cos ⁡ θ 2 + cos ⁡ θ 1 sin ⁡ θ 2 ) i ) = ( a b cos ⁡ ( θ 1 + θ 2 ) , a b sin ⁡ ( θ 1 + θ 2 ) i ) A \times B = (ab(\cos \theta_1 \cos \theta_2 - \sin \theta_1 \sin \theta_2), ab(\sin \theta_1 \cos \theta_2 + \cos \theta_1 \sin \theta_2) i) = (ab \cos (\theta_1 + \theta_2), ab \sin (\theta_1 + \theta_2) i) A×B=(ab(cosθ1​cosθ2​−sinθ1​sinθ2​),ab(sinθ1​cosθ2​+cosθ1​sinθ2​)i)=(abcos(θ1​+θ2​),absin(θ1​+θ2​)i)。

引入 n n n次复根，即 x n = 1 x ^ n = 1 xn=1，这样的解显然有 n n n个，设 w n i = e 2 π n i w_{n} ^{i} = e ^{\frac{2 \pi}{n} i} wni​=en2π​i，在复平面内即是把一个圆分成了 n n n等份。

我们取 n n n等分的第一个交所对应的向量 w n = cos ⁡ ( 2 π n ) + sin ⁡ ( 2 π n ) i w_n = \cos(\frac{2 \pi}{n}) + \sin(\frac{2 \pi}{n}) i wn​=cos(n2π​)+sin(n2π​)i，则其他复根都可用 w n w_n wn​的 i i i次幂来表示。

快速傅里叶变换

考虑如何分治求解：

f ( x ) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 + a 4 x 4 + a 5 x 5 + a 6 x 6 a 7 x 7 f(x) = a_0 + a_1 x + a_2 x ^ 2 + a_3 x ^ 3 + a_4 x ^ 4 + a_5 x ^ 5 + a_6 x ^ 6 a_ 7 x ^ 7 f(x)=a0​+a1​x+a2​x2+a3​x3+a4​x4+a5​x5+a6​x6a7​x7

按照 x x x的次幂，分奇偶，再在右边提出一个 x x x，

f ( x ) = ( a 0 + a 2 x 2 + a 4 x 4 + a 6 x 6 ) + x ( a 1 + a 3 x 2 + a 5 x 4 + a 7 x 6 ) f(x) = (a_0 + a_2 x ^ 2 + a_4 x ^ 4 + a_6 x ^ 6) + x (a_1 + a_3 x ^ 2 + a_5 x ^ 4 + a_7 x ^ 6) f(x)=(a0​+a2​x2+a4​x4+a6​x6)+x(a1​+a3​x2+a5​x4+a7​x6)

G ( x ) = a 0 + a 2 x + a 4 x 2 + a 6 x 3 G(x) = a_0 + a_2 x + a_4 x ^ 2 + a_6 x ^ 3 G(x)=a0​+a2​x+a4​x2+a6​x3

H ( x ) = a 1 + a 3 x + a 5 x 2 + a 7 x 3 H(x) = a_1 + a_3 x + a_5 x ^ 2 + a_7 x ^ 3 H(x)=a1​+a3​x+a5​x2+a7​x3

有 f ( x ) = G ( x 2 ) + x H ( x 2 ) f(x) = G(x ^ 2) + x H(x ^ 2) f(x)=G(x2)+xH(x2)

由单位复根有

D F T ( f ( w n k ) ) = D F T ( G ( w n 2 k ) ) + w n k D F T ( H ( w n 2 k ) ) = D F T ( G ( w n 2 k ) ) + w n k D F T ( H ( w n 2 k ) ) DFT(f(w _n ^ k)) = DFT(G(w _n ^{2k})) + w_n ^ k DFT(H(w _n ^{2k})) = DFT(G(w _{\frac{n}{2}} ^ k)) + w_n ^ k DFT(H(w_{\frac{n}{2}} ^ k)) DFT(f(wnk​))=DFT(G(wn2k​))+wnk​DFT(H(wn2k​))=DFT(G(w2n​k​))+wnk​DFT(H(w2n​k​))

D F T ( f ( w n k + n 2 ) ) = D F T ( G ( w n 2 k ) ) − w n k D F T ( H ( w n 2 k ) ) DFT(f(w _n ^{k + \frac{n}{2}})) = DFT(G(w_{\frac{n}{2}} ^ k)) - w_n ^ k DFT(H(w _{\frac{n}{2}} ^ k)) DFT(f(wnk+2n​​))=DFT(G(w2n​k​))−wnk​DFT(H(w2n​k​))

由此，求出 D F T ( G ( w n 2 k ) ) DFT(G(w _{\frac{n}{2}} ^k)) DFT(G(w2n​k​))和 D F T ( H ( w n 2 k ) ) DFT(H(w_{\frac{n}{2}} ^ k)) DFT(H(w2n​k​))即可知 D F T ( f ( w n k ) ) , D F T ( f ( w n k + n 2 ) ) DFT(f(w _n ^ k)), DFT(f(w _n ^{k + \frac{n}{2}})) DFT(f(wnk​)),DFT(f(wnk+2n​​))然后对 G , H G, H G,H再分别递归求解即可。

快速傅里叶逆变换

把单位复根值代入多项式，得到的是如下结果：

[ y 0 y 1 y 2 ⋮ y n − 2 y n − 1 ] \left[ \begin{matrix} y_0\\ y_1\\ y_2\\ \vdots\\ y_{n - 2}\\ y_{n - 1}\\ \end{matrix} \right] ⎣⎢⎢⎢⎢⎢⎢⎢⎡​y0​y1​y2​⋮yn−2​yn−1​​⎦⎥⎥⎥⎥⎥⎥⎥⎤​ = [ 1 1 1 ⋯ 1 1 1 w n 1 w n 2 ⋯ w n n − 2 w n n − 1 1 w n 2 w n 4 ⋯ w n 2 ( n − 2 ) w n 2 ( n − 1 ) ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ 1 w n n − 2 w n 2 ( n − 2 ) ⋯ w n ( n − 2 ) ( n − 2 ) w n ( n − 1 ) ( n − 2 ) 1 w n n − 1 w n 2 ( n − 1 ) ⋯ w n ( n − 2 ) ( n − 1 ) w n ( n − 1 ) ( n − 1 ) ] \left[ \begin{matrix} 1 & 1 & 1 & \cdots & 1 & 1\\ 1 & w_n ^ 1 & w_n ^ 2 & \cdots & w_n ^ {n - 2} & w_n ^{n - 1}\\ 1 & w_n ^ 2 & w_n ^ 4 & \cdots & w_n ^ {2(n - 2)} & w_n ^{2(n - 1)}\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ 1 & w_n ^{n - 2} & w_n ^ {2(n - 2)} & \cdots & w_n ^{(n - 2)(n - 2)} & w_n ^{(n - 1)(n - 2)}\\ 1 & w_n ^{n - 1} & w_n ^ {2(n - 1)} & \cdots & w_n ^{(n - 2)(n - 1)} & w_n ^{(n - 1)(n - 1)}\\ \end{matrix} \right] ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​111⋮11​1wn1​wn2​⋮wnn−2​wnn−1​​1wn2​wn4​⋮wn2(n−2)​wn2(n−1)​​⋯⋯⋯⋱⋯⋯​1wnn−2​wn2(n−2)​⋮wn(n−2)(n−2)​wn(n−2)(n−1)​​1wnn−1​wn2(n−1)​⋮wn(n−1)(n−2)​wn(n−1)(n−1)​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​ [ a 0 a 1 a 2 ⋮ a n − 2 a n − 1 ] \left[ \begin{matrix} a_0\\ a_1\\ a_2\\ \vdots\\ a_{n - 2}\\ a_{n - 1}\\ \end{matrix} \right] ⎣⎢⎢⎢⎢⎢⎢⎢⎡​a0​a1​a2​⋮an−2​an−1​​⎦⎥⎥⎥⎥⎥⎥⎥⎤​

经过 D F T DFT DFT我们已经得到了左边的矩阵，考虑如何变换得到右边的系数矩阵，线性代数我们知道，只要在左边乘上一个中间大矩阵的逆，我们即可得到右边的系数矩阵。

由于这个矩阵的元素非常特殊，他的逆矩阵也有特殊的性质，就是每一项取倒数，再除以 n n n，就能得到他的逆矩阵。

每一项取倒数有 1 w n = w n − 1 = e − 2 π i n = cos ⁡ ( 2 π n ) + i sin ⁡ ( − 2 π n ) \frac{1}{w_n} = w_{n} ^{-1} = e ^{-\frac{2 \pi i}{n}} = \cos(\frac{2 \pi}{n}) + i \sin (- \frac{2 \pi}{n}) wn​1​=wn−1​=e−n2πi​=cos(n2π​)+isin(−n2π​)，所以我们只要将这个代入做一次 D F T DFT DFT，也就是 I D F T IDFT IDFT，最后再对整体除以 n n n即可得到系数矩阵。

对以上进行证明

f ( x ) = ∑ i = 0 n − 1 a i x i f(x) = \sum\limits_{i = 0} ^{n - 1} a_i x ^ i f(x)=i=0∑n−1​ai​xi， y i = f ( w n i ) y_i = f(w_n ^ i) yi​=f(wni​)，构造 A ( x ) = ∑ i = 0 n − 1 y i x i A(x) = \sum\limits_{i = 0} ^{n - 1}y_i x ^ i A(x)=i=0∑n−1​yi​xi，将 b i = w n − i b_i = w_{n} ^{-i} bi​=wn−i​代入多项式 A ( x ) A(x) A(x)

有 A ( b k ) = ∑ i = 0 n − 1 y i w n − i k = ∑ i = 0 n 1 w n − i k ∑ j = 0 n − 1 a j w n i j = ∑ j = 0 n − 1 a j ∑ i = 0 n − 1 ( w n j − k ) i A(b_k) = \sum\limits_{i = 0} ^{n - 1} y_i w_n ^{-ik} = \sum\limits_{i = 0} ^{n 1}w_n ^{-ik} \sum\limits_{j = 0} ^{n - 1} a_j w_{n} ^{ij} = \sum\limits_{j = 0} ^{n - 1} a_j\sum\limits_{i = 0} ^{n - 1} (w_{n} ^{j - k}) ^ i A(bk​)=i=0∑n−1​yi​wn−ik​=i=0∑n1​wn−ik​j=0∑n−1​aj​wnij​=j=0∑n−1​aj​i=0∑n−1​(wnj−k​)i

令 S ( w n a ) = ∑ i = 0 n − 1 ( w n a ) i S(w_{n} ^ a) = \sum\limits _{i = 0} ^{n - 1} (w_{n} ^{a}) ^ i S(wna​)=i=0∑n−1​(wna​)i

显然有 a = 0 a = 0 a=0， S ( w n a ) = n S(w_n ^ a) = n S(wna​)=n

a ≠ 0 a \neq 0 a​=0，时我们取 S ( w n a ) , w n a S ( w n a ) S(w_n ^ a),w_n ^ a S(w_n ^ a) S(wna​),wna​S(wna​)，两者相减，除以一个系数有 S ( w n a ) = ∑ i = 1 n ( w n a ) i − ∑ i = 0 n − 1 ( w n a ) i w n a − 1 = ( w n a ) n − ( w n a ) 0 w n a − 1 = 0 S(w_n ^ a) = \frac{\sum\limits_{i = 1} ^{n} (w_n ^ a) ^ i - \sum\limits_{i = 0} ^{n - 1} (w_n ^ a) ^ i}{w_n ^ a - 1} = \frac{(w_n ^ a) ^ n - (w_n ^ a) ^ 0}{w_n ^ a - 1} = 0 S(wna​)=wna​−1i=1∑n​(wna​)i−i=0∑n−1​(wna​)i​=wna​−1(wna​)n−(wna​)0​=0

所以有 S ( w n a ) = [ a = 1 ] S(w_n ^ a) = [a = 1] S(wna​)=[a=1]

A ( b k ) = a k × n A(b_k) = a_k \times n A(bk​)=ak​×n

仔细想想，这个证明，就是把我们 D F T DFT DFT过程中得到的点值作为系数去做一遍 D F T DFT DFT，得到的也就是 A ( x ) A(x) A(x)的点值表达式，同时对其除以 n n n，也就是 f ( x ) f(x) f(x)的系数表达式了。

如何优化（蝴蝶变换）

分治过程中考虑系数如何变换
{ a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 } { a 0 , a 2 , a 4 , a 6 } { a 1 , a 3 , a 5 , a 7 } { a 0 , a 4 } { a 2 , a 6 } { a 1 , a 5 } { a 3 , a 7 } { a 0 } { a 4 } { a 2 } { a 6 } { a 1 } { a 5 } { a 3 } { a 7 } \{a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7\}\\ \{a_0, a_2, a_4, a_6\}\{a_1, a_3, a_5, a_7\}\\ \{a_0, a_4\}\{a_2, a_6\}\{a_1, a_5\}\{a_3, a_7\}\\ \{a_0\}\{a_4\}\{a_2\}\{a_6\}\{a_1\}\{a_5\}\{a_3\}\{a_7\}\\ {a0​,a1​,a2​,a3​,a4​,a5​,a6​,a7​}{a0​,a2​,a4​,a6​}{a1​,a3​,a5​,a7​}{a0​,a4​}{a2​,a6​}{a1​,a5​}{a3​,a7​}{a0​}{a4​}{a2​}{a6​}{a1​}{a5​}{a3​}{a7​}

这个过程中有一个规律，例如 1 = 001 1 = 001 1=001，倒置后变成了 100 100 100， 4 4 4，也即是最后 a 1 a_1 a1​所在的位置。

r [ i ] r[i] r[i]表示 i i i翻转之后的数字，考虑如何从小到大递推得到 r [ i ] r[i] r[i]，有 r [ 0 ] = 0 r[0] = 0 r[0]=0，当我们在求 x x x时，先考虑除个位数以外的数，就是 r [ x > > 1 ] > > 1 r[x >> 1] >> 1 r[x>>1]>>1了，如果个位是 1 1 1则加上 l i m > > 1 lim >> 1 lim>>1，就有了如下代码

void change(int lim) {
    for (int i = 0; i < lim; i++) {
        r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
    }
}

真正可用的 F F T FFT FFT代码

P3803 【模板】多项式乘法（FFT）

#include <bits/stdc++.h>

using namespace std;

struct Complex {
  double r, i;

  Complex(double _r = 0, double _i = 0) : r(_r), i(_i) {}
};

Complex operator + (const Complex &a, const Complex &b) {
  return Complex(a.r + b.r, a.i + b.i);
}

Complex operator - (const Complex &a, const Complex &b) {
  return Complex(a.r - b.r, a.i - b.i);
}

Complex operator * (const Complex &a, const Complex &b) {
  return Complex(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r);
}

Complex operator / (const Complex &a, const Complex &b) {
  return Complex((a.r * b.r + a.i * b.i) / (b.r * b.r + b.i * b.i), (a.i * b.r - a.r * b.i) / (b.r * b.r + b.i * b.i));
}

typedef long long ll;

const int N = 5e6 + 10;

int r[N], n, m;

Complex a[N], b[N];

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  }
}

void FFT(Complex *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  const double pi = acos(-1.0);
  for (int mid = 1; mid < lim; mid <<= 1) {
    Complex wn = Complex(cos(pi / mid), rev * sin(pi / mid));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      Complex w = Complex(1, 0);
      for (int k = 0; k < mid; k++, w = w * wn) {
        Complex x = f[cur + k], y = w * f[cur + mid + k];
        f[cur + k] = x + y, f[cur + mid + k] = x - y;
      }
    }
  }
  if (rev == -1) {
    for (int i = 0; i < lim; i++) {
      a[i].r /= lim;
    }
  }
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  scanf("%d %d", &n, &m);
  n += 1, m += 1;
  for (int i = 0; i < n; i++) {
    scanf("%lf", &a[i].r);
  }
  for (int i = 0; i < m; i++) {
    scanf("%lf", &b[i].r);
  }
  int lim = 1;
  while (lim <= n + m) {
    lim <<= 1;
  }
  get_r(lim);
  FFT(a, lim, 1);
  FFT(b, lim, 1);
  for (int i = 0; i < lim; i++) {
    a[i] = a[i] * b[i];
  }
  FFT(a, lim, -1);
  for (int i = 0; i < n + m - 1; i++) {
    printf("%lld ", ll(a[i].r + 0.5));
  }
  puts("");
  return 0;
}

快速数论变换(NTT)

#include <bits/stdc++.h>

using namespace std;

const int N = 5e6 + 10, mod = 998244353;

int a[N], b[N], r[N], n, m;

int quick_pow(int a, int n) {
  int ans = 1;
  while (n) {
    if (n & 1) {
      ans = 1ll * ans * n % mod;
    }
    a = 1ll * a * a % mod;
    n >>= 1;
  }
  return ans;
}

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  }
}

void NTT(int *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  for (int mid = 1; mid < lim; mid <<= 1) {
    int wn = quick_pow(3, (mod - 1) / (mid << 1));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      int w = 1;
      for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {
        int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;
        f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if (rev == -1) {
    int inv = quick_pow(lim, mod - 2);
    reverse(f + 1, f + lim);
    for (int i = 0; i < lim; i++) {
      f[i] = 1ll * f[i] * inv % mod;
    }
  }
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  return 0;
}

多项式求逆

f ( x ) g ( x ) ≡ 1 ( m o d x n ) f(x) g(x) \equiv 1 \pmod {x ^ n} f(x)g(x)≡1(modxn)，称 f ( x ) f(x) f(x)为 g ( x ) g(x) g(x)或者 g ( x ) g(x) g(x)为 f ( x ) f(x) f(x)膜 x n x ^ n xn意义下的逆元。

下面我们讨论给定 f ( x ) f(x) f(x)，求其逆 f − 1 ( x ) f ^{-1}(x) f−1(x)。

倍增求解

假设我们已经求得 f ( x ) f(x) f(x)膜 x ⌈ n 2 ⌉ x ^{\lceil \frac{n}{2}} \rceil x⌈2n​⌉下的逆元 f 0 − 1 ( x ) f_0 ^{-1} (x) f0−1​(x)，要求 f − 1 ( x ) f ^{-1}(x) f−1(x)，即膜 x n x ^{n} xn下的逆元，则

f ( x ) f 0 − 1 ( x ) ≡ 1 ( m o d x ⌈ n 2 ⌉ ) f(x) f_0 ^{-1}(x) \equiv 1 \pmod{x ^{\lceil\frac{n}{2}\rceil} } f(x)f0−1​(x)≡1(modx⌈2n​⌉)

显然 f ( x ) f − 1 ( x ) ≡ 1 ( m o d x ⌈ n 2 ⌉ ) f(x) f^{-1}(x) \equiv 1 \pmod {x ^{ \lceil\frac{n}{2}\rceil}} f(x)f−1(x)≡1(modx⌈2n​⌉)也是成立的

对两边同时乘以 f 0 − 1 ( x ) f_0 ^{-1}(x) f0−1​(x)并移项有

f − 1 ( x ) − f 0 − 1 ( x ) ≡ 0 ( m o d x ⌈ n 2 ⌉ ) f ^{-1}(x) - f_0 ^{-1}(x) \equiv 0 \pmod{x ^{\lceil\frac{n}{2}\rceil}} f−1(x)−f0−1​(x)≡0(modx⌈2n​⌉)

对两边同时开方得到

f − 2 ( x ) − 2 f − 1 f 0 − 1 ( x ) + f 0 − 2 ( x ) ≡ 0 ( m o d x n ) f ^{-2}(x) - 2 f^{-1} f_0 ^{-1}(x) + f_0 ^{-2}(x) \equiv 0 \pmod {x ^n} f−2(x)−2f−1f0−1​(x)+f0−2​(x)≡0(modxn)

我们再对两边乘上一个 f ( x ) f(x) f(x)，则有

f − 1 ( x ) − 2 f 0 − 1 + f ( x ) f 0 − 2 ( x ) ≡ 0 ( m o d x n ) f ^{-1}(x) - 2 f_0 ^{-1} + f(x) f_0 ^{-2}(x) \equiv 0 \pmod{x ^n} f−1(x)−2f0−1​+f(x)f0−2​(x)≡0(modxn)

再对其进行移项可得

f − 1 ( x ) ≡ f 0 − 1 ( x ) ( 2 − f ( x ) f 0 − 1 ( x ) ) ( m o d x n ) f ^{-1}(x) \equiv f_0 ^{-1}(x)\left( 2 - f(x) f_0 ^{-1}(x) \right) \pmod {x ^n} f−1(x)≡f0−1​(x)(2−f(x)f0−1​(x))(modxn)

由此我们递归求解即可。

#include <bits/stdc++.h>

using namespace std;

const int N = 1e6 + 10, mod = 998244353;

int a[N], b[N], c[N], r[N], n;

int quick_pow(int a, int n) {
  int ans = 1;
  while (n) {
    if (n & 1) {
      ans = 1ll * ans * a % mod;
    }
    a = 1ll * a * a % mod;
    n >>= 1;
  }
  return ans;
}

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = ((i & 1) * (lim >> 1)) + (r[i >> 1] >> 1);
  }
}

void NTT(int *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  for (int mid = 1; mid < lim; mid <<= 1) {
    int wn = quick_pow(3, (mod - 1) / (mid << 1));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      int w = 1;
      for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {
        int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;
        f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if (rev == -1) {
    int inv = quick_pow(lim, mod - 2);
    reverse(f + 1, f + lim);
    for (int i = 0; i < lim; i++) {
      f[i] = 1ll * inv * f[i] % mod;
    }
  }
}

void polyinv(int *a, int *b, int n) {
  if (n == 1) {
    b[0] = quick_pow(a[0], mod - 2);
    return ;
  }
  polyinv(a, b, n + 1 >> 1);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  for (int i = 0; i < n; i++) {
    c[i] = a[i];
  }
  for (int i = n; i < lim; i++) {
    c[i] = 0;
  }
  NTT(b, lim, 1);
  NTT(c, lim, 1);
  for (int i = 0; i < lim; i++) {
    int cur = (2 - 1ll * c[i] * b[i] % mod + mod) % mod;
    b[i] = 1ll * b[i] * cur % mod;
  }
  NTT(b, lim, -1);
  for (int i = n; i < lim; i++) {
    b[i] = 0;
  }
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  scanf("%d", &n);
  for (int i = 0; i < n; i++) {
    scanf("%d", &a[i]);
  }
  polyinv(a, b, n);
  for (int i = 0; i < n; i++) {
    printf("%d%c", b[i], i + 1 == n ? '\n' : ' ');
  }
  return 0;
}

多项式开根

给定多项式 g ( x ) g(x) g(x)，求 f ( x ) f(x) f(x)，满足 f 2 ( x ) = g ( x ) f ^ 2(x) = g(x) f2(x)=g(x)。

假设我们已经得到了 g ( x ) g(x) g(x)，膜 x ⌈ n 2 ⌉ x ^{\lceil \frac{n}{2} \rceil} x⌈2n​⌉下的根 f 0 ( x ) f_0 (x) f0​(x)，要求膜 x n x ^ n xn下的根 f ( x ) f(x) f(x)

有 f 0 2 ( x ) ≡ g ( x ) ( m o d x ⌈ n 2 ⌉ ) f_0 ^2(x) \equiv g(x) \pmod {x ^{\lceil \frac{n}{2} \rceil}} f02​(x)≡g(x)(modx⌈2n​⌉)

移项再开方有 ( f 0 2 ( x ) − g ( x ) ) 2 ≡ 0 ( m o d x n ) \left(f_0 ^2(x) - g(x) \right) ^ 2 \equiv 0 \pmod {x ^ n} (f02​(x)−g(x))2≡0(modxn)

则， ( f 0 2 ( x ) + g ( x ) ) 2 ≡ 4 f 0 2 ( x ) g ( x ) ( m o d x n ) \left( f_0 ^ 2(x) + g(x) \right) ^ 2 \equiv 4 f_0 ^ 2(x) g(x) \pmod {x ^ n} (f02​(x)+g(x))2≡4f02​(x)g(x)(modxn)

g ( x ) ≡ ( f 0 2 ( x ) + g ( x ) 2 f 0 ( x ) ) 2 ( m o d x n ) g(x) \equiv \left(\frac{f_0 ^ 2(x) + g(x)}{2f_0 (x)} \right) ^ 2 \pmod {x ^ n} g(x)≡(2f0​(x)f02​(x)+g(x)​)2(modxn)

所以 f ( x ) ≡ f 0 2 ( x ) + g ( x ) 2 f 0 ( x ) ( m o d x n ) f(x) \equiv \frac{f_0 ^ 2(x) + g(x)}{2f_0(x)} \pmod {x ^ n} f(x)≡2f0​(x)f02​(x)+g(x)​(modxn)。

所以有 f ( x ) ≡ 2 − 1 f 0 ( x ) + 2 − 1 f 0 − 1 ( x ) g ( x ) ( m o d x n ) f(x) \equiv 2 ^{-1} f_0 (x) + 2 ^{-1} f_0 ^{-1}(x) g(x) \pmod {x ^ n} f(x)≡2−1f0​(x)+2−1f0−1​(x)g(x)(modxn)，

对于 g ( 0 ) = 1 g(0) = 1 g(0)=1的特殊情况

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int N = 5e6 + 10, mod = 998244353, inv2 = mod + 1 >> 1;

int a[N], b[N], c[N], d[N], r[N];

int quick_pow(int a, int n) {
  int ans = 1;
  while (n) {
    if (n & 1) {
      ans = 1ll * ans * a % mod;
    }
    a = 1ll *  a * a % mod;
    n >>= 1;
  }
  return ans;
}

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  }
}

void NTT(int *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  for (int mid = 1; mid < lim; mid <<= 1) {
    int wn = quick_pow(3, (mod - 1) / (mid << 1));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      int w = 1;
      for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {
        int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;
        f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if (rev == -1) {
    int inv = quick_pow(lim, mod - 2);
    reverse(f + 1, f + lim);
    for (int i = 0; i < lim; i++) {
      f[i] = 1ll * f[i] * inv % mod;
    }
  }
}

void polyinv1(int *a, int *b, int n) {
  if (n == 1) {
    b[0] = quick_pow(a[0], mod - 2);
    return ;
  }
  polyinv1(a, b, n + 1 >> 1);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  for (int i = 0; i < n; i++) {
    c[i] = a[i];
  }
  for (int i = n; i < lim; i++) {
    c[i] = 0;
  }
  NTT(b, lim, 1);
  NTT(c, lim, 1);
  for (int i = 0; i < lim; i++) {
    int cur = (2 - 1ll * c[i] * b[i] % mod + mod) % mod;
    b[i] = 1ll * b[i] * cur % mod;
  }
  NTT(b, lim, -1);
  for (int i = n; i < lim; i++) {
    b[i] = 0;
  }
}

void polysqrt(int *a, int *b, int n) {
  if (n == 1) {
    b[0] = 1;
    return ;
  }
  polysqrt(a, b, n + 1 >> 1);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  for (int i = 0; i < lim; i++) {
    d[i] = 0;
  }
  polyinv1(b, d, n);
  for (int i = 0; i < n; i++) {
    c[i] = a[i];
  }
  for (int i = n; i < lim; i++) {
    c[i] = 0;
  }
  get_r(lim);
  NTT(b, lim, 1);
  NTT(c, lim, 1);
  NTT(d, lim, 1);
  for (int i = 0; i < lim; i++) {
    b[i] = (1ll * inv2 * b[i] % mod + 1ll * inv2 * d[i] % mod * c[i] % mod) % mod; 
  }
  NTT(b, lim, -1);
  for (int i = n; i < lim; i++) {
    b[i] = 0;
  }
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  int n;
  scanf("%d", &n);
  for (int i = 0; i < n; i++) {
    scanf("%d", &a[i]);
  }
  polysqrt(a, b, n);
  for (int i = 0; i < n; i++) {
    printf("%d%c", b[i], i + 1 == n ? '\n' : ' ');
  }
  return 0;
}

二次剩余解一般情况

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int N = 5e6 + 10, mod = 998244353, inv2 = mod + 1 >> 1;

int a[N], b[N], c[N], d[N], r[N];

namespace Quadratic_residue {
  struct Complex {
    int r, i;

    Complex(int _r = 0, int _i = 0) : r(_r), i(_i) {}
  };

  int I2;

  Complex operator * (const Complex &a, Complex &b) {
    return Complex((1ll * a.r * b.r % mod  + 1ll * a.i * b.i % mod * I2 % mod) % mod, (1ll * a.r * b.i % mod + 1ll * a.i * b.r % mod) % mod);
  }

  Complex quick_pow(Complex a, int n) {
    Complex ans = Complex(1, 0);
    while (n) {
      if (n & 1) {
        ans = ans * a;
      }
      a = a * a;
      n >>= 1;
    }
    return ans;
  }

  int get_residue(int n) {
    mt19937 e(233);
    if (n == 0) {
      return 0;
    }
    if(quick_pow(n, (mod - 1) >> 1).r == mod - 1) {
      return -1;
    }
    uniform_int_distribution<int> r(0, mod - 1);
    int a = r(e);
    while(quick_pow((1ll * a * a % mod - n + mod) % mod, (mod - 1) >> 1).r == 1) {
      a = r(e);
    }
    I2 = (1ll * a * a % mod - n + mod) % mod;
    int x = quick_pow(Complex(a, 1), (mod + 1) >> 1).r, y = mod - x;
    if(x > y) swap(x, y);
    return x;
  }
}

int quick_pow(int a, int n) {
  int ans = 1;
  while (n) {
    if (n & 1) {
      ans = 1ll * ans * a % mod;
    }
    a = 1ll *  a * a % mod;
    n >>= 1;
  }
  return ans;
}

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  }
}

void NTT(int *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  for (int mid = 1; mid < lim; mid <<= 1) {
    int wn = quick_pow(3, (mod - 1) / (mid << 1));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      int w = 1;
      for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {
        int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;
        f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if (rev == -1) {
    int inv = quick_pow(lim, mod - 2);
    reverse(f + 1, f + lim);
    for (int i = 0; i < lim; i++) {
      f[i] = 1ll * f[i] * inv % mod;
    }
  }
}

void polyinv1(int *a, int *b, int n) {
  if (n == 1) {
    b[0] = quick_pow(a[0], mod - 2);
    return ;
  }
  polyinv1(a, b, n + 1 >> 1);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  for (int i = 0; i < n; i++) {
    c[i] = a[i];
  }
  for (int i = n; i < lim; i++) {
    c[i] = 0;
  }
  NTT(b, lim, 1);
  NTT(c, lim, 1);
  for (int i = 0; i < lim; i++) {
    int cur = (2 - 1ll * c[i] * b[i] % mod + mod) % mod;
    b[i] = 1ll * b[i] * cur % mod;
  }
  NTT(b, lim, -1);
  for (int i = n; i < lim; i++) {
    b[i] = 0;
  }
}

void polysqrt(int *a, int *b, int n) {
  if (n == 1) {
    b[0] = Quadratic_residue::get_residue(a[0]);
    return ;
  }
  polysqrt(a, b, n + 1 >> 1);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  for (int i = 0; i < lim; i++) {
    d[i] = 0;
  }
  polyinv1(b, d, n);
  for (int i = 0; i < n; i++) {
    c[i] = a[i];
  }
  for (int i = n; i < lim; i++) {
    c[i] = 0;
  }
  get_r(lim);
  NTT(b, lim, 1);
  NTT(c, lim, 1);
  NTT(d, lim, 1);
  for (int i = 0; i < lim; i++) {
    b[i] = (1ll * inv2 * b[i] % mod + 1ll * inv2 * d[i] % mod * c[i] % mod) % mod; 
  }
  NTT(b, lim, -1);
  for (int i = n; i < lim; i++) {
    b[i] = 0;
  }
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  int n;
  scanf("%d", &n);
  for (int i = 0; i < n; i++) {
    scanf("%d", &a[i]);
  }
  polysqrt(a, b, n);
  for (int i = 0; i < n; i++) {
    printf("%d%c", b[i], i + 1 == n ? '\n' : ' ');
  }
  return 0;
}

分治FFT

考虑计算这么一个式子 f ( i ) = ∑ j = 1 i f i − j g ( j ) f(i) = \sum\limits_{j = 1} ^{i} f_{i - j}g(j) f(i)=j=1∑i​fi−j​g(j)，给定 g ( x ) g(x) g(x)，求 f ( x ) f(x) f(x)，边界条件 f ( 0 ) = 1 f(0) = 1 f(0)=1。

假设我们已经算出 [ l , m i d ] [l, mid] [l,mid]，考虑计算其对 [ m i d + 1 , r ] [mid + 1, r] [mid+1,r]的贡献 w ( i ) w(i) w(i)，

有 w ( x ) = ∑ i = l m i d f ( i ) g ( x − i ) w(x) = \sum\limits_{i = l} ^{mid} f(i) g(x - i) w(x)=i=l∑mid​f(i)g(x−i)，因为 [ m i d + 1 , r ] [mid + 1, r] [mid+1,r]区间还没开始计算，所以 f ( i ) = 0 , i ∈ [ m i d , x − 1 ] f(i) = 0, i \in [mid, x - 1] f(i)=0,i∈[mid,x−1]，则 w ( x ) = ∑ i = l x − 1 f ( i ) g ( x − i ) w(x) = \sum\limits_{i = l} ^{x - 1} f(i) g(x - i) w(x)=i=l∑x−1​f(i)g(x−i)

我们设 a ( i ) = f ( i + l ) , b ( i ) = g ( i + 1 ) a(i) = f(i + l), b(i) = g(i + 1) a(i)=f(i+l),b(i)=g(i+1)，上式 w ( x ) = ∑ i = 0 x − l − 1 a ( i ) b ( x − l − i + 1 ) w(x) = \sum\limits_{i = 0} ^{x - l - 1} a(i)b(x - l - i + 1) w(x)=i=0∑x−l−1​a(i)b(x−l−i+1)

有 w ( x ) = c ( x − l + 1 ) = ∑ i = 0 x − l − 1 a ( i ) b ( x − l − i + 1 ) w(x) = c(x - l + 1) = \sum\limits_{i = 0} ^{x - l - 1}a(i) b(x - l - i + 1) w(x)=c(x−l+1)=i=0∑x−l−1​a(i)b(x−l−i+1)

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int N = 5e6 + 10, mod = 998244353, inv2 = mod + 1 >> 1;

int a[N], b[N], c[N], d[N], r[N];

int quick_pow(int a, int n) {
  int ans = 1;
  while (n) {
    if (n & 1) {
      ans = 1ll * ans * a % mod;
    }
    a = 1ll *  a * a % mod;
    n >>= 1;
  }
  return ans;
}

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  }
}

void NTT(int *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  for (int mid = 1; mid < lim; mid <<= 1) {
    int wn = quick_pow(3, (mod - 1) / (mid << 1));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      int w = 1;
      for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {
        int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;
        f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if (rev == -1) {
    int inv = quick_pow(lim, mod - 2);
    reverse(f + 1, f + lim);
    for (int i = 0; i < lim; i++) {
      f[i] = 1ll * f[i] * inv % mod;
    }
  }
}

void DACFFT(int l, int r) {
  if (l == r) {
    return ;
  }
  int mid = l + r >> 1;
  DACFFT(l, mid);
  for (int i = 0; i <= mid - l; i++) {
    c[i] = b[i + l];
  }
  for (int i = 0; i < r - l; i++) {
    d[i] = a[i + 1];
  }
  int lim = 1;
  while (lim <= r - l + mid - l + 1) {
    lim <<= 1;
  }
  get_r(lim);
  NTT(c, lim, 1);
  NTT(d, lim, 1);
  for (int i = 0; i < lim; i++) {
    c[i] = 1ll * c[i] * d[i] % mod;
  }
  NTT(c, lim, -1);
  for (int i = mid - l; i <= r - l - 1; i++) {
    b[i + l + 1] = (b[i + l + 1] + c[i]) % mod;
  }
  for (int i = 0; i < lim; i++) {
    c[i] = d[i] = 0;
  }
  DACFFT(mid + 1, r);
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  int n;
  scanf("%d", &n);
  for (int i = 1; i < n; i++) {
    scanf("%d", &a[i]);
  }
  b[0] = 1;
  DACFFT(0, n - 1);
  for (int i = 0; i < n; i++) {
    printf("%d%c", b[i], i == n ? '\n' : ' ');
  }
  return 0;
}

牛顿迭代

迭代求函数零点

对于 f ( x ) f(x) f(x)，任意选取一个点 ( x 0 , f ( x 0 ) ) (x_0, f(x_0)) (x0​,f(x0​))作为当前我们预估的零点，取他的泰勒展开前两项 g ( x ) = f ( 0 ) + f ′ ( x ) ( x − x 0 ) g(x) = f(0) + f'(x)(x - x_0) g(x)=f(0)+f′(x)(x−x0​)，

解出方程 g ( x ) = 0 g(x) = 0 g(x)=0，得到 x 1 x_1 x1​，然后重复上述操作，最后 x n x_n xn​会趋近于我们所要的正解。

考虑如何应用到多项式上

边界条件 n = 1 n = 1 n=1时， [ x 0 ] g ( f ( x ) ) = 0 [x ^0]g(f(x)) = 0 [x0]g(f(x))=0，的解单独求出。

假设我们已经求得膜 x ⌈ n 2 ⌉ x ^{\lceil\frac{n}{2} \rceil} x⌈2n​⌉下的解 f 0 ( x ) f_0(x) f0​(x)，要求膜 x n x ^n xn下的解 f ( x ) f(x) f(x)，得到该点的泰勒展开：
∑ i = 0 ∞ g ( i ) ( f 0 ( x ) ) i ! ( f ( x ) − f 0 ( x ) ) n f ( x ) − f 0 ( x ) 前 面 的 项 已 经 被 截 了 , 所 以 最 低 次 幂 是 大 于 ⌈ n 2 ⌉ 的 有 ( f ( x ) − f 0 ( x ) ) i ≡ 0 ( m o d x n ) , i > = 2 有 g ( f 0 ( x ) ) + g ′ ( f 0 ( x ) ) ( f ( x ) − f 0 ( x ) ) ≡ 0 ( m o d x n ) f ( x ) ≡ f 0 ( x ) − g ( f 0 ( x ) ) g ′ ( f 0 ( x ) ) ( m o d x n ) \sum_{i = 0} ^{\infty} \frac{g ^{(i)}(f_0(x))}{i !}(f(x) - f_0(x)) ^n\\ f(x) - f_0(x)前面的项已经被截了,所以最低次幂是大于\lceil \frac{n}{2} \rceil的\\ 有(f(x) - f_0 (x)) ^ i \equiv 0 \pmod {x ^ n}, i >= 2\\ 有g(f_0(x)) + g'(f_0(x))(f(x) - f_0(x)) \equiv 0 \pmod {x ^ n}\\ f(x) \equiv f_0(x) - \frac{g(f_0(x))}{g'(f_0(x))} \pmod {x ^ n}\\ i=0∑∞​i!g(i)(f0​(x))​(f(x)−f0​(x))nf(x)−f0​(x)前面的项已经被截了,所以最低次幂是大于⌈2n​⌉的有(f(x)−f0​(x))i≡0(modxn),i>=2有g(f0​(x))+g′(f0​(x))(f(x)−f0​(x))≡0(modxn)f(x)≡f0​(x)−g′(f0​(x))g(f0​(x))​(modxn)

应用

多项式求逆

对 于 给 定 的 h ( x ) ， 有 g ( f ( x ) ) = 1 f ( x ) − h ( x ) ≡ 0 ( m o d x n ) f ( x ) ≡ f 0 ( x ) − 1 f 0 ( x ) − h ( x ) − 1 f 0 2 ( x ) ( m o d x n ) f ( x ) ≡ 2 f 0 ( x ) − h ( x ) f 0 2 ( x ) ( m o d x n ) f ( x ) ≡ f 0 ( x ) ( 2 − h ( x ) f 0 ( x ) ) ( m o d x n ) 对于给定的h(x)，\\ 有g(f(x)) = \frac{1}{f(x)} - h(x) \equiv 0 \pmod {x ^ n}\\ f(x) \equiv f_0(x) - \frac{\frac{1}{f_0(x)} - h(x)}{-\frac{1}{f_0 ^ 2(x)}} \pmod{x ^n}\\ f(x) \equiv 2f_0(x) - h(x) f_0 ^ 2(x) \pmod {x ^ n}\\ f(x) \equiv f_0(x)(2 - h(x)f_0(x)) \pmod {x ^ n}\\ 对于给定的h(x)，有g(f(x))=f(x)1​−h(x)≡0(modxn)f(x)≡f0​(x)−−f02​(x)1​f0​(x)1​−h(x)​(modxn)f(x)≡2f0​(x)−h(x)f02​(x)(modxn)f(x)≡f0​(x)(2−h(x)f0​(x))(modxn)

多项式开根

对 于 给 定 的 h ( x ) 有 g ( f ( x ) ) = f 2 ( x ) − h ( x ) ≡ 0 ( m o d x n ) f ( x ) ≡ f 0 ( x ) − f 0 2 ( x ) − h ( x ) 2 f 0 ( x ) ( m o d ( ) x n ) f ( x ) ≡ f 0 2 ( x ) + h ( x ) 2 f 0 ( x ) ( m o d x n ) f ( x ) ≡ 2 − 1 f 0 ( x ) + 2 − 1 f 0 − 1 ( x ) h ( x ) ( m o d x n ) 对于给定的h(x)\\ 有g(f(x)) = f ^ 2(x) - h(x) \equiv 0 \pmod {x ^ n}\\ f(x) \equiv f_0(x) - \frac{f_0 ^2(x) - h(x)}{2f_0(x)} \pmod (x ^n)\\ f(x) \equiv \frac{f_0 ^ 2(x) + h(x)}{2f_0(x)} \pmod {x ^ n}\\ f(x) \equiv 2 ^{-1} f_0 (x) + 2 ^{-1} f_0 ^{-1}(x) h(x) \pmod {x ^n}\\ 对于给定的h(x)有g(f(x))=f2(x)−h(x)≡0(modxn)f(x)≡f0​(x)−2f0​(x)f02​(x)−h(x)​(mod()xn)f(x)≡2f0​(x)f02​(x)+h(x)​(modxn)f(x)≡2−1f0​(x)+2−1f0−1​(x)h(x)(modxn)

多项式 exp ⁡ \exp exp

对 于 给 定 的 h ( x ) 有 g ( f ( x ) ) = ln ⁡ f ( x ) − h ( x ) ≡ 0 ( m o d x n ) f ( x ) ≡ f 0 ( x ) − ln ⁡ f 0 ( x ) − h ( x ) 1 f 0 ( x ) ( m o d x n ) f ( x ) ≡ f 0 ( x ) ( 1 − ln ⁡ f 0 ( x ) + h ( x ) ) ( m o d x n ) 对于给定的h(x)\\ 有g(f(x)) = \ln f(x) - h(x) \equiv 0 \pmod {x ^ n}\\ f(x) \equiv f_0(x) - \frac{\ln f_0(x) - h(x)}{\frac{1}{f_0(x)}} \pmod {x ^ n}\\ f(x) \equiv f_0(x)(1 - \ln f_0(x) + h(x)) \pmod {x ^ n}\\ 对于给定的h(x)有g(f(x))=lnf(x)−h(x)≡0(modxn)f(x)≡f0​(x)−f0​(x)1​lnf0​(x)−h(x)​(modxn)f(x)≡f0​(x)(1−lnf0​(x)+h(x))(modxn)

多项式对数函数 ln ⁡ f ( x ) \ln f(x) lnf(x)

如果存在解必然有 [ x 0 ] f ( x ) = 1 [ x ^ 0]f(x) = 1 [x0]f(x)=1，

对 ln ⁡ f ( x ) \ln f(x) lnf(x)求导，有 d ln ⁡ f ( x ) d x ≡ f ′ ( x ) f ( x ) ( m o d x n ) \frac{d \ln f(x)}{dx} \equiv \frac{f'(x)}{f(x)} \pmod {x ^ n} dxdlnf(x)​≡f(x)f′(x)​(modxn)，

d x dx dx乘到右边，再求积分有：
∫ d ln ⁡ f ( x ) ≡ ∫ f ′ ( x ) ) f ( x ) d x ( m o d x n ) ln ⁡ f ( x ) ≡ ∫ f ′ ( x ) f ( x ) ( m o d x n ) \int d \ln f(x) \equiv \int \frac{f'(x))}{f(x)} dx \pmod {x ^ n}\\ \ln f(x) \equiv \int \frac{f'(x)}{f(x)} \pmod {x ^ n}\\ ∫dlnf(x)≡∫f(x)f′(x))​dx(modxn)lnf(x)≡∫f(x)f′(x)​(modxn)
然后只要对先对 f ( x ) f(x) f(x)求个导，求个逆，最后求一次积分即可，整体复杂度 O ( n log ⁡ n ) O(n \log n) O(nlogn)。

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int N = 5e6 + 10, mod = 998244353, inv2 = mod + 1 >> 1;

int a[N], b[N], c[N], d[N], r[N], inv[N];

int quick_pow(int a, int n) {
  int ans = 1;
  while (n) {
    if (n & 1) {
      ans = 1ll * ans * a % mod;
    }
    a = 1ll *  a * a % mod;
    n >>= 1;
  }
  return ans;
}

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  }
}

void get_inv(int n) {
  inv[1] = 1;
  for (int i = 2; i < n; i++) {
    inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
  }
}

void NTT(int *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  for (int mid = 1; mid < lim; mid <<= 1) {
    int wn = quick_pow(3, (mod - 1) / (mid << 1));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      int w = 1;
      for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {
        int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;
        f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if (rev == -1) {
    int inv = quick_pow(lim, mod - 2);
    reverse(f + 1, f + lim);
    for (int i = 0; i < lim; i++) {
      f[i] = 1ll * f[i] * inv % mod;
    }
  }
}

void polyinv(int *a, int *b, int n) {
  if (n == 1) {
    b[0] = quick_pow(a[0], mod - 2);
    return ;
  }
  polyinv(a, b, n + 1 >> 1);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  for (int i = 0; i < n; i++) {
    c[i] = a[i];
  }
  for (int i = n; i < lim; i++) {
    c[i] = 0;
  }
  NTT(b, lim, 1);
  NTT(c, lim, 1);
  for (int i = 0; i < lim; i++) {
    int cur = (2 - 1ll * c[i] * b[i] % mod + mod) % mod;
    b[i] = 1ll * b[i] * cur % mod;
  }
  NTT(b, lim, -1);
  for (int i = n; i < lim; i++) {
    b[i] = 0;
  }
}

void derivative(int *a, int *b, int n) {
  for (int i = 0; i < n; i++) {
    b[i] = 1ll * a[i + 1] * (i + 1) % mod;
  }
}

void integrate(int *a, int n) {
  get_inv(n);
  for (int i = n - 1; i >= 1; i--) {
    a[i] = 1ll * a[i - 1] * inv[i] % mod;
  }
  a[0] = 0;
}

void polyln(int *a, int *b, int n) {
  derivative(a, b, n);
  polyinv(a, d, n);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  NTT(b, lim, 1);
  NTT(d, lim, 1);
  for (int i = 0; i < lim; i++) {
    b[i] = 1ll * b[i] * d[i] % mod;
  }
  NTT(b, lim, -1);
  integrate(b, n);
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  int n;
  scanf("%d", &n);
  for (int i = 0; i < n; i++) {
    scanf("%d", &a[i]);
  }
  polyln(a, b, n);
  for (int i = 0; i < n; i++) {
    printf("%d%c", b[i], i + 1 == n ? '\n' : ' ');
  }
  return 0;
}

多项式 exp ⁡ \exp exp

牛顿迭代推导式子有$f(x) \equiv f_0(x)(1 - \ln f_0(x) + h(x)) \pmod {x ^ n}\$

#include <bits/stdc++.h>

using namespace std;

const int mod = 998244353, inv2 = mod + 1 >> 1;

namespace Quadratic_residue {
  struct Complex {
    int r, i;

    Complex(int _r = 0, int _i = 0) : r(_r), i(_i) {}
  };

  int I2;

  Complex operator * (const Complex &a, Complex &b) {
    return Complex((1ll * a.r * b.r % mod  + 1ll * a.i * b.i % mod * I2 % mod) % mod, (1ll * a.r * b.i % mod + 1ll * a.i * b.r % mod) % mod);
  }

  Complex quick_pow(Complex a, int n) {
    Complex ans = Complex(1, 0);
    while (n) {
      if (n & 1) {
        ans = ans * a;
      }
      a = a * a;
      n >>= 1;
    }
    return ans;
  }

  int get_residue(int n) {
    mt19937 e(233);
    if (n == 0) {
      return 0;
    }
    if(quick_pow(n, (mod - 1) >> 1).r == mod - 1) {
      return -1;
    }
    uniform_int_distribution<int> r(0, mod - 1);
    int a = r(e);
    while(quick_pow((1ll * a * a % mod - n + mod) % mod, (mod - 1) >> 1).r == 1) {
      a = r(e);
    }
    I2 = (1ll * a * a % mod - n + mod) % mod;
    int x = quick_pow(Complex(a, 1), (mod + 1) >> 1).r, y = mod - x;
    if(x > y) swap(x, y);
    return x;
  }
}

const int N = 1e6 + 10;

int r[N], inv[N], a[N], b[N], c[N], d[N], e[N], t[N], n;

//a是输入数组，b存放多项式逆，c存放多项式开根，d存放多项式对数ln，e存放多项式指数exp，t作为中间转移数组

int quick_pow(int a, int n) {
  int ans = 1;
  while (n) {
    if (n & 1) {
      ans = 1ll * a * ans % mod;
    }
    a = 1ll * a * a % mod;
    n >>= 1;
  }
  return ans;
}

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  }
}

void get_inv(int n) {
  inv[1] = 1;
  for (int i = 2; i <= n; i++) {
    inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
  }
}

void NTT(int *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  for (int mid = 1; mid < lim; mid <<= 1) {
    int wn = quick_pow(3, (mod - 1) / (mid << 1));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      int w = 1;
      for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {
        int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;
        f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if (rev == -1) {
    int inv = quick_pow(lim, mod - 2);
    reverse(f + 1, f + lim);
    for (int i = 0; i < lim; i++) {
      f[i] = 1ll * f[i] * inv % mod;
    }
  }
}

void polyinv(int *f, int *g, int n) {
  if (n == 1) {
    g[0] = quick_pow(f[0], mod - 2);
    return ;
  }
  polyinv(f, g, n + 1 >> 1);
  for (int i = 0; i < n; i++) {
    t[i] = f[i];
  }
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  NTT(t, lim, 1);
  NTT(g, lim, 1);
  for (int i = 0; i < lim; i++) {
    int cur = (2 - 1ll * g[i] * t[i] % mod + mod) % mod;
    g[i] = 1ll * g[i] * cur % mod;
    t[i] = 0;
  }
  NTT(g, lim, -1);
  for (int i = n; i < lim; i++) {
    g[i] = 0;
  }
}

void polysqrt(int *f, int *g, int n) {
  if (n == 1) {
    g[0] = Quadratic_residue::get_residue(f[0]);
    return ;
  }
  polysqrt(f, g, n + 1 >> 1);
  polyinv(g, b, n);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  for (int i = 0; i < n; i++) {
    t[i] = f[i];
  }
  NTT(g, lim, 1);
  NTT(b, lim, 1);
  NTT(t, lim, 1);
  for (int i = 0; i < lim; i++) {
    g[i] = (1ll * inv2 * g[i] % mod + 1ll * inv2 * b[i] % mod * t[i] % mod) % mod;
    b[i] = t[i] = 0;
  }
  NTT(g, lim, -1);
  for (int i = n; i < lim; i++) {
    g[i] = 0;
  }
}

void derivative(int *a, int *b, int n) {
  for (int i = 0; i < n; i++) {
    b[i] = 1ll * a[i + 1] * (i + 1) % mod;
  }
}

void integrate(int *a, int n) {
  for (int i = n - 1; i >= 1; i--) {
    a[i] = 1ll * a[i - 1] * inv[i] % mod;
  }
  a[0] = 0;
}

void polyln(int *f, int *g, int n) {
  polyinv(f, b, n);
  derivative(f, g, n);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  NTT(g, lim, 1);
  NTT(b, lim, 1);
  for (int i = 0; i < lim; i++) {
    g[i] = 1ll * g[i] * b[i] % mod;
    b[i] = 0;
  }
  NTT(g, lim, -1);
  for (int i = n; i < lim; i++) {
    g[i] = 0;
  }
  integrate(g, n);
}

void polyexp(int *f, int *g, int n) {
  if (n == 1) {
    g[0] = 1;
    return ;
  }
  polyexp(f, g, n + 1 >> 1);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  polyln(g, d, n);
  for (int i = 0; i < n; i++) {
    t[i] = (f[i] - d[i] + mod) % mod;
  }
  t[0] = (t[0] + 1) % mod;
  get_r(lim);
  NTT(g, lim, 1);
  NTT(t, lim, 1);
  for (int i = 0; i < lim; i++) {
    g[i] = 1ll * g[i] * t[i] % mod;
    t[i] = d[i] =  0;
  }
  NTT(g, lim, -1);
  for (int i = n; i < lim; i++) {
    g[i] = 0;
  }
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  scanf("%d", &n);
  for (int i = 0; i < n; i++) {
    scanf("%d", &a[i]);
  }
  get_inv(4 * n);
  polyexp(a, e, n);
  for (int i = 0; i < n; i++) {
    printf("%d%c", e[i], i + 1 == n ? '\n' : ' ');
  }
  return 0;
}

多项式快速幂

给定 f ( x ) f(x) f(x)，要求 g ( x ) ≡ f k ( x ) ( m o d x n ) g(x) \equiv f ^{k}(x) \pmod{x ^ n} g(x)≡fk(x)(modxn)。
g ( x ) = exp ⁡ ln ⁡ f k ( x ) = exp ⁡ k ln ⁡ f ( x ) g(x) = \exp ^{\ln f^{k}(x)} = \exp ^{k \ln f(x)}\\ g(x)=explnfk(x)=expklnf(x)

#include <bits/stdc++.h>

using namespace std;

const int mod = 998244353, inv2 = mod + 1 >> 1;

namespace Quadratic_residue {
  struct Complex {
    int r, i;

    Complex(int _r = 0, int _i = 0) : r(_r), i(_i) {}
  };

  int I2;

  Complex operator * (const Complex &a, Complex &b) {
    return Complex((1ll * a.r * b.r % mod  + 1ll * a.i * b.i % mod * I2 % mod) % mod, (1ll * a.r * b.i % mod + 1ll * a.i * b.r % mod) % mod);
  }

  Complex quick_pow(Complex a, int n) {
    Complex ans = Complex(1, 0);
    while (n) {
      if (n & 1) {
        ans = ans * a;
      }
      a = a * a;
      n >>= 1;
    }
    return ans;
  }

  int get_residue(int n) {
    mt19937 e(233);
    if (n == 0) {
      return 0;
    }
    if(quick_pow(n, (mod - 1) >> 1).r == mod - 1) {
      return -1;
    }
    uniform_int_distribution<int> r(0, mod - 1);
    int a = r(e);
    while(quick_pow((1ll * a * a % mod - n + mod) % mod, (mod - 1) >> 1).r == 1) {
      a = r(e);
    }
    I2 = (1ll * a * a % mod - n + mod) % mod;
    int x = quick_pow(Complex(a, 1), (mod + 1) >> 1).r, y = mod - x;
    if(x > y) swap(x, y);
    return x;
  }
}

const int N = 1e6 + 10;

int r[N], inv[N], a[N], b[N], c[N], d[N], e[N], t[N], n;

int quick_pow(int a, int n) {
  int ans = 1;
  while (n) {
    if (n & 1) {
      ans = 1ll * a * ans % mod;
    }
    a = 1ll * a * a % mod;
    n >>= 1;
  }
  return ans;
}

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  }
}

void get_inv(int n) {
  inv[1] = 1;
  for (int i = 2; i <= n; i++) {
    inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
  }
}

void NTT(int *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  for (int mid = 1; mid < lim; mid <<= 1) {
    int wn = quick_pow(3, (mod - 1) / (mid << 1));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      int w = 1;
      for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {
        int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;
        f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if (rev == -1) {
    int inv = quick_pow(lim, mod - 2);
    reverse(f + 1, f + lim);
    for (int i = 0; i < lim; i++) {
      f[i] = 1ll * f[i] * inv % mod;
    }
  }
}

void polyinv(int *f, int *g, int n) {
  if (n == 1) {
    g[0] = quick_pow(f[0], mod - 2);
    return ;
  }
  polyinv(f, g, n + 1 >> 1);
  for (int i = 0; i < n; i++) {
    t[i] = f[i];
  }
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  NTT(t, lim, 1);
  NTT(g, lim, 1);
  for (int i = 0; i < lim; i++) {
    int cur = (2 - 1ll * g[i] * t[i] % mod + mod) % mod;
    g[i] = 1ll * g[i] * cur % mod;
    t[i] = 0;
  }
  NTT(g, lim, -1);
  for (int i = n; i < lim; i++) {
    g[i] = 0;
  }
}

void polysqrt(int *f, int *g, int n) {
  if (n == 1) {
    g[0] = Quadratic_residue::get_residue(f[0]);
    return ;
  }
  polysqrt(f, g, n + 1 >> 1);
  polyinv(g, b, n);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  for (int i = 0; i < n; i++) {
    t[i] = f[i];
  }
  NTT(g, lim, 1);
  NTT(b, lim, 1);
  NTT(t, lim, 1);
  for (int i = 0; i < lim; i++) {
    g[i] = (1ll * inv2 * g[i] % mod + 1ll * inv2 * b[i] % mod * t[i] % mod) % mod;
    b[i] = t[i] = 0;
  }
  NTT(g, lim, -1);
  for (int i = n; i < lim; i++) {
    g[i] = 0;
  }
}

void derivative(int *a, int *b, int n) {
  for (int i = 0; i < n; i++) {
    b[i] = 1ll * a[i + 1] * (i + 1) % mod;
  }
}

void integrate(int *a, int n) {
  for (int i = n - 1; i >= 1; i--) {
    a[i] = 1ll * a[i - 1] * inv[i] % mod;
  }
  a[0] = 0;
}

void polyln(int *f, int *g, int n) {
  polyinv(f, b, n);
  derivative(f, g, n);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  NTT(g, lim, 1);
  NTT(b, lim, 1);
  for (int i = 0; i < lim; i++) {
    g[i] = 1ll * g[i] * b[i] % mod;
    b[i] = 0;
  }
  NTT(g, lim, -1);
  for (int i = n; i < lim; i++) {
    g[i] = 0;
  }
  integrate(g, n);
}

void polyexp(int *f, int *g, int n) {
  if (n == 1) {
    g[0] = 1;
    return ;
  }
  polyexp(f, g, n + 1 >> 1);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  polyln(g, d, n);
  for (int i = 0; i < n; i++) {
    t[i] = (f[i] - d[i] + mod) % mod;
  }
  t[0] = (t[0] + 1) % mod;
  get_r(lim);
  NTT(g, lim, 1);
  NTT(t, lim, 1);
  for (int i = 0; i < lim; i++) {
    g[i] = 1ll * g[i] * t[i] % mod;
    t[i] = d[i] =  0;
  }
  NTT(g, lim, -1);
  for (int i = n; i < lim; i++) {
    g[i] = 0;
  }
}

/*
  a是输入数组，
  b存放多项式逆，
  c存放多项式开根，
  d存放多项式对数ln，
  e存放多项式指数exp，
  t作为中间转移数组,
  如果要用到polyinv，得提前调用get_inv(n)先预先得到我们想要得到的逆元范围。
*/

char str[N];

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  scanf("%d %s", &n, str + 1);
  for (int i = 0; i < n; i++) {
    scanf("%d", &a[i]);
  }
  int k = 0;
  for (int i = 1; str[i]; i++) {
    k = (1ll * k * 10 + (str[i] - '0')) % mod;
  }
  get_inv(4 * n);
  polyln(a, d, n);
  for (int i = 0; i < n; i++) {
    a[i] = 1ll * d[i] * k % mod;
    d[i] = 0;
  }
  polyexp(a, e, n);
  for (int i = 0; i < n; i++) {
    printf("%d%c", e[i], i + 1 == n ? '\n' : ' ');
  }
  return 0;
}

多项式除法

给定一个 n n n次多项式 F ( x ) F(x) F(x)和 m m m次多项式 G ( x ) G(x) G(x)，要求 R ( x ) , Q ( x ) R(x), Q(x) R(x),Q(x)，满足 F ( x ) = R ( x ) G ( x ) + Q ( x ) F(x) = R(x)G(x) + Q(x) F(x)=R(x)G(x)+Q(x)。

R ( x ) R(x) R(x)是一个 n − m n - m n−m阶多项式， Q ( x ) Q(x) Q(x)是一个小于 m m m阶的多项式。
有 F ( x ) ≡ R ( x ) G ( x ) + Q ( x ) ( m o d x n + 1 ) F ( 1 x ) ≡ R ( 1 x ) G ( 1 x ) + Q ( 1 x ) ( m o d x n + 1 ) 同 时 乘 上 一 个 x n , F r e v ( x ) ≡ ( x m R r e v ( x ) ) ( x n − m G r e v ( x ) ) + x n − d e g Q Q r e v ( x ) ( m o d x n + 1 ) F r e v ( x ) ≡ R r e v ( x ) G r e v ( x ) + Q r e v ( x ) x n − d e g Q ( m o d x n + 1 ) 有 d e g Q < m , n − d e g Q > = n − m + 1 , 所 以 有 F r e v ( x ) ≡ R r e v ( x ) G r e v ( x ) ( m o d x n − m + 1 ) 有F(x) \equiv R(x) G(x) + Q(x) \pmod{x ^ {n + 1}}\\ F(\frac{1}{x}) \equiv R(\frac{1}{x})G(\frac{1}{x}) + Q(\frac{1}{x}) \pmod {x ^{n + 1}}\\ 同时乘上一个x ^ n, F^{rev}(x) \equiv \left(x ^m R ^{rev}(x)\right) \left(x ^{n - m}G ^{rev}(x)\right) + x ^{n - deg_Q} Q ^{rev} (x) \pmod{x ^{n + 1}} \\ F^{rev}(x) \equiv R^{rev}(x) G^{rev}(x) + Q^{rev}(x) x ^{n - deg_Q}\pmod{x ^{n + 1}}\\ 有deg_Q < m, n - deg_Q >= n - m + 1,所以有F^{rev}(x) \equiv R^{rev}(x) G^{rev}(x) \pmod{x ^{n - m + 1}}\\ 有F(x)≡R(x)G(x)+Q(x)(modxn+1)F(x1​)≡R(x1​)G(x1​)+Q(x1​)(modxn+1)同时乘上一个xn,Frev(x)≡(xmRrev(x))(xn−mGrev(x))+xn−degQ​Qrev(x)(modxn+1)Frev(x)≡Rrev(x)Grev(x)+Qrev(x)xn−degQ​(modxn+1)有degQ​<m,n−degQ​>=n−m+1,所以有Frev(x)≡Rrev(x)Grev(x)(modxn−m+1)
只要多项式求逆，即可得到 R ( x ) R(x) R(x)，然后代入原式求得 Q ( x ) Q(x) Q(x)。

#include <bits/stdc++.h>

using namespace std;

const int mod = 998244353, inv2 = mod + 1 >> 1;

namespace Quadratic_residue {
  struct Complex {
    int r, i;

    Complex(int _r = 0, int _i = 0) : r(_r), i(_i) {}
  };

  int I2;

  Complex operator * (const Complex &a, Complex &b) {
    return Complex((1ll * a.r * b.r % mod  + 1ll * a.i * b.i % mod * I2 % mod) % mod, (1ll * a.r * b.i % mod + 1ll * a.i * b.r % mod) % mod);
  }

  Complex quick_pow(Complex a, int n) {
    Complex ans = Complex(1, 0);
    while (n) {
      if (n & 1) {
        ans = ans * a;
      }
      a = a * a;
      n >>= 1;
    }
    return ans;
  }

  int get_residue(int n) {
    mt19937 e(233);
    if (n == 0) {
      return 0;
    }
    if(quick_pow(n, (mod - 1) >> 1).r == mod - 1) {
      return -1;
    }
    uniform_int_distribution<int> r(0, mod - 1);
    int a = r(e);
    while(quick_pow((1ll * a * a % mod - n + mod) % mod, (mod - 1) >> 1).r == 1) {
      a = r(e);
    }
    I2 = (1ll * a * a % mod - n + mod) % mod;
    int x = quick_pow(Complex(a, 1), (mod + 1) >> 1).r, y = mod - x;
    if(x > y) swap(x, y);
    return x;
  }
}

const int N = 6e5 + 10;

int r[N], inv[N], b[N], c[N], d[N], e[N], t[N];

int quick_pow(int a, int n) {
  int ans = 1;
  while (n) {
    if (n & 1) {
      ans = 1ll * a * ans % mod;
    }
    a = 1ll * a * a % mod;
    n >>= 1;
  }
  return ans;
}

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  }
}

void get_inv(int n) {
  inv[1] = 1;
  for (int i = 2; i <= n; i++) {
    inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
  }
}

void NTT(int *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  for (int mid = 1; mid < lim; mid <<= 1) {
    int wn = quick_pow(3, (mod - 1) / (mid << 1));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      int w = 1;
      for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {
        int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;
        f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if (rev == -1) {
    int inv = quick_pow(lim, mod - 2);
    reverse(f + 1, f + lim);
    for (int i = 0; i < lim; i++) {
      f[i] = 1ll * f[i] * inv % mod;
    }
  }
}

void polyinv(int *f, int *g, int n) {
  if (n == 1) {
    g[0] = quick_pow(f[0], mod - 2);
    return ;
  }
  polyinv(f, g, n + 1 >> 1);
  for (int i = 0; i < n; i++) {
    t[i] = f[i];
  }
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  NTT(t, lim, 1);
  NTT(g, lim, 1);
  for (int i = 0; i < lim; i++) {
    int cur = (2 - 1ll * g[i] * t[i] % mod + mod) % mod;
    g[i] = 1ll * g[i] * cur % mod;
    t[i] = 0;
  }
  NTT(g, lim, -1);
  for (int i = n; i < lim; i++) {
    g[i] = 0;
  }
}

void polysqrt(int *f, int *g, int n) {
  if (n == 1) {
    g[0] = Quadratic_residue::get_residue(f[0]);
    return ;
  }
  polysqrt(f, g, n + 1 >> 1);
  polyinv(g, b, n);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  for (int i = 0; i < n; i++) {
    t[i] = f[i];
  }
  NTT(g, lim, 1);
  NTT(b, lim, 1);
  NTT(t, lim, 1);
  for (int i = 0; i < lim; i++) {
    g[i] = (1ll * inv2 * g[i] % mod + 1ll * inv2 * b[i] % mod * t[i] % mod) % mod;
    b[i] = t[i] = 0;
  }
  NTT(g, lim, -1);
  for (int i = n; i < lim; i++) {
    g[i] = 0;
  }
}

void derivative(int *a, int *b, int n) {
  for (int i = 0; i < n; i++) {
    b[i] = 1ll * a[i + 1] * (i + 1) % mod;
  }
}

void integrate(int *a, int n) {
  for (int i = n - 1; i >= 1; i--) {
    a[i] = 1ll * a[i - 1] * inv[i] % mod;
  }
  a[0] = 0;
}

void polyln(int *f, int *g, int n) {
  polyinv(f, b, n);
  derivative(f, g, n);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  NTT(g, lim, 1);
  NTT(b, lim, 1);
  for (int i = 0; i < lim; i++) {
    g[i] = 1ll * g[i] * b[i] % mod;
    b[i] = 0;
  }
  NTT(g, lim, -1);
  for (int i = n; i < lim; i++) {
    g[i] = 0;
  }
  integrate(g, n);
}

void polyexp(int *f, int *g, int n) {
  if (n == 1) {
    g[0] = 1;
    return ;
  }
  polyexp(f, g, n + 1 >> 1);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  polyln(g, d, n);
  for (int i = 0; i < n; i++) {
    t[i] = (f[i] - d[i] + mod) % mod;
  }
  t[0] = (t[0] + 1) % mod;
  get_r(lim);
  NTT(g, lim, 1);
  NTT(t, lim, 1);
  for (int i = 0; i < lim; i++) {
    g[i] = 1ll * g[i] * t[i] % mod;
    t[i] = d[i] =  0;
  }
  NTT(g, lim, -1);
  for (int i = n; i < lim; i++) {
    g[i] = 0;
  }
}

/*
  b存放多项式逆，
  c存放多项式开根，
  d存放多项式对数ln，
  e存放多项式指数exp，
  t作为中间转移数组,
  如果要用到polyinv，得提前调用get_inv(n)先预先得到我们想要得到的逆元范围。
*/

int f[N], fr[N], g[N], gr[N], rr[N], n, m;

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  scanf("%d %d", &n, &m);
  for (int i = 0; i <= n; i++) {
    scanf("%d", &f[i]);
    fr[n - i] = f[i];
  }
  for (int i = 0; i <= m; i++) {
    scanf("%d", &g[i]);
    gr[m - i] = g[i];
  }
  for (int i = n - m + 1; i <= n; i++) {
    fr[i] = gr[i] = 0;
  }
  polyinv(gr, b, n - m + 1);
  for (int i = 0; i < n - m + 1; i++) {
    gr[i] = b[i];
    b[i] = 0;
  }
  int lim = 1;
  while (lim < 2 * (n - m + 1)) {
    lim <<= 1;
  }
  get_r(lim);
  NTT(fr, lim, 1);
  NTT(gr, lim, 1);
  for (int i = 0; i < lim; i++) {
    fr[i] = 1ll * fr[i] * gr[i] % mod;
  }
  NTT(fr, lim, -1);
  for (int i = 0; i <= n - m; i++) {
    rr[i] = fr[n - m - i];
  }
  for (int i = 0; rr[i]; i++) {
    printf("%d ", rr[i]);
  }
  puts("");
  lim = 1;
  while (lim <= 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  NTT(rr, lim, 1);
  NTT(g, lim, 1);
  for (int i = 0; i < lim; i++) {
    g[i] = 1ll * g[i] * rr[i] % mod;
  }
  NTT(g, lim, -1);
  for (int i = 0; i < m; i++) {
    f[i] = (f[i] - g[i] + mod) % mod;
  }
  for (int i = 0; i < m; i++) {
    printf("%d ", f[i]);
  }
  puts("");
  return 0;
}
_lifehappy_
关注
  • 0
    点赞
  • 踩
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
多项式全家桶
Recalling_Clouds的博客
01-23 724
下面的所有知识均属多项式范畴。 前言 本文略去多项式的定义，毕竟初中数学都讲过。 本文会出现一些可能没有见过的表达： 对于一个关于 xxx 的多项式 f(x)f(x)f(x)，有时也称其最高次项次数为 f(x)f(x)f(x) 的度数。 多项式同样定义了四则运算，包括常数与多项式的四则运算和多项式之间的四则运算，同时也衍生出了同余符号 ≡\equiv≡、对数函数 ln⁡(f(x))\ln(f(x))ln(f(x))，指数函数 exp⁡(f(x))\exp(f(x))exp(f(x)) 等等。 下文一般采
多项式拟合平滑.docx
最新发布
06-27
多项式拟合平滑是一种数据分析技术，用于揭示数据集中的潜在趋势和模式。该方法通过构建一个多项式函数，使这个函数尽可能地贴近原始数据点，从而实现对数据的平滑表示。这种方法在处理含有噪声或者不规则变化的数据...
切比雪夫多项式(Chebyshev Polynomials)
deye1979的专栏
12-28 6097
切比雪夫多项式在逼近理论中有重要的应用。这是因为第一类切比雪夫多项式的根（被称为切比雪夫节点）可以用于多项式插值。相应的插值多项式能最大限度地降低龙格现象，并且提供多项式在连续函数的最佳一致逼近。 参考资料：https://wenku.baidu.com/view/ba41a20f767f5acfa1c7cd9c.html https://wenku.baidu.com/vie...
【学习笔记】多项式全家桶（包含全套证明）
热门推荐
繁凡さん的博客
12-29 1万+
目录0x00.多项式基本概述0x01.多项式的基础概念0x02. 多项式的度0x03. 多项式的乘法0x04. 多项式的逆元0x10.拉格朗日插值0x20.多项式乘法0x21.快速傅里叶变换（FFT）0x22.快速数论变换（NTT）0x23.快速沃尔什变换（FWT）0x24.快速莫比乌斯变换（FMT）0x30.多项式求逆0x40.多项式开方0x50.多项式除法 / 取模0x60.多项式牛顿迭代0x70.多项式多点求值 / 快速插值0x80.多项式三角函数0x90.多项式反...
多项式全家桶板子
lahlah的博客
04-07 141
文章目录多项式FWTFFTNTT多项式求逆分治FFT(假，这是求逆)任意模数NTT任意模数多项式求逆多项式对数函数（多项式 ln）多项式指数函数（多项式 exp）多项式开根多项式快速幂 多项式 基本常识： 不会了就拿这个推 FWT #include<bits/stdc++.h> #define mod 998244353 #define ll long long using namespace std; const int N = (1 << 17) + 5; void fwt_
多项式拟合_多项式拟合_VBa_
10-01
在Excel VBA（Visual Basic for Applications）环境中，我们可以利用VBA强大的编程能力来实现复杂的数学运算，例如多项式拟合。多项式拟合是一种通过拟合数据点来创建一个多项式函数的过程，常用于数据分析、预测和...
多项式拟合c++程序
01-18
对数据进行多项式拟合，可得拟合多项式的系数，法方程的系数矩阵，平方误差。
legendre_勒让德多项式拟合_多项式拟合_
09-29
勒让德多项式是一种在数学和工程领域广泛应用的特殊函数，尤其在数值分析和插值问题中占有重要地位。勒让德多项式是由勒让德（Pierre-Simon Laplace）提出的，它是一组正交于[ -1, 1 ]区间上的多项式序列，具有良好...
多项式.zip_多项式
07-15
用链表实现多项式的加减乘除以及多项式的打印显示
快速傅里叶变换（完整推导过程 + 模板）
lifehappy
01-28 3317
快速傅里叶变换 多项式表示 系数表示法： 一个nnn次多项式可以用n+1n + 1n+1个系数表示出来：f(x)=a0+a1x+a2x2+⋯+an−1xn−1+anxnf(x) = a_0 + a_1 x + a_2 x ^ 2 + \dots + a_{n - 1} x ^{n- 1} + a_n x ^nf(x)=a0​+a1​x+a2​x2+⋯+an−1​xn−1+an​xn。 点值表示法： 通过线性代数，高斯消元我们可以知道，一个nnn次多项式可以通过n+1n + 1n+1个点联立方程组解得： f(
伯努利数（详解 + 例题 ：P3711 仓鼠的数学题）
lifehappy
03-05 1632
伯努利数 定义Sk(n)=∑i=0n−1ikS_k(n) = \sum\limits_{i = 0} ^{n - 1} i ^ kSk​(n)=i=0∑n−1​ik。 从二项式出发 (0+1)k+1=∑i=0kCk+1i0i+0k+1⋮(n−1+1)k+1=∑i=0kCk+1i(n−1)i+(n−1)k+1把次方k+1的移项，再整体相加，得nk+1=∑i=0kCk+1iSi(n)nk+1=∑i=0k−1Ck+1iSi(n)+(k+1)Sk(n)Sk(n)=1k+1(nk+1−∑i=0k−1Ck+1iSi(n
多项式对数函数ln f(x)
lifehappy
01-30 1581
多项式对数函数ln⁡f(x)\ln f(x)lnf(x) 如果存在解必然有[x0]f(x)=1[ x ^ 0]f(x) = 1[x0]f(x)=1， 对ln⁡f(x)\ln f(x)lnf(x)求导，有dln⁡f(x)dx≡f′(x)f(x)(modxn)\frac{d \ln f(x)}{dx} \equiv \frac{f'(x)}{f(x)} \pmod {x ^ n}dxdlnf(x)​≡f(x)f′(x)​(modxn)， dxdxdx乘到右边，再求积分有： ∫dln⁡f(x)≡∫f′(x))f(
多项式除法，多项式取模
lifehappy
02-05 1389
多项式除法 给定一个nnn次多项式F(x)F(x)F(x)和mmm次多项式G(x)G(x)G(x)，要求R(x),Q(x)R(x), Q(x)R(x),Q(x)，满足F(x)=R(x)G(x)+Q(x)F(x) = R(x)G(x) + Q(x)F(x)=R(x)G(x)+Q(x)。 R(x)R(x)R(x)是一个n−mn - mn−m阶多项式，Q(x)Q(x)Q(x)是一个小于mmm阶的多项式。 有F(x)≡R(x)G(x)+Q(x)(modxn+1)F(1x)≡R(1x)G(1x)+Q(1x)(modx
多项式模板（inv, sqrt, ln, exp）
lifehappy
02-04 1364
#include <bits/stdc++.h> using namespace std; const int mod = 998244353, inv2 = mod + 1 >> 1; namespace Quadratic_residue { struct Complex { int r, i; Complex(int _r = 0, int _i = 0) : r(_r), i(_i) {} }; int I2; Complex op
拉格朗日插值
lifehappy
09-03 1064
拉格朗日插值 简要阐述 结论给定n+1n + 1n+1个点最多可以得到一个nnn次多项式的表达式，并且f(x)=∑i=1nyi∏j∤ix−xjxi−yjf(x) = \sum_{i = 1} ^{n} y_i \prod\limits_{j \nmid i}\frac{x - x_j}{x_i - y_j}f(x)=∑i=1n​yi​j∤i∏​xi​−yj​x−xj​​我们随便往里带入一个xi,yjx_i, y_jxi​,yj​可以得到等式是成立的，所以我们只要套模板即可求得某个函数的值。 特殊情况 当xx
导数卷积 （NTT）
lifehappy
03-04 902
导数卷积 有f(x)=∑i=0n−1ai×xi求g(x)=∑i=0n−1f(i)(x)f(n−i−1)(x)∑i=0n−1∑j=0n−1ai+j(i+j)!j!∑k=0n−1ak+(n−i−1)(k+(n−i−1))!k!设F(n)=an×n!我们单独求解g(x)的第m项，有gm=∑i=0n−1∑j=1d(F(i+j)∗F(n−i−1+d−j))×(1j!∗1(d−j)!)设H(n)=1n! 有f(x) = \sum_{i = 0} ^{n - 1}a_i \times x ^ i\\ 求g(x) = \s
生成函数简单入门
lifehappy
01-31 785
生成函数 可表示为F(x)=∑nankn(x)F(x) = \sum\limits_{n} a_n k_n(x)F(x)=n∑​an​kn​(x)，对于不同类型的生成函数，有不同的核函数kn(x)k_n(x)kn​(x)。 普通生成函数：kn(x)=xnk_n(x) = x ^ nkn​(x)=xn。 指数生成函数：kn(x)=xnn!k_n(x) = \frac{x ^ n}{n !}kn​(x)=n!xn​。 迪利克雷生成函数：kn(x)=1nxk_n(x) = \frac{1}{n ^ x}kn​(x
C语言课程设计一元多项式计算系统
首先提示用户输入两个多项式A和B，输入的形式是先输入多项式A的项数，然后回车，接着输入多项式A的第一项的系数，每一项按照先输入系数，然后输入空格或回车，再输入多项式A的第一项的指数，然后继续输入多项式A的...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值