对一些简单数论算法的介绍。

UPD 22.05.25 : 这篇博文感觉写得云里雾里的，重构应该提上日程。

BSGS (Baby Step Giant Step)
二次剩余 (Cipolla)
Miller rabin
- 引理
- 步骤
- 代码
Pollard Rho
- 原理
- 优化
DLP
- p is composite
- Pohlig-Hellman algorithm
  - 素数阶
  - 从素数阶到素数积阶
Factor n = p * q (close-prime)
- 代码
B-M Algorithm
莫比乌斯反演
- 前置知识
  - 积性函数
    - 定义
    - 性质
- 莫比乌斯函数
  - 定义
  - 性质 1
    - 证明
  - 性质 2
    - 证明
  - 线性筛
- 莫比乌斯反演
  - 约数形式
    - 证明
  - 倍数形式
- 例题
- 参考链接

算法，是彩色的。

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$ 再对应相减去：

![BM](https://coinc1dens.github.io/pics/algorithm/BM.png)

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} )\]

证明略。

例题

四元组统计

参考链接

OI-Wiki

https://blog.csdn.net/tomandjake_/article/details/81083703

https://www.cnblogs.com/ticmis/p/13210767.html