快速数论变换（ntt）

CaF3

已于 2025-01-20 15:22:57 修改

阅读量1.1k

点赞数 14

分类专栏：笔记算法文章标签：算法

于 2025-01-20 15:18:35 首次发布

本文链接：https://blog.youkuaiyun.com/qq_54180412/article/details/145260080

版权

笔记同时被 2 个专栏收录

13 篇文章

订阅专栏

算法

6 篇文章

订阅专栏

一、快速傅里叶变换的缺陷

1.单位根的计算

由于 $w_n=e^{\frac{2\pi i}{n}}=\cos\frac{2\pi}{n}+i\cdot\sin\frac{2\pi}{n}$ 中带有三角函数，计算会产生浮点数，带来误差。

2.单位根的性质

fft之所以能在 $O (n l o g n)$ 的时间内进行dft，是因为它通过蝴蝶变换将递归层数优化为了 $O (l o g n)$ 层，而要实现蝴蝶变换，就要满足以下三个性质：

$w_n^0=w_n^n=1$
$w_n^k=-w_n^{k+n/2}$
$w_n^k=w_{2n}^{2k}$

所以，如果我们能找到一个可以避免浮点数运算且具有这三个性质的东西，那么我们就能用它来替换单位根，以此来避免误差。

二、原根

1.阶

若a与p互素且p>1，对于满足 $a^n\equiv 1\pmod{p}$ 的最小的n，我们称之为a模p的阶，记作 $\delta_p(a)$ 。

例如： $2^1\equiv 2\pmod{7},2^2\equiv 4\pmod{7},2^3\equiv 1\pmod{7}$ ，故 $\delta_7(2)=3$ 。

2.原根

假设p是正整数，a是整数，a与p互素。若 $\delta_p(a)=\phi(p)$ ，则称a为模p的一个原根。原根有很多性质，这里不做赘述，具体可以参考文章原根的概念、性质及其存在性。

例如： $\delta_7(3)=6=\phi(7)$ ，故3是模7的一个原根。

定理：若p为素数，g为p的原根，那么 $g^i\bmod p(1<g<p,0<i<p)$ 的结果两两不同。

由定理可得， $\because (g^\frac{p-1}{2})^2\equiv g^{p-1}\equiv 1\pmod{p}$ 且 $g^\frac{p-1}{2}\not\equiv g^{p-1}\pmod{p}$ ， $\therefore g^\frac{p-1}{2}\equiv -1\pmod{p}$ 。

3.替换

假设p为素数，g为p的一个原根。
记 $w_n\equiv g^{\frac{p-1}{n}}\pmod{p}$

证明：

$w_n^0\equiv g^0\equiv 1\pmod{p},w_n^n\equiv g^{p-1}\equiv 1\pmod{p}$
$w_n^{k+n/2}\equiv g^{\frac{(p-1)(k+n/2)}{n}}\equiv g^\frac{(p-1)k}{n}\cdot g^\frac{p-1}{2}\equiv-w_n^k\pmod{p}$
显然。

由此可得，原根具有蝴蝶变换所需的三个性质，故可以替换单位根。

三、代码

实际应用时，一般取p=998244353，g=3。

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

// Problem: P3803 【模板】多项式乘法（FFT）
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P3803
// Memory Limit: 500 MB
// Time Limit: 2000 ms
// 
// Powered by CP Editor (https://cpeditor.org)

#include <bits/stdc++.h>
using namespace std;

#define rep(n) for (int _i_i_ = (n); _i_i_; --_i_i_)
#define mem(arr, x) memset((arr), x, sizeof (arr))
#define mkp make_pair
#define edge(y, x, ver) for (int _i_ = h[x], y = ver[_i_]; _i_; y = ver[_i_ = nxt[_i_]])
#define import mt19937
#define random rnd(chrono::system_clock::now().time_since_epoch().count())
#define setdb(v) setprecision(v)
using ll = long long;
using ull = unsigned long long;
using pii = pair<int, int>;
using pll = pair<ll, ll>;

const int maxn = 1e6 + 5;
const int mod = 998244353, g = 3, gi = 332748118;

int n, m;

int pw(int x, int n) {
	int res = 1;
	for (; n; n >>= 1) {
		if (n & 1) res = 1ll * res * x % mod;
		x = 1ll * x * x % mod;
	}
	return res;
}

void ntt(vector<int> &x, int t) {
	int n = x.size(), l = 31 - __builtin_clz(n);
	vector<int> r(n);
	for (int i = 0; i < n; ++i) {
		r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
		if (i < r[i]) swap(x[i], x[r[i]]);
	}
	for (int k = 1; k < n; k <<= 1) {
		int wn = pw(t == 1 ? g : gi, (mod - 1) / (k << 1));
		for (int i = 0; i < n; i += k << 1) {
			int w = 1;
			for (int j = 0; j < k; ++j, w = 1ll * w * wn % mod) {
				int x1 = x[i + j], y1 = 1ll * x[i + j + k] * w % mod;
				x[i + j] = (x1 + y1) % mod;
				x[i + j + k] = (x1 - y1 + mod) % mod;
			}
		}
	}
}

void solve() {
	cin >> n >> m;
	vector<int> a(n + 1), b(m + 1);
	for (auto &x : a) cin >> x;
	for (auto &y : b) cin >> y;
	int p = 1;
	while (p <= m + n) p <<= 1;
	rep (p - n - 1) a.emplace_back(0);
	rep (p - m - 1) b.emplace_back(0);
	ntt(a, 1);
	ntt(b, 1);
	for (int i = 0; i < p; ++i) a[i] = 1ll * a[i] * b[i] % mod;
	ntt(a, -1);
	int inv = pw(p, mod - 2);
	for (int i = 0; i <= n + m; ++i) cout << 1ll * a[i] * inv % mod << ' ';
}

int main() {
	ios_base::sync_with_stdio(false);
    cin.tie(0);
    cout.tie(0);
	int t = 1;
	//cin >> t;
	rep (t) {
		solve();
	}
	return 0;
}

/*
+---------------------------------------------+
|GREY|GREEN|CYAN|BLUE|PURPLE|YELLOW|ORANGE|RED|
+---------------------------------------------+
*/