素数
即 质数，巨大多内容板块
全文 绝大多数 内容是对 [0] 中讲述的 粗略抄写 和 胡乱加工
1. 基本概念
定义
如果数 \(p\) 有且仅有 \(1\) 和 \(p\) 两个因子，这个数就是 素数（质数）
反之即为 合数，显然，每个 \(>1\) 的整数要么是 素数，要么是 合数
特别的，规定 \(1\) 既不是素数，也不是合数
性质
-
任何 正整数 \(n\) 均可表示为 若干素数 的 乘积
根据定义容易证明，若一个数 有多个因数，可将其 分解为若干因数的乘积
直到分解到 每个因数均不可再分解，此时的每个因数 即为素数
-
算术基本定理：正整数 \(n\) 有且仅有 一种 分解成 若干素数乘积 的方式（不记顺序）
考虑 反证法，若存在数 \(n\) 有 至少两种 分解方式，且 \(n\) 是 最小的 这样的数
即 \(n = \prod p ^ {a} = \prod q ^ b ~ (\exist ~p \neq q)\)
显然，\(p_1 \mid \prod q ^ b\)，则 \(q_1, q_2, ..., q_m\) 中至少存在一个 \(q_i\) 使得 \(p_1 \mid q_i\)
不妨设 \(q_i\) 就是 \(q_1\)，于是 \(p_1 \mid q_1\)，又知道 \(p, q\) 均为 素数，故 \(p = q\)
而若 \(a_1 > b_1\)，则 \(p_1 ^ {a_1 - b_1} \times p_2 ^ {a_2} \times ... \times p_n ^ {a_n} = q_2 ^ {b_2} \times q_3 ^ {b_3} \times ... \times q_n ^ {b_n}\)
于是 \(q_2, q_3, ..., q_m\) 中至少存在一个 \(q_i\) 使得 \(p_1 \mid q_i\)，但 \(p_1 = q_1 \neq q_i ~ (i \in [2, m])\)
故矛盾，同理当 \(a_1 < b_1\) 时 亦矛盾，则仅有 \(a_1 = b_1\) 这种可能
于是 \(n' = \prod \limits _ {i = 1} ^ n p_i ^ {a_i} = \prod \limits _ {i = 1} ^ m q_i ^ {b_i}\)，我们可以将上述过程 推广到 \(n, m\) 来证明
但其实由于我们 认为 \(n\) 是 最小的 这样的数，但显然 \(n' < n\)，故 矛盾
则 原命题已经得证
-
第 \(n\) 个 素数 \(P_n\) 大致大小为 自然对数 的 \(n\) 倍，即 \(P_n \sim n \ln n\)
-
大小不超过 \(x\) 的 素数个数 \(\pi (x) \sim \dfrac x {\ln x}\)
\(3./4.\) 的 证明十分的硬，此处不表
梅森数
和 素数 有些关系，梅森数 是 形如
\[2 ^ p - 1 \]的数，其中 \(p\) 是 素数
显然，当 \(p\) 为 合数 时，\(2 ^ p - 1\) 一定 不是素数，因为其至少有 \(2 ^ m - 1\) 这个因子
\[2 ^ {km} - 1 = (2 ^ m - 1) (2 ^ {m (k - 1)} + 2 ^ {m (k - 2)} + ... + 1) \]
但当 \(p\) 为 素数 时，\(2 ^ p - 1\) 并 不一定是素数，例如 \(2 ^ {11} - 1 = 2047 = 23 \times 89\) 就假了
2. 互素关系
定义
当 \(\gcd (n, m) = 1 ~ (n, m \in \Z)\) 时，我们就称它们是 互素（互质）的
\[n, m \in \Z, \gcd (n, m) = 1 \Longleftrightarrow n \perp m \]\(Stern-Brocot\) 树
可以用于构造所有 最简分数（\(\dfrac m n ~ (m \perp n)\)）的集合
其核心操作 十分简单，考虑从 \((\dfrac 0 1, \dfrac 1 0)\) 两个分数出发，重复若干次 下面这个操作
- 向 相邻两分数 \(\dfrac m n, \dfrac {m'} {n'} ~ (\dfrac m n < \dfrac {m'} {n'})\) 之间插入 \(\dfrac {m + m'} {n + n'}\) 这个 “中位分数”
这里就 忽略分母为 \(0\) 这个问题 了
图片来源于 [1]
这相当于 一棵无尽的二叉树，我们 构造的序列 即是其 中序遍历 得到的序列
其中 每一个结点都代表一个分数，每一个分数都可以表示成 \(\dfrac {m + m'} {n + n'}\) 的 形式
其中 \(\dfrac m n\) 是 在它左边离它最近的祖先，\(\dfrac {m'} {n'}\) 是 在它右边离它最近的祖先
如何证明这种构造能 不重不漏的 得到 所有的 最简分数 呢？
首先，根据我们 小学奥数 中的 糖水原理，容易证明 两分数 生成的 中位分数 满足
\[\dfrac m n < \dfrac {m + m'} {n + n'} < \dfrac {m'} {n'} \]故 不可能重复，而 由于我们每次生成的数 均严格在两数之间，故 每次递归区间变小
容易理解，当 无尽的递归下去，区间将变得 充分小，不会漏掉
同理，它应当能得到 所有的 符合条件的数，我们只需论证为什么其得到的一定是 最简分数
先考虑来证明一个 引理
若 \(ad - bc = 1\)，则 \(a \perp b, c \perp d\)（\(a, b, c, d \in \Z\)）
仍然考虑反证，若存在 \(\gcd (a, b) = k > 1\)，则有 \(a = ka', b = kb'\)
于是 \(k a'd - k b'c = 1\)，即 \(k (a'd - b'c) = 1\)，又 \(k > 1, a, b, c, d \in \Z\)，显然矛盾
故 \(\gcd (a, b) = 1\)，即 \(a \perp b\)，而 \(c, d\) 与 \(a, b\) 显然对称，故原命题得证
根据上述引理，我们容易推导出，若 \(mn' - m'n = 1\)，则 \(\dfrac m n, \dfrac {m'} {n'}\) 为 最简分数
即 \(n \perp m, n' \perp m'\)
于是我们只需论证 对于序列上任意相邻两数，均满足 \(mn' - m'n = 1\) 即可
显然 \(\dfrac 0 1, \dfrac 1 0\) 满足该性质，于是我们可以考虑 向下递推
设有 \(\dfrac m n, \dfrac {m'} {n'}\) 已满足条件，我们只需证明 \(\dfrac m n, \dfrac {m' + m} {n' + n}\) 与 \(\dfrac {m' + m} {n' + n}, \dfrac {m'} {n'}\) 均满足条件即可
于是有 \(m (n' + n) - n (m' + m) = (m' + m) n' - (n' + n) m' = mn' - m'n = 1\)
故原命题得证，即 核心操作 得到的一定是 最简分数
\(Farey\) 序列
与 \(Stern-Brocot\) 树 相辅相成，互相统一，我们将第 \(i\) 个 \(Farey\) 序列记作 \(F_i\)
\(Farey\) 序列，即 法里序列，也有称 \(F_i\) 为 阶为 \(i\) 的法里级数 的说法，本质相同
其表示 分母小于等于 \(i\) 的 最简真分数 按 从小到大 顺序排列形成的集合
图片来源于 [1]
容易发现，每个 \(Farey\) 序列就是 \(Stern-Brocot\) 的 一棵子树，故 构造是简单的
且显然的性质是，对于序列中任意 连续的三个数 \(\dfrac a b, \dfrac c d, \dfrac e f\)，其满足 \(c = a + e, d = b + f\)
Luogu P8058 [BalkanOI2003] Farey 序列
求第 \(n\) 个 \(Farey\) 序列中的第 \(k\) 项
考虑直接构造复杂度是 \(O (N ^ 2)\) 的，不可承受
注意到 \(Farey\) 序列具有 单调性，容易想到 二分答案 \(x\)
于是我们只需要求在第 \(N\) 个 \(Farey\) 序列上，小于 \(x\) 的分数有多少个即可，记这个个数为 \(W\)
由 \(Farey\) 序列的定义，容易知道二分 上下界 \(x \in [0, 1], x \in \Q\)
我们定义 \(f_i\) 为 小于 \(x\) 的 分数 个数，当 不要求最简 时，显然 \(f_i = \lfloor ix \rfloor\)
若我们 要求最简，则容易得出 \(f_i = \lfloor ix \rfloor - \sum \limits _ {d \mid i} f_d\)，而 \(W = \sum \limits _ {i = 1} ^ n f_i\)
由于每次求和均从 \(1\) 到 \(n\)，故后面这部分可以 \(O (N \log N)\) 预处理系数（筛一下）
然后就结束力，预处理后 二分单次 \(Check\) 可以做到 \(O (N)\)，故总复杂度 \(O (N (\log V + \log N))\)
#include <bits/stdc++.h>
const int MAXN = 40005;
const double EPS = 1e-10;
using namespace std;
int F[MAXN], C[MAXN];
int N, K;
inline void Init () {
for (int i = 1; i <= N; ++ i) C[i] = 1;
for (int i = N; i >= 1; -- i)
for (int j = i * 2; j <= N; j += i)
C[i] -= C[j];
}
inline bool Check (const double M) {
int Cnt = 0;
for (int i = 1; i <= N; ++ i) Cnt += int (M * i) * C[i];
return Cnt < K;
}
inline void Solve (const double x) {
int Ans;
for (int i = 1; i <= N; ++ i) {
Ans = ceil (i * x);
if (fabs (1.0 * Ans / i - x) <= EPS)
return cout << Ans << ' ' << i << '\n', void ();
}
}
double L = 0, R = 1, M, A;
int main () {
cin >> N >> K;
Init ();
while (R - L >= EPS) {
M = (L + R) / 2;
if (Check (M)) L = A = M;
else R = M;
}
Solve (A);
return 0;
}
[ABC333G] Nearest Fraction
双倍经验 Luogu P1298 最接近的分数 但是可以暴力过掉，注意多 判断分子
调了半天，就不应该尝试 卡精度，老老实实用 __int128 和 分数类 就行了
这个题又 完全不卡时间
思路是 简单的，就是在 \(Stern-Brocot\) 树上二分
但是如果直接 暴力一层一层地跳，就会 \(TLE / RE\) 掉（\(RE\) 是因为 递归层数过多）
比如 \(1e-18\) 这个数，就是 一直向左跳，而每次分母只增加 \(1\)，就会挂掉
故考虑 倍增，每次尝试往 一个方向跳到底 再换方向跳
考虑 \(Stern-Brocot\) 树的性质，一个结点 由 左右离他最近的祖先 分子分母 分别求和得到
容易发现，每次换向，分母就会变成 其父亲和其父亲的父亲 的和，成 \(Fib\) 数列
故可以证明，这样的换向应当是 \(O (\log N)\) 次的（分母最大为 \(N\)）
加上倍增，总时间复杂度就是 \(O (\log ^ 2 N)\)，甚至 可以多组询问
#include <bits/stdc++.h>
const __int128 INF = 10000000000000000000ull;
using namespace std;
struct Fract {
__int128 x, y;
inline Fract operator + (const Fract &a) const {
return {x + a.x, y + a.y};
}
inline Fract operator - (const Fract &a) const {
return {x * a.y - a.x * y, y * a.y};
}
inline bool operator < (const Fract &a) const {
return x * a.y < y * a.x;
}
inline bool operator > (const Fract &a) const {
return x * a.y > y * a.x;
}
inline __float128 Val () const {
return (__float128) 1.0 * x / y;
}
};
long long N;
inline Fract Mul (const Fract &a, const long long x) {
return {a.x * x, a.y * x};
}
inline Fract ABS (const Fract &a) {
if ((a.x > 0 && a.y < 0) || (a.x < 0 && a.y > 0)) return {- a.x, a.y};
return {a.x, a.y};
}
inline Fract Que (const Fract val, Fract L = {0, 1}, Fract R = {1, 0}) {
Fract M = L + R, l = L, r = R, x, y;
while (1) {
M = L + R, x = R;
if (M.y > N) return ABS (L - val) > ABS (R - val) ? R : L;
for (int i = 33; ~ i; -- i) {
y = x, x = x + Mul (L, 1ll << i);
if (x.y > N || x < val) x = y;
}
r = x, x = L;
for (int i = 33; ~ i; -- i) {
y = x, x = x + Mul (R, 1ll << i);
if (x.y > N || x > val) x = y;
}
l = x, L = l, R = r;
}
}
Fract val;
long double v;
int main () {
cin >> v >> N;
val = {(unsigned long long) (v * INF), INF};
Fract Ans = Que (val);
cout << (unsigned long long) Ans.x << ' ' << (unsigned long long) Ans.y << '\n';
return 0;
}
SP26073 DIVCNT1 - Counting Divisors
定义 \(d (n)\) 为 \(n\) 的 正整数因数个数，多次询问 \(\sum \limits _ {i = 1} ^ {n} d (i)\)
下文中部分 图片 / 思路 来源于 [3]
先考虑 对式子转化，令 \(B = \left \lfloor \sqrt N \right \rfloor\)，于是有
\[\begin {aligned} & \sum \limits _ {i = 1} ^ {n} { d(i) } \\ =& \sum \limits _ {i = 1} ^ {n} { \left \lfloor \dfrac n i \right \rfloor } \\ =& \sum \limits _ {i = 1} ^ {n} { \sum \limits _ {j = 1} ^ {n} { \left [ i \times j \le n \right ] } } \\ =& \sum \limits _ {i = 1} ^ {B} { \sum \limits _ {j = 1} ^ {n} { \left [ i \times j \le n \right ] } } + \sum \limits _ {i = 1} ^ {n} { \sum \limits _ {j = 1} ^ {B} { \left [ i \times j \le n \right ] } } + \sum \limits _ {i = B + 1} ^ {n} { \sum \limits _ {j = B + 1} ^ {n} { \left [ i \times j \le n \right ] } } - \sum \limits _ {i = 1} ^ {B} { \sum \limits _ {j = 1} ^ {B} { \left [ i \times j \le n \right ] } } \\ =& \sum \limits _ {i = 1} ^ {B} { \sum \limits _ {j = 1} ^ {n} { \left [ i \times j \le n \right ] } } + \sum \limits _ {i = 1} ^ {n} { \sum \limits _ {j = 1} ^ {B} { \left [ i \times j \le n \right ] } } - \sum \limits _ {i = 1} ^ {B} { \sum \limits _ {j = 1} ^ {B} { \left [ i \times j \le n \right ] } } \\ =& 2 \sum \limits _ {i = 1} ^ {B} { \left \lfloor \dfrac n i \right \rfloor } - B ^ 2 \end {aligned} \]我们要求的就是 \(\sum \limits _ {i = 1} ^ B {\left \lfloor \dfrac n i \right \rfloor}\)，相当于 反比例函数 \(y = \dfrac n x\) 内侧，在 \(y \in [1, B]\) 段的 整点数
直接做是 \(O (B)\) 的，多组询问，似得很透彻
考虑用一些 很新的方式 来数这个整点，比如尝试用 一些向量 来 拟合 这个 反比例函数
容易发现，对于反比例函数上一点，斜率最接近该点导数 的 整向量 \((u, v)\)（向量两端点均为整点）
若保证 一端点为该点，另一端点在 反比例函数外侧
靠反比例函数内侧 包含且仅包含了 反比例函数内侧 所有整点（见上图）
于是我们不断地 构造这样最贴近的向量 即可
先考虑 一个这样的向量 带来的贡献，显然就是 一个矩形 和 一个三角形
于是贡献为 \(x v + \dfrac {(u - 1) (v + 1)} {2}\)
可以认为 最上方一列 算成 上面一个向量的贡献，不能重复计算，故第一部分不是 \(x (v + 1)\)
现在只需要考虑 构造 的方法了，显然，我们应当按照如下的 大致思路
- 选定 起始点 和 起始斜率
- 不断 递增斜率，找到 最贴近当前点 \((x, y)\) 的向量
- 将当前点 移动到 选中向量的右端
- 重复上述步骤直到 终点
分步考虑，选定起点 显然 十分的简单，对于一个 反比例函数，我们需求 \(y \in [1, B]\) 的部分
于是起始 \(x = y = B\) 即可，斜率 显然就是 \(- B\)
为了 不漏算，显然最后我们拟合的图像上的整点 均不在原反比例函数内侧
故具体而言，此处起点 \(x = \left \lfloor \dfrac n B \right \rfloor, y = B + 1\)，而 为了方便，斜率均采用 绝对值后的值
中间两步是 重点，对于 当前点 和 当前向量，为了 贴近原函数
显然我们先 执行多次 \(3.\) 操作，直到 下一次移动会使当前点移动到原函数内侧
这时候开始找 新的向量，我们可以使用一个 单调栈 来维护 斜率递减的向量集
这里的 斜率 是 绝对值后的，故称为 递减
我们每次 取出栈顶，若该线段 仍会使得下一次移动后当前点到原函数内侧，则 将其弹出
于是我们最终能找到 两条向量，使得其另一端 恰好跨过原函数两侧，如下图
接着我们就可以 依据这两个向量，来找到最贴近 **原函数 **的 新向量 了
考虑用 \(Stern-Brocot\) 树 的 构造方式，得到向量 \(M = (L_x + R_x, L_y + R_y)\)
容易发现，若从 当前点 起始，向量 \(M\) 的 另一端 在 原函数外侧
则我们可以 弹出 \(R\)，加入 \(M\)，重复找 \(M\) 的上述操作
反之我们 弹出 \(L\)，加入 \(M\)，重复找 \(M\) 的上述操作
考虑 什么时候终止 呢？如果 \(R\) 的斜率 大于 原函数在 \(x + x_L\) 点的斜率
显然无论如何分下去，我们都 不可能再得到 一个 \(M\) 使得其另一端在原函数外侧了
显然 \(L\) 另一端在内侧，而 \(R\) 斜率比 原函数 在 \(L\) 另一端斜率 大
故加上 \(R\) 后得到的 \(M\) 只会 在原函数内侧，离原函数更远（相较于 \(L\)）
这时候我们就 停止二分，而 最后一次合法的 \(M\) 就是我们想要的
这一部分可以 配合上图理解
最后我们需要来考虑一下 终点 的问题
考虑到 后期函数斜率很小，基本都会是 \(\dfrac 1 k\) 的形式，相当于 \(Stern-Brocot\) 上 最靠右的几列
容易发现，此时 层数 与 \(k\) 十分接近，即 二分 这一步的 时间复杂度 将退化到 \(O (k)\)
或说 \(O (y)\)
不太可承受，于是考虑 这一部分可以直接暴力做掉，而 这一部分 我们取 \(O (\sqrt [3] n)\)
据可靠消息，这样询问一次时间复杂度就是 \(O (\sqrt [3] {n} \log n)\) 的
据可靠消息 的意思就是 笔者不会证明
然后再看一下这个 原题 的数据范围
完全能过
注意各种
long long和大量__int128的使用，记得重载输出
#include <bits/stdc++.h>
const int MAXN = 100005;
using namespace std;
int T;
struct Writer {
template <typename T> Writer& operator << (T x) {
if (x == 0) { putchar('0'); return *this; }
if (x < 0) putchar('-'), x = -x;
static short sta[40], top = 0;
while (x > 0) sta[++ top] = x % 10, x /= 10;
while (top > 0) putchar(sta[top] + '0'), top --;
return *this;
}
inline Writer& operator << (char c) {putchar(c); return *this;}
Writer() {}
} out;
struct Vector {
int x, y;
inline Vector operator + (const Vector &a) const {
return {x + a.x, y + a.y};
}
};
long long N;
inline bool Checkin (const long long x, const long long y) {
return (__int128) x * y <= N;
}
inline bool Checked (const long long x, const Vector v) {
return (__int128) N * v.x <= (__int128) v.y * x * x; // Abs
}
vector <Vector> S;
inline void Solve () {
cin >> N;
long long b = sqrt (N), c = cbrt (N), x = N / b, y = b + 1;
__int128 Ans = 0;
Vector L, R, M;
S.push_back ({1, 0}), S.push_back ({1, 1});
while (1) {
for (L = S.back (), S.pop_back (); !Checkin (x + L.x, y - L.y); x += L.x, y -= L.y)
Ans += (__int128) x * L.y + (__int128) (L.y + 1) * (L.x - 1) / 2;
if (y <= c) break ;
for (R = S.back (); Checkin (x + R.x, y - R.y); S.pop_back (), R = S.back ()) L = R;
while (1) {
M = L + R;
if (!Checkin (x + M.x, y - M.y)) S.push_back (R = M);
else if (Checked (x + (L = M).x, R)) break ;
}
}
for (int i = 1; i < y; ++ i) Ans += N / i;
out << ((Ans << 1) - (b * b)) << '\n';
}
int main () {
cin >> T;
while (T --) Solve ();
return 0;
}
3. 素数筛法
\(Sieve ~ of ~ Prime\)
十分重要，从 普及 - 到 NOI，许多数论题 都要求处理 \([1, N]\) 的 所有素数
埃式筛法
\(Sieve ~ of ~ Eratosthenes\)
全称 埃拉托斯特尼筛法，十分经典，时间复杂度 \(O (N \ln \ln N)\)
其核心在于 用已经筛出的素数 的整倍标记 合数，代码就很简单
bool Vis[];
inline void Prime (const int N, vector <int> &Pri) {
for (int i = 2; i <= N; ++ i)
if (!Vis[i]) {
Pri.push_back (i);
for (int k = i; k <= N; k += i)
Vis[k] = 1;
}
}
比较 精妙 的点在于其 复杂度的证明，即 \(\sum \limits _ {p \in \mathbb P, p \le n} \dfrac n p = O (n \ln \ln n)\)
注意到 每个数会被它 所有质因数 筛到一遍，故可以有
\[\begin {aligned} \sum \limits _ {p \in \mathbb P, p \le n} \dfrac n p =& \sum \limits _ {i = 1} ^ {n} \omega (i) \end {aligned} \]以下证明 来自 djq 老师，并不保证经过笔者理解或验证
根据 \(d (i)\) 的计算式有
\[\begin {aligned} \sum \limits _ {i = 1} ^ {n} { 2 ^ {\omega (i)} } \le& \sum \limits _ {i = 1} ^ {n} { d (i) } = O (n \ln n) \end {aligned} \]根据 \(2 ^ x\) 的 凸性 和 琴生不等式 有
\[\begin {aligned} \sum \limits _ {i = 1} ^ {n} { 2 ^ {\omega (i)} } \ge & 2 ^ { \frac { \sum _ {i = 1} ^ {n} { \omega (i) } } {n} } n \\ \therefore { 2 ^ { \frac { \sum _ {i = 1} ^ {n} { \omega (i) } } {n} } } \le & O (\ln n) \end {aligned} \]于是两边 同时取对数，有
\[\begin {aligned} \frac { \sum _ {i = 1} ^ {n} { \omega (i) } } {n} \le & O (\log \ln n) \\ \therefore \sum \limits _ {i = 1} ^ {n} { \omega (i) } =& O (n \log \ln n) = O (n \log \log n) = O (n \ln \ln n) \end {aligned} \]这已经很快了，但是 不够快，它甚至 卡不过去板子（真的不是因为板子是 线性筛 吗？
于是注意到下列优化（
-
我们其实只需要 筛到根号，于是
bool Vis[]; inline void Prime (const int N, vector <int> &Pri) { for (int i = 2; i * i <= N; ++ i) if (!Vis[i]) for (int k = (i << 1); k <= N; k += i) Vis[k] = 1; for (int i = 2; i <= N; ++ i) if (!Vis[i]) Pri.push_back (i); } -
第二重循环，可以从 \(i ^ 2\) 开始，因为更小的前面已经筛过了
bool Vis[]; inline void Prime (const int N, vector <int> &Pri) { for (int i = 2; i * i <= N; ++ i) if (!Vis[i]) for (int k = i * i; k <= N; k += i) Vis[k] = 1; for (int i = 2; i <= N; ++ i) if (!Vis[i]) Pri.push_back (i); } -
提前筛出 偶数，然后第二重循环可以 成倍加
bool Vis[];
inline void Prime (const int N, vector <int> &Pri) {
for (int i = 4; i <= N; i += 2) Vis[i] = 1;
for (int i = 3; i * i <= N; ++ i)
if (!Vis[i])
for (int k = i * i; k <= N; k += (i << 1))
Vis[k] = 1;
for (int i = 2; i <= N; ++ i)
if (!Vis[i]) Pri.push_back (i);
}
上述方法任选，均可通过板子
欧拉筛法
\(Sieve ~ of ~ Euler\)
其本质仍然是 用质数标合数，但是在 埃式筛 中，一个合数可能会被 多个质数 筛到
这就造成了 非线性的复杂度，而在 欧拉筛 中，我们可以保证 一个合数只被筛到一次
这样就保证了复杂度
具体而言，我们尝试 只用每个合数的最小质因数 筛掉它
于是我们对于每个 \(i\)，从小到大 遍历 已知质数 \(p_j\)，并标记合数 \(i \times p_j\)
当 \(p_j \mid i\) 时，停止遍历，由于我们 从小到大 遍历的质数，故此时 \(p_j\) 一定是 \(i\) 的 最小质因子
复杂度 \(O(n)\)，十分正确
bool Vis[];
inline void Prime (const int N, vector <int> &Pri) {
for (int i = 2; i <= N; ++ i) {
if (!Vis[i]) Pri.push_back (i);
for (auto k : Pri) {
if (i * k > N) break ;
Vis[i * k] = 1;
if (i % k == 0) break ;
}
}
}
4. 素性测试
试除法
容易发现，我们尝试 \(2 \sim \sqrt N\) 之间有没有整数 整除 \(N\) 即可
费马素性测试
基于 费马小定理，若 \(n \in \mathbb P, a \perp N\)，则有 \(a ^ {n - 1} \equiv 1 \pmod n\)
于是 随机选一些 \(a \in [1, n)\)，但 无需保证 \(a \perp n\)，根据上述定理判断
- 若 \(a ^ {n - 1} \equiv 1 \pmod n\) 不成立，则 \(n\) 一定不是素数
- 若 \(a ^ {n - 1} \equiv 1 \pmod n\) 成立，则 \(n\) 大概率 是素数
显然，尝试的 \(a\) 越多，\(n\) 是素数的概率越高，但其始终 不是完全正确的
为什么不用保证 \(a \perp n\)？显然，当 \(n\) 为 素数 是，始终有 \(a \perp n\)，不影响判断
更令人难受的是，有一些合数 \(d\)，无论取什么 \(a \in [1, d)\) 都不能把它判掉
我们称这种数是 \(Carmicgeal\) （卡迈克尔）数，前面几个是 \(561, 1105, 1729, ...\)
它在 \(10 ^ 7\) 的范围内就有 \(255\) 个，根本 特判不掉，这就成了这个素性测试的 致命问题
[4] \(Wiki - pedia\) 上指明，当 \(n\) 充分大时，\([1, n]\) 间有 至少 \(n ^ {\frac 2 7}\) 个 \(Carmicgeal\) 数
\(Miller - Rabin\)
这才是 重点，其在 保证时间复杂度 的情况下 部分解决了 费马素性测试的不完美
此算法（及其优化变种）是 已知最快的（非确定性）素性测试算法
仍然基于 费马素性测试，然后多次进行 二次探测定理逆否命题判定，排除 \(Carmicgeal\) 数
二次探测定理 可以在 同余关系 一部分中找到
具体而言，对于数 \(n\)，我们将 \(n - 1\) 表示成 \(2 ^ t u\) 的形式，其中 \(2 \nmid u\)
我们任取一基数 \(a\)，计算 \(a ^ u \pmod n\) 的值，考虑每次将其 乘方，即得到 \(a ^ {2u} \pmod n\)
根据 二次探测定理，若 \(n \in \mathbb P\)，则 \(x ^ 2 \equiv 1 \pmod n\) 有且仅有 \(x \equiv 1\) 与 \(x \equiv n - 1\) 两解
当每次 \(a ^ {2u} \equiv 1 \pmod n\) 时，去判断是否 \(a ^ u \equiv 1 \pmod n\) 或 \(a ^ u \equiv n - 1\pmod n\) 即可
容易发现，若 \(n \in \mathbb P\)，或 \(n\) 为 \(Carmicgeal\) 数，则至少 \(a ^ {n - 1} \equiv 1 \pmod n\)，且 \(2 \mid n - 1\)
也就是说 对于一个基数 \(a\)，我们 至少可以做一次检测
据理论证明，其每次测试 错误概率 \(< \dfrac 1 4\)，于是 \(50\) 次测试后错误概率只有 \(4 ^ {-50}\)，十分安全
其优势在于 没有 类似于 \(Carmicgeal\) 数 一样的 无法得出正确解 的数
显然，我们对于 选取的基数 \(a\)，需要做 一次快速幂，同时做 \(t\) 次倍增，时间复杂度 \(O (\log N)\)
故总时间复杂度 \(O (k \log N)\)，其中 \(k\) 是 选择的基数个数，这东西 快的很啊
但是吧，非确定性 的算法 总是不那么让人放心，即使其错误率极低
而事实上，该算法初始实质上是一个 基于广义黎曼猜想 的 确定性算法
虽然 猜想没有被证明，但是我们仍能在 一定范围内验证它
从而得到一个 在一定范围内的确定性算法
具体而言，我们给出以下范围选择
在 \(n \le 2 ^ {32}\) 的范围内，可以使用 \(2, 7, 61\) 三个底数，或 前四个质数 作为底数
在 \(n \le 2 ^ {61}\) 的范围内，可以使用 前九个质数 作为底数
在 \(n \le 2 ^ {64}\) 的范围内，可以使用 \(2, 325, 9375, 28189, 450775, 9780504, 1795265022\) 七个底数
在 \(n \le 2 ^ {78}\) 的范围内，可以使用 前十二个质数 作为底数
在 \(n < 2 ^ {81}\) 的范围内，可以使用 前十三个质数 作为底数
于是这下就可以 放心使用了
namespace Miller_Rabin {
int P[15] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
int MAX_32 = 4, MAX_61 = 9, MAX_64 = 12;
inline long long Qpow (__int128 a, long long b, const long long p) {
__int128 Ret = 1;
while (b) {
if (b & 1) Ret = Ret * a % p;
a = a * a % p, b >>= 1;
}
return Ret;
}
inline bool Check (const long long x) {
if (x < 3 || !(x & 1)) return x == 2;
long long d = x - 1, t = 0;
__int128 v;
while (!(d & 1)) d >>= 1, ++ t;
for (int i = 0; i < MAX_64 && P[i] < x; ++ i) {
v = Qpow (P[i], d, x);
if (v == 1) continue ;
for (int k = 0; k < t; ++ k) {
if (v == x - 1) break ;
v = v * v % x;
if (k == t - 1) return 0;
}
}
return 1;
}
}
以上代码有一些 小小的优化，会比较快
- 对质数 \(2\) 的特判
- 对于底数 \(a\)，有 \(a ^ u \equiv 1 \pmod n\) 时，显然 \(a ^ {u 2 ^ t} \equiv 1 \pmod n\)，可以提前判定通过
- 同时当出现 \(a ^ {u 2 ^ t} \equiv n - 1 \pmod n\) 时，当且仅当 \(t = 0\) 时我们 不能直接判定通过
显然若 \(u = n - 1\) 时...
\(AKS\)
抽象算法，目前已经证明的 渐近时间复杂度 \(O (\log ^ {10.5} n)\)
而 变体 可以优化到 \(O (\log ^ 6 n) \sim O (\log ^ {7.5} n)\)（仍然毫无帮助
但其具有 独特地位，因为
它是第一个被发表的 一般的、多项式的、确定性的 和 无仰赖的 素数判定算法
—— [5] 百度百科
实现和证明 可能对于 \(OI\) 并无帮助，过于困难？了，这里只 简单介绍 基本定理
若 \(a \perp p, a \in \Z\)，则 \(p \in \mathbb P\) 的 充分必要 条件 是
\[\begin {aligned} (x - a) ^ p & \equiv (x ^ p - a) \pmod p \end {aligned} \]但是由于我们要验证 所有与 \(p\) 互质的 \(a\)，上式不可做，于是有以下的 同余多项式
\[\begin {aligned} (x + a) ^ n & \equiv (x ^ n + 1) \pmod {x ^ r - 1, n} \end {aligned} \]这等价于存在 多项式 \(f, g\)，使得
\[\begin {aligned} (x + a) ^ n - (x ^ n + a) & = (x ^ r - 1) g + nf \end {aligned} \]而有关 \(AKS\) 算法的正确性证明中，推导出了一个 足够小的 \(r\)，以及 整数集 \(a \in A\)
只要验证当 \(r\) 为该值时，上式对所有的 \(a \in A\) 均成立，则 \(n\) 为 素数
其中 \(r\) 应当是 最小的 使得 \(n \bmod r\) 的 阶* \(> \log N\) 的 正整数
* : 若 \(n \perp r\)，则称使得 \(n ^ k \equiv 1 \pmod r\) 的 最小正整数 \(k\) 为 \(n \bmod r\) 的 阶
而 \(a \in [1, \left \lfloor \sqrt {\varphi (r)} \log n \right \rfloor]\) 即可，可以证明， \(r\) 与 \(|A|\) 均是 \(poly ~ \log n\) 级别的
5. 分解质因数
试除法
显然的 \(O (\sqrt n)\) 算法，我们在 判定 \(n\) 的素性 的 过程中
每遇到一个 \(a\)，满足 \(a \mid n\)，就 将 \(n\) 除尽 \(a\)，并记录
容易证明 每次找到的 \(a\) 一定是 当前 \(n\) 的 最小质因数
\(Pollard - rho\)
核心思想仍然是 随机
我们知道 生日悖论
即在 \([1, n]\) 的范围内 均匀随机 取数，那么某两次的 \(a, b\) 满足 \(|a - b| \le O (\sqrt n)\) 的概率 大于 \(50 \%\)
也就是 \(O (\sqrt n)\) 次随机内取到至少一组相同的数 的概率，经典的 \(Hash\) 碰撞估计
由此我们知道 组合随机采样 比 单次随机采样 满足条件的概率 高得多
这曾被用来快速有效的 估计鱼塘里鱼的条数（[6] 标记重捕法）
同样，这个原理也可以用来 分解质因数
\(Pollard\) 老师提出了这么一个 随机数生成方法
\[f (x) = (x ^ 2 + c) \bmod c \]其中 \(c\) 是 随机常数，我们 随一个 \(x_1\)，然后我们 不停迭代，即 \(x_2 = f (x_1), x_3 = f (x_2)\)
然后每次取 \(k_i = f (x_i) - f (x_{i - 1})\)，求 \(\gcd (k_i, n)\)
若 该值 不为 \(1\)，我们就找到了一个 \(n\) 的 非平凡因子，返回
于是 递归求解 即可（即将 \(n\) 分为 该非平凡因子 和 \(n\) 除以该因子 两部分）
最终通过 \(Miller-Rabin\) 判定 当前因子 是否为素因子，若是则 终止并返回（无法再分解）
下面尝试证明其 时间复杂度，即 得到每个非平凡因子 的时间为 \(O (n ^ {\frac 1 4} \log n)\)
考虑这个 \(f (x)\) 的 性质， \(\forall x \equiv y \pmod p, f (x) \equiv f (y) \pmod p\)，其中 \(p \mid n\)
证明显然，\(x ^ 2 \equiv y ^ 2 \pmod p\) 即可
这等价于告诉我们，在 按照上述方法生成的序列 \(x_1, x_2, x_3, ..., x_n\) 中
若 \(x_i = x_j\) 则 \(x_{i + 1} = x_{j + 1}, ..., x_{i + k} = x_{j + k}\)
容易发现 当找到这样的一组 \(i, j\) 后，继续迭代的 \(f (x)\) 将 成环
显然 \(x_{i + (j - i)} = x_j = x_{j + (j - i)} = ... = x_{i + k (j - i)}\)，环长为 \(j - i\)
根据 生日悖论，若 \(x_k\) 在 \([0, n)\) 中 均匀随机，则只需要 \(O (\sqrt n)\) 次就可以找到相等的 \(x_i, x_j\)
故也就是说（在符合预期的情况下）\(i, j < O (\sqrt n)\)，则 序列中不同的数为 \(O(\sqrt n)\) 个种
而我们设 \(n\) 的 最小非平凡因子 为 \(m\)，显然有 \(m \le \sqrt n\)
若设 \(y_i = x_i \bmod m\)，同上述，容易说明 \(y_i\) 中有 \(O (\sqrt m) \le O (n ^ {\frac 1 4})\) 种 不同的数
注意到 \(f (x)\) 的性质不只对 \(n\) 成立，而是对任意 \(p \mid n\) 均成立
即 我们取相邻的两项 \(x_i, x_{i + 1}\) 做差，其结果为 \(m\) 整倍数 的概率约 \(\dfrac 1 {n ^ {\frac 1 4}}\)
换言之，按照算法流程，我们期望 \(O (n ^ {\frac 1 4})\) 次就可以找到一个 \(n\) 的 非平凡因子
其中一次 \(\gcd\) 需要 \(O (\log N)\) 的时间
根据上述证明，我们知道了 \(f (x)\) 会在 迭代一定次数（\(O (\sqrt n\)）后 产生环
于是 生成的整个序列 在坐标系上就形如 \(\rho\)，这也是 算法后半部分名字 的由来
显然此时的 \(k_i\) 产生重复，再去判断是 无意义的，故我们要 找到环后退出
若此时 仍未找到非平凡因子，则 换一个 \(c\) 与 \(x_1\) 重新开始即可
这里我们可以使用 \(Floyd\) 的方法，我们维护 两个 \(f (x)\) 生成的值 \(x, x'\)
初始时 \(x = x' = x_1\)，之后 \(x\) 每迭代一次，\(x'\) 就 迭代两次
即当 \(x = x_p\) 时 \(x' = x_{2p}\)，显然，当（初始之后）首次 \(x = x'\) 时，我们就找到了环
容易得出，此时有 $2p - p = $ 环长，也就是上述证明过程中的 \(j - i\)，期望 \(O (\sqrt n)\)
于是我们 迭代 \(O (\sqrt n)\) 次即可 正常退出，不会有大量浪费
可能会发现，当退出时 \(x\) 不一定遍历完整个环
但这是 无关紧要的，因为本身这只是一个 随机序列，遍历完整并不会有什么额外的优势
但一般而言 不要尝试跳过前半部分的遍历（即不要什么初始时 \(x = x_{10}\) 这样的事情）
因为显然，到环之前 出现有效答案的效率很高，不应浪费大量时间在判环上
于是我们可以写出 如下代码
namespace Pollard_Rho {
namespace Miller_Rabin {}
mt19937 r (time (0));
inline long long Nxt (const long long x, const long long c, const long long n) {
return ((__int128) x * x + c) % n;
}
inline long long solve (const long long n) {
long long c = r () % n, x = r () % n, t = Nxt (x, c, n), pre = 0, tmp = 0, gcd_ = 0, cnt = 0, sqt = sqrt (n);
while (x != t) {
tmp = abs (x - pre);
if ((gcd_ = __gcd (tmp, n)) > 1) return gcd_;
pre = x, x = Nxt (x, c, n), t = Nxt (Nxt (t, c, n), c, n);
if (++ cnt > sqt) break ;
}
return n;
}
inline void Solve (const long long n, vector <long long> &a) {
if (Miller_Rabin::Check (n)) return a.push_back (n), void ();
long long d = n;
while (d == n) d = solve (n);
Solve (d, a), Solve (n / d, a);
}
}
其中
Miller_Rabin部分就是上面部分，这里就省略了
但这东西 不太行，甚至过不去板子题 Luogu P4718 【模板】Pollard-Rho
考虑优化，看上去这个 \(\gcd\) 就 很令人火大，多了个 \(\log\)，和 \(\color {black} \textsf {H} \color {red} \textsf {anghang}\) 一样
显然，我们发现 \(\gcd\) 这个东西存在 这么一个性质
若 \(\gcd (a, n) > 1\)，则一定有 \(\gcd (ab, n) > 1\)
于是在计算的时候，我们可以把 每次的 \(pre\) 相乘并累积成一块
一定次数后再将这个 一块的 \(pre'\) 与 \(n\) 求一次 \(\gcd\)，若 \(> 1\)，直接返回即可
若存在一次使得 \(pre' = 0\)，为防止 丢失答案，我们 退出循环并暴力判断
玄学部分
- 最佳 ”块长“ 约为 \(128\)（\(OI\) 范围）
- 仍然注意前面 有效信息 密集，故实际块长按 \(1, 2, 4, ..., 128 , 128, 128 ...\) 这样取
- 小质因数（\(2\)）之类的 较难被判断出来，性价比很低，可以直接特判
- \(c\) 一般在 \([3, n)\) 范围内 比较容易出答案
于是我们就得到了 最终代码，也就是 能通过板子题 的版本
#include <bits/stdc++.h>
const int MAXN = 1000005;
using namespace std;
namespace Pollard_Rho {
namespace Miller_Rabin {
int P[15] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
int MAX_64 = 12;
inline long long Qpow (__int128 a, long long b, const long long p) {
__int128 Ret = 1;
while (b) {
if (b & 1) Ret = Ret * a % p;
a = a * a % p, b >>= 1;
}
return Ret;
}
inline bool Check (const long long x) {
if (x < 3 || !(x & 1)) return x == 2;
long long d = x - 1, t = 0;
__int128 v;
while (!(d & 1)) d >>= 1, ++ t;
for (int i = 0; i < MAX_64 && P[i] < x; ++ i) {
v = Qpow (P[i], d, x);
if (v == 1) continue ;
for (int k = 0; k < t; ++ k) {
if (v == x - 1) break ;
v = v * v % x;
if (k == t - 1) return 0;
}
}
return 1;
}
}
mt19937 r (time (0));
inline long long Mul (const long long x, const long long y, const long long n) {
long long d = x * (long double) y / n, ret = x * y - d * n;
if (ret < 0) return ret + n;
if (ret > n) return ret - n;
return ret;
}
inline long long add (const long long a, const long long b, const long long n) {
long long d = (a + (long double) b) / n + 0.5, ret = a + b - d * n;
if (ret < 0) return ret + n;
return ret;
}
inline long long Add (const __int128 a, const __int128 b, const long long n) {
__int128 c = a + b;
if (c >= n) return c - n;
else return c;
}
inline long long Nxt (const long long x, const long long c, const long long n) {
return Add (Mul (x, x, n), c, n);
// return ((__int128) x * x + c) % n;
}
inline long long solve (const long long n) {
if ((n & 1) == 0) return 2;
long long c = r () % (n - 3) + 3, x = r () % n, t = Nxt (x, c, n);
long long pre = 0, tmp = 0, gcd_ = 0;
for (int lim = 1; x != t; lim = min (1 << 7, lim << 1)) {
pre = 1, tmp = 0;
for (int i = 1; i < lim; ++ i, pre = tmp) {
tmp = Mul (pre, abs (Add (x, - t, n)), n);
// tmp = (__int128) pre * abs (t - x) % n;
if (!tmp) break ;
x = Nxt (x, c, n), t = Nxt (Nxt (t, c, n), c, n);
}
gcd_ = __gcd (pre, n);
if (gcd_ > 1) return gcd_;
}
return n;
}
long long Ans = 0;
inline void Solve (const long long n) {
if (Miller_Rabin::Check (n)) return Ans = max (Ans, n), void ();
long long d = n;
while (d == n) d = solve (n);
Solve (d), Solve (n / d);
}
}
int T;
long long A;
int main () {
cin >> T;
for (int i = 1; i <= T; ++ i) {
cin >> A, Pollard_Rho::Ans = 0;
Pollard_Rho::Solve (A);
if (Pollard_Rho::Ans == A) cout << "Prime" << '\n';
else cout << Pollard_Rho::Ans << '\n';
}
return 0;
}
注意到我们每次乘 \(t - x\) 与 乘上 相邻两项 在 随机序列 的情况下
显然等价（出现因数的概率相等）
只是这样 好写，并没有什么 特殊性质
使用 快速乘 和 快速加 替换掉
__int128后，效率 极大的提升，\(2.27 s \rightarrow 722 ms\)
Luogu P4358 [CQOI2016] 密钥破解
就是板子
给出 \(N = p \times q (p, q \in \mathbb P), e, c\)
而有 \(r = (p - 1) (q - 1), de \equiv 1 \pmod r, c ^ d \equiv n \pmod N\)
试求 \(d, n\)
这实质上就是 \(RSA\) 加密算法，这也是其需要 大质数 来 防止破译的原因（有 多项式算法）
容易发现 \(d\) 即是 \(e \bmod r\) 的 逆元，\(\operatorname {exgcd}\) 即可（\(r \notin \mathbb P\)，不可以用 费马小定理）
故重点在于求 \(r\)，显然知道 \(p, q\) 即可，于是 \(Pollard - rho\) 将 \(N\) 分解就行了
时间复杂度 \(O (N ^ {\frac 1 4})\)，只有一组数据，很轻松，甚至优化的 \(O (\sqrt N)\) 都能过
#include <bits/stdc++.h>
const int MAXN = 100005;
using namespace std;
namespace Pollard_Rho {
namespace Miller_Rabin {
int P[15] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
int MAX_64 = 12;
inline long long Qpow (__int128 a, long long b, const long long p) {
__int128 Ret = 1;
while (b) {
if (b & 1) Ret = Ret * a % p;
a = a * a % p, b >>= 1;
}
return Ret;
}
inline bool Check (const long long x) {
if (x < 3 || !(x & 1)) return x == 2;
long long d = x - 1, t = 0;
__int128 v;
while (!(d & 1)) d >>= 1, ++ t;
for (int i = 0; i < MAX_64 && P[i] < x; ++ i) {
v = Qpow (P[i], d, x);
if (v == 1) continue ;
for (int k = 0; k < t; ++ k) {
if (v == x - 1) break ;
v = v * v % x;
if (k == t - 1) return 0;
}
}
return 1;
}
}
mt19937 r (time (0));
inline long long Mul (const long long x, const long long y, const long long n) {
long long d = x * (long double) y / n, ret = x * y - d * n;
if (ret < 0) return ret + n;
if (ret > n) return ret - n;
return ret;
}
inline long long add (const long long a, const long long b, const long long n) {
long long d = (a + (long double) b) / n + 0.5, ret = a + b - d * n;
if (ret < 0) return ret + n;
return ret;
}
inline long long Add (const __int128 a, const __int128 b, const long long n) {
__int128 c = a + b;
if (c >= n) return c - n;
else return c;
}
inline long long Nxt (const long long x, const long long c, const long long n) {
return Add (Mul (x, x, n), c, n);
// return ((__int128) x * x + c) % n;
}
inline long long solve (const long long n) {
if ((n & 1) == 0) return 2;
long long c = r () % (n - 3) + 3, x = r () % n, t = Nxt (x, c, n);
long long pre = 0, tmp = 0, gcd_ = 0;
for (int lim = 1; x != t; lim = min (1 << 7, lim << 1)) {
pre = 1, tmp = 0;
for (int i = 1; i < lim; ++ i, pre = tmp) {
tmp = Mul (pre, abs (Add (x, - t, n)), n);
// tmp = (__int128) pre * abs (t - x) % n;
if (!tmp) break ;
x = Nxt (x, c, n), t = Nxt (Nxt (t, c, n), c, n);
}
gcd_ = __gcd (pre, n);
if (gcd_ > 1) return gcd_;
}
return n;
}
inline void Solve (const long long n, vector <long long> &a) {
if (Miller_Rabin::Check (n)) return a.push_back (n), void ();
long long d = n;
while (d == n) d = solve (n);
Solve (d, a), Solve (n / d, a);
}
}
using namespace Pollard_Rho::Miller_Rabin;
inline long long EXGCD (const long long a, const long long b, long long &x, long long &y) {
if (b == 0) return x = 1, y = 0, a;
long long d = EXGCD (b, a % b, y, x);
y -= (a / b) * x; return d;
}
vector <long long> p;
long long N, e, c;
long long k, r, d, n;
int main () {
cin >> e >> N >> c;
Pollard_Rho::Solve (N, p);
r = N - p[0] - p[1] + 1;
EXGCD (e, r, d, k);
d = (d % r + r) % r;
n = Qpow (c, d, N);
cout << d << ' ' << n << '\n';
return 0;
}
\(Pollard-rho\) 的板子 过长了，后续考虑省略，则以此版本为准
Luogu P4714 「数学」约数个数和
需要用到一些 数论函数基础 里的内容
给定整数 \(N, K\)，试求 \(N\) 的（（所有约数的）\(\times K\)）的约数和
输出对 \(998244353\) 取模，\(N, K < 10 ^ {18}\)
考虑 \(K = 0\)，有 \(f_0 (x) = \sigma (x)\) 为 积性函数*
这里我们利用 积性函数 存在 \(f (x) f(y) = f (xy), x, y \in N ^ *\) 这个性质
此处 暂不给出 \(\sigma (x)\) （约数和）为积性函数的证明
然后有 \(K > 0\) 时，\(f_i (x) = \sum \limits _ {d \mid x} f_{i - 1} (d)\) ，这是 狄利克雷卷积 的形式
根据 狄利克雷卷积 的性质，可以知道 \(f_i (x)\) 仍为 积性函数
则我们可以将 \(n\) 质因数分解，将最后 结果的 \(f\) 相乘即可，形式化表达则是
\[f_K (n) = \prod (f_K (p_i ^ {q_i})) \]其中 \(n = \prod p_i ^ {q_i}\) 是 \(n\) 的 唯一分解式
于是我们只需要考虑如何求 \(f_K (p ^ q)\) 即可，带回 \(f_i (x) = \sum \limits _ {d \mid x} f_{i - 1} (d)\)，可将其变形成
\[f_i (x) = \sum \limits _ {i = 1} ^ q f_{i - 1} (p ^ i) \]可以发现，这等价于将 \(q\) 个位置分成 \(K + 1\) 份，每份可以为 \(0\)，插板法很好理解
故我们就是求 \({q + K + 1} \choose {K + 1}\)，容易变形成 \({q + K + 1} \choose {q}\)，显然 \(q \le 62\)，暴力做就行
然后其中 分解质因数 的部分需求 \(Pollard-rho\)，此处不放板子，参照上一道题
注意 \(N = 1\)！
#include <bits/stdc++.h>
using namespace std;
const int MOD = 998244353;
namespace Pollard_Rho {}
using namespace Pollard_Rho::Miller_Rabin;
struct Prime {
long long v;
int cnt;
} Pri[36];
int Cnt = 0;
inline void Add (const int p) {
for (int i = 1; i <= Cnt; ++ i)
if (Pri[i].v == p) return Pri[i].cnt ++, void ();
Pri[++ Cnt].v = p, Pri[Cnt].cnt = 1;
}
vector <long long> Frac;
long long N, K;
long long Ans = 1;
inline long long F (const int q) {
long long Ret = 1;
for (int i = q; i >= 1; -- i)
Ret = (__int128) Ret * (i + K + 1) % MOD * Qpow (i, MOD - 2, MOD) % MOD;
return Ret;
}
int main () {
cin >> N >> K;
if (N == 1) return cout << 1 << '\n', 0;
Pollard_Rho::Solve (N, Frac);
for (auto p : Frac) Add (p);
for (int i = 1; i <= Cnt; ++ i)
Ans = Ans * F (Pri[i].cnt) % MOD;
cout << Ans << '\n';
return 0;
}
总算结束力！
6. 引用资料
[1] Stern–Brocot 树与 Farey 序列 —— OI - Wiki
[2] 「P8058 [BalkanOI2003] Farey 序列」—— _Fontainebleau_
[3] DIVCNT1 - Counting Divisors —— yhx-12243
[5] AKS素数测试 —— 百度百科
[6] 标记重捕法 —— 百度百科
标签：return,int,dfrac,long,素数,inline,const From： https://www.cnblogs.com/FAKUMARER/p/18316762