首页 > 其他分享 >多项式小寄

多项式小寄

时间：2024-02-24 14:55:47浏览次数：27  
标签：int 多项式 sum Len Complex varphi equiv

多项式小记

目录

目录

序

十分的基础，还有很多东西 没有学习，准备 鸽了 学了 \(DS\) 之后再来补充

这篇博客主讲 \(FFT\)，\(NTT\) 的 基本原理 以及 多项式简单运算，而 题目涉及较少（没做


什么是多项式

即形如 \(F(x) = a_1 x ^ {k_1} + a_2 x ^ {k_2} + ... + a_{k_1} x ^ 1 + a_{k_1 + 1} x ^ 0\) 的 式子，也可以理解为一个 函数

\(OI\) 中一般讨论的 多项式 默认 只有一元，也就是 一个变量 \(x\)

上式可以直接写作 左式 \(F(x)\) 的形式，也可以直接称作 多项式 \(F\)

这都不影响其 实际指代（也就是 右式）

（上述 右式 也可以写为 \(\sum a_k x^{k}\)）

\(eg. F(x) = x ^ 2 + 5x + 2\) 故有 \(F(5) = 52\)


多项式四则运算

加法

存在多项式 \(F(x) = \sum a_k x ^ k, G(x) = \sum b_k x ^ k\)

则有 \(F(x) + G(x) = \sum (a_k + b_k) x ^ k\)

\(eg. F(x) = 2 x ^ 2 + 5x + 4\)

\(G (x) = 4x ^ 3 + 3x + 6\)

则 \(F(x) + G(x) = (0 + 4) x ^ 3 + (2 + 0) x ^ 2 + (5 + 3) x ^ 1 + (4 + 6) x ^ 0 = 4 x ^ 3 + 2 x ^ 2 + 8 x + 10\)


减法

和加法 显然互逆，不多言


乘法

存在多项式 \(F(x) = \sum_i a_i x ^ i, G(x) = \sum_j b_j x ^ j\)

则有

\[F(x) * G (x) = \sum_i\sum_j a_ib_j x^{i + j} \]

\[eg. \\ F(x) = 2x ^ 2 + 5x + 4, \\ G(x) = 4x ^ 3 + 3x + 6 \\ F(x) * G(x) = (2 * 4) x ^ {2 + 3} + (2 * 0) x ^ {2 + 2} + (2 * 3) x ^ {2 + 1} + (2 * 6) x ^ {2 + 0} + (5 * 4) x ^ {1 + 3} + (5 * 0) x ^ {1 + 2} \\ + (5 * 3) x ^ {1 + 1} + (5 * 6) x ^ {1 + 0} + (4 * 4) x ^ {0 + 3} + (4 * 0) x ^ {0 + 2} + (4 * 3) x ^ {0 + 1} + (4 * 6) x ^ {0 + 0} \\ = 8 x ^ 5 + 20 x ^ 4 + 14 x ^ 3 + 27 x ^ 2 + 42 x + 24 \]


除法

和乘法 显然互逆，注意 余式


\(FFT\)

以下内容多引用自 [1] 对应的文章

多项式的四则运算 中，加减法都可以在 线性时间 内很快的解决

而如果按照定义暴力计算，多项式乘法（加法卷积）则需要 \(O(N^2)\) 的时间来解决

非常的慢！

而 \(FFT~(Fast ~ Fourier ~ Transformation)\) ，即 快速傅里叶离散变换，

就是在 \(O(N \log N)\) 时间下解决 多项式乘法 的一种算法


这玩意儿十分的精妙，他利用到一个非常重要的性质，即


多项式的点值表达

就像 二次函数 有 两点式（加上二次项系数确定这个二次函数） 一样

显然，一个 \(N\) 次多项式 可以由其式子 对应函数 上的 \(N + 1\) 个点来确认

所以我们可以把一个 \(N\) 次多项式从 系数表达 转化成 \(N + 1\) 个点来表达，这一步也叫 插值

这有什么好处呢？

可以发现，多项式的四则运算 可以 直接映射到 点坐标的四则运算 上

也就是若 \(F(x)\) 过点 \((x_1, y_1)\) ，而 \(G (x)\) 过点 \((x_2, y_2)\)

则 \(F(x) + G(x)\) 一定过点 \((x_1 + x_2, y_1 + y_2)\)，多项式的乘法 也满足这个性质

（即 \(F(x) * G(x)\) 过点 \((x_1 x_2, y_1 y_2)\)）

而对于 \(N\) 次多项式 \(F(x)\)，\(M\) 次多项式 \(G(x)\)，它们的积最高只有 \(N + M\) 次

也就是我们至多在两个函数上 各找 \(N + M + 1\) 个点运算后就可以确定 \(F(x) * G (x)\)


显然 知道点值表达 后的运算只需要 \(O(N + M)\) 的时间 （\(N, M\) 同阶时认为是 \(O(N)\)）

于是重点就变成了 用系数表达求点值表达 和 用点值表达反推系数表达

亦即 插值 和 求值


\(DFT\)


单位根

按理说 插值 代入 任意点 都没有问题，可是 求取任意点 就只能 代入多项式 做到

这玩意儿求一个点就是 \(O(N)\) 的时间，而我们需要 \(N + M\) 个点，完蛋

于是我们需要一些 有性质的数 来带入求值，而 单位根 就是一种好东西

单位根 / \(n\) 次单位根 ：即 \(n\) 次幂 为 \(1\) 的 复数，用 \(w_n\) 表示

显然 \(n\) 次单位根 共 \(n\) 个，写作 \(w_n ^ k\) 的形式

根据 复数乘法 模长相乘，幅角相加 的法则，发现 单位根 只存在于 半径为 \(1\) 的 单位圆 上

且 \(n\) 次单位根 将此 单位圆 \(n\) 等分

故可以得到 单位根 的如下性质（之后会用到的）

  1. \(w_n ^ 0 = 1 + 0i = 1\)

  2. \(w_n ^ k = (w_n ^ 1) ^ k\)

  3. \(w_n ^ k = w_n ^ {k \% n}\)

  4. \((w_n ^ k) ^ p = w_n ^ {kp} = w_n ^ {kp \% n}\)

  5. \(w_n ^ i * w_n ^ j = w_n ^ {i + j}\)

  6. \(w_{an} ^ {ak} = w_n ^ k \to w_{2n} ^ {2k} = w_n ^ k\)

  7. \(w_n ^ {k + n/2} = - w_n ^ {k} ~ (2 | n)\)

（通过 复数乘法 法则都可以简单证明）


算法原理

现在我们可以回到 \(DFT\) 的原理了

考虑一个 \(n - 1\) 次（\(n\) 项）多项式 \(F(x) = \sum_{k = 0} ^ {n - 1} a_k x^k\)

我们按 次数奇偶 分成两部分，这里保证 \(N\) 是 \(2\) 的整数次幂，从而可以 平均分开

得到

\[F(x) = \sum_{k = 0} ^ {(n - 1) / 2} a_{2k} x ^ {2k} + a_{2k + 1} x ^ {2k + 1} \]

