使用外部库和接口¶
为Sage编写代码时，使用Python作为基本结构和接口。为了提高速度、效率或方便起见，您可以使用以下任何语言实现部分代码： Cython ，C/C++，FORTRAN 95，GAP，通用LISP，单数，和PARI/GP。还可以使用所有包含SCAGE的C/C++库 [SageComponents]. 如果你对你的代码依赖于可选的Sage包，你可以使用Octave，甚至Magma，Mathematica或Maple。
在本章中，我们将讨论Sage和 PARI ， GAP 和 单数 .
并行库接口¶
这里是一个逐步的指南，以添加新的巴黎函数到Sage。我们以矩阵的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) 召唤造就了一位新的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(flag) 方法到 matrix_integer 我们称之为 matfrobenius() 方法对与矩阵关联的PARI对象执行一些健全性检查。然后我们将输出从PARI转换为Sage对象：
def frobenius(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 matfrobenius()
EXAMPLES::
sage: A = MatrixSpace(ZZ, 3)(range(9))
sage: A.frobenius(0)
[ 0 0 0]
[ 1 0 18]
[ 0 1 12]
sage: A.frobenius(1)
[x^3 - 12*x^2 - 18*x]
sage: A.frobenius(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矩阵做一个包装器。卡坦矩阵 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矩阵——请参阅目录中的代码） 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))
定义“容易”和“难”是主观的，但这里有一个定义。如果Python或Sage中已经有一个与您试图包装的GAP函数的输出数据类型对应的类，那么包装一个GAP函数是“容易的”。例如，包装任何GUAVA（GAP的纠错码包）函数都是“容易的”，因为纠错码是有限域上的向量空间，而GUAVA函数返回以下数据类型之一：
有限域上的向量，
有限域上的多项式，
有限域上的矩阵，
置换群或其元素，
整数。
Sage已经为每一个都设置了课程。
“练习”是一个很难的例子！这里有一些想法。
写一个盖普的包装
FreeLieAlgebra函数（或者，更广泛地说，所有有限表示的李代数函数在GAP中）。这需要创建新的Python对象。写一个盖普的包装
FreeGroup函数（或者，更广泛地说，所有有限呈现的群函数在GAP中）。这需要编写一些新的Python对象。为GAP的字符表编写一个包装器。虽然这可以在不创建新的Python对象的情况下完成，但是为了最大限度地利用这些表，最好使用新的Python对象。
LibGAP¶
通过接口使用其他程序的缺点是在发送输入和接收结果时有一个不可避免的延迟（大约10毫秒）。如果你必须在一个紧循环中调用函数，这可能会非常慢。调用共享库的延迟要低得多，而且还可以避免将所有内容转换为字符串。这就是Sage包含GAP内核的共享库版本的原因 libgap 在萨奇。中第一个示例的libgap模拟 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函数在概念上没有太大区别。与GAP一样，这可以从简单到困难不等，这取决于Sage中已有多少奇异函数输出的数据结构。
首先，一些术语。对我们来说 曲线 X 在有限域上 F 是一个形式方程 f(x,y) = 0 在哪里 f in F[x,y] 是多项式。它可能是单数，也可能不是单数。A 学位授予地 d 是伽罗瓦轨道 d 点 X(E) 在哪里 E/F 有学位 d . 例如，一个学位的地方 1 也是一个获得学位的地方 3 不过是一个学位的地方 2 不是因为没有学位 3 扩展 F 包含学位 2 扩展。学位授予地 1 也叫 F -理性点。
作为Sage/Singular接口的示例，我们将解释如何包装Singular的 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接口中对Singular进行相同计算的另一种方法：
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)
[1]:
[1]:
0
[2]:
1
[3]:
0
[2]:
[1]:
-2
[2]:
-1
[3]:
1
...
通过查看输出，请注意，我们的包装器函数将需要解析 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 .
还有一个示例（除了docstring中的示例）：
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))
单数：另一种方法¶
对于Singular，还有一个更像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)]
接下来，我们实现通用函数（为了简洁起见，我们省略了docstring，与上面相同）。注意 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/Singular接口中实现，而上一节中的代码只使用了该接口中最基本的部分。
创建新的伪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 也。在这之后是一个docstring，我们在这里省略了它（有关详细信息，请参阅文件）。接下来是：
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 要设置倍频程接口：
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:
An 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/ 其他软件包的接口示例。
- SageComponents
看到了吗网址：www.sagemath.org/一份名单