对一些简单数论算法的介绍。
UPD 22.05.25 : 这篇博文感觉写得云里雾里的,重构应该提上日程。
- BSGS (Baby Step Giant Step)
- 二次剩余 (Cipolla)
- Miller rabin
- Pollard Rho
- DLP
- Factor n = p * q (close-prime)
- B-M Algorithm
- 莫比乌斯反演
算法,是彩色的。
BSGS (Baby Step Giant Step)
拔山盖世 / 北上广深算法
概述
大步小步算法是一种解决 离散对数问题 (DLP) 的常用方法,也即求解方程:$a^x\equiv b\ (mod\ n)$ 中的 $x$。
其时间复杂度为 $O(\sqrt{n}log{\sqrt n})$,本质上就是一种 (以空间换时间/中间相遇攻击) 暴力方法。
(某些情况下,也可以用 Pohlig-Hellman Algorithm 等方法解决 DLP)
基础知识
欧拉函数
对于正整数 $n$,其欧拉函数是小于等于 $n$ 的正整数中与 $n$ 互质的数的数目,用 $\phi(n)$ 表示。
欧拉定理 (数论四大定理之一)
若 $n,a$ 为正整数,且 $gcd(n,a) = 1$,则 $a^{\phi(n)}\equiv 1\ (mod\ n)$
由欧拉定理知:若 $n,a$ 互素,那么 $a^{x+i\phi(n)}\equiv a^x\ (mod\ n)$,也就是说我们要求出的最小自然数 $x<\phi(n)<n$ 。
过程
前面说过这是一种暴力方法,也就是说,只要枚举出 $i\in [0,\phi(n))$ 的每一种情况,我们就能得到答案。这个时候采用中途相遇攻击,用空间换时间,就可有效降低时间复杂度。
素数模
指 $gcd(a,n)=1$
- 首先我们设 $a^x\equiv b\ (mod\ n)$ 中的 $x=im-j\in [0,\phi(n))$。这样一来,问题就转化为求解方程 $(a^m)^i\equiv a^jb\ (mod\ n)$。
- 枚举 $i$,将所有的 $(a^m)^i\ (mod\ n)$ 都存入容器内 (map/hash/sorted_array + binary search/…)。
- 枚举 $j$,在容器中查询 $a^jb\ (mod\ n)$:若匹配到,则 $im-j$ 就为该 DLP 的一个解。
- 时间复杂度:$O(i+j*O(query))$
- 为什么要用 map/others:查询复杂度低,仅为 $log$ 级别。
- $m$ 的取值:一般为 $\sqrt{n} + 1$。(平衡 $i,j$ 的枚举范围)
- 在不知道 $\phi(n) $ 的前提下,$(im-j)$ 需要覆盖 $[0,n)$;在以上取值下,$i,j$ 各枚举 $m$ 次即可。
- 这里也可以存储 $a^jb\ (mod\ n)$,枚举 $(a^m)^i\ (mod\ n)$;或者令 $x=im+j$ 等等。其核心思想只在于用空间换时间。
任意模
如果 $gcd(a,n)>1$,$x$ 的范围就更小了,虽然普通 $BSGS$ 也可以求出 $x$,但该如何优化呢?
我们设 $c=gcd(a,n)$,则原式可化为 ${a^k\over c^k}a^{x-k} \equiv {b\over c^k}\ (mod\ {n\over c^k})$ 直到 $gcd({n\over c^k},a)=1$。
(若 $c^k \nmid b$,那么该式一定无解)
令 $a’={a^k\over c^k},b’={b\over c^k},n’={n\over c^k},x’=x-k$,那么原式化为 $a’·a^{x’}\equiv b’\ (mod\ n’)$,那么既可以把它当成一个多了个系数的普通 BSGS 来解,也可以求 $a’$ 的逆元,化为 $a^{x’}\equiv b’·{a’}^{-1} (mod\ n’)$ 来做,区别不大。最后的结果是 $x=x’+k$。
代码实现
#include<bits/stdc++.h>
#define ll long long
using namespace std;
map<ll, ll>mp;
ll a, p, b, cnt, coefficient;
ll gcd(ll a,ll b){
return b == 0 ? a : gcd(b, a % b);
}
ll pow(ll a, ll b, ll c){
ll ans = 1;
while(b){
if(b & 1) ans = ans * a % c;
a = a * a % c;
b >>= 1;
}
return ans;
}
bool pre(){ // 处理底数与模数不互素的情况
if(b >= p) { puts("No Solution");return 0; }
if(!a){
if(!b) puts("1");
else puts("No Solution");
return 0;
}
ll c = gcd(a, p);
cnt = 0, coefficient = 1;
while(c > 1 and p % c == 0){
if(b % c){
if(b == coefficient) cout<<cnt<<endl;
else puts("No Solution");
return 0;
}
b /= c, p /= c, coefficient = coefficient * a/c % p;
c = gcd(a, p);
++ cnt;
}
return 1;
}
void bsgs(){ // 朴素bsgs
mp.clear();
ll m = (ll)sqrt(p) + 1, base = a, now = b;
for(int i = 0;i < m; ++ i){
mp[now] = i + 1; // 打标记,同时记录幂数
now = now * base % p;
}
base = pow(a, m, p), now = coefficient;
for(int i = 1;i <= m; ++ i){
if(mp[now = now * base % p]){
cout << i * m - (mp[now] - 1) + cnt << endl;
// x = x' + k = im - j + k
return;
}
}
puts("No Solution");
}
int main(){
while(cin>>a>>p>>b){
if(a + p + b == 0) return 0;
if(pre()) bsgs();
}
}
例题
- Luogu P4195 (BSGS 板子题,记得开 O2 优化) / SPOJ MOD
- SDOI 2011 计算器 (pow + exgcd + bsgs)
二次剩余 (Cipolla)
当存在某个 $x$,使 \(x^2\equiv a\ (mod\ p)\) 成立时,称 $a$ 是模 $p$ 的二次剩余
本方法对于 $p$ 为奇素数成立
引理
$(a+b)^p\equiv a^p + b^p\ (mod\ p)$
证明
\[(a+b)^p=\sum_{i=0}^p\binom{p}{i}a^ib^{p-i}\\ \binom{p}{i}\equiv 0\ (mod\ p)\\(i\ne p,i\ne 0,a,b\in \Bbb Z)\]勒让德符号
定义
\[\left(\frac{a}{p}\right)= 1,& \text{$a$ 是模 $p$ 的二次剩余} \\ -1,& \text{$a$ 是模 $p$ 的非二次剩余} \\ 0,& a\equiv0 \pmod p \\ ({a\over p})\equiv a^{ {p-1}\over 2}\ (mod\ p)\]- 其中 $p$ 为奇素数
证明
参考 https://www.cnblogs.com/3200Pheathon/p/10800065.html
设对 $p$,有原根 $g$,则知 ${g^1,…,g^{p-1}\%p}$ 两两不同,又 $g^{2i}={(g^i)}^2$,故这 $p-1\over 2$ 个数 ${g^i}$ 恰对应表示 $p$ 的所有正二次剩余; 又因为两两不同,故剩下的 ${g^{2k+1} }$ 即为所有的非二次剩余。
由费马小定理:$ (a^{ {p-1}\over 2}-1 )(a^{ {p-1}\over 2}+1)\equiv 0\ (mod\ p)$
故:$a^{ {p-1}\over 2}\equiv ±1\ (mod\ p)$
-
若为二次剩余:显然;
-
若为非二次剩余, 只能 $a^{ {p-1}\over 2}\equiv g^1\equiv -1\ (mod\ p)$,证毕。
二次剩余的数量
在 $[0,p)$ 中有 ${ {p-1}\over 2}+1$ 个(非二次剩余的数量是 ${ {p-1}\over 2}$)
证明
从正向考虑,两个不同的数 $x,y \in (0,p)$,若 $x^2\equiv y^2\equiv a\ (mod\ p)$,则 $(x-y)(x+y)\equiv 0\ (mod\ p)$,故知 $x+y\equiv 0\ (mod\ p)$,即 $y=p-x$,故这样的互不相同的 $a$ 恰有 ${ {p-1}\over 2}$ 对,加上 0 的存在得到以上。
则模域内剩下的数无法得到,即为非二次剩余。
过程
- 随机选取 $t$,使得 $t^2 - a$ 为非二次剩余(选取 $t$ 的期望次数极小)
- 令 $b=\sqrt {t^2 - a}$(实数域下,其与 $t$ 的运算类比复数运算进行书写),则 $x=(t+b)^{ {p+1}\over 2}$
证明
\[\begin{align} x^2&=(t+b)^{p+1}\\ &\equiv (t^{p+1}+b^{p+1})\\ &\equiv 1·t^2+(-1)·(t^2 - a)\\ &\equiv a\ (mod\ p) \end{align}\]代码(二次剩余)
期望复杂度$O(log n)$
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
long long T,t,n,p,vi;
struct node{
long long x,y;
};
node mul(node a,node b,long long c){
node ans = {0,0};
ans.x = (a.x * b.x % c + a.y * b.y % c * vi) % c;
ans.y = (a.x * b.y % c + a.y * b.x % c) % c;
return ans;
}
long long pow(long long a,long long b,long long c){
long long ans = 1;
while(b){
if(b & 1) ans = ans * a % c;
b >>= 1;
a = a * a % c;
}
return ans;
}
long long powp(node a,long long b,long long c){
node ans = {1,0};
while(b){
if(b & 1)ans = mul(ans,a,c);
b >>= 1;
a = mul(a,a,c);
}
return ans.x;
}
int main(){
scanf("%d",&T);
while(T--){
scanf("%d%d",&n,&p);
if(!n) puts("0");
else if(pow(n,p - 1 >> 1,p) != 1){
puts("Hola!");
continue;
}
else{
t = sqrt(n)+1;
while(pow(t * t - n,p - 1 >> 1,p) == 1) t++;
vi = t * t - n;
t = powp((node){t,1},p + 1 >> 1,p);
long long g = p - t;
if(!g or !t) puts("0");
else printf("%d %d\n",min(g,t),max(g,t));
}
}
}
代码 (Python 版)
def powp(a, b, e, mod, base):
ansa = 1
ansb = 0
while e:
if e & 1:
ansa, ansb = (ansa * a + ansb * b * base) % mod, (ansa * b + ansb * a) % mod
e >>= 1
a, b = (a * a + b * b * base) % mod, (a * b + b * a) % mod
return ansa
def Cipolla(a, q):
# assert isPrime(q) and q&1
# get m s.t. m**2 % q == a
t = random.randint(2, q - 1)
while pow(t**2 - a, q - 1 >> 1, q) == 1:
t = random.randint(2, q - 1)
base = t**2 - a
m = powp(t, 1, q + 1 >> 1, q, base)
return m
Miller rabin
int范围内使用2,3,5,7验证即可;
在 $\phi$ 已知的情况下,也可以用来分解 multi-prime RSA。
引理
-
费马小定理:
\[if\ p\ is\ a\ prime,\\ a^{p-1}\equiv 1\ (mod\ p)\] -
二次探测定理:
\[if\ p\ is\ a\ prime,\\and\ a^2\equiv 1\ (mod\ p),\\ a\equiv ±1\ (mod\ p)\]
步骤
- 设 $p-1 = {2^k}·t $($t$ 为奇数)
- 随机选取 $a\in (1,p-1)$,从 $a^{2^0·t}$ 开始,直到 $a^{2^k·t}$ 为止,套用二次探测定理及费马小定理验证
- 若验证失败,则为合数;否则多取几个随机 $a$ 进行验证
代码
期望复杂度$O(klog^2 n)$
#include<iostream>
using namespace std;
int p[5]={0,2,3,5,7};
long long pow(long long a,int b,int c){
long long ans = 1;
while(b){
if(b & 1)ans = ans * a % c;
b >>= 1;
a = a * a % c;
}
return ans;
}
bool check(long long a,int k,int x){
long long temp = a;
while(k--){
a = a * a % x;
if(a == 1 and temp - 1 != 0 and temp + 1 != x) return 0;
temp = a;
}
return temp == 1;
}
bool mb(int x){
if(x < 10) return (x == 2 or x == 3 or x == 5 or x == 7);
if(x % 2 == 0) return 0;
int k = 0,t = x - 1;
while(t % 2 == 0)
t >>= 1,
k ++;
for(int i = 1;i <= 4;++i)
if(!check(pow(p[i],t,x),k,x)) return 0;
return 1;
}
int main(){
int n,x;
cin>>n;
for(int i = 1;i <= n;++i){
cin>>x;
cout << (mb(x)?"Yes":"No") << endl;
}
}
Pollard Rho
反正我觉得挺玄学的
原理
随机,生日悖论
需要选取数的个数约在 $O(N^{1\over 4})$
${a_i}$ 这个函数是玄学函数
- 随机选取 $a_0\in (1,N-1)$,不断计算得到 $a_i\equiv a^2_{i-1}+k\ (mod\ N)$($k$ 为常数)
- 模域下,根据生日悖论,数列 ${a_i}$ 形成环且期望环长为 $\sqrt N$
- 判断是否满足 $gcd(abs(a_i-a_{i-1}),N) \in (1,N)$,若满足则找到 $N$ 的一个因子
- 成环后仍未找到因子则表示分解失败,可调整 $k$ 重新分解
- 找到后通过 Miller rabin 算法判断因子是否为质数,若不是则使用该算法继续分解
优化
-
路径倍长
某题解认为 127 是个不错的 loop
在对 $N$ 模域下,$gcd(\prod abs(a_i-a_{i-1}),N)$ 与非模域下相同,故每相隔一段时间后计算 $gcd(\prod abs(a_i-a_{i-1}),N)$ 以优化对于 $gcd$ 的时间复杂度
DLP
$g^x=h\ (mod\ p)$, Given $p(prime), g, h$, how to get $x$ ?
p is composite
直接分解 p, 对于每个方程各自枚举 $x_i$,最后对 $x\equiv x_i\ (mod\ (p_i-1))$ 进行 CRT 。
Pohlig-Hellman algorithm
素数阶
条件:若 $ord(g) = q^e$, 我们可以在约 $O(eq)$ 内得到 $x$
- 设 $x=x_0+x_1q+x_2q^2+…+x_{e-1}q^{e-1}$
-
令 $g_{e-1}=g^{q^{e-1} },h_{e-1}=h^{q^{e-1} }$,则 $g_{e-1} ^ x\equiv g_{e-1} ^ {x_0}\equiv h_{e-1}\ (mod\ p)$,枚举求出 $x_0$
- 令 $h_{e-2}=h^{q^{e-2} }·(g^{q^{e-2} x_0})^{-1}\ (mod\ p)$, 则解决 $g_{e-1}^{x_1}\equiv h_{e-2}\ (mod\ p)$ 即可;如此反复求出 $x$
从素数阶到素数积阶
条件:若 $ord(g)=N=q_1^{e_1}q_2^{e_2}…q_t^{e_t}$, 我们可以在约 $O(q_1^{e_1}+q_2^{e_2}+…+q_t^{e_t}+log_N)$ 内得到 $x$
- 令 $g_i=g^{N\over {q_i^{e_i } } },h_i=h^{N\over {q_i^{e_i } } }$,设 $x_i\equiv x\ (mod\ q_i^{e_i})$ 也即 $x=x_i+z_iq_i^{e_i}$
- 则 $g_i^{x_i}=h_i$, 枚举得到 $x_i$
- 对 $x_i$ 使用 CRT 即可得到 $x$
这里为了节约时间,也可以将第二步更换为求素数阶使用的方法,时间复杂度在 $O(e_1q_1+e_2q_2+…+e_tq_t+log_N)$ 左右,在 ECC 中就介绍的是这种攻击方法
参考链接:离散对数问题之 Pohlig-Hellman algorithm
Factor n = p * q (close-prime)
BSGS 的应用
这里的推导过程见:
https://grocid.net/2017/09/16/finding-close-prime-factorizations/
http://www.soreatu.com/essay/An%20Interesting%20Factorization%20Method%20on%20Close-prime.html
时空复杂度约为 $O(\phi_{approx} - \phi)$
代码
from gmpy2 import *
from Crypto.Util.number import *
def close_factor(n, b):
# approximate phi
phi_approx = int(n - 2 * sqrt(n) + 1)
# create a look-up table
look_up = {}
z = 1
zz = 1
for i in range(0, b + 1):
look_up[zz] = -i
z = (z * 2) % n
zz= invert(z, n)
# check the table
mu = invert(pow(2, phi_approx, n), n)
fac = pow(2, b, n)
j = 0
while True:
mu = (mu * fac) % n
j += b
if mu in look_up:
phi = phi_approx + (look_up[mu] - j)
break
if j > b * b:
return
m = n - phi + 1
roots = (m - int(sqrt(m ** 2 - 4 * n))) // 2, \
(m + int(sqrt(m ** 2 - 4 * n))) // 2
return roots
def Test():
while 1:
p = getPrime(40)
q = getPrime(40)
n = p * q
b = 2**39
# O(sqrt(b) + X)
# b is o(phi_approx - phi)
print(close_factor(n, int(sqrt(b) + 1000000)))
Test()
B-M Algorithm
用于生成满足条件的(级数最小的)线性递推式,密码学中常用于还原 LFSR。
参考链接:https://www.mina.moe/archives/10758、https://www.cnblogs.com/zzqsblog/p/6877339.html、https://www.cnblogs.com/zhouzhendong/p/Berlekamp-Massey.html、https://cdcq.github.io/2022/08/28/20220828a/
不待续了,有机会看原论文..
## 原理
设误差 $d_n=\Sigma^m_{i=1} C_iS_{(n-m-1)+i} - S_n$,其中 $m$ 是指阶数(长度)。那么对于一个给定的序列 $S$,我们的目标是找到最短的递推式 $C(m)=<C_1,C_2,...,C_m>$,其应当满足 $d_{all} = 0,all < n$,特别地,$C(0)=\{\}$。
这里直捣黄龙,给出 key point:
设修正后的递推式为 $C'(m)$,而在某一次失配时,递推式为 $B(m')$,$b=d_k$;当前失配时,递推式为 $C(m)$。
那么,新的递推式就为 $C'(m)=C(m)-{d\over b}B(m')$。
这里某一次失配的 $B(m')$ 的选择规则是:取决于 $min(n-Fail\_Place_B+m')$,它在最近一次引起递推式长度变长的位置。
用说的不太清楚,这里举个例子,也即将 $B(m^{'})$ 的部分偏移 $n-m$ 再对应相减去:

tips:
+ 归纳法证明该方法下的序列最短,其中需要用到引理 $$
+ 特别地,对于 LFSR,失配时直接 $C'(m)=C(m)+B(m') mod\ 2$
## 代码
> 时间复杂度 $O(n^2)$
代码 DEBUG 没看出来
莫比乌斯反演
前置知识
已知有 $k=\lfloor {n\over i}\rfloor$,求 $i_{max}\ s.t.\ k=\lfloor {n\over i_{max} }\rfloor$:
显然 $n=ki+r,r<i$,$i=i_{max}$ 时等价于 $r’=n\%k$,故 $i_{max}=\lfloor {n\over k}\rfloor$。
积性函数
定义
$f(x)$,满足 $\forall x,y \in N_+,gcd(x,y)=1$ 都有 $f(xy)=f(x)f(y)\ (f(1)=1)$。
若对 $gcd$ 无要求,则称之为完全积性函数。
性质
若 $f(x),g(x)$ 均为积性函数,则下列函数也是积性函数: \(h(x)=f(x^p)\\ h(x)=f^p(x)\\ h(x)=f(x)g(x)\\ h(x)=\sum_{d\mid x}f(d)\) 有关最后一条的证明:
设 $gcd(m,n)=1$,则 $d\mid mn$ 可唯一表示为 $d=d_1d_2,d_1\mid m,d_2\mid n$。
故 $h(mn)=\sum\limits_{d\mid mn}f(d)=\sum\limits_{d_1\mid m}\sum\limits_{d_2\mid n}f(d_1d_2)=(\sum\limits_{d_1\mid m}f(d_1))(\sum\limits_{d_2\mid n}f(d_2))=h(m)h(n)$
莫比乌斯函数
定义
\[\mu(n)=\left\{ \begin{array}{l} 1,& \mbox {n=1} \\ (-1)^k,& \mbox {$n=p_1p_2...p_k(e_i=1)$}\\ 0,& \mbox {else}\\ \end{array} \right.\]显然,$\mu(n)$ 为一积性函数。
性质 1
\[\sum_{d\mid n}\mu(d)=\left\{ \begin{array}{l} 1, & \mbox {$n=1$}\\ 0, & \mbox {$n\ne 1$} \end{array} \right.\]证明
设 $n=\prod\limits_{i=1}^k p_i^{e_i},n’=\prod\limits_{i=1}^k p_i$,显然有: \(\begin{aligned} \sum_{d\mid n}\mu(d) &= \sum_{d\mid n'}\mu(d)\\ &= \sum_{d\mid n'}C_k^i·(-1)^i\\ &= (1 + (-1))^k\\ &= 0 \end{aligned}\)
性质 2
若 $\theta(x)$ 为一积性函数,则:$\sum_{d\mid n}\mu(d)\theta(d)=(1-\theta(p_1))(1-\theta(p_2))…(1-\theta(p_k))$。
特别地,有:
性质 1、 $\theta(d)={1\over d},\sum\limits_{d\mid n}\mu(d)\theta (d)=\prod(1-{1\over p_i} )$
证明
由,$f(n)=\sum_{d\mid n}\mu(d)\theta(d)$ 为积性函数,则:$f(n)=\prod\limits_{i=1}^k f(p_i^{e_i} )=\prod (1+\mu(p_i)\theta(p_i)+…+\mu(p_i^{e_i} )\theta(p_i^{e_i} ) )=\prod (1-\theta(p_i))$。
线性筛
void getMu() {
mu[1] = 1;
for (int i = 2; i <= n; ++i) {
if (!flg[i]) p[++tot] = i, mu[i] = -1;
for (int j = 1; j <= tot && i * p[j] <= n; ++j) {
flg[i * p[j]] = 1;
if (i % p[j] == 0) {
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
}
莫比乌斯反演
约数形式
\[f(n)=\sum\limits_{d|n} g(d) \Rightarrow g(n)=\sum\limits_{d|n} \mu(d)· f(\frac{n}{d})\]证明
$\sum\limits_{d | n} \mu(d)· f(\frac{n}{d})=\sum\limits_{d | n}( \mu(d)\sum\limits_{k\mid {n\over d} }g(k))=\sum\limits_{k | n}( g(k)\sum\limits_{d\mid {n\over k} }\mu(d))=g(n)$,最后一步由 性质 1。 |
倍数形式
\[g(n)=\sum_{n|i} f(i) \Rightarrow f(n)=\sum_{n|i}g(i)\times \mu(\frac{i}{n} )\]证明略。
例题
参考链接
https://blog.csdn.net/tomandjake_/article/details/81083703
https://www.cnblogs.com/ticmis/p/13210767.html