于是设

\[FL(x) = \sum_{k = 0} ^ {(n - 1) / 2} a_{2k} x ^ k \\ FR(x) = \sum_{k = 0} ^ {(n - 1) / 2} a_{2k + 1} x ^ k \]

显然我们可以发现

\[F(x) = FL(x ^ 2) + x FR(x ^ 2) \]

我们代入 \(x = w_n ^ k\)，因为

\[(w_n ^ k) ^ 2 = w_n ^ k * w_n ^ k = w_n ^ {2 * k} = w_{n / 2} ^ k \]

故上式变为

\[① ~ F(w_n ^ k) = FL (w_{n / 2} ^ k) + w_n ^ k FR (w_{n / 2} ^ k) \]

显然 \(k < n / 2\)，故我们再代入 \(x = w_n ^ {k + n / 2}\)，有

\[F(w_n ^ {k + n / 2}) = FL ((w_n ^ {k + n / 2}) ^ 2) + w_n ^ {k + n / 2} FR ((w_n ^ {k + n / 2}) ^ 2) \]

又 \((w_n ^ k) ^ 2 = w_n ^ {2k}\)，有

\[F(w_n ^ {k + n / 2}) = FL (w_n ^ {2k + n}) + w_n ^ {k + n / 2} FR (w_n ^ {2k + n}) \]

又 \(w_n ^ k = w_n ^ {k \% n}\)，有

\[F(w_n ^ {k + n / 2}) = FL (w_n ^ {2k}) + w_n ^ {k + n / 2} FR (w_n ^ {2k}) \]

又 \(w_{2n} ^ {2k} = w_n ^ k\)，有

\[F(w_n ^ {k + n / 2}) = FL (w_{n / 2} ^ k) + w_n ^ {k + n / 2} FR (w_{n / 2} ^ k) \]

又 \(w_n ^ {k + n/2} = - w_n ^ {k} ~ (2 | n)\)，有

\[② ~ F(w_n ^ {k + n / 2}) = FL (w_{n / 2} ^ k) - w_n ^ k FR (w_{n / 2} ^ k) \]

核心流程

把 刚刚的玩意儿 \(②\) 和 上面推出的一个式子 \(①\) 放在一起

\[① ~ F(w_n ^ k) = FL (w_{n / 2} ^ k) + w_n ^ k FR (w_{n / 2} ^ k) \\ ② ~ F(w_n ^ {k + n / 2}) = FL (w_{n / 2} ^ k) - w_n ^ k FR (w_{n / 2} ^ k) \]

艹？长得差不多？只有符号的区别？

非常的好，于是现在我们只需要知道 \(FL, FR\) 在 \(w_{n / 2} ^ {0}, w_{n / 2} ^ {1}, ..., w_{n / 2} ^ {n / 2 - 1}\) 处的 点值表达 即可

而 \(FL, FR\) 的性质与原式 \(F\) 完全相同（右式 形式一致）仅有最高次的区别

可以把它们看作次数 减小为一半 的一般多项式，我们很容易可以想到 分治

通过上式将 原问题 不断划分为规模只有 一半 的 子问题，递归求解，项数为 \(1\) 时 直接返回

最后按上述两式子，即可分别得出 \(F\) 在 \(w_n ^ 0 ... w_n ^ {n / 2 - 1}\) 与 \(w_n ^ {n / 2} ... w_n ^ {n - 1}\) 处的 点值表达 了

时间复杂度 \(T(N) = 2 T(N / 2) + O(N) = O(N \log N)\)

注意到我们需要 多项式项数 为 \(2\) 的 整次幂 从而保证 划分始终均匀

所以当 项数不足 时直接在 最高次补足项 即可，显然 补的项数不超过原项数，复杂度不变


代码

我们可以 重载复数运算 方便后续编写，这里提供一个 重载好的复数结构体

struct Complex {
    double R, I;

    inline Complex (double r = 0, double i = 0) {
        R = r, I = i;
    }
    inline Complex operator + (const Complex &a) const {
        return Complex (R + a.R, I + a.I);
    }
    inline Complex operator - (const Complex &a) const {
        return Complex (R - a.R, I - a.I);
    }
    inline Complex operator * (const Complex &a) const {
        return Complex (R * a.R - I * a.I, I * a.R + R * a.I);
    }
    // 除法是不必要的
};

这里给出一个 按照上述流程 的 \(DFT\) 递归实现，后面会讲优化

注意到由于存在 \(0\) 次项，以及为了后续分治均匀，下标最好从 \(0\) 开始

const double PI = acos (-1);
Complex S[MAXN];

inline void DFT (Complex* p, const int Len) {
	if (Len == 1) return ;
	Complex *l = p, *r = (p + (Len >> 1));
	for (int i = 0; i < Len; ++ i) S[i] = p[i];
	for (int i = 0; i < (Len >> 1); ++ i)
		l[i] = S[i << 1], r[i] = S[i << 1 | 1];
	DFT (l, Len >> 1), DFT (r, Len >> 1);
	Complex U = Complex (cos (2 * PI / Len), sin (2 * PI / Len)); // 第一单位根
	Complex Now = Complex (1, 0); // 当前单位根
	for (int k = 0; k < (Len >> 1); ++ k, Now = Now * U)
		S[k] = l[k] + Now * r[k], S[k + (Len >> 1)] = l[k] - Now * r[k];
	for (int i = 0; i < Len; ++ i) p[i] = S[i];
}

这个函数将多项式 系数数列 转化为 点值数列，初始传入 系数数列 \(p\)

做完后的 \(p\) 数组就是 积多项式 分别在 \(w_n ^ 0 ... w_n ^ {n - 1}\) 处的 点值


\(IDFT\)

插值 是 插完了，我们得到了 一堆 表示 积多项式 的 点

但是这玩意儿它没用啊，我们要的是 系数表达，所以还需要 求值

也就是还得要个东西 把刚刚的流程反过来，得到 最终的多项式（系数表达）

结论

这东西流程上是很简单的，可以先来个 结论

把 \(DFT\) 中的 \(w_n ^ 1\) 换成 \(w_n ^ {-1}\)，做完之后 除以 \(n\) 即可

（可以做题了


证明

以下给出证明，设有 \(G = DFT (F)\)，也就是 \(G\) 数列 是 \(F\) 数列 \(DFT\) 变换后的 结果

模拟暴力代入，（设将 \(w_n ^ k\) 代入以求得第 \(k\) 个点 的 点值表达），则有

（\(G_k\) 表示 \(G\) 多项式的第 \(k\) 项，同形式同理）

\[① ~ G_k = \sum_{i = 0} ^ {n - 1} F_i (w_n ^ k) ^ i \]

我们需要证得 的结论是

\[② ~ F_k = \dfrac {\sum_{i = 0} ^ {n - 1} G_i (w_n ^ {-k}) ^ i} {n} \]

于是将 \(①\) 代入 \(②\) 有

