sage_basic
Introduction¶
[!QUOTE]
Sage is free, open-source math software that supports research and teaching in algebra, geometry, number theory, cryptography, numerical computation, and related areas. Both the Sage development model and the technology in Sage itself are distinguished by an extremely strong emphasis on openness, community, cooperation, and collaboration: we are building the car, not reinventing the wheel. The overall goal of Sage is to create a viable, free, open-source alternative to Maple, Mathematica, Magma, and MATLAB.
Basic use¶
在简单语法上和 python 很类似，函数的定义和使用是一样的，就不讲了；如果想要在 python 中使用，就需要先导入 sage 库 from sage.all import * 即可。
下面的内容基于 ipynb/jupyter 形式实现，使用 sageMath 10.3 。
In [1]:
Copied!
from sage.all import *
2+2, 2**3, 5*5, 5**2, 2*2**3, 2**2**3**0
from sage.all import *
2+2, 2**3, 5*5, 5**2, 2*2**3, 2**2**3**0
Out[1]:
(4, 8, 25, 25, 16, 4)
In [2]:
Copied!
factor(2260305277112127131723251710487571168300606058081777971297059135638227918807397016980532281858292026960397175871392336597693440)
factor(2260305277112127131723251710487571168300606058081777971297059135638227918807397016980532281858292026960397175871392336597693440)
Out[2]:
2^10 * 5 * 7^2 * 37 * 193 * 331 * 569 * 5153 * 65381 * 146317 * 155943343 * 9206678941 * 18623046139 * 1290117019611958829 * 6846725164570923688379767 * 575390540597317088017221077
[!QUOTE]
As the last example shows, some mathematical expressions return ‘exact’ values, rather than numerical approximations. To get a numerical approximation, use either the function N or the method n (and both of these have a longer name, numerical_approx, and the function N is the same as n). These take optional arguments prec, which is the requested number of bits of precision, and digits, which is the requested number of decimal digits of precision; the default is 53 bits of precision.
In [3]:
Copied!
exp(2), n(exp(2)), numerical_approx(pi, prec=200)
exp(2), n(exp(2)), numerical_approx(pi, prec=200)
Out[3]:
(e^2, 7.38905609893065, 3.1415926535897932384626433832795028841971693993751058209749)
In [4]:
Copied!
x = var('x')
f = 1/((1+x)*(x-1))
f.partial_fraction(x)
x = var('x')
f = 1/((1+x)*(x-1))
f.partial_fraction(x)
Out[4]:
-1/2/(x + 1) + 1/2/(x - 1)
get help¶
[!QUOTE]
Sage has extensive built-in documentation, accessible by typing the name of a function or a constant (for example), followed by a question mark (?); alse, we can use the
help()function.
data structure¶
和 python 一样，sage 也有列表，字典，元组，集合等数据结构。
In [6]:
Copied!
list(range(2, 20, 3)), vector(range(2, 20, 3))
list(range(2, 20, 3)), vector(range(2, 20, 3))
Out[6]:
([2, 5, 8, 11, 14, 17], (2, 5, 8, 11, 14, 17))
不同的是，Sage math 的列表支持不同类型的元素，而 Python 的列表只支持一种类型；同样的，使用 len() 获取列表长度，使用 append() 添加元素，使用 remove() 删除元素。
In [7]:
Copied!
l = [1, 3.2, "hello", (2, 3), [1, 2, 3]]
for i in l:
print(i, end=" ")
l = [1, 3.2, "hello", (2, 3), [1, 2, 3]]
for i in l:
print(i, end=" ")
1 3.20000000000000 hello (2, 3) [1, 2, 3]
字典等也是同理，支持不同类型的元素。
In [8]:
Copied!
d = {1: "one", 'pi': 3.1415 , "admin": "darstib"}
for k in d.keys():
print(k, d[k], end=" # ")
d = {1: "one", 'pi': 3.1415 , "admin": "darstib"}
for k in d.keys():
print(k, d[k], end=" # ")
1 one # pi 3.14150000000000 # admin darstib #
In [9]:
Copied!
x = var('x')
y = sin(x)
diff(y, x), integrate(y, x), limit(y, x=0), limit(y, x=oo), solve(y == 0, x)
x = var('x')
y = sin(x)
diff(y, x), integrate(y, x), limit(y, x=0), limit(y, x=oo), solve(y == 0, x)
Out[9]:
(cos(x), -cos(x), 0, ind, [x == 0])
求解器 solve ¶
当然也可以创建多个变量，这和 z3 求解器是类似的。
In [10]:
Copied!
a, b, c = var('a, b, c')
solve([a + b + c == 1, a + 2*b + 4*c == 2, a + 3*b + 9*c == 3], a, b, c)
a, b, c = var('a, b, c')
solve([a + b + c == 1, a + 2*b + 4*c == 2, a + 3*b + 9*c == 3], a, b, c)
Out[10]:
[[a == 0, b == 1, c == 0]]
我们可以用这个来求线性方程组等等；对于多解情况，可以直接输出，也可以存起来以便控制输出精度。
In [11]:
Copied!
var('x y p q')
eq1 = p+q==9
eq2 = q*y+p*x==-6
eq3 = q*y**2+p*x**2==24
solns = solve([eq1,eq2,eq3,p==1],p,q,x,y, solution_dict=True)
[[s[p].n(10), s[q].n(10), s[x].n(10), s[y].n(10)] for s in solns]
var('x y p q')
eq1 = p+q==9
eq2 = q*y+p*x==-6
eq3 = q*y**2+p*x**2==24
solns = solve([eq1,eq2,eq3,p==1],p,q,x,y, solution_dict=True)
[[s[p].n(10), s[q].n(10), s[x].n(10), s[y].n(10)] for s in solns]
Out[11]:
[[1.0, 8.0, -4.9, -0.14], [1.0, 8.0, 3.6, -1.2]]
求根 find_root ¶
有时候，solve 能够找到一个解，但是无法给出粗略值，这时候就可以使用 find_root 函数。
In [12]:
Copied!
var('x')
solve(sin(x) == cos(x), x), find_root(sin(x) == cos(x), 0, pi) # 限制范围在 0 到 pi 之间
var('x')
solve(sin(x) == cos(x), x), find_root(sin(x) == cos(x), 0, pi) # 限制范围在 0 到 pi 之间
Out[12]:
([sin(x) == cos(x)], 0.7853981633974484)
微积分 ¶
在前面我们也提到过，在密码学中，我们经常需要对函数进行求导和积分，在 SageMath 中，我们也可以很方便地进行求导和积分。 直接使用 diff 共可接收三个参数，按照我们习惯的说法就是 (f(x), x, n) 其中 n 表示微分阶数。
In [14]:
Copied!
diff(sin(x^2), x, 4)
diff(sin(x^2), x, 4)
Out[14]:
16*x^4*sin(x^2) - 48*x^2*cos(x^2) - 12*sin(x^2)
多元微分咋办？
In [15]:
Copied!
x, y = var('x,y')
f = x^2 + 17*y^2
f.diff(x), f.diff(y)
x, y = var('x,y')
f = x^2 + 17*y^2
f.diff(x), f.diff(y)
Out[15]:
(2*x, 34*y)
（不）定积分怎么办？
In [16]:
Copied!
integral(x*sin(x^2), x), integral(x/(x^2+1), x, 0, 1) # 后两个参数表示积分上下限
integral(x*sin(x^2), x), integral(x/(x^2+1), x, 0, 1) # 后两个参数表示积分上下限
Out[16]:
(-1/2*cos(x^2), 1/2*log(2))
矩阵 ¶
matrix 和 Matrix 相似，但是不同，具体参考 https://devv.ai/search?threadId=e2zhxp04r6yo 。
简而言之，
matrix是一个函数，在创建一个矩阵实例。Matrix是一个类，在定义一个矩阵空间，并可以在该空间中创建矩阵。
但是在简单地使用下，两者并没有什么明显的区别。
创建 ¶
In [17]:
Copied!
# 新建一个3 x 3的矩阵A
from sage.all import *
# 类似于 np
A = matrix(3,3)
"""
[0 0 0]
[0 0 0]
[0 0 0]
"""
A = matrix(3,range(9))
"""
[0 1 2]
[3 4 5]
[6 7 8]
"""
# 或者
A = matrix([[1,2,3],[3,2,1],[1,1,1]])
"""
[1 2 3]
[3 2 1]
[1 1 1]
"""
# 新建一个 3 x 3 的矩阵 A
from sage.all import *
# 类似于 np
A = matrix(3,3)
"""
[0 0 0]
[0 0 0]
[0 0 0]
"""
A = matrix(3,range(9))
"""
[0 1 2]
[3 4 5]
[6 7 8]
"""
# 或者
A = matrix([[1,2,3],[3,2,1],[1,1,1]])
"""
[1 2 3]
[3 2 1]
[1 1 1]
"""
Out[17]:
'\n[1 2 3]\n[3 2 1]\n[1 1 1]\n'
方法 ¶
In [18]:
Copied!
A, A.T, A.transpose(), A.rows(), A.columns()
A, A.T, A.transpose(), A.rows(), A.columns()
Out[18]:
( [1 2 3] [1 3 1] [1 3 1] [3 2 1] [2 2 1] [2 2 1] [1 1 1], [3 1 1], [3 1 1], [(1, 2, 3), (3, 2, 1), (1, 1, 1)], [(1, 3, 1), (2, 2, 1), (3, 1, 1)] )
In [19]:
Copied!
# matrix 第一个参数限制矩阵所在的环
AZ = matrix(ZZ, [[2,0], [0,1]])
AQ = matrix(QQ, [[2,0], [0,1]])
AR = matrix(RR, [[2,0], [0,1]])
# matrix 第一个参数限制矩阵所在的环
AZ = matrix(ZZ, [[2,0], [0,1]])
AQ = matrix(QQ, [[2,0], [0,1]])
AR = matrix(RR, [[2,0], [0,1]])
In [20]:
Copied!
matrix(Zmod(3),A), matrix(GF(2),A)
matrix(Zmod(3),A), matrix(GF(2),A)
Out[20]:
( [1 2 0] [1 0 1] [0 2 1] [1 0 1] [1 1 1], [1 1 1] )
向量 ¶
其实和矩阵基本一致，将 matrix 改为 vector 就差不多了，略过。
线性方程组 ¶
In [21]:
Copied!
A = matrix([[1,2,3],[3,2,1],[1,1,1]])
B = vector([0,-4,-1])
A.solve_right(B)
# 有教程说可以用 A\B 来求解，但是 `the backslash operator has been deprecated See https://github.com/sagemath/sage/issues/36394 for details.`
A = matrix([[1,2,3],[3,2,1],[1,1,1]])
B = vector([0,-4,-1])
A.solve_right(B)
# 有教程说可以用 A\B 来求解，但是 `the backslash operator has been deprecated See https://github.com/sagemath/sage/issues/36394 for details.`
Out[21]:
(-2, 1, 0)
In [22]:
Copied!
# mod(a, b) == a%b
a, b = 34525, 17
mod(a,b), a % b
# mod(a, b) == a%b
a, b = 34525, 17
mod(a,b), a % b
Out[22]:
(15, 15)
In [23]:
Copied!
# g.powermod(n, m) == power_mod(g, n, m) == g^n % m
g, n, m = 32546, 2456, 113
g.powermod(n, m), power_mod(g, n, m), g^n % m
# g.powermod(n, m) == power_mod(g, n, m) == g^n % m
g, n, m = 32546, 2456, 113
g.powermod(n, m), power_mod(g, n, m), g^n % m
Out[23]:
(49, 49, 49)
factorial(n)求 n 的阶乘gcd(a,b)求 a 和 b 的最大公约数xgcd(a,b)求 a 和 b 的最大公约数和扩展欧几里得算法euler_phi(n)欧拉函数 $\varphi(x)=x\prod_{i=1}^n(1-\frac1{p_i})$inverse_mod(a, p)求 a 在模 p 意义下的逆元
In [24]:
Copied!
# quadratic_residues(p) == \exit a => x=a^2(mod p) 二次剩余
quadratic_residues(23) # [0, 1, 2, 3, 4, 6, 8, 9, 12, 13, 16, 18]
# kronecker(a,p) == pow(a, (p-1)//2, p) legendre symbol # 判断是否二次剩余
a , p = 3, 11
kronecker(a, p), pow(a, (p-1)//2, p)
# quadratic_residues(p) == \exit a => x=a^2(mod p) 二次剩余
quadratic_residues(23) # [0, 1, 2, 3, 4, 6, 8, 9, 12, 13, 16, 18]
# kronecker(a,p) == pow(a, (p-1)//2, p) legendre symbol # 判断是否二次剩余
a , p = 3, 11
kronecker(a, p), pow(a, (p-1)//2, p)
Out[24]:
(1, 1)
In [25]:
Copied!
a = 2334
p = 0x10001
print(inverse_mod(a, p))
assert(inverse_mod(a, p) * a % p == 1)
a = 2334
p = 0x10001
print(inverse_mod(a, p))
assert(inverse_mod(a, p) * a % p == 1)
35071
x.nth_root(n)求 x 的 n 次方根nth_root_mod(a, n, p)求 $x^n\equiv a\pmod{p}$ 中的 x ( 对应sympy.nthroot_mod(a, n, p, all_roots=True)，sage 中似乎没有，手动实现 )
In [26]:
Copied!
x = 3**8-1
x.nth_root(3, truncate_mode='ceil'), x.nth_root(3, truncate_mode='floor'), x.nth_root(3, truncate_mode=1), x.nth_root(3, truncate_mode=True)
# ((18, False), (18, False), (18, False), (18, False))
def nth_root_mod(a, n, p):
roots = []
for x in range(p):
if power_mod(x, n, p) == a % p:
roots.append(x)
return roots
nth_root_mod(18, 2, 23)
x = 3**8-1
x.nth_root(3, truncate_mode='ceil'), x.nth_root(3, truncate_mode='floor'), x.nth_root(3, truncate_mode=1), x.nth_root(3, truncate_mode=True)
# ((18, False), (18, False), (18, False), (18, False))
def nth_root_mod(a, n, p):
roots = []
for x in range(p):
if power_mod(x, n, p) == a % p:
roots.append(x)
return roots
nth_root_mod(18, 2, 23)
Out[26]:
[8, 15]
平方和分解 ¶
is_square(n)判断 n 是否为完全平方数two_squares(n)将 n 分解为两个平方和three_squares(n)将 n 分解为三个平方和four_squares(n)将 n 分解为四个平方和
1. 两平方和定理（Sum of Two Squares Theorem）¶
- 定理内容：A positive integer n is equal to the sum of two perfect squares if and only if the prime factorization of n contains no odd exponent on any prime that is congruent to 3 modulo 4.
- 具体来说，如果 n 的质因数分解中包含形如 $p^k$ 的项，其中 p 是一个质数， k 是一个奇数，且 p 模 4 余 3，那么 n 就不能表示为两个完全平方数之和。
- 也可以说，一个正整数 nn 可以表示为两个平方的和 $n=a^2+b^2$ 当且仅当在其素因数分解中，所有形如 $4^k(4m+3)$ 的素数的指数都是偶数。
- 例子：例如，$5=1^2+2^2$ 可以表示为两个平方的和，而 71 不能。
2. 三平方和定理（Legendre's Three Square Theorem）¶
- 定理内容：一个正整数 nn 可以表示为三个平方的和 n=a2+b2+c2n=a2+b2+c2 当且仅当 n 不是形如 $4^k(8m+7)$ 的数。
- 例子：例如，$13=3^2+2^2+0^2$ 可以表示为三个平方的和，而 71 不能。
3. 四平方和定理（Lagrange's Four Square Theorem）¶
- 定理内容：每个正整数 nn 都可以表示为四个平方的和 $n=a^2+b^2+c^2+d^2$
- 例子：这下 71 终于能分解了：$71=2^2+3^2+3^2+7^2$
In [27]:
Copied!
big_num = 173178061442550241596295506150572803829268102881297542445649200353047297914764783385643705889370567071577408829104128703765633248277722687055281420899564198724968491216409225857070531370724352556864154450614891750313803499101686782558259953244119778256806332589612663957000269869144555485216828399422391672121
two_squares(big_num), four_squares(71)
big_num = 173178061442550241596295506150572803829268102881297542445649200353047297914764783385643705889370567071577408829104128703765633248277722687055281420899564198724968491216409225857070531370724352556864154450614891750313803499101686782558259953244119778256806332589612663957000269869144555485216828399422391672121
two_squares(big_num), four_squares(71)
Out[27]:
((2124080185874205807105261594884049506071023711463225511583343733907861918114771262515207146329816119356450361177578637163739714518747645791197729692842795, 12987160767716970127272498846249389966243604043045970877384806399031879779694511338759174534221141662365315103597141408214751747758739021312936535825747936), (2, 3, 3, 7))
素数与分解 ¶
is_prime(n)n in Primes()判断是否为质数random_prime(n)生成一个 n 的随机质数previous_prime(n)小于 n 的最大质数next_prime(n)大于 n 的最小质数prime_pi(n)pari(n).primepi()小于 n 的质数个数，根据素数定理有 $prime\_pi(n) \approx \frac{n}{ln(n)}$factor(n)大数质分解，除 factordb 查表外，这个最快；测试 big_num 用时 53.3sprime_factors(n)同上，去除重复因子，保留为列表divisor(n)n 的因子；测试 big_num 用时 49.7snumber_of_divisors(n)n 的因子个数；测试 big_num 用时 49.0ssigma(n, k)n 的因子的 k 次幂之和
In [28]:
Copied!
big_num = 8699621268124163273600280057569065643071518478496234908779966583664908604557271908267773859706827828901385412151814796018448555312901260592
factor(big_num)
big_num = 8699621268124163273600280057569065643071518478496234908779966583664908604557271908267773859706827828901385412151814796018448555312901260592
factor(big_num)
Out[28]:
2^4 * 3^2 * 31 * 61 * 223 * 4013 * 281317 * 4151351 * 339386329 * 370523737 * 5404604441993 * 26798471753993 * 25866088332911027256931479223 * 64889106213996537255229963986303510188999911
In [ ]:
Copied!
divisors(big_num) # 因数分解
divisors(big_num) # 因数分解
In [30]:
Copied!
big_num = 8699621268124163273600280057569065643071518478496234908779966583664908604557271908267773859706827828901385412151814796018448555312901260592
number_of_divisors(big_num)
big_num = 8699621268124163273600280057569065643071518478496234908779966583664908604557271908267773859706827828901385412151814796018448555312901260592
number_of_divisors(big_num)
Out[30]:
61440
CRT（中国剩余定理）¶
$$ \begin{cases}x\equiv2\pmod3\\x\equiv3\pmod5\\x\equiv2\pmod7 \end{cases} $$
求解 x ；可以使用 crt 函数解决。
In [31]:
Copied!
crt_problem = ([2,3,2],[3,5,7])
# 仅适用模两两互素 https://lazzzaro.github.io/2020/05/10/crypto-crypto%E5%B8%B8%E7%94%A8%E5%B7%A5%E5%85%B7/#Sage
def chinese_remainder(problem: tuple) -> int:
remainders = problem[0]
modulus = problem[1]
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(crt_problem), crt(crt_problem[0], crt_problem[1])
crt_problem = ([2,3,2],[3,5,7])
# 仅适用模两两互素 https://lazzzaro.github.io/2020/05/10/crypto-crypto%E5%B8%B8%E7%94%A8%E5%B7%A5%E5%85%B7/#Sage
def chinese_remainder(problem: tuple) -> int:
remainders = problem[0]
modulus = problem[1]
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(crt_problem), crt(crt_problem[0], crt_problem[1])
Out[31]:
(23, 23)