快速數論變換

簡介

數論變換（number-theoretic transform, NTT）是離散傅里葉變換（DFT）在數論基礎上的實現；快速數論變換（fast number-theoretic transform, FNTT）是快速傅里葉變換（FFT）在數論基礎上的實現。

數論變換 是一種計算卷積（convolution）的快速算法。最常用算法就包括了前文提到的快速傅里葉變換。然而快速傅立葉變換具有一些實現上的缺點，舉例來説，資料向量必須乘上覆數係數的矩陣加以處理，而且每個複數係數的實部和虛部是一個正弦及餘弦函數，因此大部分的係數都是浮點數，也就是説，必須做複數而且是浮點數的運算，因此計算量會比較大，而且浮點數運算產生的誤差會比較大。

NTT 解決的是多項式乘法帶模數的情況，可以説有些受模數的限制，數也比較大。目前最常見的模數是 998244353。

前置知識

學習數論變換需要前置知識：離散傅里葉變換、生成子羣、原根、離散對數。相關知識可以在對應頁面中學習，此處不再贅述。

定義

數論變換

在數學中，NTT 是關於任意環上的離散傅立葉變換（DFT）。在有限域的情況下，通常稱為數論變換（NTT）。

數論變換（NTT）是通過將離散傅立葉變換化為 \(F={\mathbb {Z}/p}\)，整數模質數 \(p\)。這是一個 有限域，只要 \(n\) 可除 \(p-1\)，就存在本原 \(n\) 次方根，所以我們有 \(p=\xi n+1\) 對於正整數 \(ξ\)。具體來説，對於質數 \(p=qn+1, (n=2^m)\)，原根 \(g\) 滿足 \(g^{qn} \equiv 1 \pmod p\), 將 \(g_n=g^q\pmod p\) 看做 \(\omega_n\) 的等價，則其滿足相似的性質，比如 \(g_n^n \equiv 1 \pmod p, g_n^{n/2} \equiv -1 \pmod p\)。

因為這裏涉及到數論變化，所以 \(N\)（為了區分 FFT 中的 \(n\)，我們把這裏的 \(n\) 稱為 \(N\)）可以比 FFT 中的 \(n\) 大，但是隻要把 \(\frac{qN}{n}\) 看做這裏的 \(q\) 就行了，能夠避免大小問題。

常見的有：

\[ p = 167772161 = 5 \times 2^{25}+1, g=3 \]

\[ p = 469762049 = 7 \times 2^{26}+1, g=3 \]

\[ p = 754974721 = 3^2 \times 5 \times 2^{24}+1, g=11 \]

\[ p = 998244353 = 7 \times 17 \times 2^{23}+1, g=3 \]

\[ p = 1004535809 = 479 \times 2^{21}+1, g=3 \]

就是 \(g^{qn}\) 的等價 \(\mathrm{e}^{2\pi n}\)。

迭代到長度 \(l\) 時 \(g_l = g^{\frac{p-1}{l}}\)，或者 \(\omega_n = g_l = g_N^{\frac{N}{l}} = g_N^{\frac{p-1}{l}}\)。

快速數論變換

快速數論變換（FNTT）是數論變換（NTT）增加分治操作之後的快速算法。

快速數論變換使用的分治辦法，與快速傅里葉變換使用的分治辦法完全一致。這意味着，只需在快速傅里葉變換的代碼基礎上進行簡單修改，即可得到快速數論變換的代碼。

在算法競賽中常提到的 NTT 一詞，往往實際指的是快速數論變換，一般默認「數論變換」是指「快速數論變換」。

這樣簡寫的邏輯與快速傅里葉變換相似。事實上，「快速傅里葉變換」（FFT）一詞指的是「快速離散傅里葉變換」（FDFT），但由於「快速」只能作用於離散，甚至是本原單位根階數為 \(2\) 的冪的特殊情形，不能作用於連續，因此「離散」一詞被省略掉，FDFT 變為 FFT，即 FFT 永遠指的是特殊的離散情形。

數論變換或快速數論變換是在取模意義下進行的操作，不存在連續的情形，永遠是離散的，自然也無需提到離散一詞。