\[F_k = \dfrac {\sum_{i = 0} ^ {n - 1} (w_n ^ {-k}) ^ i (\sum_{j = 0} ^ {n - 1} F_j (w_n ^ i) ^ j)} {n} \\ = \dfrac {\sum_{i = 0} ^ {n - 1} \sum_{j = 0} ^ {n - 1} w_n ^ {-ik} w_n ^ {ij} F_j} {n} \\ = \dfrac {\sum_{i = 0} ^ {n - 1} \sum_{j = 0} ^ {n - 1} w_n ^ {i (j - k)} F_j} {n} \]

分两部分贡献讨论，（最后 两部分相加）若 \(j = k\)，则（有 \(w_n ^ 0 = 1\)）

\[R_1 = \dfrac {\sum_{i = 0} ^ {n - 1} w_n ^ 0 F_j} {n} \\ = \dfrac {\sum_{i = 0} ^ {n - 1} F_j} {n} \\ = F_j = F_k \]

若 \(j \neq k\)，设 \(p = j - k\)，则

\[R_2 = \dfrac {\sum_{i = 0} ^ {n - 1} \sum_{j = 0} ^ {n - 1} w_n ^ {ip} F_j} {n} ~ (j \neq k) \]

即需要证得 \(F_k = R_1 + R_2\)，而我们注意到

\[\sum_{i = 0} ^ {n - 1} w_n ^ {ip} F_j = \sum_{i = 0} ^ {n - 1} w_n ^ {ip} F_{k + p} \\ w_n ^ p \sum_{i = 0} ^ {n - 1} w_n ^ i F_{k + p} \]

等比数列求和有

\[\sum_{i = 0} ^ {n - 1} w_n ^ i F_{k + p} \\ = \dfrac {\sum_{i = 1} ^ {n} w_n ^ i F_{k + p} - \sum_{i = 0} ^ {n - 1} w_n ^ i F_{k + p}} {w_n ^ 1 - 1} \\ = \dfrac {F_{k + p} * (w_n ^ n - w_n ^ 0)} {w_n ^ 1 - 1} \\ = \dfrac {F_{k + p} * 0} {w_n ^ 1 - 1} \\ = 0 \]

所以

\[R_2 = \dfrac {\sum_{i = 0} ^ {n - 1} \sum_{j = 0} ^ {n - 1} w_n ^ {ip} F_j} {n} ~ (j \neq k) \\ = \dfrac {\sum_{j = 0} ^ {n - 1} (\sum_{i = 0} ^ {n - 1} w_n ^ {ip} F_j)} {n} \\ = \dfrac {\sum_{j = 0} ^ {n - 1} w_n ^ p \sum_{i = 0} ^ {n - 1} w_n ^ i F_{k + p}} {n} \\ = \dfrac {\sum_{j = 0} ^ {n - 1} w_n ^ p * 0} {n} \\ = 0 \]

于是最终有 \(R_1 + R_2 = F_k + 0 = F_k\)，即 结论得证


代码

就，和 \(DFT\) 不能说 一模一样，只能说 长得很像

const double PI = acos (-1);
Complex S[MAXN];

inline void DFT (Complex* p, const int Len) {
	if (Len == 1) return ;
	Complex *l = p, *r = (p + (Len >> 1));
	for (int i = 0; i < Len; ++ i) S[i] = p[i];
	for (int i = 0; i < (Len >> 1); ++ i)
		l[i] = S[i << 1], r[i] = S[i << 1 | 1];
	DFT (l, Len >> 1), DFT (r, Len >> 1);
	Complex U = Complex (cos (2 * PI / Len), - sin (2 * PI / Len)); // 第一单位根 (注意取负)
	Complex Now = Complex (1, 0); // 当前单位根
	for (int k = 0; k < (Len >> 1); ++ k, Now = Now * U)
		S[k] = l[k] + Now * r[k], S[k + (Len >> 1)] = l[k] - Now * r[k];
	for (int i = 0; i < Len; ++ i) p[i] = S[i];
}

最后的实现

显然，有了 \(DFT / IDFT\)，最后的乘法是简单的

我们只需要 补齐项数，对最开始的两个多项式分别做 \(DFT\)

然后将求出的 点值数列 相乘，对乘积做 \(IDFT\)

最后输出结果前 \(N + M\) 项即可，这部分代码如下

inline Mul () {

    cin >> N >> M;

    for (P = 1; P <= N + M; P <<= 1);

    for (int i = 0; i <= N; ++ i) cin >> F[i].R;

    for (int i = 0; i <= M; ++ i) cin >> G[i].R;

    FFT (F, P, 1), FFT (G, P, 1);

    for (int i = 0; i < P; ++ i) F[i] = F[i] * G[i];

    FFT (F, P, - 1);

    for (int i = 0; i <= N + M; ++ i)
        cout << (int) (F[i].R / P + 0.49) << ' ';
}

融融

用 \(EPS\) 个 \(zhicheng\) 的注意力 可以发现

\(DFT\) 和 \(IDFT\) 的 代码差距 和 它们的 名字差距 一样巨大

据完全统计，有整整 一个字符

所以为了节省码量，我们可以把它们 绒绒 起来

const double PI = acos (-1);
Complex S[MAXN];

inline void FFT (Complex* p, const int Len, const int Typ) {
	if (Len == 1) return ;
	Complex *l = p, *r = (p + (Len >> 1));
	for (int i = 0; i < Len; ++ i) S[i] = p[i];
	for (int i = 0; i < (Len >> 1); ++ i)
		l[i] = S[i << 1], r[i] = S[i << 1 | 1];
	FFT (l, Len >> 1, Typ), FFT (r, Len >> 1, Typ);
	Complex U = Complex (cos (2 * PI / Len), Typ * sin (2 * PI / Len));
	Complex Now = Complex (1, 0);
	for (int k = 0; k < (Len >> 1); ++ k, Now = Now * U)
		S[k] = l[k] + Now * r[k], S[k + (Len >> 1)] = l[k] - Now * r[k];
	for (int i = 0; i < Len; ++ i) p[i] = S[i];
}
// Use DFT -> FFT (p, Len, 1), Use IDFT -> FFT (p, Len, -1)

感觉可以了？交一波 P3803 【模板】多项式乘法（FFT）

可以看到因为常数过大而 \(TLE\) 了,看来我们需要一个常数更小的写法 —— \(command\_block\)

（显然，\(command\_block\) 老师 没有体验到 洛 谷 4 . 0 的 强 大）

（以上代码可以 轻 松 通 过 ）


于是我们继续来做


没用的 常数优化

观察 如下代码

for (int k = 0; k < (Len >> 1); ++ k, Now = Now * U)
	S[k] = l[k] + Now * r[k], S[k + (Len >> 1)] = l[k] - Now * r[k];

显然 Now * r[k] 这种 复数乘法 常数极大，可以 尝试省去，如下

Complex tmp = (1, 0);
for (int k = 0; k < (Len >> 1); ++ k, Now = Now * U)
	tmp = Now * r[k], S[k] = l[k] + tmp, S[k + (Len >> 1)] = l[k] - tmp;

实测 会快 \(EPS\) ？笔者更倾向于 编译器比较聪明

