Sagemath 教程

N1Key2024/2/6大约 15 分钟基本数论函数线性代数相关函数多项式环及其函数格相关函数椭圆曲线函数

原文链接： https://tl2cents.github.io/2022/03/27/SageMath简明教程/

启动命令： sage

打开 jupyter： sage -n jupyter

小技巧：

1、tab键可以自动补全函数

2、查看函数源码使用 “??” ，如：discrete_log??

3、查看函数用法使用 “?” ，如：discrete_log?

4、sage shell 中可以执行 Python 第三方库的安装命令，如：pip install pwntools

5、sage 中也可以使用 !pip install pwntools 命令安装

jupyter notebook使用

编辑与命令模式的切换

编辑模式 -> 命令模式，按下Esc键即可；
命令模式 -> 编辑模式，按下Enter键即可。

添加单元格： 在命令模式下，按 A 在当前单元格的前面添加单元格，按 B 是在当前单元格的后面添加单元格。

删除单元格： 在命令模式下，连续按下两次D。

添加行号： 在命令模式下，按下L。

1、执行单元格：ctrl + enter
2、补全代码：在编辑模式下，按下Tab键可以补全在前面已经出现过的字符。

sage输出markdown格式

sage 中的 latex() 函数是解决排版很好的方案，要通过xdvi在屏幕上显示内容，使用 view() 函数

sage 基本函数使用

sage与python的不同

5       # 5的类型是整数环上的元素，而不是int
sage: type(5)
<class 'sage.rings.integer.Integer'>

# 注意int(5)和5的方法是不一样的，如5的比特长方法
sage: int(5).bit_length()
3
sage: 5.nbits()
3

5**3	#  **代表幂
5^3     #  ^也是幂，python中是异或
5^^3    #  ^^表示sage中的异或
10/3    # 注意10/3是一个有理数，而在python中它会变成小数
sage: type(10/3)
<class 'sage.rings.rational.Rational'>

另外在编程中特别需要注意的一点是 $p o w (a, x, N)$ 运算，它的返回结果是 Zmod(N)上的元素，你之后不在在这个有限环上运算，最好转成正常整数 $Z Z (p o w (a, x, N))$ 再运算:

sage: type(pow(2,4,7))
<class 'sage.rings.finite_rings.integer_mod.IntegerMod_int'>

sage: type(ZZ(pow(2,4,7)))
<class 'sage.rings.integer.Integer'>

基本的环、域

整数环 ：ZZ，也等价于Integers

有理数环： QQ

实数域：RR

复数域：CC

多项式环：PolynomialRing()

# 定义在有理数环上的多项式环
sage: PR.<x> = PolynomialRing(QQ)

PR：多项式环的名字，自定义

x：变量的名字，自定义

其他有限域、环：

伽罗瓦域，又称有限域(N必须是某个素数或者某个素数的k次方)：GF(N)、GF(2^8,modulus=[1,0,0,1,1,1,0,0,1])

一般有限环：Zmod(N)

一个元素在各种域上的转换示例如下：

sage: e=12
sage: CC(e)
12.0000000000000
sage: ZZ(e)
12
sage: G=GF(7)
sage: G(e)
5
sage: Z7=Zmod(7)
sage: Z7(e)
5

基本数论函数

1.求逆

e 模 n 的逆： 定义在Zmod(N)上的元素e，直接e^-1或者1/e

否则直接使用 inverse_mod(e,n) 函数获得 e 模 n 的逆

2.最大公因数

（a,b)的最大公因数：gcd(a,b)

3.最大公倍数

(a,b)的最大公倍数：lcm(a,b)

4.模幂

$e^{x} m o d n$ ：定义在 Zmod(N) 上的元素e，直接 e^x

否则直接使用pow(e,x,n)

5.素数判断

sage: x=123721984710347123123
sage: x.is_prime()
False

6.阶乘

factorial(x)

7.欧拉函数：

euler_phi(n) 							
#  求n的欧拉函数 phi = n*prod([1 - 1/p for p in prime divisors(n)])

8.中国剩余定理求解

