算法,是彩色的。

介绍一些简单的数论算法,以后可能会补充。

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$

  1. 首先我们设 $a^x\equiv b\ (mod\ n)$ 中的 $x=im-j\in [0,\phi(n))$。这样一来,问题就转化为求解方程 $(a^m)^i\equiv a^jb\ (mod\ n)$。
  2. 枚举 $i$,将所有的 $(a^m)^i\ (mod\ n)$ 都存入容器内 (map/hash/sorted_array + binary search/…)。
  3. 枚举 $j$,在容器中查询 $a^jb\ (mod\ n)$:若匹配到,则 $im-j$ 就为该 DLP 的一个解。

任意模

如果 $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();
	}
}

例题

二次剩余 (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)\]

证明

参考 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)$

二次剩余的数量

在 $[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 的存在得到以上。

则模域内剩下的数无法得到,即为非二次剩余。

过程

  1. 随机选取 $t$,使得 $t^2 - a$ 为非二次剩余(选取 $t$ 的期望次数极小)
  2. 令 $b=\sqrt {t^2 - a}$,则 $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验证即可

引理

步骤

  1. 设 $p-1 = {2^k}·t $($t$ 为奇数)
  2. 随机选取 $a\in (1,p-1)$,从 $a^{2^0·t}$ 开始,直到 $a^{2^k·t}$ 为止,套用二次探测定理费马小定理验证
  3. 若验证失败,则为合数;否则多取几个随机 $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}$ 这个函数是玄学函数

  1. 随机选取 $a_0\in (1,N-1)$,不断计算得到 $a_i\equiv a^2_{i-1}+k\ (mod\ N)$($k$ 为常数)
    • 模域下,根据生日悖论,数列 ${a_i}$ 形成环且期望环长为 $\sqrt N$
  2. 判断是否满足 $gcd(abs(a_i-a_{i-1}),N) \in (1,N)$,若满足则找到 $N$ 的一个因子
    • 成环后仍未找到因子则表示分解失败,可调整 $k$ 重新分解
  3. 找到后通过 Miller rabin 算法判断因子是否为质数,若不是则使用该算法继续分解

优化

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$

  1. 设 $x=x_0+x_1q+x_2q^2+…+x_{e-1}q^{e-1}$
  2. 令 $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$

  3. 令 $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$

  1. 令 $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}$
  2. 则 $g_i^{x_i}=h_i$, 枚举得到 $x_i$
  3. 对 $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(the greatest one)

原理

设 $d_n=\Sigma^m_{i=1} C_iS_{(n-m-1)+i}$,其中 $m$ 是指阶数(长度)。那么对于一个给定的序列 $S$,我们的目标是找到最短的递推式 $C(m)=<C_1,C_2,…,C_m>$,其应当满足 $d_{all} = 0$。

特别地,$C(0)={}$,$d$ 第一次不为 0 时,$C(i)={0,…, 0}$。

这里直捣黄龙,给出key point:

设某一次失配时,递推式为 $B(m’)$,$b=d_k$;而当前失配时,递推式为 $C(m)$,修正后的递推式为 $C’(m)$。

那么,新的递推式就为 $C’(m)=C(m)-{d\over b}B(m’)$

这里 $B(m’)$ 的选择规则是:取决于 $min(n-Fail_Place_B+m’)$。这里取不取等号尚待探究,稍微有点小问题。

用说的不太清楚,这里举个例子:

BM

待补充

tips:

代码

时间复杂度 $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