还是得来点 有用的常数优化


蝴蝶变换

\(S\) 这个临时数组 是 大 常 数 的 一 大 助 力，减少数组复制总能 显著减小常数

我们可以发现，递归下去的过程中 \(p\) 数组 值的位置在变化，具体来讲

0   1   2   3   4   5   6   7 - Flr1
0   2   4   6 | 1   3   5   7 - Flr2
0   4 | 2   6 | 1   5 | 3   7 - Flr3
0 | 4 | 2 | 6 | 1 | 5 | 3 | 7 - Flr4

可以发现，最终位置 是 原位置 的 二进制反转 的值，而这个值可以 \(O(N)\) 递推得到

for (int i = 0; i < N; ++ i) U[i] = (U[i >> 1] >> 1) | ((i & 1) ? N >> 1 : 0);

数组复制，没了，效 果 显 著！


递归变递推（迭代）

显然，递归 令人不喜，如果有能够 递推 的形式那再好不过

我们可以直接尝试把 递归的过程反过来

也就是 最开始 就把数放到 递归的时候对应数最终的位置上，然后 往上合并

inline void FFT (Complex *p, const int Len, const int Typ) {
	for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
	int len = 1;
	Complex Now = Complex (1, 0), tmp, U;
	for (int k = 2; k <= Len; len = k, k <<= 1) {
		U = Complex (cos (2 * PI / k), Typ * sin (2 * PI / k));
		for (int i = 0; i < Len; i += k) {
			Now = Complex (1, 0);
			for (int l = i; l < i + len; ++ l, Now = Now * U)
				tmp = Now * p[l + len], p[l + len] = p[l] - tmp, p[l] = p[l] + tmp;			
		}	
	}
}
// k 即是当前 合并后 的区间长, i 是区间起始处，我们每次将 [l, l + len) 与 [l + len, l + k) 这两个区间合并

进 步 十 分 巨 大，这个东西得配合上面的 蝴蝶变换 来用才行


三步变两步

简单来说，这基于 下面的式子，设存在多项式 \(F, G\)，则

\[P = F + Gi \\ P = F ^ 2 - G ^ 2 + 2 FGi \]

这个东西不很显然，但是可以用 多项式四则运算的性质 与 实数 基本相同 来 感性理解

也就是 多项式乘法 本身具有一般意义上的 分配律

所以我们可以把多项式 \(F, G\) 分别 放到 \(p\) 数组的 实部，虚部

做一遍 \(DFT\)，然后自身乘方后做 \(IDFT\)，我们需要的是 \(F*G\)

于是最后输出 虚部 的 \(\dfrac {1}{2}\) 即可，这东西就是 不严格的 \(\dfrac {1}{3}\) 巨量提升

   inline void Mul () {
    cin >> N >> M;
       
    for (P = 1; P <= N + M; P <<= 1);
       
    for (int i = 0; i < P; ++ i) U[i] = (U[i >> 1] >> 1) | ((i & 1) ? P >> 1 : 0);

    for (int i = 0; i <= N; ++ i) cin >> F[i].R;

    for (int i = 0; i <= M; ++ i) cin >> F[i].I;

    FFT (F, P, 1);

    for (int i = 0; i < P; ++ i) F[i] = F[i] * F[i];

    FFT (F, P, - 1);

    for (int i = 0; i <= N + M; ++ i)
        cout << (int) (F[i].I / P / 2 + 0.49) << ' ';
}

破除封建迷信

谁说自己重载 \(Complex\) 会快的？标准库的轮子 还是 有实力的

所以我们可以把

struct Complex {
	...
};

改成（当然直接修改 每个地方的定义 也行）

#define Complex complex <double>

然后修改对应的 输入输出（这玩意儿好像不支持直接 cin）

// Input Part

int x;

for (int i = 0; i <= N; ++ i) cin >> x, F[i] += x;
	
for (int i = 0; i <= M; ++ i) cin >> x, F[i] += Complex (0, x);

// Output Part

for (int i = 0; i <= N + M; ++ i) cout << (int) (F[i].imag () / P / 2 + 0.499) << ' '; // 注意 .imag ()

好了，现在再交交看，可以发现 std::complex 十 分 强 大


最终代码

2024-02-23 17:33 - C++20 O2 - 1.24 KB / 952 ms / 40.54 MB | Record

#include <bits/stdc++.h>
#define Complex complex <double>

const int MAXN = 2100005;
const double PI = acos (-1);

using namespace std;

int U[MAXN];

inline void FFT (Complex *p, const int Len, const int Typ) {
	for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
	int len = 1;
	Complex Now = Complex (1, 0), tmp, U;
	for (int k = 2; k <= Len; len = k, k <<= 1) {
		U = Complex (cos (2 * PI / k), Typ * sin (2 * PI / k));
		for (int i = 0; i < Len; i += k) {
			Now = Complex (1, 0);
			for (int l = i; l < i + len; ++ l, Now = Now * U)
				tmp = Now * p[l + len], p[l + len] = p[l] - tmp, p[l] = p[l] + tmp;			
		}	
	}
}

Complex F[MAXN];

int N, M, P, x;

int main () {
	
	ios::sync_with_stdio(0);
	cin.tie(0), cout.tie(0);
	
	cin >> N >> M;
	
	for (P = 1; P <= N + M; P <<= 1);
	
	for (int i = 0; i < P; ++ i) 
        U[i] = (U[i >> 1] >> 1) | ((i & 1) ? P >> 1 : 0);
	
	for (int i = 0; i <= N; ++ i) cin >> x, F[i] += x;
	
	for (int i = 0; i <= M; ++ i) cin >> x, F[i] += Complex (0, x);
	
	FFT (F, P, 1);
	
	for (int i = 0; i < P; ++ i) F[i] = F[i] * F[i];
	
	FFT (F, P, - 1);
	
	for (int i = 0; i <= N + M; ++ i)	
        cout << (int) (F[i].imag () / P / 2 + 0.499) << ' ';
	
	return 0;
}

\(NTT\)

这东西还真只能接着 \(FFT\) 讲，没看 \(FFT\) 先看 \(FFT\)


咋来的

\(FFT\) 这玩意儿好吧？但是要用到 \(double\) 这 破东西

这就不可避免地带来了 精度 的问题，数据范围一大就 令人恼火

但是 单位根 这东西吧，它不好替代，至少

数学家们已经证明，在 复数域 中，单位根是 唯一一类 满足要求的数 —— \(command\_block\)

但是吧，大部分 \(OI\) 中的计数题都是在 模意义下的，于是咱对 ”这类数“ 的限制放宽了

在 模意义同余系 中，我们找到了 单位根 的 替代品，原根 \(g\)


原根

得讲讲这玩意儿是啥...

以下内容多引用自 [3] 对应的文章


定义

如果 \(a\), \(m\) 互质，那么我们把 \(a ^ k \equiv 1 ~ (mod ~ m)\) 的 最小正整数 \(k\) 称为 \(a\) 的模 \(m\) 的阶

在 原文 中老师将前提 写成了 \(a, m\) 互质

