Number Theory（2）

5 素数

接下来我们来讨论素数。

5.1 素数的例子

如果一个正整数 $p$ 恰好只有两个因子，即 $1$ 和 $p$ ，那么这个数就称为素数。

而其他那些还有非平凡因子的数都称为合数，每一个大于 $1$ 的整数要么是素数，要么是合数。

任何正整数 $n$ 都可以表示成素数的乘积：

\[n = p1 \cdots p_m=\prod_{k=1}^{m} p_k\quad, p1 \le \dots \le p_m \]

而 算数基本定理 表示：仅有一种方法将 $n$ 按照素数非减的次序写成素数的乘积。（证明略）

关于素数的其他定理：

第 $n$ 个素数 $P_n$ 大约是 $n$ 的自然对数的 $n$ 倍：

\[P_n \sim n \ln n \]
对于不超过 $x$ 的素数个数 $\pi(x)$:

\[\pi(x) \sim \frac{x}{\ln x} \]

证明过于复杂，这里不给出证明，大家可以自行了解。

还有一种与素数相关的数：梅森数。

形如

\[2^p-1 \]

的数，其中 $p$ 总表示素数。

如果 $n$ 是合数，那么 $2^n -1$ 不可能是素数，因为 $2^{km}-1$ 以 $2^m-1$ 作为一个因子：

\[2^{km}-1=(2^m-1)(2^{m(k-1)}+2^{m(k-2)}+\cdots + 1) \]

但是当 $p$ 是素数时，$2^p-1$ 并不总是素数，$2^{11}-1=2047=23 \times 89$ 是最小的这类非素数。

5.2 互素关系

当 $\gcd(n,m)= 1$ 时，整数 $m$ 和 $n$ 没有公共的素因子，我们就称它们是 互素的。

我们写成

\[m \bot n \Longleftrightarrow m,n 是整数，且 \gcd(i,j)=1 \]

而一个分数是最简分数，当且仅当 $m \bot n$，通常我们通过消去分子分母的最大公约数来简化成最简分数，于是我们有

\[\frac m {\gcd(n,m)} \bot \frac n {\gcd(n,m)} \]

那我们有没有什么好的方法来构造满足 $m \bot n$ 的全部非负的分数 $\frac m n$ 组成的集合呢？

答案是肯定的，他被成为 Stern-Brocot 树，其思想是从两个分数 $\left( \frac 0 1 ,\frac 1 0 \right)$ 出发，然后依照你希望的次数重复下面的操作

在两个相邻的分数 $\frac m n$ 和 $\frac {m'} {n'}$ 之间插入 $\frac {m + m'} {n + n'}$。

新的分数就是 $\frac {m+m'}{n+n'}$ 称为 $\frac m n$ 和 $\frac {m'} {n'}$ 的 中位分数，例如

\[\frac 0 1 ,\frac 1 1,\frac 10 \]

接下来就是

\[\frac 01,\frac 1 2,\frac 11,\frac 21,\frac 10 \]

于是我们把这个分数阵列可以看成是一棵无线的二叉树构造，就是这样

每一个分数都是 $\frac {m+m'} {n+n'}$，其中 $\frac m n$ 是位于左上方且离它最近的祖先，而 $\frac{m'} {n'}$ 则是右上方离它最近的祖先。

为什么这种构造是对的呢？为什么每次的 $\frac {m+m'} {n+n'}$ 是最简分数？而又只会出现一次？

因为我们发现，在这个构造中任意一个阶段的相邻分数都满足

\[m'n-mn'=1 \]

首先这个关系在 $1 \times 1 - 0 \times 0=1$ 是成立的，而每一次插入一个新的中位分数 $\frac {m+m'} {n+n'}$ 时，我们需要检验的就是

\[\begin{aligned} (m+m')n - m(n+n') & =1 \\ m'(n+n')-(m+m')n' & =1 \end{aligned} \]

括号展开就一样了，而这也让我们发现

\[\frac m n \lt \frac {m+m'} {n+n'} \lt \frac {m'} {n'} \]

一个中位分数不一定在原先两个的中间，但它的确位于它们两个中间的某个地方，于是我们就不可能得到相同的分数。

但是我们现在仍然存在一个问题——会遗漏么？

这个问题时简单的，由于它是无限的，所以我们会不断用上面的单调性去逼近它，这样就不会漏。

这里我们再引入一个概念，阶为 $N$ 的法里级数 $\mathscr F_n$，它是介于 $0$ 和 $1$ 之间的分母不超过 $N$ 的所有最简分数组成的集合，且按照递增的次序排列，那么当 $N=6$ 时就是

\[\mathscr F_n = \frac 0 1,\frac 1 6,\frac 1 5,\frac 1 4 ,\frac 1 3,\frac 2 5,\frac 1 2,\frac 3 5,\frac 2 3,\frac 3 4,\frac 4 5,\frac 5 6,\frac 1 1 \]

容易发现它是 S-B 树 的子树，好像也没什么好讲的。

事实上，我们可以把 S-B 树 看成一个表示有理数的数系，因为每一个正的最简分数都恰好出现一次。

我们用字母 $L,R$ 表示它在这一层往哪边走，就可以唯一确定一个数。

我们记当前的字符串为 $S$，那么就有

\[f(S) = 与 S 对应的分数 \]

例如 $f(LRRL)= \frac 5 7$。

这个东西和矩阵乘法的构造是比较像的，也就是我们尝试去构造矩阵

\[M(s) = \left( \begin{matrix} n & n' \\ m & m' \end{matrix} \right) \]

