欧几里得算法（辗转相除法）
def gcd(a, b):
if b == 0:
return a
else:
return gcd(b, a % b)扩展欧几里得算法
def ext_euclid(a, b):
if b == 0:
return 1, 0, a
else:
x, y, q = ext_euclid(b, a % b)
x, y = y, (x - (a // b) * y)
return x, y, q大步小步算法（BSGS算法）
//C++
#include <iostream>
#include <cstdio>
#include <map>
using namespace std;
typedef long long LL;
LL quick_mod(LL a, LL b, LL c)//费马小定理+快速幂求逆元
{
LL ans = 1;
while (b)
{
if (b % 2 == 1)
ans = (ans*a) % c;
b /= 2;
a = (a*a) % c;
}
return ans;
}
int log_mod(int a, int b, int n)
{
int m, v, e = 1, i;
m = (int)sqrt(n + 0.5);
v = quick_mod(quick_mod(a, m, n),n-2, n);
map<int, int> x;
x[1] = 0;
for (int i = 1; i < m; i++) {
e = e*a%n;
if (!x.count(e)) x[e] = i;
}
for (i = 0; i < m; i++) {
if (x.count(b)) return i*m + x[b];
b = b*v%n;
}
return -1;
}#python
def BSGS(g, y, p):
m = int(sqrt(p))
if not is_square(p):
m += 1
S = {pow(g, j, p): j for j in range(m)}
gs = pow(g, inverse(m, p), p)
for i in range(m):
if y in S:
return i * m + S[y]
y = y * gs % p
return None扩展大步小步算法（扩展BSGS算法）
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
inline int gcd(int a,int b){
if(!b)
return a;
else{
while(int i=a%b){
a=b;
b=i;
}
return b;
}
}
inline int qpow(ll a,int n,int m) {
//这个快速幂保证p不是1，少模一次是一次
ll s=1;
while(n) {
if(n&1)
s=s*a%m;
a=a*a%m;
n>>=1;
}
return s;
}
unordered_map<int,int> M;
//要求a,n互质 a^x=b mod n .k,t是留给exbsgs调用的
int bsgs(int a,int b,int n,int k=1,int t=0) {
if(b==1)
return 0;
M.clear();
int m=ceil(sqrt(n));
ll s=b;//BS
for(int i=0; i<m; i++,s=s*a%n)
M[s]=i;
s=k;//GS
k=qpow(a,m,n);
for(ll i=1; i<=m; i++) {
s=s*k%n;
if(M.count(s))
return i*m-M[s]+t; //这样就保证找到的是最小解了
}
return -1;
}
//a^x=b mod n
int exbsgs(int a,int b,int n) {
if(b==1) {
return 0;
}
int d=gcd(a,n),k=1,t=0;
while(d^1) {
if(b%d) {
return -1;
}
++t;
b/=d;
n/=d;
k=(ll)k*(a/d)%n;
if(b==k) {
return t;
}
d=gcd(a,n);
}
return bsgs(a,b,n,k,t);
}
int main() {
int a,b,n;
while(1) {
scanf("%d%d%d",&a,&n,&b);
if(!a&&!n&&!b)
break;
a%=n;
b%=n;
int ans=exbsgs(a,b,n);
if(ans==-1)
puts("No Solution");
else
printf("%d\n",ans);
}
return 0;
}威尔逊定理
欧拉定理
费马小定理
费马商
中国剩余定理（孙子定理 / CRT）
#sage
def chinese_remainder(modulus, remainders):
Sum = 0
prod = reduce(lambda a, b: a*b, modulus)
for m_i, r_i in zip(modulus, remainders):
p = prod // m_i
Sum += r_i * (inverse_mod(p,m_i)*p)
return Sum % prod
chinese_remainder([3,5,7],[2,3,2]) #23扩展中国剩余定理（扩展CRT / ExCRT）
#互质与不互质两种情况下都能工作良好的中国剩余定理（解同余方程组）
def GCRT(mi, ai):
# mi,ai分别表示模数和取模后的值,都为列表结构
assert (isinstance(mi, list) and isinstance(ai, list))
curm, cura = mi[0], ai[0]
for (m, a) in zip(mi[1:], ai[1:]):
d = gmpy2.gcd(curm, m)
c = a - cura
assert (c % d == 0) #不成立则不存在解
K = c / d * gmpy2.invert(curm / d, m / d)
cura += curm * K
curm = curm * m / d
return (cura % curm, curm) #(解,最小公倍数)变种1：Noisy CRT
参考：
Noisy Polynomial Interpolation and Noisy Chinese Remaindering
DeadSec CTF 2023 - Loud系列
费马因式分解
import gmpy2
from Crypto.Util.number import *
def factor(n):
a = gmpy2.iroot(n, 2)[0]
while 1:
B2 = pow(a, 2) - n
if gmpy2.is_square(B2):
b = gmpy2.iroot(B2, 2)[0]
p = a + b
q = a - b
return p, q
a += 1from gmpy2 import *
def fermat_factorization(n):
factor_list = []
get_context().precision = 2048
x = int(sqrt(n))
print(x)
while True:
x += 1
y2 = x ** 2 - n
if is_square(y2):
print('x = ',x)
y2 = mpz(y2)
get_context().precision = 2048
y = int(sqrt(y2))
factor_list.append([x+y, x-y])
if len(factor_list) == 2:
break
return factor_list高斯整数
z =
GI = GaussianIntegers()
GI(z).factor()
#因子组合CopperSmith攻击
#Sage
#单元
PR.<x> = PolynomialRing(Zmod(n))
f = (a + x)^e - c
root = f.small_roots(X=2^256, beta=1)[0] # find root < 2^256 with factor = n
#调参，增大格
#beta=0.48, epsilon=0.02
#多元
import itertools
def small_roots(f, bounds, m=1, d=None):
if not d:
d = f.degree()
R = f.base_ring()
N = R.cardinality()
f /= f.coefficients().pop(0)
f = f.change_ring(ZZ)
G = Sequence([], f.parent())
for i in range(m+1):
base = N^(m-i) * f^i
for shifts in itertools.product(range(d), repeat=f.nvariables()):
g = base * prod(map(power, f.variables(), shifts))
G.append(g)
B, monomials = G.coefficient_matrix()
monomials = vector(monomials)
factors = [monomial(*bounds) for monomial in monomials]
for i, factor in enumerate(factors):
B.rescale_col(i, factor)
B = B.dense_matrix().LLL()
B = B.change_ring(QQ)
for i, factor in enumerate(factors):
B.rescale_col(i, 1/factor)
H = Sequence([], f.parent().change_ring(QQ))
for h in filter(None, B*monomials):
H.append(h)
I = H.ideal()
if I.dimension() == -1:
H.pop()
elif I.dimension() == 0:
roots = []
for root in I.variety(ring=ZZ):
root = tuple(R(root[var]) for var in f.variables())
roots.append(root)
return roots
return []
PR.<a, b> = PolynomialRing(Zmod(n))
f = 4*r^2*a*b + 2*r*(a+b) + 1 - n
roots = small_roots(f, (2^256, 2^256), m=3)
a, b = roots[0]
PR.<x, y> = PolynomialRing(Zmod(q)) # n = x*y
# PR.<x, y> = Polygen(RealField(1000)) # n ≈ x*y
f = (2^256 * a + x) * s - (2^256 + 1) * y * b - c
roots = small_roots(f, [2^256, 2^256], m=4, d=4)
#其他多元
load('coppersmith.sage')
P.<x, y> = PolynomialRing(GF(p))
f = 2^170 * a^2 + 2^86 * a * x + x^2 - 2^85 * b + c - y
roots = coron(f, X=2^85, Y=2^85, k=1, debug=True)Gröbner基
空间和域
理想
Gröbner基定义
#Sage
###ZZ/QQ/RR
#Example-1
P.<x, y> = PolynomialRing(QQ)
f1 = x^2 + x*y - 10
f2 = x^3 + x^2*y - 20
f3 = x^4 + x*y^3 - 70
G = Ideal([f1, f2, f3]).groebner_basis()
print(G)
#Example-2
PR = PolynomialRing(Zmod(N), 'x', len(Cs))
x = PR.gens()
f1 = (65537*x[0] - 66666*x[1] + 12345*x[2] - x[3])
f2 = x[0] + x[1] + x[2] - s
Fs = [f1, f2]
Fs.extend([(x[i]**e - Cs[i]) for i in range(l)])
I = Ideal(Fs)
B = I.groebner_basis()
print(B)
m = ''
for b in B:
assert b.degree() == 1
mi = ZZ(-b(0, 0, 0, 0))
print(mi)
###Zmod(p)
from sage.matrix.matrix2 import Matrix
def resultant(f1, f2, var):
return Matrix.determinant(f1.sylvester_matrix(f2, var))
P.<Rx, Ry, Qx, Qy> = PolynomialRing(Zmod(p))
f1 = Ry^2 - Rx^3 - a*Rx - b
f2 = Qy^2 - Qx^3 - a*Qx - b
f3 = (Qy + Ry)^2 + (Qx - Rx)^2 * (- Rx - Qx - Px)
f4 = (- Qy - Ry) * (Rx - Px) + (Qx - Rx) * (- Ry - Py)
f5 = Rx * Qx - N
G = Ideal([f1, f2, f3, f4, f5]).groebner_basis()
#结式+矩阵子式（西尔维斯特矩阵）
print('[!] computing resultant h1...')
h1 = resultant(G[0], G[1], Rx) # Ry, Qx, Qy
print('[!] computing resultant h2...')
h2 = resultant(G[0], G[2], Rx) # Ry, Qx, Qy
print('[!] computing resultant h3...')
h3 = resultant(G[3], G[4], Rx) # Ry, Qx, Qy
print('[!] computing resultant h4...')
h4 = resultant(G[3], G[5], Rx) # Ry, Qx, Qy
print('[!] computing resultant h5...')
h5 = resultant(h1, h2, Ry) # Qx, Qy
print('[!] computing resultant h6...')
h6 = resultant(h3, h4, Ry) # Qx, Qy
print('[!] computing resultant h7...')
h7 = resultant(h5, h6, Qy) # Qx
print('[!] computing resultant h8...')
h8 = resultant(h7, f5, Qx) # Rx
roots = h8.univariate_polynomial().roots()
p, q = [ZZ(t[0]) for t in roots if ZZ(t[0]).is_prime()]Diffie-Hellman密钥交换（DH密钥交换 / DHKE）
矩阵快速幂
int qpow(int x, int n, int m) {
int res = 1;
while (n) {
if (n & 1)
res = res * x % m;
x = x * x % m;
n >>= 1;
}
return res;
}根据矩阵乘法运算改写为矩阵快速幂：
#include <iostream>
using namespace std;
#define N 2
struct matrix {
int m[N][N];
matrix() {
memset(m, 0, sizeof(m));
}
void prt();
};
void matrix::prt() {
for (int i = 0; i < N; ++ i) {
for (int j = 0; j < N; ++ j) {
cout << this -> m[i][j] << " ";
}
cout << endl;
}
}
matrix operator * (const matrix a, const matrix b) {
matrix ans;
for (int i = 0; i < N; ++ i) {
for (int j = 0; j < N; ++ j) {
for(int k = 0; k < N; ++ k) {
ans.m[i][j] += a.m[i][k] * b.m[k][j];
}
}
}
return ans;
}
matrix qpow(matrix x, int n) {
matrix res;
for (int i = 0; i < N; ++ i) {
res.m[i][i] = 1;
}
while (n) {
if (n & 1) res = res * x;
x = x * x;
n >>= 1;
}
return res;
}
int fib(int n) {
matrix a;
a.m[0][0] = a.m[1][0] = a.m[0][1] = 1;
matrix base;
base.m[0][0] = 1;
matrix ans = qpow(a, n - 1);
ans = ans * base;
return ans.m[0][0];
}
int main() {
cout << fib(1) << endl; // 1
cout << fib(2) << endl; // 1
cout << fib(3) << endl; // 2
cout << fib(4) << endl; // 3
cout << fib(5) << endl; // 5
cout << fib(6) << endl; // 8
cout << fib(7) << endl; // 13
}
#以2x2矩阵相乘为例
m = [[1 for i in range(2)]for j in range(2)]
m[1][1] = 0
n = int(input())
def mulMatrix(x,y): #定义二阶矩阵相乘的函数
ans = [[0 for i in range(2)]for j in range(2)]
for i in range(2):
for j in range(2):
for k in range(2):
ans[i][j] += x[i][k] * y[k][j]
return ans
def quickMatrix(m,n):
E = [[0 for i in range(2)]for j in range(2)] #先定义一个单位矩阵
for i in range(2):
E[i][i] = 1
while(n):
if n % 2 != 0:
E = mulMatrix(E,m)
m = mulMatrix(m,m)
n >>= 1
return E
print(quickMatrix(m,n))佩尔方程 / Pell方程
def solve_pell(N, numTry = 100):
cf = continued_fraction(sqrt(N))
for i in range(numTry):
denom = cf.denominator(i)
numer = cf.numerator(i)
if numer^2 - N * denom^2 == 1:
return numer, denom
return None, None
N =
solve_pell(N)