虽然这篇博客对笔者以下内容编写有很大作用，但...火大

于是显然，如果 \(a ^ n \equiv 1 ~ (mod ~ m)\)，则 \(k | n\)（余数的可乘性）

由 欧拉定理 可知，当 \(a, m\) 互质时，$a ^ {\varphi (m)} \equiv 1 ~ (mod ~ m) $

故 \(a\) 模 \(m\) 的阶 \(k\) 为 \(\varphi (m)\) 的因数，而

我们把满足 \(k = \varphi (m)\) 对应的 \(a\) 称为 模 \(m\) 的原根

\(eg. m = 3\)，有 \(\varphi (3) = 2\)，且 \(2\) 模 \(3\) 的阶 \(k = 2\)，故 \(2\) 是 \(3\) 的 原根


称之为 根，因为其是 \(x ^ {\varphi (m)} \equiv 1 ~ (mod ~ m)\) 的 一个解

称之为 原，因为 原神 因为


性质

先来点用不到的
  1. （上面提到过） \(a\) 模 \(m\) 的阶 \(k\) 为 \(\varphi (m)\) 的因数

  1. \(g\) 是模 \(m\) 的 原根 当且仅当 对于任何与 \(m\) 互质的数 \(a\) 都能找到 \(k\)，使得 \(a \equiv g ^ k ~ (mod ~ m)\)

显然由于 1. 如果 \(g\) 不是 模 \(m\) 的 原根，则 \(g\) 的 阶 \(k\) 小于 \(\varphi (m)\)

又根据阶的定义，则 \(g ^ k \equiv 1 ~ (mod ~ m)\)，故 \(g ^ {t} \equiv g ^ {t + ik} ~ (mod ~ m)\)

即由 \(g ^ i\) 生成的 模 \(m\) 意义下互不同余 的数至多 \(k\) 个，小于 \(\varphi (m)\) 个，必要性得证

而若 \(g\) 为 模 \(m\) 的 原根，则 \(g, g ^ 2, ..., g ^ {\varphi (m)}\) 除以 \(m\) 的余数 互不相同

反之若存在 \(g ^ i \equiv g ^ j ~ (mod ~ m)\)，不妨设 \(i < j\)，由于 \(g, m\) 互质，\(g ^ i\) 与 \(m\) 必然互质

根据 同余的除法性质，即 同余式两边可以同时除以 \(a\) 的条件是 \(a, m\) 互质

我们得到 \(g ^ {j - i} \equiv 1 ~ (mod ~ m)\)，与 原根 的 阶 \(k = \varphi (m)\) 矛盾

故 \(g, g ^ 2, ..., g ^ {\varphi (m)}\) 除以 \(m\) 的余数 互不相同 为 真命题，同上可知

\(g, g ^ 2, ... , g ^ {\varphi (m)}\) 均与 \(m\) 互质，而与 \(m\) 互质且小于 \(m\) 的数一共 \(\varphi (m)\) 个

又知 \(g, g ^ 2, ..., g ^ {\varphi (m)}\) 互不相同，故其一定为 所有小于 \(m\) 且与 \(m\) 互质的数 的 一种排列


  1. 设 \(g\) 是模 \(m\) 的一个原根，则 \(g ^ k\) 是原根 当且仅当 \(k, \varphi (m)\) 互质

这也是 原根的判定定理，可由此得出 原根的个数 以及延伸出 原根的存在性定理

对 必要性 考虑 反证法，因为 \(g ^ k\) 为 原根

则根据 2. \((g ^ k) ^ 1, (g ^ k) ^ 2, ... , (g ^ k) ^ {\varphi (m)}\) 模 \(m\) 互不相同，而若 \(gcd (k, \varphi (m)) = t ~ (t \neq 1)\)

则 \((g ^ {k}) ^ {\varphi (m) / t} \equiv (g ^ {\varphi (m)}) ^ {k / t} \equiv 1 ~ (mod ~ m)\)，显然 \(\varphi (m) / t < \varphi (m)\)

此处与 原根 的 阶 为 \(\varphi (m)\) 矛盾，故 必要性得证

充分性则考虑到若 \(k, \varphi (m)\) 互质，则总能找到 \(s, t\) 使得 \(sk + t\varphi (m) = 1\)

故可以有 \(g ^ i \equiv g ^ {i (sk + t \varphi (m))} ~ (mod ~ m)\)

则 \(g ^ i \equiv g ^ {isk} g ^ {it\varphi (m)} ~ (mod ~ m)\)，显然 \(g ^ {it \varphi (m)} \equiv g ^ {\varphi (m)} \equiv 1 ~ (mod ~ m)\)

故 \(g ^ i \equiv g ^ {isk} \equiv (g ^ k) ^ {is}\)，也就是 \(g ^ i\) 可以表示成 \(g ^ k\) 的 \(is\) 次方

而根据 2. \(g ^ i\) 可以生成所有与 \(m\) 互质的数，故 \(g ^ k\) 亦可生成，故 \(g ^ k\) 也为原根

故 如果模数 \(m\) 有原根存在，则模数 \(m\) 有 \(\varphi (\varphi (m))\) 个原根

（对于原根 \(g\)，满足 \(i, \varphi (m)\) 互质的 \(g ^ i\) 都是原根）


  1. 对于质数 \(p\)，模 \(p\) 的 原根 \(g\) 存在，且有 \(\varphi (\varphi (p)) = \varphi (p - 1)\) 个

这是 3. 的推广，同时会用到下面性质

4.1. 如果整数 \(a\) 模 \(p\) 的 阶 是 \(t\)，则 当且仅当 \(k, t\) 互质时，\(a ^ k\) 的 阶 也是 \(t\)

对 必要性 同样考虑 反证法，因为 \(a\) 模 \(p\) 的 阶 为 \(t\)，故 \(a ^ t \equiv 1 ~ (mod ~ p)\)

显然 \(a ^ {it} \equiv 1 ~ (mod ~ p)\)，又设若 \(k, t\) 不互质，且存在 公因数 \(c > 1\)

则有 \(t | kt / c\)，也就是 存在 \(i\) 使得 \(it = kt / c\)，故 \(a^{ {kt} / {c}} \equiv 1 ~ (mod ~ p)\)

即 \((a ^ k) ^ {t/c} \equiv 1 ~ (mod ~ p)\)，显然 \(t/c < t\)，又 \(t\) 为 \(a ^ k\) 的 阶，矛盾，故 必要性得证

充分性 也可以利用上述结论 证得，考虑因为 \(a\) 模 \(p\) 的 阶 是 \(t\)

则 根据阶的定义，显然，对于一切 \(a ^ i \equiv 1 ~ (mod ~ p)\)，存在 \(t | i\)

我们设 \(a ^ k\) 存在 阶 \(c\)，即有 \((a ^ k) ^ c \equiv a ^ {kc} \equiv 1 ~ (mod ~ p)\)

\(k, t\) 互质，即 \(lcm (k, t) = kt\)，即不存在数 \(c\) 使得 \(c < t\) 且 \(t | ck\)