那么我们知道

\[M(SL)= \left(\begin{matrix} n & n+n' \\ m & m+m' \end{matrix}\right) = \left(\begin{matrix} n & n' \\ m & m' \end{matrix}\right) \left(\begin{matrix} 1 & 1 \\ 0 & 1 \end{matrix}\right) \]

然后我们就发现其实

\[L= \left(\begin{matrix} 1 & 1 \\ 0 & 1 \end{matrix}\right) ,R= \left(\begin{matrix} 1 & 0 \\ 1 & 1 \end{matrix}\right) \]

于是我们就可以通过一个 $LR$ 字符串定位出它所定义的值了。

还有一种方法就是改变 $m,n$，而不是保持状态 $S$，我们发现

\[f(RS) = f(S) +1 \]

因为 $RS$ 像 $S$，只是将上面一行加到了下面一行上，也就是

\[S= \left(\begin{matrix} n & n' \\ m & m' \end{matrix}\right) ,RS= \left(\begin{matrix} n & n' \\ m+n & m'+n' \end{matrix}\right) \]

类似的性质对 $L$ 也有，于是我们就有

\[\begin{cases} \frac m n = f(RS) \Leftrightarrow & \frac {m-n} n = f(S) \\ \frac m n = f(LS) \Leftrightarrow & \frac {m} {n-m} = f(S) \end{cases} \]

这样的定位就可以不用矩阵乘法。

关于互素的讨论到此结束，更多可以参见《具体数学》！

P8058 [BalkanOI2003] Farey 序列

把所有分子和分母都 $\le n$ 的 最简真分数 从小到大排成一行，形成的序列称为 Farey 序列。

求出 $n$ 所对应的 Farey 序列中第 $k$ 小的数。

$2 \le n \le 4 \times 10^4,1 \le k \le$ 符合条件分数的个数。

sto smb orz

给出一个 smb 的亚线性做法。

有了上面的一些简单铺垫，我们不难想到在 S-B 树 上面进行二分。

那么对于一个最简分数 $\frac x y$，有多少个数比它小呢？

容易列出式子

\[\begin{aligned} & \sum_{i=1}^n \sum_{j=1}^i [\gcd(i,j)=1] \left[ \frac j i \le \frac x y \right] \end{aligned} \]

对这个式子进行莫反，可以得到

\[\newcommand {\lf} {\left \lfloor} \newcommand {\rf} {\right \rfloor} \begin{aligned} & \sum_{i=1}^n \sum_{j=1}^i [\gcd(i,j)=1] \left[ \frac j i \le \frac x y \right] \\ = & \sum_{i=1}^n \sum_{j=1}^i \sum_{d \mid i,d \mid j} \mu(d) \left[ \frac j i \le \frac x y \right] \\ = & \sum_{d=1}^n \mu(d) \sum_{i=1}^{\lf \frac n d \rf} \sum_{j=1}^i \left[ \frac j i \le \frac x y \right] \\ = & \sum_{d=1}^n \mu(d) \sum_{i=1}^{\lf \frac n d \rf} \sum_{j=1}^i \left[ j \le \frac {ix} y \right] \\ = & \sum_{d=1}^n \mu(d) \sum_{i=1}^{\lf \frac n d \rf} \lf \frac {ix} y \rf \end{aligned} \]

我们发现前面是可以用杜教筛处理的（其实这道题中不用），而后面这一块是类欧几里得的板子，于是这个值我们是可以在 $\mathcal O(n^{\frac 2 3} + \sqrt n \log n)$ 的时间内完成的。

而我们如何去二分到这个值呢？

不难想到直接在 S-B 树 上面走，每次枚举走到左子树还是走到右子树，但是在 $k=1$ 的情况我们会走 $n$ 步，很明显是不满足条件的。

我们考虑把往一个方向走的连续段缩起来，也就是只枚举在哪些地方拐（eg：原先一直往左子树走，现在变成往右子树走），设 $F(n)$ 表示拐 $n$ 次能到达的分母最小值，那么容易得到 $F(n) = F(n-1)+F(n-2)$，也就是一直拐 $n$ 次，那么左边就是 $F(n-1)$，右边是 $F(n-2)$。

那么这样我们只会拐 $\log n$ 次，每次走多长可以用倍增实现，于是我们就做完了。

总时间复杂度 $\mathcal O(n^{\frac 2 3} + \sqrt n \log^3 n)$。代码。

ABC333G Nearest Fraction

给定一个小于 $1$ 的正实数 $r$ 和一个正整数 $n$。要求在满足 $0 \le p \le q \le n$ 和 $\gcd(p,q)=1$ 的前提下，找到使 $|r-\frac{q}{p}|$ 最小的二元组 $(p,q)$ 。如果存在多个这样的二元组 $(p,q)$，输出 $\frac{q}{p}$ 值最小的那个。

数据范围：$1\le n \le 10^{10}$，$r\in (0,1)$ 且最多包含 18 位有效数字。

首先 $r$ 容易被转化成分数形式，把它设为 $m$，我们考虑用 S-B 树 来解决这个问题。

对于当前的区间 $[l,r]$，我们考虑在树上二分来找到最近的位置，那么每一次一定是形如先从 $r$ 一直往左跳到第一个 $\lt m$ 的地方，再把这个点设成左端点，又不断往右跳。

这样我们就可以不断减小这个区间，从而不断逼近 $m$。

每一次搜索都用两个答案 $ans1,ans2$ 分别记录在 $m$ 两侧得到的最近点，这样就做完了。

用上一道题类似的倍增优化可以做到 $\mathcal O(\log^2 n)$。代码。

双倍经验：P1298 最接近的分数

SP26073 DIVCNT1 - Counting Divisors

$\sigma_0(i)$ 表示 $i$ 的约数个数。

求

\[S(n) = \sum_{i=1}^n \sigma_0(i) \]
多测，$T \le 10^5,n \le 2^{63}$。

比较难的应用。

首先我们可以对式子进行转化

\[\newcommand {\lf} {\left \lfloor} \newcommand {\rf} {\right \rfloor} \begin{aligned} S(n) = & \sum_{i=1}^n \sigma_0(i) \\ = & \sum_{i=1}^n \lf \frac n i \rf \\ = & \sum_{i=1}^n \sum_{j=1}^n [i \times j \le n] \\ = & \sum_{i=1}^{\lf \sqrt n \rf} \sum_{j=1}^n [i \times j \le n] + \sum_{i=1}^n \sum_{j=1}^{\lf \sqrt n \rf} [i \times j \le n] - \sum_{i=1}^{\lf \sqrt n \rf} \sum_{j=1}^{\lf \sqrt n \rf} [i \times j \le n] \\ = & 2 \sum_{i=1}^{\lf \sqrt n \rf} \lf \frac n i \rf - {\lf \sqrt n \rf}^2 \end{aligned} \]

那么问题就转化成求

\[\newcommand {\lf} {\left \lfloor} \newcommand {\rf} {\right \rfloor} \sum_{i=1}^{\lf \sqrt n \rf} \lf \frac n i \rf \]

我们尝试把这个反比例函数画出图来，那么我们要求的其实就是下图中红色部分中的点数（from x义x）：

不妨将其称作 R 区域。

我们希望用一些向量去不断拟合这个反比例函数，于是就想到了这个有理数系：S-B 树。

于是我们需要一个单调栈，里面维护的向量斜率单调递减，进行这样的操作：

每次从栈顶取出一个元素，不断加到当前点 $(x,y)$ 上面，直到它差一步就走进 R 区域。

由于我们保证了斜率递减，且向量的斜率表示出来已经是最简分数，所以每一步都是对答案有贡献的，也就是图中棕色区域中的点（from x义x）：

容易发现这样的点数一定不会多算，可以做到不重不漏。（在斜线下的点一定在反比例函数中，并且不会多）

如何计算点数？（我们假设每次上都要计算最下面的一行，从而不能计算上面的一行）

画出图来就是这个样子：

把它分成蓝色和橙色的两个部分，而橙色部分一定是上下对称的，所以假设当前向量为 $L$，当前起点在 $(x,y)$，终点在 $(x+L.x,y-L.y)$，则贡献是

\[ x \cdot L.y + \frac {(L.x-1)(L.y+1)} 2 \]

不断弹出栈顶，使得攻栈顶刚好走不进 R 区域，称它为 $r$。设最近一次弹的是 $l$（有可能是上一步的栈顶向量）。

现在就是，加上 $l$ 刚好走进 R，加上 $r$ 刚好走不进 R，就是这个样子（from x义x）：

开始二分，用 S-B 树 得到中间向量 $(l_x+r_x,l_y+r_y)$，我们称之为 $mid$。

如果 $mid$ 走一步不会进 R，那么 $mid$ 进栈并把 $r = mid$ 继续二分。
否则
- 如果 $r \le f'(x+x_{mid})$，那么继续二分肯定都走不出 R 了，直接停止二分，就是这个样子（from x义x）：

可以发现如果 $r \le f'(x+x_mid)$，$mid,r$ 的中间向量 $midmid$ 显然也一定会进 R，而且因为 $f'$ 单调增，$r \le f'(x+x_{midmid})$ 也一定满足。

- 否则 $l=mid$ 继续二分。

当 $y \le \sqrt[3] n$ 时可以直接暴力。

这样的时间复杂度是 $\mathcal O(\sqrt[3] n \log n)$。（笔者不会证）

注意到其他的一些东西：

最开始的点我们选择 $\left(\left \lfloor \frac n B \right \rfloor,B+1\right)$，因为这样才能做到不重不漏（细节会少很多）。
初始栈中的向量是 $(1,0)$ 和 $(1,1)$。

这样我们就做完了。代码。

5.3 素数筛法

最重要的知识点，几乎所有题目都需要筛出 $1 \sim n$ 的所有素数。

首先我们存在一个非常显然的做法——埃氏筛法。

埃氏筛用所有已经筛出的素数倍数标记合数：从 $2$ 到 $n$ 枚举所有 $i$，若 $i$ 未被标记，则 $i$ 是素数，将 $i$ 除 $i$ 以外的倍数标记为合数。

埃氏筛的精髓在于复杂度的证明：不超过 $n$ 的指数的倒数之和为 $\mathcal O(\ln \ln n)$ 级别。

也就是

\[\sum_{p \in \mathbb P ,p \le n} \frac 1 p = \mathcal O(n \ln\ln n) \]

这说明埃氏筛的时间复杂度是 $\mathcal O(n \ln \ln n)$。

这里是来自 djq 的证明：

因为每个数只会被其素因子筛到，所以

\[\sum_{p \in \mathbb P,p \le n} \frac 1 p = \sum_{i=1}^n \omega (i) \]

根据 $d(i)$ 的计算式

\[\sum_{i=1}^n 2^{\omega (i)} \le \sum_{i=1}^n d(i) = \mathcal O(n \ln n) \]

根据 $2^x$ 的凸性和琴生不等式得

\[\begin{aligned} \sum_{i=1}^n 2^{\omega (i)} & \ge n 2^{\frac {\sum_{i=1}^n \omega(i)} n} \\ \therefore 2^{\frac {\sum_{i=1}^n \omega(i)} n} & \le \mathcal O(\ln n) \end{aligned} \]

两边同时取对数

\[\frac {\sum_{i=1}^n \omega(i)} n \le \mathcal O(\ln \ln n) \]

因此

\[\sum_{i=1}^n \omega(i) \le \mathcal O(n \ln \ln n) \]

埃氏筛已经非常接近与线性了，但是它也还不是最优的，我们存在一个 线性筛素数 的做法。

它的思想和埃氏筛类似，埃氏筛的优秀之处在于仅用质数的倍数筛去合数，但合数会被多个质因子筛到。让每个合数仅被晒一次就能做到线性。

考虑用每一个合数 最小质因数 筛它：从 $2$ 到 $n$ 枚举所有数 $i$。对于每个 $i$，令其最小的质因子为 $p$，则对于不大于 $p$ 的质数 $q$，$iq$ 的最小质因子为 $q$。将所有 $iq$ 标记为合数，则每个合数 $c$ 仅在 $i = \frac c m$ 时以 $im$ 的形式晒去，其中 $m$ 是 $c$ 的最小质因子。

综上，有如下步骤：

从小到大遍历当前筛出的所有素数 $pr_j$，将 $i \times pr_j$ 标记为合数。
若 $pr_j \mid i$，退出循环，因为 $pr_j \mid i \times pr_k(k \gt j)$，所以 $i \times pr_k$ 的最小质因子为 $pr_j$ 不是 $pr_k$，再筛就重复了。

时间复杂度线性，于是就可以通过 P3383 【模板】线性筛素数。

  for(int i=2;i<N;i++){
  	if(!vis[i]) prm[++cnt]=i;
  	for(int j=1;j<=cnt&&prm[j]*i<N;j++){
  	  vis[prm[j]*i]=true;
  	  if(i%prm[j]==0) break;
  	}
  }

5.4 素数判定（Miller Rabin 素性测试）

如何判断一个数是否是素数？

小素数的判定是容易的：直接试除法即可，时间复杂度 $O(\sqrt{n})$。

至于大素数判定……

费马曾提出一种方法：费马素性测试。

它基于费马小定理：设 $n$ 是素数， $a$ 是正整数且与 $n$ 互素，那么有 $a^{n-1} \equiv 1 \pmod n$。

为了测试 $n$ 是否是素数，我们在 $1 \sim n$ 中随机选择一个基数 $a$ ，而 $a$ 不需要与 $n$ 互素。

如果 $a^{n-1} \equiv 1 \pmod n$ 不成立，那么 $n$ 肯定不是素数。
如果 $a^{n-1} \equiv 1 \pmod n$ 成立，那么 $n$ 很大概率是素数，尝试的 $a$ 越多， $n$ 是素数的概率就越大。

但是从第 $2$ 种情况看出来费马素数测试不是完全正确的。

对于某个 $a$ 值，总有一些合数被误判通过了测试；不同的 $a$ 值，被误判的合数不太一样。

特别地，有一些合数，不管选什么 $a$ 值，都能通过测试，这种数叫 Carmicheal 数 $(561,1105,1729,……)$。

不过这种数很少，$1$ 亿以内只有 $255$ 个，而越到大的越稀疏，所以费马素性测试几乎不会出错，编码也很简单。

由于上面 费马素性测试 的不完美，我们就引入了 Miller-Rabin 素性测试。

把费马素性测试改进一下，它是已知最快的随机素数测试算法。

原理也非常简单，就是费马素性测试排除 Carmicheael 数。

因为并且 Carmicheael 数在进行次数较多的二次探测定理逆否命题判定后，都会被判定出来。

如果 $n > 2$，且 $n$ 是奇数，测试它是否为素数。

根据费马测试，如果 $a^{n-1}\equiv 1 \pmod n$ 不成立，一定不是素数。

而由于 $n-1$ 过大，我们考虑一个技巧：

把 $n-1$ 表示成幂的形式，令 $n-1=2^tu$ 其中 $u$ 是奇数，$t$ 是正整数。

于是有：

\[a^{n-1}\equiv (a^u)^{2^t}\pmod n \]

所以计算出 $a^u$ 依次乘二，这样每次平方同时使用二次探测定理进行判定。

容易发现，对于多个随机的基值 $a$，假设随机了 $s$ 个，做了 $s$ 次测试，出错的概率是 $2^{-s}$ 。

当 $s=50$ 时，出错的概率已经可以忽略不计了。

考虑做了 $s$ 次运算，每次做 $t$ 次快速幂，在 $n$ 很大的时候速度极快！

ll mul(ll a,ll b,ll m){ll z=(long double)a/m*b,res=(a*b-z*m)%m;return (res+m)%m;}
ll add(ll a,ll b,ll n){return (a+b>=n)?a+b-n:a+b;} 
ll qpow(ll a,ll b,ll m){
  ll res=1ll;
  while(b){
  	if(b&1) res=mul(res,a,m);
  	a=mul(a,a,m);
  	b>>=1;
  }
  return res;
}

const ll prm[12]={2,3,5,7,11,13,17,19,23,29,31};
bool miller_rabin(ll n){
  if(n<2) return false;
  for(int i=0;i<=10;i++){
  	if(n==prm[i]) return true;
  	if(n%prm[i]==0) return false;
  } 
  ll u=n-1,t=0;
  while(!(u&1)) u>>=1,t++;
  for(int i=0;i<=10;i++){
    ll a=prm[i],x1=qpow(a,u,n),x2;
    for(int j=1;j<=t;j++){
      x2=mul(x1,x1,n);
      if(x2==1&&x1!=1&&x1!=n-1) return false;//二次探测定理
      x1=x2;
	}
	if(x1!=1) return false;
  }
  return true;
}

5.5 分解质因数（Pollard-rho）

判断完了素数，我们来尝试分解质因数。

众所周知，可以用试除法做到 $O(\sqrt{n})$，但是同样 $n$ 一旦大了就不行了。

所以我们考虑一种更高消的方法——Pollard-rho 启发式方法 分解质因数。

事实上，我们还有一种奇妙的方法来寻找一个合数的因子（from TQX）:

scanf("%d",&n);
vector<int> fac;
for(int i=1;i<=???;++i){
	int a=rand()%(n-1)+1;
	if(n%a==0) fac.push_back(a);
}

这个代码试图通过随即一些数并判断它们是否是 $n$ 的因数，很明显这是一个非常蠢的方法，因为单次操作找到 $n$ 的因数的概率实在是太小了。

这个算法非常的劣质，但我们的 Pollard_rho 正是基于这个算法而来的，它的核心思想就是随机。

而我们知道 生日悖论：

如果现在有 $1000$ 个数，在里面随机到 $37$ 的概率非常小，但如果我们找两个数 $i,j$ 使 $|i-j|=37$ 的概率就要大一倍；

而当 $k=30$ 时，概率已经超过 $50\%$ ，当 $k=100$ 时，概率已经到了 $99.99\%$。

这个悖论告诉我们：组合随机采样比单个个体满足条件要高得多，这样可以提高正确率。

于是，Pollard 就提出了一个随机数的生成方法：

\[f(x)=(x^2+c) \bmod n \]

其中 $c$ 可以随机生成，在随机生成了一个 $x_1$ 后， $x2=f(x1),x3=f(x2),\cdots$。

这个函数的随机性非常好，但是进行了数次生成之后就会出现一个环，在坐标系中表示出来就变成了这个样子：

发现这个图像很像 $\rho$，于是这就是它的名字由来。

那么考虑如何判环？

可以用 Floyd 的方法：

假设我们有 $2$ 个函数生成的值 $a,b$，我们每次让 $a$ 走一步， $b$ 走两步，当 $a,b$ 相等时就出现了环。

于是我们可以用这样随机出来的数相邻两个的差与 $n$ 取 gcd 来分解 $n$，即如果得到的 gcd $\neq 1$ ，那么就直接返回。

但是 gcd 的时间复杂度时 $O(\log n)$ ，每一次都去调用一定会很慢。

在计算时，我们将每次产生的 $abc(a-b)$ 相乘并积累下来，最后直接判断这个乘积与 $n$ 的 $\gcd$。

但是，如果某一时刻累积下来的样本的乘积为 $0$ 了，为例不让样本丢失，我们直接退出循环进行判断即可。

每次计算的阈值可以倍增，并且加上一个上限，使用 $128$ 是一种不错的选择（玄学）。

至此，你就可以通过 P4718 【模板】Pollard-Rho 了。代码。

ll mul(ll a,ll b,ll m){ll z=(long double)a/m*b,res=(a*b-z*m)%m;return (res+m)%m;}
ll add(ll a,ll b,ll n){return (a+b>=n)?a+b-n:a+b;} 
ll qpow(ll a,ll b,ll m){
  ll res=1ll;
  while(b){
  	if(b&1) res=mul(res,a,m);
  	a=mul(a,a,m);
  	b>>=1;
  }
  return res;
}

const ll prm[12]={2,3,5,7,11,13,17,19,23,29,31};
bool miller_rabin(ll n){
  if(n<2) return false;
  for(int i=0;i<=10;i++){
  	if(n==prm[i]) return true;
  	if(n%prm[i]==0) return false;
  } 
  ll u=n-1,t=0;
  while(!(u&1)) u>>=1,t++;
  for(int i=0;i<=10;i++){
    ll a=prm[i],x1=qpow(a,u,n),x2;
    for(int j=1;j<=t;j++){
      x2=mul(x1,x1,n);
      if(x2==1&&x1!=1&&x1!=n-1) return false;
      x1=x2;
	}
	if(x1!=1) return false;
  }
  return true;
}

ll gcd(ll a,ll b){return (b==0)?a:gcd(b,a%b);}
ll f(ll a,ll c,ll n){return add(mul(a,a,n),c,n);} 
ll rad(ll n){return 1ll*rand()*rand()%n*rand()%n+1ll;}

ll pollard_rho(ll n){
  if(n==4) return 2;
  ll x=rad(n-1),c=rad(n-1),y=x;
  x=f(x,c,n);y=f(x,c,n);
  for(int lim=1;x!=y;lim=min(128,lim<<1)){
  	ll cur=1ll,nxt=0;
  	for(int i=1;i<lim;i++,cur=nxt){
  	  nxt=mul(cur,abs(add(x,n-y,n)),n);
  	  if(!nxt) break;//x=y	
	  x=f(x,c,n);y=f(f(y,c,n),c,n);
	}
	ll d=gcd(cur,n);
	if(d!=1) return d; 
  }
  return n;
}

void sol(ll n){
  ll d=pollard_rho(n);
  while(d==n) d=pollard_rho(n);
  if(miller_rabin(d)) ans=max(d,ans);
  else sol(d);
  if(miller_rabin(n/d)) ans=max(ans,n/d);
  else sol(n/d);
}

void wrk(){
  n=read();
  if(miller_rabin(n)){puts("Prime");return;}
  ans=0;sol(n);print(ans);
}

P4358 [CQOI2016] 密钥破解

一种非对称加密算法的密钥生成过程如下：

1.任选两个不同的质数$p,q$

2.计算$N=p \times q$，$r=(p-1)(q-1)$

3.选取小于$r$，且与$r$互质的整数$e$

4.计算整数$d$，使得$ed≡1(mod r)$

5.二元组$(N,e)$称为公钥，二元组$(N,d)$称为私钥。

当需要加密消息$n$时，（假设$n$是一个小于$N$的整数，因为任何格式的消息都可转为整数表示），使用公钥$(N,e)$，按照

\[n^e≡c(mod N) \]
运算，可得到密文$c$

对密文$c$解密时，用私钥$(N,d)$，按照

\[c^d≡n(mod N) \]
运算，可得到原文 $n$。算法正确性证明省略。

由于用公钥加密的密文仅能用对应的私钥解密，而不能用公钥解密，因此称为非对称加密算法。通常情况下，公钥由消息的接收方公开，而私钥由消息的接收方自己持有。这样任何发送消息的人都可以用公钥对消息加密，而只有消息的接收方自己能够解密消息。

现在，你的任务是寻找一种可行的方法来破解这种加密算法，即根据公钥破解出私钥，并据此解密密文。

$e,N,c \le 2^{62},c<N$。

题意写得非常长，但是并不难。

发现只要我们找到 $p,q$，就做完了，然而这个 $p,q$ 可以直接交给 Pollard_rho 即可。代码。

P4714 「数学」约数个数和

给你一个正整数 $N$，请计算 $N$ 的 $($所有约数的$)\times K$ 约数个数和。

答案可能很大，请输出对 $998244353$ 取模的结果。

$1 \leq N \leq 10 ^ {18}, 0 \leq K \leq 10^{18} $。

首先，对于 $K=0$ 的情况，$f_0(x) = \sigma(x)$ 是积性函数。

而当 $K \gt 0$ 时，我们有

\[f_i(n) = \sum_{d \mid n} f_{i-1}(n) \]

这相当于一个狄利克雷卷积，而根据狄利克雷卷积卷积的性质，我们可以知道 $f_i(n)$ 也是积性函数。

于是现在问题就变成了如何求 $f_K(p^q)$，因为我们把 $n$ 分解之后就可以得到若干个 $f_K(p^q)$ 的乘积。

发现这是简单的，等价于组合数学中在 $q+1$ 个空位中插入 $K+1$ 个板，可以重复。这时经典问题，你把每两个板之间就加一个球就可以得到

\[\binom {q+K+1} {K+1} \]

根据对称公式，相当于

\[\binom {q+K+1} {q} \]

而 $q \le 62$，所以可以直接用下降幂暴力计算。

至于分解质因数，交给 Pollard_Rho 就可以了。代码。

6 离散对数问题

离散对数问题就是解决模 $p$ 意义下的 $\log_ab$，也就等价于解离散对数方程

\[a^x \equiv b \pmod p \]

这样的 $x$ 就是 $\log_a b$，$x$ 可能存在多个，也可能一个都不存在，我们希望找到这样的一个 $x$。

当 $a \bot p$ 的时候，我们可以用大步小步算法解决，反之则可以用扩展大步小步算法解决。

6.1 大步小步算法（BSGS）

大步小步算法全称 Baby Step Giant Step，简称 BSGS，适用于 $a \bot p$ 的情形。

整体思路用到了 根号平衡 的思想，设块长为 $M$，那么最终的 $x$ 一定可以被表示成 $x=AM-B,0 \le B \lt M$。

于是这个方程就可以被表示成

\[a^{AM} \equiv ba^{B} \pmod p \]

只要我们直接暴力枚举 $A$ 和 $B$，用哈希表（pb_ds）存一下是否出现过，就可以用 $\mathcal O(\max(A,B))$ 的时间复杂度之内计算出答案。

根据欧拉定理（稍后会讲），有解就一定在 $[0,p-1]$ 中有至少一个解，所以这里 $0 \le B \lt M,0 \le A \le \left \lceil \frac p M \right \rceil$。

那么设块长为 $\left \lceil \sqrt p \right \rceil$ 时，复杂度最优，为 $\mathcal O(\left \lfloor \sqrt p \right \rfloor)$。

于是你就可以通过 P3846 [TJOI2007] 可爱的质数/【模板】BSGS 啦！代码。

int BSGS(int a,int b,int H){
  int A=1,B=sqrt(H)+1;a%=H,b%=H;
  if(H==1||b==1) return 0;
  if(!a) return !b?1:-1;
  gp_hash_table<int,int> mp;
  for(int i=1;i<=B;i++) mp[1ll*(A=1ll*A*a%H)*b%H]=i;
  for(int i=1,AA=A;i<=B;i++,AA=1ll*AA*A%H)
    if(mp.find(AA)!=mp.end()) return i*B-mp[AA];
  return -1;
}

容易发现你从小到大枚举时就可以得到最小的非负整数解，而关于 $x$ 的循环性质我们将会在后面阶与原根中讨论（循环的长度为 $\delta_a(p)$）。

在这里还有一些板题：

P3306 [SDOI2013] 随机数生成器

最近小 W 准备读一本新书，这本书一共有 $p$ 页，页码范围为 $0 \sim p-1$。

小 W 很忙，所以每天只能读一页书。为了使事情有趣一些，他打算使用 NOI2012 上学习的线性同余法生成一个序列，来决定每天具体读哪一页。

我们用 $x_i$ 来表示通过这种方法生成出来的第 $i$ 个数，也即小 W 第 $i$ 天会读哪一页。这个方法需要设置 $3$ 个参数 $a,b,x_1$，满足 $0\leq a,b,x_1\lt p$，且 $a,b,x_1$ 都是整数。按照下面的公式生成出来一系列的整数：

\[x_{i+1} \equiv a \times x_i+b \pmod p \]

其中 $\bmod$ 表示取余操作。

但是这种方法可能导致某两天读的页码一样。

小 W 要读这本书的第 $t$ 页，所以他想知道最早在哪一天能读到第 $t$ 页，或者指出他永远不会读到第 $t$ 页。

$1 \le T \le 50,0 \le a,b,x_1,t \lt p,2 \le p \le 10^9,p$ 是质数

考虑先全部左移一维，那么得到的数其实就是 $x_0$，我们令 $s_i = \frac {x_i} {a^i}$，那么根据递推式，我们可以得到 $a^i s_i = a^i s_{i-1} + b$，所以 $s_i = s_{i-1} + \frac b {a^i}$ 且 $s_0 = x_0$，则 $s_i = x_0 + \sum_{j=1}^i \frac b {a^j}$。

于是我们就可以将 $x_i$ 表示出来：

\[\begin{aligned} x_i & = a^i s_i \\ & = a^i x_0 + b \sum_{j=0}^{i-1} a^j \\ & = a^i x_0 + b \cdot \frac {a^i-1}{a-1} \\ & = \left(x_0 + \frac b {a-1} \right)a^i - \frac b{a-1} \end{aligned} \]

这就可以用 BSGS 直接完成啦——注意特判 $a=0,a=1$ 的情况。

总时间复杂度 $O(T\sqrt p)$。代码。

P4884 多少个 1？

给定整数 $K$ 和质数 $m$，求最小的正整数 $N$，使得 $11\cdots1$($N$ 个 $1$)$\equiv K \pmod m$。

说人话：就是 $111\cdots 1111 \bmod m = K$。

$6\leq m\leq 10^{11}$，$0< K< m$，保证 $m$ 是质数。

a 如果我们直接把 $11111 \cdots$ 分解出来是 $10^0+10^1+10^2 + \cdots$，这是非常不好处理的。

我们考虑将这个数表示成一个数，也就是它乘 $9$ 之后就会变成 $\frac{10^n - 1} 9$。

那么此时这个方程就变成了

\[\begin{aligned} \frac {10^n-1} 9 & \equiv K \pmod m \\ 10^n & \equiv 9K+1 \pmod m \end{aligned} \]

于是这个东西就可以直接用 BSGS 完成了。

注意可能在乘的过程中爆 long long。代码。

CF1106F Lunar New Year and a Recursive Sequence

有一串 $n$ 个数的数列，给你 $b_1\sim b_k$。当 $i>k$ 时：

\[f_i=(\prod\limits_{j=1}^kf_{i-j}^{b_j})\bmod{998244353} \]
已知 $f_1=f_2=\cdots=f_{k-1}=1,f_n=m$，问 $f_k$ 可能是多少。

$k \leq 100$，$n \leq 10^9$。

会用到原根相关知识，若不会原根建议先看后面原根部分。

有关于 $\prod$ 的式子是非常不好求的，我们考虑把它转化成有关加法的式子从而可以通过矩阵快速幂来完成。

而这种需要让我们想到了原根，众所周知，$3$ 是 $998244353$ 的一个原根，我们记它为 $G$，并记 $f_n = G^{g_n}$，于是我们就把式子转化成了

\[G^{g_i} = \prod_{j=1}^k G^{g_{i-j}b_j} \bmod 998244353 \]

把指数提下来，同时模数变成了 $\phi (p)$，我们就可以得到

\[g_i = \sum_{j=1}^k g_{i-j}b_j \bmod 998244352 \]

那么只要我们能求得 $g_i$ 那么求得 $f_i$ 也就是非常容易的了。

发现这个东西刚好是矩阵乘法的形式，也就是我们可以构造 $k \times k$ 矩阵

\[A = \left[ \begin{matrix} 0 & 0 & \cdots & 0 & b_k \\ 1 & 0 & \cdots & 0 & b_{k-1} \\ 0 & 1 & \cdots & 0 & b_{k-2} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & b_{1} \end{matrix} \right] \]

从而得到

\[\left[ \begin{matrix} 0 & 0 & \cdots & g_k \end{matrix}\right] A^{n-k}= \left[ \begin{matrix} g_{n-k+1} & g_{n-k+2} & \cdots & g_n \end{matrix}\right] \]

于是当我们知道 $g_k$ 时，可以快速算出 $g_n$ 的值。

而由于前面都是 $0$，所以实际上我们有

\[g_k \times A^{n-k}[k][k]\equiv g_n \pmod {998244352} \]

其中 $A^{n-k}[k][k]$ 是一个可以用矩阵快速幂处理出来的常数。

由于 $998244352$ 并不是质数，所以当我们知道 $g_n$ 的时候 $g_k$ 可以直接用 exGCD 求出。

现在问题就转化成如何求出 $g_n$，由于 $G^{g_n} \equiv f_n \pmod {998244353}$，容易发现这时直接可以用 BSGS 求出的。

于是这道题就做完了，实现时注意时模 $998244352$ 还是 $998244353$ 即可。代码。

6.2 扩展大步小步算法（exBSGS）

接下来我们来讨论更为普遍的问题 $a^x \equiv b \pmod p$，如果没有 $a \bot p$ 怎么处理？

从已知到未知： 由于 $a \bot p$ 的情形我们已经可以解决了，那么我们现在考虑把 $a,p$ 从不互质凑到互质。

我们运用同余方程的运算法则，我们把两边同时除以 $d= \gcd(a,p)$ 就可以得到

\[\frac a d a^{x-1} \equiv \frac b d \pmod {\frac p d} \]

由于 $\frac a d$ 和 $\frac p d$ 一定是互质的，所以在模 $\frac p d$ 的意义下 $\frac a d$ 是又逆元的，于是我们可以把 $\frac a d$ 除过去，就是

\[a^{x-1} \equiv \left( \frac b d \right) \times \left( \frac a d \right)^{-1} \pmod {\frac pd} \]

然而这时的 $a$ 和 $\frac pd$ 并不一定互质，比如说 $a=2,p=4$ 的情况，那么我们就重复上面的操作直到 $a \bot p$，而这个操作只会最多进行 $\log p$ 次，这时容易证明的，在 $a=2,p=2^k$ 取到极值。

而注意到每除依次，$x$ 就会 $-1$。于是当操作 $k$ 次后，我们需要计算其实是 $x-k$。此时如果 $b=0$，直接返回 $k$ 就可以了，这时需要特判的。

容易发现这样的时间复杂度是 $\mathcal O( \sqrt p + \log^2 p )$，后者有每一次 exgcd 求逆元的时间复杂度。

这样就可以通过 P4195 【模板】扩展 BSGS/exBSGS 和 SP3105 MOD - Power Modulo Inverted 了。代码。

int gcd(int a,int b){return b==0?a:gcd(b,a%b);}
void exgcd(int a,int b,int &x,int &y){
  if(b==0) return x=1,y=0,void();
  exgcd(b,a%b,y,x),y-=x*(a/b);
}
int inv(int a,int p){exgcd(a,p,x,y);return (x%p+p)%p;}

int BSGS(int a,int b,int H){
  int A=1,B=sqrt(H)+1;a%=H,b%=H;
  if(H==1||b==1) return 0;
  if(!a) return !b?1:-1;
  gp_hash_table<int,int> mp;
  for(int i=1;i<=B;i++) mp[1ll*(A=1ll*A*a%H)*b%H]=i;
  for(int i=1,AA=A;i<=B;i++,AA=1ll*AA*A%H)
    if(mp.find(AA)!=mp.end()) return i*B-mp[AA];
  return -1;
}

int exBSGS(int a,int b,int H){
  int d=gcd(a,H),cnt=0;a%=H,b%=H;
  while(d>1){
	if(b==1) return cnt;
	if(b%d) return -1;
	H/=d,b=1ll*(b/d)*inv(a/d,H)%H;
	d=gcd(H,a%=H),++cnt;
  }
  int res=BSGS(a,b,H);
  return ~res?res+cnt:-1;
}

P5345 【XR-1】快乐肥宅

给定 $n$ 组 $k_i,r_i,g_i$，求关于 $x$ 的方程组

\[\begin{cases} {k_1}^x & \equiv r_1 & \pmod {g_1} \\ {k_2}^x & \equiv r_2 & \pmod {g_2} \\ \cdots \\ {k_n}^x & \equiv r_n & \pmod {g_n} \end{cases} \]
（定义 $0^0=1$）在 $[0,10^9]$ 范围内的最小整数解，或判断在这个范围内无解。

$1 \le n \le 10^3,1 \le k_i,r_i \le g_i \le 10^7$。

难绷。

容易发现，我们的解题过程分成两部分：求解每一个方程的解 和 求线性同余方程组的解。

首先考虑第一部分，如何求解每一个方程的解呢？

根据 exBSGS 的求解过程可以知道，解的情况一定分成了三种：无解、唯一解和无穷解。

画出图来，解一定是一个 $\rho$ 形的，也就是这个样子（from 小粉兔）

其中 $\rho$ 的尾巴是我们依次 $\div \gcd$ 的过程，这种是只有一个解；

而进入环上就变成了无穷解，而这个解是形如线性同余方程；

反之，没有出现过的就是无解。

尾巴的长度是可以直接求 $a^x \equiv 1 \pmod b$ 的最小正整数解（这里和上面的板子会略有区别），这里的 a 是 exBSGS 化出来的值。那么我们就可以用一次 exBSGS 的过程求出第一个解的位置 $a$ 和循环节的长度 $b$，这就转化成了方程 $x \equiv a \pmod b$。

反之对于只有一个解的情况，我们直接记它为答案，最后特判一下就可以了。

接下来考虑第二部分，求 $n$ 个线性同余方程组的解。

由于不保证模数互质，所以这个需要用 exCRT 解决，而 $10^3$ 个 $10^7$ 是容易爆 long long 的，如何判断呢？

考虑 exCRT 的过程，中间每一次我们都是 ans+=x*M，那么只要这次的 $M \gt 10^9$，就不需要继续算后面的部分了，因为下一次一定就暴了。

所以这个时候我们只需要判断当前的解对于后面的方程是否都成立即可。

这样就做完了，时间复杂度为 $n$ 次 exBSGS 的复杂度。具体实现中有细节。代码。

posted @ 2024-05-21 13:10 H_W_Y 阅读(25) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

H_W_Y

用拼搏，为自己注解；用努力，为成功奠基