多項式科技初步

時間  2019-12-11
標籤 多項式 科技 初步 简体版
原文   原文鏈接

寫在前面

爲了體現簡潔，在每一部分只會放關鍵代碼。c++

代碼具備必定的通用性，保證代碼之間函數的調用是合法的。數組

若是 MathJax 加載不出來或加載有誤，請您多刷新幾回。函數

總結可能會比較長，能夠點右邊的小火箭回到目錄。測試

有什麼問題請聯繫我，萬分感謝。優化

參考資料

Picks 的博客 Picks's Blogui

Miskcoo 的博客 Miskcoo's Spacespa

多項式乘法

問題描述

給定兩個多項式 $A(x) $ 和 \(B(x)\) ：
\[ A(x)=\sum_{i=0}^{n}a_ix^i\ ,\ B(x)=\sum_{i=0}^{n} b_i x^i \]
求卷積 \(C(x) =A(x) * B(x)\) ，知足
\[ C(x)=\sum_{i=0}^{n} c_i x^i\ ,\ c_i=\sum_{k=0}^{i} a_k\times b_{i-k} \]code

解決方法

這部分上面提到的神仙們講解的都十分詳細。orm

個人小結限於篇幅，就直接掛 pdf 版本的鏈接了：快速傅里葉變換初步blog

代碼實現

使用 [ UOJ 34 ] 多項式乘法 做爲測試題。

  • Fast Fourier Transform，使用的是 3 次變換的最基本寫法，用時 362 ms

    inline void FFT(Complex *f, int len, int o) {
    for (int i = 0; i < len; ++i)
      if (rev[i] > i) swap(f[i], f[rev[i]]);
    for (int i = 1; i < len; i <<= 1) {
      Complex wn = Complex(cos(PI / i) , o * sin(PI / i));
      for (int j = 0; j < len; j += (i << 1)) {
        Complex w = Complex(1, 0), x, y;
        for (int k = 0; k < i; ++k, w = w * wn) {
          x = f[j + k]; y = w * f[i + j + k];
          f[i + j + k] = x - y; f[j + k] = x + y;
        }
      }
    }
    if (o == -1) for (int i = 0; i < len; ++i) f[i].x /= len;
}

  • Fast Number-Theoretic Transform，模數爲 998244353，用時 403 ms

    inline void NTT(int *f, int len, int o) {
  for (int i = 1; i < len; ++i)
    if (i > rev[i]) swap(f[i], f[rev[i]]);
  for (int i = 1; i < len; i <<= 1) {
    int wn = qpow(3, (mod - 1) / (i << 1));
    if (o == -1) wn = qpow(wn, mod - 2);
    for (int j = 0; j < len; j += (i << 1)) {
      int w = 1, x, y;
      for (int k = 0; k < i; ++k, w = 1ll * w * wn % mod) {
        x = f[j + k]; y = 1ll * w * f[i + j + k] % mod;
        f[j + k] = mo(x + y); f[i + j + k] = mo(x - y + mod);
      }
    }
  }
  if (o == -1) {
    int invl = qpow(len, mod - 2);
    for (int i = 0; i < len; ++i) f[i] = 1ll * f[i] * invl % mod;
  }
}

多項式求逆

問題描述

給定一個 \(n\) 次多項式 \(A(x)\) ，求出一個多項式 \(B(x)\)， 知足
\[ A(x) * B(x) \equiv 1 \pmod{ x^{n+1}} \]
係數對 998244353 取模。

解決方法

採用倍增的思想。

考慮只有常數項的時候，\(A(x)\equiv c\pmod{x}\) ，那麼 \(A^{-1}(x)\) 即爲 \(c^{-1}\)。