故 \(c \ge t\)，又显然 \(c = t\) 时，\((a ^ k) ^ c \equiv 1 ~ (mod ~ p)\)，故此时 \(c ~ (t)\) 即是 \(a ^ k\) 模 \(p\) 的 阶

充分性得证


4.2. 对于 \(n\) 次 同余方程 \(f(x) \equiv 0 ~ (mod ~ p)\)，根 至多 有 \(n\) 个

考虑 数学归纳法，一次同余方程 \(ax + b \equiv 0 ~ (mod ~ p)\) 至多一个根

注意到这里 \(p\) 是 质数，故满足 一次同余方程唯一根 的前提，即 \(a, p\) 互质

考虑 \(n\) 次同余方程 \(f(x) \equiv 0 ~ (mod ~ p)\)

若 方程无根，命题成立，反之则 不妨设一根 \(t\)，代入则有

\(f (t) = \sum_i ^ n a_i t ^ i\)，又 原式 \(\sum _ i ^ n a_i x ^ i\)，两式相减有

\(f(x) = \sum _ i ^ n a _ i (x ^ i - t ^ i)\)，显然有 公因式 \((x - t)\)，故原式可化为

\(f (x) = (x - t) g(x)\)，其中 \(g(x)\) 为一个一般的 \(n - 1\) 次多项式

即 \(n\) 次多项式 \(f\) 至多 比 \(n - 1\) 次多项式 \(g\) 多出一个根 \(t\)，由于 一次同余方程 仅 一根

故 原命题得证


4.3. 如果 \(d | p - 1\)，则 同余方程 \(x ^ d \equiv 1 ~ (mod ~ p)\) 恰有 \(d\) 个根

考虑 \(p\) 为质数，故 费马小定理 可知 \(x ^ {p - 1} \equiv 1 ~ (mod ~ p)\)

又由于 \(4.2.\)，该 \(d\) 次 同余方程 \(x ^ d - 1 \equiv 0 ~ (mod ~ p)\) 至多 \(d\) 个根

只需证明其 下界为 \(d\)

又考虑到 同余方程 \(x ^ {p - 1} \equiv 1 ~ (mod ~ p)\) 显然存在 \(p - 1\) 个根

即是 正整数域 内所有 \(p\) 的 可能余数

设 \(p - 1 = qd\)，显然 \(q\) 为 正整数，根据 代数恒等式（等比数列求和逆向）

\(x ^ {p - 1} - 1 = (x ^ d) ^ q - 1 = (x ^ d - 1) ((x ^ d) ^ {q - 1} + (x ^ d) ^ {q - 2} + ... + (x ^ d) ^ 1 + 1)\)

注意到这里 原文老师 又 笔误 将 \(x ^ d - 1\) 写成 \(x ^ q - 1\) 了

即 \(x ^ {p - 1} - 1 = (x ^ d - 1) g (x)\)

而 \(g (x)\) 为 \(dq - d\) 次 多项式，根据 4.2.，其 至多贡献 \(dq - d\) 个根

又 \(x ^ {p - 1} - 1 \equiv 0 ~ (mod ~ p)\) 共有 \(p - 1 = dq\) 个根

故 \(x ^ d - 1 \equiv 0 ~ (mod ~ p)\) 至少贡献 \(d\) 个根，即其 根的数量下界为 \(d\)，原命题得证


4.4. 若 \(a\) 的 阶 为 \(k\)，则 同余方程 \(x ^ k \equiv 1 ~ (mod ~ p)\) 的 所有根 恰为 \(a, a ^ 2, ... a ^ k\)

\(p\) 为质数，显然 \(a, p\) 互质，根据 4.3.，同余方程 \(x ^ k \equiv 1 ~ (mod ~ p)\) 恰有 \(k\) 个根

且显然 \(a ^ k \equiv (a ^ 2) ^ k \equiv ... \equiv (a ^ k) ^ k \equiv 1 ~ (mod ~ p)\)

故只需证明 \(a \% p, (a ^ 2) \% p, ... , (a ^ k) \% p\) 互不相等 即可

考虑 反证法，若存在 \(a ^ i \equiv a ^ j ~ (mod ~ p)，1 \le i < j \le k\)

由于显然 \(a ^ i\) 与 \(p\) 互质，故 根据同余的除法性质 \(a ^ {j - i} \equiv 1 ~ (mod ~ p)\)

又考虑 \(j - i < k\)，而 \(a\) 的 阶 为 \(k\)，矛盾，故原命题得证


由以上四个引理，我们可以证明 4.

考虑我们记 \(n (t)\) 表示 满足 模 \(p\) 的 阶 为 \(t\) 且与 \(p\) 互质 的 \(a\) 的 个数

（显然与 \(p\) 互质只是摆设）

而显然，对于所有与 \(p\) 互质的数 \(a\)，\(a\) 模 \(p\) 的 阶 存在且唯一，故 \(\sum n(t) = \varphi (p) = p - 1\)

唯一性 根据 阶的定义可知，存在性则，由于 \(p\) 是质数，根据 费马小定理 可构造

由于 \(p\) 是质数，\(\varphi(p) = p - 1\)，故根据 1.，对于所有 \(t\)，存在 \(t | \varphi (p) = p - 1\)

则显然 \(\sum_{t | (p - 1)} n(t) = \sum n(t)\ = varphi (p) = p - 1\)

而 满足 \(gcd (p - 1, a) = t\) 的 \(a\) 一定共有 \(\varphi ((p - 1) / t)\) 个（\(gcd ((p - 1)/t, a/t) = 1\)）

故 \(\sum_{t | p - 1} \varphi ((p - 1) / t) = \sum_{t | p - 1} \varphi (t) = p - 1\)

即 \(\sum_{t | p - 1} n (t) = \sum_{t | p - 1} \varphi (t)\)

而对于 \(t | p - 1\) 的任意 \(t\)

① 若 \(n(t) = 0\)，则显然 \(n (t) < \varphi (t)\)

② 反之则存在 模 \(p\) 的 阶 为 \(t\) 且与 \(p\) 互质的 \(a\)（此处 \(a\) 代指其一而非整体）

根据 4.4.，\(x ^ t \equiv 1 ~ (mod ~ p)\) 恰有 \(t\) 个解，且为 \(a, a ^ 2, ..., a ^ t\)

又考虑 模 \(p\) 的 阶 为 \(t\) 且与 \(p\) 互质的数 均为同余方程的根，也就是一定在 \(a, a ^ 2, ..., a ^ t\) 中

而根据 4.1.，这样的数应当为 满足 \(i, t\) 互质的 \(a ^ i\)，显然共有 \(\varphi (t)\) 个

而 模 \(p\) 的 阶 为 \(t\) 且与 \(p\) 互质的数 即是 \(n(t)\)，则表明这种情况下 \(n(t) = \varphi(t)\)

故有 当 \(n(t) = 0\) 时，\(n(t) < \varphi(t)\)，而其余情况 \(n(t) = \varphi(t)\)

又知 \(\sum_{t | p - 1} n (t) = \sum_{t | p - 1} \varphi (t)\)

