使用外部库和接口¶
在为Sage编写代码时，请使用Python作为基本结构和接口。为提高速度、效率或方便起见，您可以使用以下任何一种语言实现部分代码： Cython 、C/C++、Fortran 95、GAP、Common Lisp、Single和Pari/GP。您还可以使用Sage附带的所有C/C++库，请参见 spkg 。如果您对依赖于外部程序的代码没有意见，您可以使用Octave，甚至Magma、MATHEMICATA或Maple。
在本章中，我们将讨论Sage和 PARI ， GAP 和 单数 。
Pari C库接口¶
下面是向Sage添加新的Pari函数的分步指南。我们以矩阵的Frobenius形式为例。使用PARI库实现了整数上矩阵的一些繁重提升。为了在Pari中计算Frobenius形式， matfrobenius 函数被使用。
有两种方式可以与Sage的Pari库交互。GP接口使用GP解释器。Pari接口使用对Pari C函数的直接调用-这是首选方式，因为它的速度要快得多。因此，本节重点介绍如何使用Pari。
我们将把一个新方法添加到 gen 班级。这是所有Pari库对象的抽象表示。这意味着，一旦我们将一个方法添加到这个类中，每个Pari对象，无论它是数字、多项式还是矩阵，都将拥有我们的新方法。所以你可以这样做 pari(1).matfrobenius() ，但由于Pari想要申请 matfrobenius 对于矩阵，而不是数字，您将收到一个 PariError 在这种情况下。
这个 gen 类在中定义 SAGE_ROOT/src/sage/libs/cypari2/gen.pyx ，这就是我们添加方法的地方 matfrobenius ：
def matfrobenius(self, flag=0):
r"""
M.matfrobenius(flag=0): Return the Frobenius form of the square
matrix M. If flag is 1, return only the elementary divisors (a list
of polynomials). If flag is 2, return a two-components vector [F,B]
where F is the Frobenius form and B is the basis change so that
`M=B^{-1} F B`.
EXAMPLES::
sage: a = pari('[1,2;3,4]')
sage: a.matfrobenius()
[0, 2; 1, 5]
sage: a.matfrobenius(flag=1)
[x^2 - 5*x - 2]
sage: a.matfrobenius(2)
[[0, 2; 1, 5], [1, -1/3; 0, 1/3]]
"""
sig_on()
return self.new_gen(matfrobenius(self.g, flag, 0))
请注意 sig_on() statement 。
这个 matfrobenius 调用只是对Pari C库函数的调用 matfrobenius 具有适当的参数。
这个 self.new_gen(GEN x) Call构建了一个新的Sage gen 来自给定参数的对象 GEN 巴黎在哪里？ GEN 被存储为 .g 属性。除此之外， self.new_gen() 叫停了 sig_off() 宏，还清除了pari堆栈，因此在 return 如上所示的语句。所以在那之后 self.new_gen() ，全部平价 GEN 未转换为Sage的 gen 已经没有了。还有 self.new_gen_noclear(GEN x) 它的作用与 self.new_gen(GEN x) 但事实并非如此 not 打电话 sig_off() 也不会清理平价协议。
有关调用哪个函数以及如何调用它的信息可以从PARI用户手册中获取(注意：SAGE包含PARI的开发版本，因此请查看该用户手册的版本)。寻找 matfrobenius 您可以找到：
库语法为
GEN matfrobenius(GEN M, long flag, long v = -1)，在哪里v是一个可变的数字。
如果您熟悉GP，请注意Pari C函数的名称可能与相应的GP函数不同(例如，请参阅 mathnf )，所以一定要检查手册。
我们还可以添加一个 frobenius_form(flag) 方法添加到 matrix_integer 类，在该类中我们调用 matfrobenius() 在执行一些健全性检查之后，对与矩阵相关联的pari对象使用。然后，我们将输出从Pari转换为Sage对象：
def frobenius_form(self, flag=0, var='x'):
"""
Return the Frobenius form (rational canonical form) of this matrix.
INPUT:
- ``flag`` -- 0 (default), 1 or 2 as follows:
- ``0`` -- (default) return the Frobenius form of this
matrix.
- ``1`` -- return only the elementary divisor
polynomials, as polynomials in var.
- ``2`` -- return a two-components vector [F,B] where F
is the Frobenius form and B is the basis change so that
`M=B^{-1}FB`.
- ``var`` -- a string (default: 'x')
ALGORITHM: uses PARI's :pari:`matfrobenius`
EXAMPLES::
sage: A = MatrixSpace(ZZ, 3)(range(9))
sage: A.frobenius_form(0)
[ 0 0 0]
[ 1 0 18]
[ 0 1 12]
sage: A.frobenius_form(1)
[x^3 - 12*x^2 - 18*x]
sage: A.frobenius_form(1, var='y')
[y^3 - 12*y^2 - 18*y]
"""
if not self.is_square():
raise ArithmeticError("frobenius matrix of non-square matrix not defined.")
v = self.__pari__().matfrobenius(flag)
if flag == 0:
return self.matrix_space()(v.python())
elif flag == 1:
r = PolynomialRing(self.base_ring(), names=var)
retr = []
for f in v:
retr.append(eval(str(f).replace("^","**"), {'x':r.gen()}, r.gens_dict()))
return retr
elif flag == 2:
F = matrix_space.MatrixSpace(QQ, self.nrows())(v[0].python())
B = matrix_space.MatrixSpace(QQ, self.nrows())(v[1].python())
return F, B
GAP¶
在Sage中包装Gap函数就是用Python语言编写一个程序，该程序使用pExpect接口将各种命令输送到Gap并将输入回读到Sage中。这有时很容易，有时很难。
例如，假设我们想要制作一个用于计算简单李代数的Cartan矩阵的包装器。的Cartan矩阵 G_2 使用以下命令在GAP中可用：
gap> L:= SimpleLieAlgebra( "G", 2, Rationals );
<Lie algebra of dimension 14 over Rationals>
gap> R:= RootSystem( L );
<root system of rank 2>
gap> CartanMatrix( R );
在Sage中，用户可以通过键入以下命令来访问这些命令：
sage: L = gap.SimpleLieAlgebra('"G"', 2, 'Rationals'); L
Algebra( Rationals, [ v.1, v.2, v.3, v.4, v.5, v.6, v.7, v.8, v.9, v.10,
v.11, v.12, v.13, v.14 ] )
sage: R = L.RootSystem(); R
<root system of rank 2>
sage: R.CartanMatrix()
[ [ 2, -1 ], [ -3, 2 ] ]
请注意 '"G"' 它在GAP中作为字符串进行计算 "G" 。
本节的目的是使用此示例来说明如何编写一个Python/Sage程序，其输入是，比如说， ('G',2) 并且其输出是上面的矩阵(但作为Sage Matrix-参见目录中的代码 SAGE_ROOT/src/sage/matrix/ 和Sage参考手册的相应部分)。
首先，必须将输入转换为由合法GAP命令组成的字符串。然后，GAP输出(也是一个字符串)必须被解析并在可能的情况下转换为相应的Sage/Python对象。
def cartan_matrix(type, rank):
"""
Return the Cartan matrix of given Chevalley type and rank.
INPUT:
- type -- a Chevalley letter name, as a string, for
a family type of simple Lie algebras
- rank -- an integer (legal for that type).
EXAMPLES::
sage: cartan_matrix("A",5)
[ 2 -1 0 0 0]
[-1 2 -1 0 0]
[ 0 -1 2 -1 0]
[ 0 0 -1 2 -1]
[ 0 0 0 -1 2]
sage: cartan_matrix("G",2)
[ 2 -1]
[-3 2]
"""
L = gap.SimpleLieAlgebra('"%s"' % type, rank, 'Rationals')
R = L.RootSystem()
sM = R.CartanMatrix()
ans = eval(str(sM))
MS = MatrixSpace(QQ, rank)
return MS(ans)
输出 ans 是一个Python列表。最后两行将该列表转换为Sage类实例 Matrix 。
或者，可以将上述函数的第一行替换为：
L = gap.new('SimpleLieAlgebra("%s", %s, Rationals);'%(type, rank))
“容易”和“难”的定义是主观的，但这里有一个定义。如果您试图包装的GAP函数的输出数据类型在Python或Sage中已经有对应的类，那么包装GAP函数是“容易的”。例如，包装任何芭乐(GAP的纠错码包)函数很容易，因为纠错码是有限域上的向量空间，并且芭乐函数返回以下数据类型之一：
有限域上的向量，
有限域上的多项式，
有限域上的矩阵，
置换群或其元素，
整数。
Sage已经为这些课程中的每一个开设了课程。
作为练习，我们留下了一个“难”的例子！以下是一些想法。
为Gap编写包装器
FreeLieAlgebra函数(或者，更一般地，间隙中的所有有限表示的李代数函数)。这将需要创建新的Python对象。为Gap编写包装器
FreeGroup函数(或者，更一般地，指GAP中所有有限呈现的群函数)。这将需要编写一些新的Python对象。为GAP的字符表编写一个包装器。虽然这可以在不创建新的Python对象的情况下完成，但为了最大限度地利用这些表，最好使用新的Python对象。
LibGAP¶
通过接口使用其他程序的缺点是，在发送输入和接收结果时，存在一定的不可避免的延迟(约为10ms)。如果您必须在紧凑的循环中调用函数，这可能会慢得令人无法接受。调用共享库的延迟要低得多，而且还避免了将所有内容转换为中间的字符串。这就是为什么Sage包含GAP内核的共享库版本，可以通过以下方式获得 libgap 在塞奇。中第一个示例的libap模拟 GAP 是：：
sage: SimpleLieAlgebra = libgap.function_factory('SimpleLieAlgebra')
sage: L = SimpleLieAlgebra('G', 2, QQ)
sage: R = L.RootSystem(); R
<root system of rank 2>
sage: R.CartanMatrix() # output is a GAP matrix
[ [ 2, -1 ], [ -3, 2 ] ]
sage: matrix(R.CartanMatrix()) # convert to Sage matrix
[ 2 -1]
[-3 2]
单数¶
使用Sage中的奇异函数与使用Sage中的间隙函数在概念上没有太大区别。与GAP一样，这可能是容易的，也可能是困难的，这取决于单一函数的输出的数据结构在Sage中已经存在了多少。
首先，是一些术语。对我们来说，一个 curve X 在有限域上 F 是以下形式的方程 f(x,y) = 0 ，在哪里 f in F[x,y] 是一个多项式。它可能是也可能不是单数。一个 place of degree d 是伽罗瓦轨道 d 积分输入 X(E) ，在哪里 E/F 是有程度的 d 。例如，一个学位的地方 1 也是一个有学位的地方 3 ，而是一个有学位的地方 2 不是因为没有学位 3 展期 F 包含一个学位 2 分机。学位课程 1 也被称为 F --理性点。
作为Sage/Single接口的一个示例，我们将解释如何包装Single的 NSplaces ，它计算有限域上的曲线上的位置。(该命令 closed_points 在某些情况下也会这样做。)这很“容易”，因为在Sage中不需要新的Python类来实现这一点。
以下是如何使用单数形式使用此命令的示例：
A Computer Algebra System for Polynomial Computations / version 3-0-0
0<
by: G.-M. Greuel, G. Pfister, H. Schoenemann \ May 2005
FB Mathematik der Universitaet, D-67653 Kaiserslautern \
> LIB "brnoeth.lib";
[...]
> ring s=5,(x,y),lp;
> poly f=y^2-x^9-x;
> list X1=Adj_div(f);
Computing affine singular points ...
Computing all points at infinity ...
Computing affine singular places ...
Computing singular places at infinity ...
Computing non-singular places at infinity ...
Adjunction divisor computed successfully
The genus of the curve is 4
> list X2=NSplaces(1,X1);
Computing non-singular affine places of degree 1 ...
> list X3=extcurve(1,X2);
Total number of rational places : 6
> def R=X3[1][5];
> setring R;
> POINTS;
[1]:
[1]:
0
[2]:
1
[3]:
0
[2]:
[1]:
-2
[2]:
1
[3]:
1
[3]:
[1]:
-2
[2]:
1
[3]:
1
[4]:
[1]:
-2
[2]:
-1
[3]:
1
[5]:
[1]:
2
[2]:
-2
[3]:
1
[6]:
[1]:
0
[2]:
0
[3]:
1
以下是在Sage界面中对Single执行相同计算的另一种方法：
sage: singular.LIB("brnoeth.lib")
sage: singular.ring(5,'(x,y)','lp')
polynomial ring, over a field, global ordering
// coefficients: ZZ/5
// number of vars : 2
// block 1 : ordering lp
// : names x y
// block 2 : ordering C
sage: f = singular('y^2-x^9-x')
sage: print(singular.eval("list X1=Adj_div(%s);"%f.name()))
Computing affine singular points ...
Computing all points at infinity ...
Computing affine singular places ...
Computing singular places at infinity ...
Computing non-singular places at infinity ...
Adjunction divisor computed successfully
<BLANKLINE>
The genus of the curve is 4
sage: print(singular.eval("list X2=NSplaces(1,X1);"))
Computing non-singular affine places of degree 1 ...
sage: print(singular.eval("list X3=extcurve(1,X2);"))
<BLANKLINE>
Total number of rational places : 6
<BLANKLINE>
sage: singular.eval("def R=X3[1][5];")
''
sage: singular.eval("setring R;")
''
sage: L = singular.eval("POINTS;")
sage: print(L) # random
[1]:
[1]:
0
[2]:
1
[3]:
0
...
从输出中可以看到，我们的包装器函数需要解析由 L 因此，让我们编写一个单独的函数来实现这一点。这需要找出如何确定点的坐标在字符串中的位置 L 。要做到这一点，Python有一些非常有用的字符串操作命令。
def points_parser(string_points, F):
"""
This function will parse a string of points
of X over a finite field F returned by Singular's NSplaces
command into a Python list of points with entries from F.
EXAMPLES::
sage: F = GF(5)
sage: points_parser(L,F)
((0, 1, 0), (3, 4, 1), (0, 0, 1), (2, 3, 1), (3, 1, 1), (2, 2, 1))
"""
Pts = []
n = len(L)
# start block to compute a pt
L1 = L
while len(L1) > 32:
idx =L1.index(" ")
pt = []
# start block1 for compute pt
idx = L1.index(" ")
idx2 = L1[idx:].index("\n")
L2 = L1[idx:idx+idx2]
pt.append(F(eval(L2)))
# end block1 to compute pt
L1 = L1[idx+8:] # repeat block 2 more times
# start block2 for compute pt
idx = L1.index(" ")
idx2 = L1[idx:].index("\n")
L2 = L1[idx:idx+idx2]
pt.append(F(eval(L2)))
# end block2 to compute pt
L1=L1[idx+8:] # repeat block 1 more time
# start block3 for compute pt
idx=L1.index(" ")
if "\n" in L1[idx:]:
idx2 = L1[idx:].index("\n")
else:
idx2 = len(L1[idx:])
L2 = L1[idx:idx+idx2]
pt.append(F(eval(L2)))
# end block3 to compute pt
# end block to compute a pt
Pts.append(tuple(pt)) # repeat until no more pts
L1 = L1[idx+8:] # repeat block 2 more times
return tuple(Pts)
Now it is an easy matter to put these ingredients together into a Sage function which takes as input a triple (f,F,d): a polynomial f in F[x,y] defining X: f(x,y)=0 (note that the variables x,y must be used), a finite field F of prime order, and the degree d. The output is the number of places in X of degree d=1 over F. At the moment, there is no "translation" between elements of GF(p^d) in Singular and Sage unless d=1. So, for this reason, we restrict ourselves to points of degree one.
def places_on_curve(f, F):
"""
INPUT:
- f -- element of F[x,y], defining X: f(x,y)=0
- F -- a finite field of *prime order*
OUTPUT:
integer -- the number of places in X of degree d=1 over F
EXAMPLES::
sage: F = GF(5)
sage: R = PolynomialRing(F,2,names=["x","y"])
sage: x,y = R.gens()
sage: f = y^2-x^9-x
sage: places_on_curve(f,F)
((0, 1, 0), (3, 4, 1), (0, 0, 1), (2, 3, 1), (3, 1, 1), (2, 2, 1))
"""
d = 1
p = F.characteristic()
singular.eval('LIB "brnoeth.lib";')
singular.eval("ring s="+str(p)+",(x,y),lp;")
singular.eval("poly f="+str(f))
singular.eval("list X1=Adj_div(f);")
singular.eval("list X2=NSplaces("+str(d)+",X1);")
singular.eval("list X3=extcurve("+str(d)+",X2);")
singular.eval("def R=X3[1][5];")
singular.eval("setring R;")
L = singular.eval("POINTS;")
return points_parser(L,F)
请注意，此Sage函数返回的排序与奇异变量中的排序完全相同 POINTS 。
再举一个例子(除了文档字符串中的那个)：
sage: F = GF(2)
sage: R = MPolynomialRing(F,2,names = ["x","y"])
sage: x,y = R.gens()
sage: f = x^3*y+y^3+x
sage: places_on_curve(f,F)
((0, 1, 0), (1, 0, 0), (0, 0, 1))
单数：另一种方法¶
对于Single，还有一个更类似于Python的接口。使用这种方法，代码变得简单得多，如下所示。首先，我们演示在特定情况下计算曲线上的位置：
sage: singular.lib('brnoeth.lib')
sage: R = singular.ring(5, '(x,y)', 'lp')
sage: f = singular.new('y^2 - x^9 - x')
sage: X1 = f.Adj_div()
sage: X2 = singular.NSplaces(1, X1)
sage: X3 = singular.extcurve(1, X2)
sage: R = X3[1][5]
sage: singular.set_ring(R)
sage: L = singular.new('POINTS')
请注意，L的这些元素是以单数5为模定义的，它们的比较与您从其印刷表示中预期的不同：
sage: sorted([(L[i][1], L[i][2], L[i][3]) for i in range(1,7)])
[(0, 0, 1), (0, 1, 0), (2, 2, 1), (2, -2, 1), (-2, 1, 1), (-2, -1, 1)]
接下来，我们实现通用函数(为简洁起见，我们省略了文档字符串，与上面相同)。请注意， point_parser 功能不是必需的：
def places_on_curve(f, F):
p = F.characteristic()
if F.degree() > 1:
raise NotImplementedError
singular.lib('brnoeth.lib')
R = singular.ring(5, '(x,y)', 'lp')
f = singular.new('y^2 - x^9 - x')
X1 = f.Adj_div()
X2 = singular.NSplaces(1, X1)
X3 = singular.extcurve(1, X2)
R = X3[1][5]
singular.setring(R)
L = singular.new('POINTS')
return [(int(L[i][1]), int(L[i][2]), int(L[i][3])) \
for i in range(1,int(L.size())+1)]
这段代码更短、更美观、更具可读性。然而，它取决于某些功能，例如 singular.setring 已在Sage/Single接口中实现，而上一节中的代码仅使用该接口的最小部分。
创建新的伪TTY接口¶
您可以创建Sage伪tty接口，这些接口允许Sage与几乎任何命令行程序一起工作，并且不需要对该程序进行任何修改或扩展。它们还出人意料地快速和灵活(考虑到它们的工作方式！)，因为所有I/O都是缓冲的，而且Sage和命令行程序之间的交互可以是非阻塞的(异步)。伪tty Sage接口是异步的，因为它派生自Sage类 Expect ，它处理Sage和外部进程之间的通信。
例如，以下是文件的一部分 SAGE_ROOT/src/sage/interfaces/octave.py ，它定义了Sage和Octave之间的接口，Octave是一个用于进行数值计算的开源程序，等等：
import os
from expect import Expect, ExpectElement
class Octave(Expect):
...
前两行导入库 os ，它包含操作系统例程和类 Expect ，这是接口的基本类。第三行定义类 Octave ；它源自 Expect 也是。在这之后是一个文档字符串，我们在这里省略它(有关详细信息，请参阅文件)。接下来是：
def __init__(self, script_subdirectory="", logfile=None,
server=None, server_tmpdir=None):
Expect.__init__(self,
name = 'octave',
prompt = '>',
command = "octave --no-line-editing --silent",
server = server,
server_tmpdir = server_tmpdir,
script_subdirectory = script_subdirectory,
restart_on_ctrlc = False,
verbose_start = False,
logfile = logfile,
eval_using_file_cutoff=100)
这使用了类 Expect 要设置Octave接口，请执行以下操作：
def set(self, var, value):
"""
Set the variable var to the given value.
"""
cmd = '%s=%s;' % (var,value)
out = self.eval(cmd)
if out.find("error") != -1:
raise TypeError("Error executing code in Octave\nCODE:\n\t%s\nOctave ERROR:\n\t%s"%(cmd, out))
def get(self, var):
"""
Get the value of the variable var.
"""
s = self.eval('%s' % var)
i = s.find('=')
return s[i+1:]
def console(self):
octave_console()
用户可以通过它们来输入 octave.set('x', 3) ，之后， octave.get('x') 退货 ' 3' 。正在运行 octave.console() 将用户放入Octave交互Shell中：
def solve_linear_system(self, A, b):
"""
Use octave to compute a solution x to A*x = b, as a list.
INPUT:
- A -- mxn matrix A with entries in QQ or RR
- b -- m-vector b entries in QQ or RR (resp)
OUTPUT:
A list x (if it exists) which solves M*x = b
EXAMPLES::
sage: M33 = MatrixSpace(QQ,3,3)
sage: A = M33([1,2,3,4,5,6,7,8,0])
sage: V3 = VectorSpace(QQ,3)
sage: b = V3([1,2,3])
sage: octave.solve_linear_system(A,b) # optional - octave
[-0.333333, 0.666667, 0]
AUTHOR: David Joyner and William Stein
"""
m = A.nrows()
n = A.ncols()
if m != len(b):
raise ValueError("dimensions of A and b must be compatible")
from sage.matrix.all import MatrixSpace
from sage.rings.all import QQ
MS = MatrixSpace(QQ, m, 1)
b = MS(list(b)) # converted b to a "column vector"
sA = self.sage2octave_matrix_string(A)
sb = self.sage2octave_matrix_string(b)
self.eval("a = " + sA )
self.eval("b = " + sb )
soln = octave.eval("c = a \\ b")
soln = soln.replace("\n\n ", "[")
soln = soln.replace("\n\n", "]")
soln = soln.replace("\n", ",")
sol = soln[3:]
return eval(sol)
此代码定义了该方法 solve_linear_system ，它的工作方式如文档所示。
这些只是节选自 octave.py ；有关更多定义和示例，请查看该文件。查看目录中的其他文件 SAGE_ROOT/src/sage/interfaces/ 获取与其他软件包的接口的示例。