對於 \(n>1\) 的時候，設 \(B(x)=A^{-1}(x)\) ，有
\[ A(x)B(x)\equiv 1\pmod{x^n} \]
由於此時模 \(x^n\) 至關於只保留多項式前 \(n\) 項，因此該同餘式在模 \(x^k,0\le k\le n\) 時都成立
\[ A(x)B(x)\equiv 1\pmod{x^{\lceil\frac{n}{2}\rceil}} \]
假設在 \(\pmod {x^{\lceil \frac n2\rceil}}\) 意義下 \(A(x)\) 的逆元是 \(B′(x)\) 而且咱們已經求出，那麼
\[ A(x)B'(x) \equiv 1 \pmod {x^{\lceil \frac{n}{2} \rceil}} \]
兩式相減，得
\[ B(x)-B'(x) \equiv 0 \pmod {x^{\lceil \frac{n}{2} \rceil}} \]
兩邊平方，得
\[ B^2(x)-2B(x)B'(x)+B'^2(x)\equiv 0\pmod{x^n} \]
模數平方的合法性在於，大於 \(n\) 的係數卷積中每組乘法必然有一項下標小於 \(n\) ，所以乘起來必然爲 \(0\) 。

兩側同稱 \(A(x)\) ，整理得
\[ B(x)\equiv2B'(x)-A(x)B'^2(x)\pmod{x^n} \]
所以只需將 \(B'(x)\) 和 \(A(x)\) 在模 \(x^n\) 意義下的插值求出，有
\[ B_i=2B'_i-A_iB'^2_i=B'_i(2-A_iB') \]
所以一遍 NTT 就能夠由 \(B'(x)\) 求出 \(B(x)\) 了。

總的時間複雜度爲
\[ T(n) = T(\frac{n}{2}) + \mathcal O(n \log n) = \mathcal O(n \log n) \]
由此過程也能夠獲得一個結論：一個多項式有沒有逆元徹底取決於其常數項是否有逆元。

代碼實現

使用 [ Luogu P4238 ] 多項式求逆 做爲測試題。

  • 遞歸版本，使用 O2 優化，用時 562 ms

    inline void Inv(int *a, int *b, int n) {
  if (n == 1) {b[0] = qpow(a[0], mod - 2); return;}
  Inv(a, b, (n + 1) >> 1);
  int len = Rev(n << 1);
  for (int i = 0; i < n; ++i) tmp[i] = a[i];
  for (int i = n; i < len; ++i) b[i] = tmp[i] = 0;
  NTT(b, len, 1); NTT(tmp, len, 1);
  for (int i = 0; i < len; ++i)
    b[i] = (2ll - 1ll * tmp[i] * b[i] % mod + mod) * b[i] % mod;
  NTT(b, len, -1);
  for (int i = 0; i < len; ++i) tmp[i] = 0;
  for (int i = n; i < len; ++i) b[i] = 0;
}

  • 迭代版本，使用 O2 優化，用時 570 ms

    inline void Inv(int *a, int *b, int n) {
  b[0] = qpow(a[0], mod - 2);
  int len;
  for (int l = 1; l < (n << 1); l <<= 1) {
    len = Rev(l << 1);
    for (int i = 0; i < l; ++i) tmp[i] = a[i];
    NTT(b, len, 1); NTT(tmp, len, 1);
    for (int i = 0; i < len; ++i)
      b[i] = (2ll - 1ll * tmp[i] * b[i] % mod + mod) * b[i] % mod;
    NTT(b, len, -1);
    for (int i = l; i < len; ++i) b[i] = 0;
  }
}

多項式開根

問題描述

給定一個 \(n\) 次多項式 \(A(x)\)，求一個在 \(\bmod{x^{n+1}}\) 意義下的多項式 \(B(x)\)，使得

\[ B^2(x) \equiv A(x) \pmod{x^{n+1}} \]

多項式的係數在模 998244353 意義下進行運算，保證常數項 \(a_0=1\)。

解決方法

一樣採用倍增的思想。

考慮只有常數項的時候，\(A(x)\equiv c\pmod x\) ，那麼 \(\sqrt {A(x)}\) 即爲 \(\sqrt c \equiv 1\pmod{x}\) （二次剩餘）。

對於 \(n>1​\) 的時候，一樣根據上一題的結論，咱們能夠把問題範圍縮小到 \(\bmod {x^{\lceil \frac{n}2\rceil}}​\) ，有
\[ B^2(x)\equiv A(x)\pmod{x^n}\ \Rightarrow\ B^2(x)\equiv A(x)\pmod{x^{\lceil \frac{n}2\rceil}} \]
不妨設咱們已經求出來了 \(\bmod{x^{\lceil \frac{n}2\rceil}} ​\) 意義下的根 \(D(x)​\)，即
\[ D^2(x)\equiv A(x)\pmod{x^{\lceil \frac{n}2\rceil}} \]
所以 \(B(x)​\) 與 \(D(x)​\) 在模 \(x^{\lceil \frac{n}2\rceil}​\) 意義下同餘，移項得
\[ B(x)-D(x)\equiv 0\pmod{x^{\lceil \frac{n}2\rceil}} \]
兩側平方，得
\[ B^2(x)+D^2(x)-2B(x)D(x)\equiv 0\pmod{x^n} \]
模數能平方的緣由與上一題相同。

咱們知道\(\bmod{x^n}\) 時 \(B^2(x)\) 即爲 \(A(x)\) ，所以
\[ A(x)+D^2(x)-2B(x)D(x)\equiv 0\pmod{x^n} \]
移項，得
\[ B(x)\equiv \frac{D^2(x) + A(x)}{2D(x)}\equiv \bigg(D(x) +\frac{A(x)}{D(x)}\bigg)\times 2^{-1}\pmod{x^n} \]
所以倍增時進行多項式求逆便可，總的時間複雜度爲
\[ T(n) = T(\frac{n}{2}) + \mathcal O(n \log n) = \mathcal O(n \log n) \]

代碼實現

使用 [ Luogu P5205 ] 多項式開根 做爲測試題，多項式求逆部分均採用遞歸版本。

  • 遞歸版本，使用 O2 優化，用時 3081 ms

    inline void Sqrt(int *a, int *b, int n) {
  if (n == 1) {b[0] = 1; return;}
  Sqrt(a, b, (n + 1) >> 1);
  Inv(b, b0, n);
  int len = Rev(n << 1);
  for (int i = 0; i < n; ++i) a0[i] = a[i];
  for (int i = n; i < len; ++i) a0[i] = 0;
  NTT(a0, len, 1); NTT(b0, len, 1);
  for (int i = 0; i < len; ++i) a0[i] = 1ll * a0[i] * b0[i] % mod;
  NTT(a0, len, -1);
  for (int i = 0; i < n; ++i) b[i] = 1ll * (b[i] + a0[i]) % mod * inv2 % mod;
  for (int i = n; i < len; ++i) b[i] = 0;
}

  • 迭代版本，使用 O2 優化，用時 3104 ms

    inline void Sqrt(int *a, int *b, int n) {
b[0] = 1;
int len;
for (int l = 1; l < (n << 1); l <<= 1) {
  Inv(b, b0, l);
  len = Rev(l << 1);
  for (int i = 0; i < l; ++i) a0[i] = a[i];
  for (int i = l; i < len; ++i) a0[i] = 0;
  NTT(a0, len, 1); NTT(b0, len, 1);
  for (int i = 0; i < len; ++i) a0[i] = 1ll * a0[i] * b0[i] % mod;
  NTT(a0, len, -1);
  for (int i = 0; i < l; ++i) b[i] = 1ll * (b[i] + a0[i]) % mod * inv2 % mod;
  for (int i = l; i < len; ++i) b[i] = 0;
}
}

多項式除法和取模

問題描述

給定一個 \(n\) 次多項式 \(A(x)\) 和一個 \(m\) 次多項式 \(B(x)\)，求出多項式 \(D(x)\), \(R(x)\)，知足

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

\(D(x)\) 次數爲 \(n-m\)，\(R(x)\) 次數小於 \(m\) ，全部的運算在模 998244353 意義下進行。

解決方法

注意到帶着 \(R(x)​\) 在這裏很麻煩，前人們想到了一個神奇的解決辦法。

設 \(A^R(x)=x^nA(\frac{1}{x})\) ，咱們將右側展開：
\[ A^R(x)=x^nA(\frac{1}{x})=x^n\sum_{i=0}^na_i\frac{1}{x^i}=\sum_{i=0}^na_ix^{n-i} \]
因此 \(A^R(x)\) 就是將 \(A(x)\) 的係數反轉 。

咱們將所求的等式中 \(x\) 所有換成 \(\frac{1}{x}\) ，而後兩側同乘 \(x^n\) ：
\[ x^n A(\frac{1}{x}) = x^{n - m}D(\frac{1}{x}) x^mB(\frac{1}{x}) + x^{n - m + 1}x^{m - 1}R(\frac{1}{x}) \]
\[ A^R(x) = D^R(x)B^R(x) + x^{n - m + 1}R^R(x) \]

注意到 \(D^R(x)\) 最高次反轉後不變，依然爲 \(n-m\)。

而右側的 \(R^R(x)\) 由於前面有 \(x^{n-m+1}\) 因此最低次爲 \(n-m+1\) 。

因此咱們能夠把多項式運算在 \(\bmod{x^{n-m+1}}\) 意義下進行，這樣 \(R(x)\) 就消失了：
\[ A^R(x) \equiv D^R(x)B^R(x) \pmod{x^{n-m+1}} \]
所以就能夠獲得 \(D^R(x)\) 的解法：
\[ \frac{A^R(x)}{B^R(x)} \equiv D^R(x) \pmod{x^{n-m+1}} \]
一個多項式求逆就能夠求出 \(D^R(x)\)了，再將 \(D^R(x)\) 進行反轉就獲得了答案。

將求出的 \(D(x)\) 回代，再進行一次減法便可求出 \(R(x)\)。

複雜度與多項式求逆同階，爲 \(\mathcal O(n \log n)​\) 。

代碼實現

使用 [ Luogu P4512 ] 多項式除法 做爲測試題，再也不區分多項式求逆部分的實現方式。

多項式求逆遞部分歸版本，使用 O2 優化，用時 710 ms

注意求 \(D(x)\) 部分的那次卷積是在 \(\bmod{x^{n-m+1}}\)意義下的，因此 \(A^R(x)\) 和 \(B^R(x)\) 中高於 \(n-m+1\) 的項須要清空（由於 FFT 卷積過程當中高次係數也會對低次係數形成影響）

inline void Div(int *a, int *b, int n, int m) {
  for (int i = 0; i <= n; ++i) ar[i] = a[n - i];
  for (int i = 0; i <= m; ++i) br[i] = b[m - i];
  for (int i = n - m + 2; i <= n; ++i) ar[i] = 0;
  for (int i = n - m + 2; i <= m; ++i) br[i] = 0;
  Inv(br, invb, n - m + 1);
  int len = Rev((n - m + 1) << 1);
  NTT(ar, len, 1); NTT(invb, len, 1);
  for (int i = 0; i < len; ++i) ar[i] = 1ll * ar[i] * invb[i] % mod;
  NTT(ar, len, -1);
  for (int i = 0; i <= n - m; ++i) tmp[i] = d[i] = ar[n - m - i];
  len = Rev(n << 1);
  for (int i = n - m + 1; i <= len; ++i) tmp[i] = 0;
  NTT(b, len, 1); NTT(tmp, len, 1);
  for (int i = 0; i < len; ++i) b[i] = 1ll * b[i] * tmp[i] % mod;
  NTT(b, len, -1);
  for (int i = 0; i <= m; ++i) r[i] = mo(a[i] - b[i] + mod);
}

分治 FFT

問題描述

給定長度爲 \(n-1\) 的數組 \(g[1],\dots,g[n-1]\)，求長度爲 \(n\) 的數組 \(f[0],\dots,f[n-1]\)，其中

\[ f[i]=\sum_{j=1}^if[i-j]g[j] \]

邊界爲 \(f[0]=1\)，運算在模 998244353 下進行。

解決方法

分治求解

暴力作是 \(\mathcal O(n^2)\) 的，考慮將相同的轉移一塊兒作以達到優化的目的。

考慮使用相似 CDQ 分治的思想，每次咱們求出 \([L, mid]\) 範圍內的 \(f\) 數組以後，把這部分 \(f\) 對 \([mid+1, R]\) 範圍內 \(f\) 的貢獻一塊兒作。

考慮對 \(x\in[mid + 1, R]\) 的 \(f[x]\) 的貢獻 \(w_x\)，有
\[ w_x=\sum_{i=l}^{mid} f[i]g[x-i] \]
所以 \(w\) 數組能夠卷積求了，注意求 \(w_x\) 時後半段的 \(f\) 須要認爲是 \(0\)，不然就存在右區間內部的貢獻了。

總的時間複雜度爲
\[ T(n) = 2T(\frac{n}{2}) + \mathcal O(n \log n) = \mathcal O(n \log^2 n) \]

多項式求逆

一階分治 FFT 是能夠看做卷積處理的。

不妨設將數組當作多項式，有

\[ F(x)=\sum_{i=0}^{n-1}f[i]x^i\ ,\ G(x)=\sum_{i=0}^{n-1}g[i]x^i \]

將兩個多項式卷積，有
\[ F(x)G(x)=\sum_{i=0}^{\infty}x^i\sum_{j}f[i-j]g[j]=F(x)-f[0] \]
後一個等式成立的緣由是，注意到後一個求和就是 \(f[i]\) 的形式，因此只有 \(f[0]\) 沒有被計數

求 \(f\) 數組能夠看做是 \(\bmod{x^n}\) 意義下進行的，所以有
\[ F(x)G(x) \equiv F(x)-f[0] \pmod x^n\ \Rightarrow\ F(x) \equiv \frac{f_0}{1-G(x)} \equiv (1-G(x))^{-1} \pmod x^n \]
因而一遍多項式求逆就能夠求出來了，複雜度爲 \(\mathcal O(n\log n)\)

代碼實現

使用 [ Luogu P4721 ] 分治 FFT 做爲測試題，多項式求逆採用遞歸版本。

  • 分治版本，使用 O2 優化，用時 974 ms

  • 多項式求逆版本，使用 O2 優化，用時 261 ms

    inline void solve(int *a, int n) {
  a[0] = 1;
  for (int i = 1; i < n; ++i) a[i] = mo(mod - a[i]);
  Inv(a, b, n);
  for (int i = 0; i < n; ++i) print(b[i], 0);
}

多項式求導和積分

問題描述

給定一個 \(n\) 次多項式 \(A(x)\) ，求一個 \(n-1\) 次多項式 \(B(x)\) ，和一個 \(n+1\) 次多項式 \(C(x)\)，知足
\[ B(x)=A'(x)\ ,\ C(x)=\int A(x) \]

解決方法

直接按照定義作便可，有
\[ B(x)=\sum_{i=1}^{n}ia_ix^{i-1}\ ,\ C(x)=\sum_{i=1}^{n}\frac{a_ix^{i+1}}{i+1} \]
咱們通常認爲 \(C(x)\) 的常數項爲 \(0\) ，複雜度顯然爲 \(\mathcal O(n)\)。

代碼實現

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

inline void Int(int *a, int n) {
  for (int i = n; i; --i) a[i] = 1ll * a[i - 1] * qpow(i, mod - 2) % mod;
  a[0] = 0;
}

多項式 ln

問題描述

給出 \(n\) 次多項式 \(A(x)\)，求一個 \(\bmod{x^{n+1}}\) 下的多項式 \(B(x)\)，知足

\[ B(x) \equiv \ln A(x)\pmod{x^{n+1}} \]

全部運算在模 998244353 下進行。

解決方法

設 \(F(x)=\ln x\)，則 \(B(x)=F(A(x))\) 。

對 \(B(x)\) 求導，根據鏈式法則，有
\[ B'(x)=F'(A(x))A'(x)=\frac{A'(x)}{A(x)} \]
所以對 \(A(x)\) 分別進行求導和求逆，卷積便可求出 \(B'(x)\) ，再對其進行積分便可。

複雜度與多項式求逆同階，爲 \(\mathcal O(n \log n)\) 。

代碼實現

使用 [ Luogu P4725 ] 多項式對數函數 做爲測試題，再也不區分多項式求逆部分的實現方式。

多項式求逆遞部分歸版本，使用 O2 優化，用時 682 ms

inline void Ln(int *a, int *b, int n) {
  Inv(a, b, n); Der(a, n);
  int len = Rev(n << 1);
  NTT(a, len, 1); NTT(b, len, 1);
  for (int i = 0; i < len; ++i) a[i] = 1ll * a[i] * b[i] % mod;
  NTT(a, len, -1); Int(a, n);
}
相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程