故显然 \(n(t) = 0\) 的情况 不存在（否则 \(\sum_{t | p - 1} n (t) < \sum_{t | p - 1} \varphi (t)\)）

显然我们知道，原根的个数 就是 模 \(p\) 的 阶 为 \(\varphi (p) = p - 1\) 的数的个数，即 \(n(p - 1)\)

由于 \(n(p - 1) \neq 0\) 故此情况下 原根一定存在

又 \(n(p - 1) = \varphi (p - 1) = \varphi (\varphi (p))\) 故此情况下 原根个数为 \(\varphi (\varphi (p))\)，命题得证


回到正题

以下内容部分引用自 [2] 对应文章


来看看在 \(FFT\) 中用到了 单位根 \(w\) 的哪些性质？

  1. \(w_n ^ k = (w_n ^ 1) ^ k\)
  2. \(w_n ^ k = w_n ^ {k \% n}\)
  3. \((w_n ^ k) ^ p = w_n ^ {kp} = w_n ^ {kp \% n}\)
  4. \(w_n ^ i * w_n ^ j = w_n ^ {i + j}\)
  5. \(w_{an} ^ {ak} = w_n ^ k \to w_{2n} ^ {2k} = w_n ^ k\)
  6. \(w_n ^ {k + n/2} = - w_n ^ {k} ~ (2 | n)\)

显然，原根可以在 \(n = p - 1\) 时满足以上性质（通过 余数的性质 易得）

考虑 \(n\) 次单位根 就是把 单位圆 平均分 \(n\) 分

原根 这玩意儿吧就是 在模意义下 把 “单位圆” 平均分 \(p - 1\) 份

然后 模数 也同样具有相似的 可乘性 啥的，就，感性理解是容易的

简单推广得到当 \(n | p - 1\) 的时候也可以满足（几块并成一块）

但是啊，这个 \(p\) 这玩意儿不是 少得可怜？咋做？

仔细想想，发现 由于要 补项，\(n\) 通常情况下可以表示为 \(2 ^ i\) 的形式，也少得可怜

好家伙，于是我们只需要找质数 \(p\) 使得 \(p = 2 ^ k * m + 1\)，然后 \(k\) 巨大就行了

这里的 \(2 ^ k\) 其实决定了 可以做的长度的最大值，然后 我褐了 \(0.125\) 个表（来源 [4]）


\(p = 2 ^ k * m + 1\) \(k\) \(g\)
\(998244353\) \(23\) \(3\)
\(1004535809\) \(21\) \(3\)
\(2013265921\) \(27\) \(31\)
\(1945555039024054273\) \(56\) \(5\)
\(4179340454199820289\) \(57\) \(3\)

对于一般题，找到 一个 合适的模数 的 原根 \(g\)，那么 \(g ^ {(p - 1) / n}\) 就是一个 单位根 的完美替代

（可以感性理解成 \(g _ n ^ {(p - 1) / n} = w _ n ^ 1\)）


实现

简单替换一下，直接开搞就行

我们先拉出刚刚 \(FFT\) 的 最终版代码...的 \(FFT\) 函数

inline void FFT (Complex *p, const int Len, const int Typ) {
	for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
	int len = 1;
	Complex Now = Complex (1, 0), tmp, U;
	for (int k = 2; k <= Len; len = k, k <<= 1) {
		U = Complex (cos (2 * PI / k), Typ * sin (2 * PI / k));
		for (int i = 0; i < Len; i += k) {
			Now = Complex (1, 0);
			for (int l = i; l < i + len; ++ l, Now = Now * U)
				tmp = Now * p[l + len], p[l + len] = p[l] - tmp, p[l] = p[l] + tmp;			
		}	
	}
}

在 \(FFT\) 的基础上稍做修改

const int MOD = 998244353, g = 3, ig = 332748118;

inline void NTT (int *p, const int Len, const int Typ) {
	for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
	int len = 1, Now = 1, tmp, u;
	for (int k = 2; k <= Len; len = k, k <<= 1) {
		u = Qpow (Typ == 1 ? g : ig, (MOD - 1) / k);
		for (int i = 0; i < Len; i += k) {
			Now = 1;
			for (int l = i; l < i + len; ++ l, Now = 1ll * Now * u % MOD)
				tmp = 1ll * Now * p[l + len] % MOD, p[l + len] = (p[l] - tmp + MOD) % MOD, p[l] = (p[l] + tmp) % MOD;			
		}	
	}
}

最后 除以 \(n\) 的时候注意 求一下 \(n\) 的逆元 即可，直 接 就 过 了

当然 这玩意儿 是 不优的，显然在 取模 上面 优化空间很大

我们可以用 \(unsigned ~ long ~ long\) 扩大存数范围

减少取模次数

这有一些效果

inline void NTT (unsigned long long *p, const int Len, const int Typ) {
	for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
	int len = 1, Now = 1, tmp, u;
	for (int k = 2; k <= Len; len = k, k <<= 1) {
		u = Qpow (Typ == 1 ? g : ig, (MOD - 1) / k);
		for (int i = 0; i < Len; i += k) {
			Now = 1;
			for (int l = i; l < i + len; ++ l, Now = 1ll * Now * u % MOD)
				tmp = 1ll * Now * p[l + len] % MOD, p[l + len] = (p[l] - tmp + MOD), p[l] = (p[l] + tmp);			
		}
		if (k == (1 << 16)) for (int i = 0; i < Len; ++ i) p[i] %= MOD;	
	}
	for (int i = 0; i < Len; ++ i) p[i] %= MOD;	
}

预处理单位根

也是 十分不错的方法

inline void NTT (unsigned long long *p, const int Len, const int Typ) {
	for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
	int len = 1, tmp, u;
	W[0] = 1;
	for (int k = 2; k <= Len; len = k, k <<= 1) {
		u = Qpow (Typ == 1 ? g : ig, (MOD - 1) / k);
		for (int i = 1; i <= len; ++ i) W[i] = W[i - 1] * u % MOD;
		for (int i = 0; i < Len; i += k) 
			for (int l = i; l < i + len; ++ l)
				tmp = 1ll * W[l - i] * p[l + len] % MOD, p[l + len] = (p[l] - tmp + MOD), p[l] = (p[l] + tmp);			
		if (k == (1 << 16)) for (int i = 0; i < Len; ++ i) p[i] %= MOD;	
	}
	for (int i = 0; i < Len; ++ i) p[i] %= MOD;	
}

最后也放一份比较完整的

板子

2024-02-24 14:39 - C++20 O2 - 1.63 KB / 981 ms / 50.55 MB | Record

#include <bits/stdc++.h>
#define Complex complex <double>

const int MAXN = 2100005;
const int MOD  = 998244353;

using namespace std;

int U[MAXN];

inline int Qpow (const int a, const int b = MOD - 2) {
	if (b == 1) return a;
	int Ret = Qpow (a, b >> 1);
	if (b & 1) return 1ll * Ret * Ret % MOD * a % MOD;
	return 1ll * Ret * Ret % MOD;
}