在算法領域，不進行提速的操作是無意義的。在快速傅里葉變換中介紹 DFT 一詞，是因為 DFT 在信號處理、圖像處理領域也有其他的具體應用，同時 DFT 也是 FFT 的原理或前置知識。

在不引起混淆的情形下，常用 NTT 來代指 FNTT。為了不引起下文進一步介紹的混淆，下文的 NTT 與 FNTT 兩個詞進行了分離。

DFT、FFT、NTT、FNTT 的具體關係是：

在 DFT 與 NTT 的基礎上，增加分治操作，得到 FFT 與 FNTT。分治操作的辦法與原理，可以參見快速傅里葉變換一文。
在 DFT 與 FFT 的基礎上，將複數加法與複數乘法替換為模 \(p\) 意義下的加法和乘法，一般大小限制在 \(0\) 到 \(p-1\) 之間；將本原單位根改為模 \(p\) 意義下的相同階數的本原單位根，階數為 \(2\) 的冪，即可得到 NTT 與 FNTT。

由於替換的運算只涉及加法和乘法，因此 DFT、FFT、NTT、FNTT 擁有相同的原理，均在滿足加法與乘法的環上進行，無需域上滿足除法運算的更加嚴格的條件。

事實上，只要擁有原根，即羣論中的生成元，該模數下的 NTT 或 FNTT 即可進行。考慮到模數為 \(1\)、\(2\) 和 \(4\) 的情形太小，不具有實際意義，對於奇素數 \(p\) 和正整數 \(\alpha\)，只要給出模數為 \(p^\alpha\) 和 \(2p^\alpha\) 的原根 \(g\)，採用同樣的辦法，則 NTT 或 FNTT 仍然可以進行。

模板

下面是一個大數相乘的模板，參考來源。

參考代碼

#include <algorithm>
#include <bitset>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <iomanip>
#include <iostream>
#include <map>
#include <queue>
#include <set>
#include <string>
#include <vector>
using namespace std;

int read() {
  int x = 0, f = 1;
  char ch = getchar();
  while (ch < '0' || ch > '9') {
    if (ch == '-') f = -1;
    ch = getchar();
  }
  while (ch <= '9' && ch >= '0') {
    x = 10 * x + ch - '0';
    ch = getchar();
  }
  return x * f;
}

void print(int x) {
  if (x < 0) putchar('-'), x = -x;
  if (x >= 10) print(x / 10);
  putchar(x % 10 + '0');
}

const int N = 300100, P = 998244353;

int qpow(int x, int y) {
  int res(1);
  while (y) {
    if (y & 1) res = 1ll * res * x % P;
    x = 1ll * x * x % P;
    y >>= 1;
  }
  return res;
}

int r[N];

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;
  }
}

int A[N], B[N], C[N];

char a[N], b[N];

int main() {
  int i, lim(1), n;
  scanf("%s", &a);
  n = strlen(a);
  for (i = 0; i < n; ++i) A[i] = a[n - i - 1] - '0';
  while (lim < (n << 1)) lim <<= 1;
  scanf("%s", &b);
  n = strlen(b);
  for (i = 0; i < n; ++i) B[i] = b[n - i - 1] - '0';
  while (lim < (n << 1)) lim <<= 1;
  for (i = 0; i < lim; ++i) r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  ntt(A, lim, 1);
  ntt(B, lim, 1);
  for (i = 0; i < lim; ++i) C[i] = 1ll * A[i] * B[i] % P;
  ntt(C, lim, -1);
  int len(0);
  for (i = 0; i < lim; ++i) {
    if (C[i] >= 10) len = i + 1, C[i + 1] += C[i] / 10, C[i] %= 10;
    if (C[i]) len = max(len, i);
  }
  while (C[len] >= 10) C[len + 1] += C[len] / 10, C[len] %= 10, len++;
  for (i = len; ~i; --i) putchar(C[i] + '0');
  puts("");
  return 0;
}

參考資料與拓展閲讀

本页面最近更新：，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：ChungZH, Yukimaikoriya, tigerruanyifan, isdanni, Saisyc, 383494, Tiphereth-A
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用