{\begin{cases} x \equiv m_{1} m o d n_{1} \\ x \equiv m_{2} m o d n_{2} \end{cases}

解为：

# crt([m1,m2],[n1,n2])
sage: crt([1,2],[3,5])
7

9.求质数

# 下一个质数：
sage: next_prime(2005)
2011

#上一个质数
sage: previous_prime(2001)
1999

10.扩展欧几里得算法

$g = g c d (a, b) = u * a + v * b$

g,u,v=xgcd(a,b)

11.素数分解

factor(1024) -> 2^10
prime_divisors(1024) -> [2]
divisors() #所有因子

12.开根(有限域开根)

整数域开根： $x^{m} = y 已知 m ， y 求 x$

sage: y=87^8
sage: y.nth_root(8)
87

有限域开根： $x^{m} \equiv y m o d N 已知 m ， y ， N 求 x$

sage: y=pow(78,888,65537)
# 此时已经定义在有限域Zmod(65537)上
sage: y.nth_root(888)
65459
sage: x=y.nth_root(888)
sage: x+78
0

注意：开根有多解，nth_root只会返回一个解，如果需要得到所有解，可以加一个参数 all=True：

sage: x=y.nth_root(888,all=True)

当然，sage的有限域开根函数实现并不好，可以手动实现AMM算法或者使用更为强大的工具mma（mma解方程的性能绝对是世界顶尖水平的）。

13.离散对数

sage上实现了多种离散对数的求解方式，包括gp接口的离散对数求解。

参考sage 手册

对于求解问题： $a \equiv b a s e^{x} (m o d m)$

参数说明：
求解以 base 为底，a 的对数；ord 为 base 的阶，可以缺省;
operation 可以是 + 与 *，默认为 *；
bounds 是一个区间 (ld,ud)，它用于指定搜索的上界。具体地说，它表示搜索的范围，即在多大的范围内搜索可能的解 x。

其中 operation 为 + 时一般是应用在椭圆曲线（ECC）的离散对数（离散对数其实也就是求指数）。

1、通用的求离散对数的方法

discrete_log(a,base,ord,operation)

2、求离散对数的 Pollard-Rho 算法

discrete_log_rho(a,base,ord,operation)

3、求离散对数的Pollard-kangaroo算法（也称为lambda算法）

discrete_log_lambda(a,base,bounds,operation)

4、小步大步法。

bsgs(base,a,bounds,operation)

from sage.groups.generic import bsgs

代码示例：

#生成32位的素数p作为模数，int为32位，超过int要在数字后加L
p=random_prime(2L**32)
 
#定义有限域GF(p)
G=GF(p)
 
#找一个模p的原根
gp ('znprimroot('+str(p)+')')
#输出Mod(rt,p)，则x是模p的原根
g=G.primitive_element()

#生成私钥
x=G(ZZ.random_element(p-1))
 
#公钥y=g^x mod p，由于已经定义在GF(p)上，因此g**x就是g^x mod p
y=g**x
 
#计算离散对数的通用方法
discrete_log(y,g)==x
#计算离散对数的lambda方法
discrete_log_lambda(y,g,(floor(ZZ(x)/2),2*ZZ(x)))==x
#小步大步法计算离散对数
bsgs(g,y,(floor(ZZ(x)/2),2*ZZ(x)))==x

#n为质数或质数幂（线性筛Index Calculus）
R = Integers(99)
a = R(4)
b = a^9
b.log(a)
n=99
#或
x = int(pari(f"znlog({int(b)},Mod({int(a)},{int(n)}))"))
x = gp.znlog(b, gp.Mod(a, n))
# 代码修改自 Lazzaro @ https://lazzzaro.github.io

输出：

Mod(2, 3223324843)
True
True
True
9

线性代数相关函数

1.一般矩阵定义

# 直接定义ZZ矩阵并且初始化
A = Matrix(ZZ,[[1,2,3],[3,2,1],[1,1,1]]) 
A = matrix(Zmod(2),[[1,2,3],[3,2,1],[1,1,1]]) 
# 模2有限域上的矩阵
B=matrix(GF(2),A)

# 向量定义
Y = vector(ZZ,[0,-4,-1]) 
Y = vector(GF(2),[0,-4,-1]) 

# 定义矩阵但不初始化：GF（2）上的10*10矩阵，默认所有元素为0
m = matrix(GF(2), 10, 10)

2.分块矩阵定义

Latex 编写矩阵相关操作语法：https://zhuanlan.zhihu.com/p/266267223

[\begin{matrix} A & B \\ C & D \end{matrix}]

A = matrix(GF(2), 10, 10)
B = matrix(GF(2), 10, 30)
C = matrix(GF(2), 5, 10)
D = matrix(GF(2), 5, 30)
M=block_matrix([[A,B],[C,D]])
#等价于
M=block_matrix([[A,B],[C,D]],nrows=2)

3.矩阵运算

m1=matrix([[1,2], [1,3]])
m2=matrix([[2,1], [2,3]])
m1+m2
m1*m2
输出
[3 3]
[3 6]
[ 6  7]
[ 8 10]

转置：M.transpose()

求逆： M.inverse()或者 $M^{- 1}$

特征值：M.eigenvalues()

特征向量（右）：M.eigenvectors_right()

行列式： M.det()

求秩：M.rank()

4.解线性方程组

对于任意线性（模）方程组，sage都可以很方便地求解：

线 性 方 程 组 ： {\begin{matrix} _{11} * x_{1} + a_{12} * x_{2} + . . . + a_{1 m} * x_{m} = b_{1} \\ a_{21} * x_{1} + a_{22} * x_{2} + . . . + a_{2 m} * x_{m} = b_{2} \\ . . . . \\ a_{n 1} * x_{1} + a_{n 2} * x_{2} + . . . + a_{n m} * x_{m} = b_{n} \end{matrix}

矩 阵 表 示 ： A * X = B ， 其 中 A \in B^{n * m}, X \in R^{m}, B \in R^{n}

其中 $R$ 可以替换成任意域： $Z m o d (N) 、 Z Z 、 Q Q$

给定有限域或者整数域上的矩阵 $A, B$ ，求解 $A * X = B$ ：

X = A.solve_right(B)
#等价于
X = A \ B

给定有限域或者整数域上的矩阵 $A, B$ ，求解 $X * A = B$ ：

X = A.solve_left(B)

特别地，sagemath还集成右（左）核空间的求解，也就是 $A * X = 0 和 X * A = 0$

X = A.left_kernel()  # X*A=0
X = A.right_kernel() # A*X=0
# 矩阵形式
X=A.right_kernel_matrix()

对于一般（模）方程的求解：

# 一般方程
x, y = var('x, y')
solve([x+y==10, x-y==0], x, y)

# 解模方程
sage: solve_mod([x+y==66, x-y==23],97)
[(93, 70)]

5.向量的计算

1、计算向量的模长

数学术语，norm 或 module，具有“长度”的概念

在 n 维欧几里德空间 $R^{n}$ 上，向量 $x = (x_{1}, x_{2}, . . ., x_{n})$ 的模长计算为：

$| x | = \sqrt{x_{1}^{2} + x_{2}^{2} + . . . + x_{n}^{2}}$

sage代码实现：

x = vector([3,6,2])
x.norm()  #square root of (9+36+4)
#7

2、计算向量的内积

# 计算两个向量的内积
sage: v1 = vector([1,2,3])
sage: v2 = vector([4,5,6])
sage: v1.dot_product(v2)	# 方法1
32
sage: v1.inner_product(v2)	# 方法2
32
sage: v1*v2		# 方法3
32

3、向量的数乘

可以使用 + 运算符执行向量加法。

# 定义向量
v1 = vector([1, 2, 3])
v2 = vector([4, 5, 6])

# 向量加法
result = v1 + v2
print(result)

在这个示例中，result 将会是向量 [5, 7, 9]。

4、向量数乘

可以使用 * 运算符执行向量数乘。

# 定义向量
v = vector([1, 2, 3])

# 向量数乘
result = 2 * v
print(result)

5、单位向量

v.normalized() 方法用于计算向量的单位向量。

v = vector([3, 4])
unit_vector = v.normalized()
print(unit_vector)

多项式环及其函数

1.多项式环定义

前面已经稍微涉及了多项式环的定义：

# 定义在有理数环上的一元多项式环
sage: PR.<x> = PolynomialRing(QQ)
# PR：多项式环的名字，自定义
# x：变量的名字，自定义

# 或者等价定义
sage: PR = PolynomialRing(QQ,'x')
sage: x=PR.gen()	# gen()函数创建一个多项式环中的生成元。
# 或者
sage: PR = QQ[‘t’]

# 或者
Zmod(p)[]  表示在模p下的多项式环
GF(p)[]   表示定义在有限域上的多项式环。

Zmod(p) 表示模 p 的整数环，其中 p 是一个素数。这个整数环中的元素是 $0$ 到 $p - 1$ 之间的整数，运算是模 $p$ 的加法和乘法。
[] 表示在这个环上构建的多项式环。也就是说，多项式的系数是取自模 $p$ 的整数环的元素。

二元多项式环可以如下定义，其他类似：

PR2.<x,y> = PolynomialRing(QQ)

2.多项式定义

在定义好多项式环后：PR.<x> = PolynomialRing(QQ)

PR.<x> = PolynomialRing(QQ)
f=4*x^3+2*x+1

一个等价的定义方式是用列表，生成随机多项式：

sage: R.<y> = PolynomialRing(GF(65537))
....: degree=7
....: S = R.random_element(degree)
sage: S
60064*y^7 + 31544*y^6 + 19047*y^5 + 19873*y^4 + 53675*y^3 + 39961*y^2 + 40289*y + 30493

特别地，对于伽罗瓦域GF(2^n)上元素的声明，如下（注意模多项式需要是本元多项式）：

F.<x> = GF(2^8,modulus=[1,0,0,1,1,1,0,0,1])
F.fetch_int(21) == F(21.bits())
F(21.bits()).integer_representation()#逆过程

3.多项式相关运算、函数

多项式分解：

sage: R.<y> = PolynomialRing(GF(65537))
sage: f=R.random_element(degree)*R.random_element(degree)
sage: f.factor()
(64343) * (y + 10474) * (y + 49729) * (y^2 + 7588*y + 41809) * (y^3 + 20885*y^2 + 60459*y + 29275) * (y^3 + 27396*y^2 + 16874*y + 56765) * (y^4 + 25258*y^3 + 12887*y^2 + 59574*y + 64190)

多项式最大公因式：

sage: g=R.random_element(degree)
sage: f1=R.random_element(degree)*g
sage: f2=R.random_element(degree)*g
sage: gcd(f1,f2)
y^7 + 10115*y^6 + 11715*y^5 + 59953*y^4 + 53514*y^3 + 2238*y^2 + 55526*y + 46381
sage: f1.gcd(f2)
y^7 + 10115*y^6 + 11715*y^5 + 59953*y^4 + 53514*y^3 + 2238*y^2 + 55526*y + 46381

首一多项式（首项系数被标准化为1）：

sage: f.monic()	# 整理多项式为最高次幂的系数为1
y^14 + 10256*y^13 + 12184*y^12 + 38005*y^11 + 59412*y^10 + 48187*y^9 + 45070*y^8 + 34871*y^7 + 50384*y^6 + 61263*y^5 + 52827*y^4 + 24231*y^3 + 4268*y^2 + 31669*y + 46586

多项式的根：

sage: f.roots()
[(55063, 1), (15808, 1)]

# small_roots,使用前必须保证f是首一的多项式
# 关于small_roots的详细参数可以参考coppersmith节
sage: f.monic().small_roots()
[]

多项式的结式：

在两个多项式有公共根时,可以求结式，之后再求解根

sage: f.resultant(g)
51209

Groebner basis：

格劳布纳基础在求解多项式方程组、理想的特征等方面发挥着重要作用，在SageMath中，可以使用groebner_basis()函数来计算一个理想的格劳布纳基础。

在已知模方程的一些关系后，我们可以通过 Groebner basis 得到一些方程的根，如下图得到的就是 $a, b, c$ 三个解

sage: G=GF(next_prime(getrandbits(512)))
sage: a=G(getrandbits(512))
sage: b=G(getrandbits(512))
sage: c=a+b+233
sage: a3=a^3
sage: b3=b^3
sage: c3=c^3
sage: x,y,z =G['x,y,z'].gens()
....: I = ideal(z-x-y-233,  x^3-a3 ,y^3-b3 , z^3 - c3)
....: B = I.groebner_basis()
verbose 0 (3837: multi_polynomial_ideal.py, groebner_basis) Warning: falling back to very slow toy implementation.
sage: B
[x + 4115365407930271672145917959064197986875879923337260019480399129522546957414196891761965008866716311717547095730704169414503027605502719768955802977966739, y + 1438196640478633988394708656097782124050629313931426401244798821904034601463695109291179031271158354211457757606425866100114378067643728009244830918500697, z + 5553562048408905660540626615161980110926509237268686420725197951426581558877892001053144040137874665929004853337130035514617405673146447778200633896467203]

格相关函数

sagemath集成了格上许多经典算法，在这里介绍crypto中常用的几个：

1.LLL算法

sage: M=matrix.random(GF(next_prime(getrandbits(10))),5,5)
sage: M=M.change_ring(ZZ)
sage: M.LLL()
[-62  55 -57  -9 -48]
[-59  58  74  22 -39]
[ 11   7 -39   8 125]
[144  60 -14 -41  24]
[ 41  25  17 176 -11]

2.coppersmith算法

也就是前文提到的 small_roots() 函数，其中对多项式有如下约束：必须是首一多项式（monic polynomial），可以使用 f=f.monic() 解决；f必须是一元多项式，暂不支持多元多项式。

# 这些参数的调整可以参加coppersmith节
x=f.small_roots(X=upperbound,beta=0.5,epsilon=1/20)

3.Babai算法

Babai算法sagemath上并没有集成的接口，这里给出一个CTF中crypto常用的实现：

# 输入矩阵 basis，目标向量，v
def approximate_closest_vector(basis, v):
    """Returns an approximate CVP solution using Babai's nearest plane algorithm.
    """
    BL = basis.LLL()
    # 施密特正交化基
    G, _ = BL.gram_schmidt()
    _, n = BL.dimensions()
    small = vector(ZZ, v)
    for i in reversed(range(n)):
        c = QQ(small * G[i]) / QQ(G[i] * G[i])
        c = c.round()
        small -= BL[i] * c
    return (v - small).coefficients()

椭圆曲线函数

1.椭圆曲线

通过命令 EllipticCurve 创建椭圆曲线

EllipticCurve([a1, a2, a3, a4, a6])：返回椭圆曲线 return the elliptic curve $y^{2} + a_{1} x y + a_{3} y = x^{3} + a_{2} x^{2} + a_{4} x + a 6$

EllipticCurve([a4, a6])：和上面的一样，但是 $a 1 = a 2 = a 3 = 0$ ，则 $y^{2} = x^{3} + a_{4} x + a_{6}$

sage: EllipticCurve([0,0,1,-1,0])
Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field

sage: EllipticCurve([3,4])
Elliptic Curve defined by y^2 = x^3 + 3*x + 4 over Rational Field

sage: EllipticCurve(GF(5),[0,0,1,-1,0])
Elliptic Curve defined by y^2 + y = x^3 + 4*x over Finite Field of size 5

$(0, 0)$ 是椭圆曲线 $E$ 上的一个点，定义为 $y^{2} + y = x^{3} - x$ 。在Sage类型 $E ([0, 0])$ 中创建这个点。Sage可以在这样的椭圆曲线上添加点(回想一下，椭圆曲线支持一种加法群结构，其中无穷多处的点是零元素，曲线上的三个共线点加为零)

sage: E = EllipticCurve([0,0,1,-1,0])
sage: E
Elliptic Curve defined by yˆ2 + y = xˆ3 - x over Rational Field

sage: P = E([0,0])
sage: P + P
(1 : 0 : 1)

sage: 10*P	# 标量乘法
(161/16 : -2065/64 : 1)

sage: 20*P
(683916417/264517696 : -18784454671297/4302115807744 : 1)

sage: E.conductor()
37

椭圆曲线ECC

# 初始化一条线
sage: ecc = EllipticCurve(GF(q), [a,b])

# 选定其上的一个点
sage: G = ecc(x,y)

# 基点G点的阶
# 生成点G的阶通常也被认为是整个椭圆曲线的阶
sage: G.order()

2.计算原根

# 得到3是模7乘法群的一个原根。
sage: gp('znprimroot(7)')
Mod(3, 7)

# 或者
sage: pari('znprimroot(7)')
Mod(3, 7)

3.无穷大

用于椭圆曲线 ECC中，定义无穷大

sage: oo
+Infinity

sage: Infinity
+Infinity

4.计算点的坐标

lift_x()：给定椭圆曲线 $E$ 上点 $x$ 的坐标，求出满足条件的所有点。

# 例题：
ECC = EllipticCurve(GF(2^255 - 19),[0,486662,0,1,0])
G = ECC.lift_x(9)
print(G)

# output (9 : 43114425171068552920764898935933967039370386198203806730763910166200978582548 : 1)

5. 具有有理系数的三元三次构造椭圆曲线

EllipticCurve_from_cubic(_F_, _P=None_, _morphism=True_)

F – 具有有理系数的三个变量的齐次立方，作为多项式环元素，定义平滑的平面三次曲线 $C$ 。
P – 3 元组 (x, y, z)定义一个投影点 $C$ ，或者 None。
morphism – 布尔值（默认值 True：）。如果 True 则返回一个双有理同构 $C$ 到 Weierstrass 椭圆曲线 $E$ ，否则只返回 $E$ 。

实现代码：

from Crypto.Util.number import *

p = 73997272456239171124655017039956026551127725934222347
d = 68212800478915688445169020404812347140341674954375635
c = 1
P = (71451574057642615329217496196104648829170714086074852, 69505051165402823276701818777271117086632959198597714)
Q = (40867727924496334272422180051448163594354522440089644, 56052452825146620306694006054673427761687498088402245)

# 曲线方程x^3+y^3+1=dxy
R.<x,y,z> = Zmod(p)[]   # 模p下的多项式环
cubic = x^3+y^3+c*z^3-d*x*y*z   # 构造有理系数的三元齐次立方
E = EllipticCurve_from_cubic(cubic,morphism=False) # 只返回一个Weierstrass形式的方程
EF = EllipticCurve_from_cubic(cubic,morphism=True) # 返回一个双有理同构$C$ 到 Weierstrass 椭圆曲线E

参考链接： https://doc.sagemath.org/html/en/reference/arithmetic_curves/sage/schemes/elliptic_curves/constructor.html#sage.schemes.elliptic_curves.constructor.EllipticCurve_from_cubic

6.Curve()函数

曲线构造函数，通过多项式方程构造一个有限域上的曲线。

sagemath 9.0 版本

Signature:      Curve(F, A=None)
Docstring:     
   输入：
	F – 多元多项式，或多项式列表或元组，或代数方案。
	A –（默认值：无）用于创建曲线的环境空间。

   EXAMPLES: A projective plane curve.

      sage: x,y,z = QQ['x,y,z'].gens()
      sage: C = Curve(x^3 + y^3 + z^3); C
      Projective Plane Curve over Rational Field defined by x^3 + y^3 + z^3
      sage: C.genus()
      1

7.singular_points()函数

返回该曲线的奇异点集。

举例：

A.<x,y,z> = AffineSpace(QQ, 3)
C = Curve([y^2 - x^5, x - z], A)
C.singular_points()

8.subs()函数

subs() 函数用于替换表达式中的变量或子表达式。这个函数通常用于代数计算中，可以帮助我们在表达式中替换特定的变量或子表达式，从而进行进一步的计算或简化。

subs() 函数接受一个字典作为输入，字典中包含了要替换的变量或子表达式及其对应的替换值。函数会将表达式中的变量或子表达式替换为指定的值，并返回一个新的表达式。

举例：

sage: R.<x,y> = QQbar[]
sage: f = x^2 + y + x^2*y^2 + 5
sage: f((5,y))
25*y^2 + y + 30
sage: f.subs({x:5})
25*y^2 + y + 30