const int g = 3, ig = Qpow (g);

unsigned long long W[MAXN];

inline void NTT (unsigned long long *p, const int Len, const int Typ) {
	for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
	int len = 1, tmp, u;
	W[0] = 1;
	for (int k = 2; k <= Len; len = k, k <<= 1) {
		u = Qpow (Typ == 1 ? g : ig, (MOD - 1) / k);
		for (int i = 1; i <= len; ++ i) W[i] = W[i - 1] * u % MOD;
		for (int i = 0; i < Len; i += k) 
			for (int l = i; l < i + len; ++ l)
				tmp = 1ll * W[l - i] * p[l + len] % MOD, p[l + len] = (p[l] - tmp + MOD), p[l] = (p[l] + tmp);			
		if (k == (1 << 16)) for (int i = 0; i < Len; ++ i) p[i] %= MOD;	
	}
	for (int i = 0; i < Len; ++ i) p[i] %= MOD;	
}

unsigned long long F[MAXN], G[MAXN];

int N, M, P, x, in;

int main () {
	
	ios::sync_with_stdio(0);
	cin.tie(0), cout.tie(0);
	
	cin >> N >> M;
	
	for (P = 1; P <= N + M; P <<= 1);
	
	for (int i = 0; i < P; ++ i) U[i] = (U[i >> 1] >> 1) | ((i & 1) ? P >> 1 : 0);
	
	for (int i = 0; i <= N; ++ i) cin >> F[i];
	
	for (int i = 0; i <= M; ++ i) cin >> G[i];
	
	NTT (F, P, 1), NTT (G, P, 1);
	
	for (int i = 0; i < P; ++ i) F[i] = 1ll * F[i] * G[i] % MOD;
	
	NTT (F, P, - 1), in = Qpow (P);
	
	for (int i = 0; i <= N + M; ++ i) cout << (1ll * F[i] * in) % MOD << ' ';
	
	return 0;
}

多项式基本运算

先鸽一会儿...


引用 & 鸣谢

这些文章都使笔者学到 褐了 很多内容，而且 质量较高，十分推荐


[1] 傅里叶变换(FFT)学习笔记 - command_block 的博客 - 洛谷博客 (luogu.com.cn)

[2] NTT与多项式全家桶 - command_block 的博客 - 洛谷博客 (luogu.com.cn)

[3] 原根的概念、性质及其存在性 - 知乎 (zhihu.com) - Mathis Wang（笔 误 の 神）

[4] FFT用到的各种素数 | Miskcoo’s Space

标签：int,多项式,sum,Len,Complex,varphi,equiv
From： https://www.cnblogs.com/FAKUMARER/p/18031086

相关文章

  • 多项式初步
    多项式初步目录多项式初步自己写的分治FFT/NTTPart1分治FFT/NTTPart2①.多项式求逆：②.多项式带余除法：③.多项式开根：④.多项式对数：⑤.多项式exp：⑥.多项式快速幂：模板基础操作MTT自己写的分治FFT/NTTPart1给定序列\(g_{1\dotsn-1}\)，求序列\(f_{0\dotsn-1}\)。其中......
  • 数学：多项式
    拉格朗日插值快速傅里叶变换（FFT）已经不知道被这玩意劝退了多少次了。但理解之后其实不算很阴间的东西？前置知识多项式负数单位根快速傅里叶变换快速傅里叶逆变换快速数论变换（NTT/FNTT）前置知识FFT注意到FFT的运算到处都是又\(\cos\)又\(\sin\)的，有不小的精度问题......
  • 【多项式】任意模数 NTT/FTT
    现在有两个整数多项式\(A\)，\(B\)，\(0\lea_i,b_i\le10^9\)，\(n\le10^5\)，求它们的卷积，同时系数对\(p\le10^9\)取模。我们会发现，最终的系数可能会达到\(10^5\times10^9\times10^9=10^{23}\)级别，FFT会爆longdouble的精度；NTT，由于模数无特殊性质，完全不能使用。接......
  • 【多项式】【模版】FFT NTT
    多项式若\(A(x)=a_0+a_1x+a_2x^2+\dotsa_nx^n(a_n\ne0)\)，则称\(A\)是\(n\)次多项式。初中概念，但在OI中可以玩出花来。多项式的表示方式像上面的介绍一样，用系数\(a_0,a_1,\dotsa_n\)来表示多项式的方法称为系数表示法。还有一种表示多项式的方法，就是对于\(n\)......
  • 【模板】多项式全家桶（多项式初等函数（部分））
    【模板】多项式初等函数同时作为https://github.com/caijianhong/template-poly的document。杂项数域为\(\mathbbF_{998244353}\)，所以定义了mint为modint<998244353>。poly是多项式的类型，从std::vector<mint>继承而来。poly的构造函数如下：poly();explicitpoly(......
  • 多项式从入门到进门
    多项式全家福多项式一个以\(x\)为变量的多项式定义在一个代数域\(F\)上，将函数\(A(x)\)表示为形式和：\[A(x)=\sum\limits_{i=0}^{n-1}a_ix^i\]我们称\(a_0,a_1,\cdots,a_n\)为多项式的系数，所有系数都属于数域\(\mathbbF\)。如果一个多项式的最高次的非零系数是\(k......
  • P4512 【模板】多项式除法
    为什么不能直接用\(F(x)\timesG(x)^{-1}\)做？\(G(x)^{-1}\)是模\(x^{m+1}\)意义下的，\(n-m>m\)时得到的\(Q(x)\)就是错的。记\(F'(x)\)为\(F(x)\)系数翻转后的多项式，即若\(F(x)=\sum\limits_{i=0}^nf_ix^i\)，则\(F'(x)=\sum\limits_{i=0}^nf_{n......
  • 多项式初等函数
    前置知识：多项式基本操作（求导、积分、FFT）和一些数学知识。ReferenceOI-wiki多项式牛顿迭代给定多项式\(g(x)\)，已知有多项式\(f(x)\)满足\(g(f(x))\equiv0\pmod{x^n}\)，求出模\(x^n\)意义下的\(x^n\)。假设我们已经求得了模\(x^{\lceil\frac{n}{2}\rceil}\)意义......
  • 特征多项式学习笔记
    介绍了方阵的特征多项式以及利用上Hessenberg矩阵的\(\mathcal{O}(n^3)\)求法。ReferenceOI-wiki特征多项式：Hessenberg法及加速矩阵幂特征值与特征向量给定\(n\timesn\)的方阵\(\mathbf{T}\)，若存在一个非零列向量\(\mathbf{v}\)和数\(\lambda\)满足\(\mathbf{T}......
  • 多项式求逆
    类似于逆原，多项式也可以求逆，具体地，定义多项式\(f(x)\)的乘法逆为能使得\(f(x)*g(x)\equiv1\pmod{x^n}\)的多项式\(g(x)\)，也可记作\(f(x)^{-1}\)。假设我们现在已经求出了\(f(x)\)在模\(x^{\lceil\fracn2\rceil}\)意义下的逆原\(f_0(x)^{-1}\)。有：\[\left.\begi......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99