孤独的猫

  博客园 :: 首页 :: 博问 :: 闪存 :: 新随笔 :: 联系 :: 订阅 订阅 :: 管理 ::

载入和附加Sage文件


下面我们说明如何将写在单独文件中的程序载入到Sage中。新建一个名为example.sage的文件,包含下列内容:

print "Hello World"
print 2^3您可以使用load命令读入和执行example.sage文件。

sage: load "example.sage"
Hello World
8您也可以使用attach命令将Sage文件附加到一个运行的会话中:

sage: attach "example.sage"
Hello World
8现在如果您修改example.sage,并在Sage中输入一个空行(即按return1)),那么example.sage的内容将自动重新载入到Sage中。

特别的,attach会在文件改动时自动重载文件,这在调试代码时很方便,而load只能载入一次。

当Sage载入example.sage时,它将其转为Python,然后由Python解释器执行。这种转换是很小的;它主要包括将字面上的整数放入Integer(),浮点数放入RealNumber(),用**代替^,以及将例如R.2这样的替换成R.gen(2)。转换后的example.sage版本和example.sage包含在同一个目录中,名为example.sage.py。这个文件包含以下代码:

print "Hello World"
print Integer(2)**Integer(3)字面上的整数被转换,并且^被**代替。(在Python中,^表示异或并且**表示求幂。)

这种预解析由sage/misc/interpreter.py实现。

您可以把多行缩进代码贴进Sage,只要有新行来构成新代码块(这在文件中并不是必须的)。然而,输入这种代码的最好方法是将它保存到一个文件中,并用attach,如上所述。

创建编译型代码
速度在数学计算中是至关重要的。尽管Python是一个方便、非常高级的语言,某些计算如果使用编译语言的静态类型实现,可以比Python快上几个数量级。如果完全用Python写成,Sage的某些方面可能会非常缓慢。为解决这个问题,Sage支持一个Python的编译型“版本”,称为Cython([Cyt] 和[Pyt] )。Cython即像Python又像C。多数Python构造,包括列表解析、条件表达式、类似+=的代码都是允许的;您也可以将在其他Python模块中写的代码导入。此外,您可以声明任意C变量,并直接调用C语言库。最后的代码转换成C代码并用C编译器进行编译。

要生成您自己的编译型Sage代码,给文件一个.spyx扩展(而不是.sage)。如果您使用命令行界面工作,您可以完全像对解释型代码那样附加和载入编译型代码(目前还不支持在notebook界面附加和载入Cython代码)。实际的编译在“幕后”完成,不需要您明确做什么。参见$SAGE_ROOT/examples/programming/sagex/factorial.spyx,来得到一个直接使用GMP C语言库的编译型阶乘函数的例子。要亲自尝试,切换到$SAGE_ROOT/examples/programming/sagex/,然后执行下列操作:

sage: load "factorial.spyx"
***************************************************
                Recompiling factorial.spyx
***************************************************
sage: factorial(50)
30414093201713378043612608166064768844377641568960512000000000000L
sage: time n = factorial(10000)
CPU times: user 0.03 s, sys: 0.00 s, total: 0.03 s
Wall time: 0.03这里,末尾的L表示一个Python长整型(参见预解析器:Sage和Python的不同)。

注意,如果您退出并重启Sage,Sage将重新编译factorial.spyx。编译好的公共对象库保存在$HOME/.sage/temp/hostname/pid/spyx。这些文件在您退出Sage时被删除。

Sage预解析不会应用到spyx文件,例如,在spyx文件中1/3将返回0,而不是有理数1/3。如果foo是Sage语言库中的函数,要在spyx文件中使用它则导入sage.all并使用sage.all.foo。

import sage.all
def foo(n):
    return sage.all.factorial(n)访问单独文件中的C函数
访问在单独的*.c文件中定义的C函数也很容易。这里是一个例子。在同一个目录中创建文件test.c和test.spyx,包含以下内容:

单纯的C代码:test.c

int add_one(int n) {
  return n + 1;
}Cython代码:test.spyx2):

cdef extern from "test.c":
    int add_one(int n)
 
def test(n):
    return add_one(n)然后进行下面操作:

sage: attach "test.spyx"
Compiling (...)/test.spyx...
sage: test(10)
11如果编译Cython文件中产生的C代码需要额外的库foo,将clib foo添加到Cython源代码中。类似的,用声明cfile bar可以将额外的C文件bar包含到编译中。

独立的Python/Sage脚本
下面的独立Sage脚本对整数和多项式等进行因式分解3):

#!/usr/bin/env sage -python
 
import sys
from sage.all import *
 
if len(sys.argv) != 2:
    print "Usage: %s <n>"%sys.argv[0]
    print "Outputs the prime factorization of n."
    sys.exit(1)
 
print factor(sage_eval(sys.argv[1]))要使用这个脚本,您的SAGE_ROOT必须在PATH路径中。假设上面的脚本名为factor,这里是一个示例用法4):

bash $ ./factor 2006
2 * 17 * 59
bash $ ./factor "32*x^5-1"
(2*x - 1) * (16*x^4 + 8*x^3 + 4*x^2 + 2*x + 1)数据类型
Sage中的每个对象都有良好定义的类型。Python有广泛的基本内置类型,并且Sage库添加了更多。Python的一些内置类型包括字符串、列表、元组和浮点数,说明如下:

sage: s = "sage"; type(s)
<type 'str'>
sage: s = 'sage'; type(s)      # 您可以使用单引号或者双引号
<type 'str'>
sage: s = [1,2,3,4]; type(s)
<type 'list'>
sage: s = (1,2,3,4); type(s)
<type 'tuple'>
sage: s = int(2006); type(s)
<type 'int'>
sage: s = float(2006); type(s)
<type 'float'>对于这些,Sage增加了很多其他类型,例如向量空间:

sage: V = VectorSpace(QQ, 1000000); V
Vector space of dimension 1000000 over Rational Field
sage: type(V)
<class 'sage.modules.free_module.FreeModule_ambient_field_with_category'>在V上只能调用特定函数。在其他数学软件系统中,这些会通过“函数”符号foo(V,…)调用。在Sage中,特定函数被附加在V的类型(或类)上,并且使用类似JAVA或C++中的面向对象的语法,例如V.foo(…)。这有助于保护全局命名空间不被无数函数弄乱,并意味着很多不同作用的函数都可以称作foo,不需要用变量类型检查(或者case语句)来决定调用哪个。此外,如果您重用了一个函数的名称,那个函数仍然有效(例如,如果您将什么命名为zeta,然后想要计算黎曼-Zeta函数在0.5的值,您仍然可以输入s=.5; s.zeta())。

sage: zeta = -1
sage: s=.5; s.zeta()
-1.46035450880959在一些非常普遍的情况下,为了方便,同时由于使用面向对象符号的数学表达式可能看着令人困惑,通常的函数符号也是支持的。这里有一些例子。

sage: n = 2; n.sqrt()
sqrt(2)
sage: sqrt(2)
sqrt(2)
sage: V = VectorSpace(QQ,2)
sage: V.basis()
    [
    (1, 0),
    (0, 1)
    ]
sage: basis(V)
    [
    (1, 0),
    (0, 1)
    ]
sage: M = MatrixSpace(GF(7), 2); M
Full MatrixSpace of 2 by 2 dense matrices over Finite Field of size 7
sage: A = M([1,2,3,4]); A
[1 2]
[3 4]
sage: A.charpoly('x')
x^2 + 2*x + 5
sage: charpoly(A, 'x')
x^2 + 2*x + 5要列出A的所有成员函数,使用tab补全。只要输入A,然后按您键盘上的[tab]键,如反向搜索和TAB补全中所述。

列表,元组和序列
列表数据类型存储任意类型的元素。类似C、C++等(但和很多标准计算机代数系统不同),列表元素的索引从0开始:

sage: v = [2, 3, 5, 'x', SymmetricGroup(3)]; v
[2, 3, 5, 'x', Symmetric group of order 3! as a permutation group]
sage: type(v)
<type 'list'>
sage: v[0]
2
sage: v[2]
5(索引列表时,如果索引不是Python整数类型也可以!)Sage整型(或者有理数,或带有__index__方法的任何东西)也可以很好工作。

sage: v = [1,2,3]
sage: v[2]
3
sage: n = 2      # SAGE整型
sage: v[n]       # 非常好!
3
sage: v[int(n)]  # 也可以。
3range函数创建一个Python整数(不是Sage整数)的列表:

sage: range(1, 15)
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]这在使用列表解析5)构造列表时很有用:

sage: L = [factor(n) for n in range(1, 15)]
sage: print L
[1, 2, 3, 2^2, 5, 2 * 3, 7, 2^3, 3^2, 2 * 5, 11, 2^2 * 3, 13, 2 * 7]
sage: L[12]
13
sage: type(L[12])
<class 'sage.structure.factorization.Factorization'>
sage: [factor(n) for n in range(1, 15) if is_odd(n)]
[1, 3, 5, 7, 3^2, 11, 13]更多如何使用列表解析创建列表的信息,参见[PyT] 。

列表切片是一个精彩的特性。如果L是一个列表,那么L[m:n]返回L从第m个元素到第n-1个元素的子列表,如下所述。

sage: L = [factor(n) for n in range(1, 20)]
sage: L[4:9]
[5, 2 * 3, 7, 2^3, 3^2]
sage: print L[:4]
[1, 2, 3, 2^2]
sage: L[14:4]
[]
sage: L[14:]
[3 * 5, 2^4, 17, 2 * 3^2, 19]元组与列表类似,除了它们是不可变的,这意味着一旦它们被创建则不能改变。

sage: v = (1,2,3,4); v
(1, 2, 3, 4)
sage: type(v)
<type 'tuple'>
sage: v[1] = 5
...
TypeError: 'tuple' object does not support item assignment序列是第三种面向列表的Sage类型。与列表和序列不同,序列不是内置的Python类型。序列默认是可变的,但用Sequence的类方法set_immutable可以将它设为不可变,如下面例子所示。序列的所有元素都有共同的parent,称之为序列的universe。

sage: v = Sequence([1,2,3,4/5])
sage: v
[1, 2, 3, 4/5]
sage: type(v)
<class 'sage.structure.sequence.Sequence'>
sage: type(v[1])
<type 'sage.rings.rational.Rational'>
sage: v.universe()
Rational Field
sage: v.is_immutable()
False
sage: v.set_immutable()
sage: v[0] = 3
...
ValueError: object is immutable; please change a copy instead.序列从列表中导出,并能在任何可以用列表的地方使用:

sage: v = Sequence([1,2,3,4/5])
sage: isinstance(v, list)
True
sage: list(v)
[1, 2, 3, 4/5]
sage: type(list(v))
<type 'list'>作为另一个列子,向量空间的基是不可变序列,因为您能不改变它们这一点是很重要的。

sage: V = QQ^3; B = V.basis(); B
[
(1, 0, 0),
(0, 1, 0),
(0, 0, 1)
]
sage: type(B)
<class 'sage.structure.sequence.Sequence'>
sage: B[0] = B[1]
...
ValueError: object is immutable; please change a copy instead.
sage: B.universe()
Vector space of dimension 3 over Rational Field字典
字典(有时也称作关联数组6))是从“可哈希7)”对象(例如,字符串、数字和此类元组;参见Python文档http://docs.python.org/tut/node7.htmlhttp://docs.python.org/lib/typesmapping.html获取细节)到任意对象的映射。

sage: d = {1:5, 'sage':17, ZZ:GF(7)}
sage: type(d)
<type 'dict'>
sage: d.keys()
 [1, 'sage', Integer Ring]
sage: d['sage']
17
sage: d[ZZ]
Finite Field of size 7
sage: d[1]
5第三个键说明字典的索引可以是复杂的,例如整数环。

您可以将字典转换为有同样数据的列表:

sage: d.items()
[(1, 5), ('sage', 17), (Integer Ring, Finite Field of size 7)]一个常见用法是遍历字典中的对:

sage: d = {2:4, 3:9, 4:16}
sage: [a*b for a, b in d.iteritems()]
[8, 27, 64]字典是无序的,如最后的输出所示。

集合
Python有内置的集合类型。它提供的主要特性是快速检查某个元素是否在集合中,并有标准的集合论运算。

sage: X = set([1,19,'a']);   Y = set([1,1,1, 2/3])
sage: X
set(['a', 1, 19])
sage: Y
set([1, 2/3])
sage: 'a' in X
True
sage: 'a' in Y
False
sage: X.intersection(Y)
set([1])Sage也有自己的集合类型,通过内置的Python集合类型实现,不过只有很少额外的Sage相关功能。用Set(…)来创建一个Sage集合。例如,

sage: X = Set([1,19,'a']);   Y = Set([1,1,1, 2/3])
sage: X
{'a', 1, 19}
sage: Y
{1, 2/3}
sage: X.intersection(Y)
{1}
sage: print latex(Y)
\left\{1, \frac{2}{3}\right\}
sage: Set(ZZ)
Set of elements of Integer Ring迭代器
迭代器是最近加入到Python中的,对于数学应用尤其有用。这里有几个例子;参见[PyT] 获取更多细节。我们在非负整数的平方上创建一个到10000000的迭代器。

sage: v = (n^2 for n in xrange(10000000))
sage: v.next()
0
sage: v.next()
1
sage: v.next()
4我们在4p+1(p也为素数)形式的素数上创建一个迭代器,并查看开始的一些值。

sage: w = (4*p + 1 for p in Primes() if is_prime(4*p+1))
sage: w         # 在下一行,0xb0853d6c是一个随机的十六进制数
<generator object at 0xb0853d6c>
sage: w.next()
13
sage: w.next()
29
sage: w.next()
53某些环,例如有限域和整数,有迭代器和它们相连:

sage: [x for x in GF(7)]
[0, 1, 2, 3, 4, 5, 6]
sage: W = ((x,y) for x in ZZ for y in ZZ)
sage: W.next()
(0, 0)
sage: W.next()
(0, 1)
sage: W.next()
(0, -1)循环、函数、控制语句和比较
我们已经看到一些例子,其中有for循环的常见用法。在Python中,for循环有缩进结构,例如

>>> for i in range(5):
       print(i)
 
0
1
2
3
4注意for语句末尾的冒号(这里不像GAP或Maple,没有“do”或者“od”),以及循环“主体”(即print(i))前面的缩进。这个缩进是很重要的。在Sage中,当您在“:”后按enter时,缩进将自动添加上去,如下所示。

sage: for i in range(5):
...       print(i)  # 现在按两次回车
0
1
2
3
4符号=用来赋值。符号==用来检查是否相等:

sage: for i in range(15):
...       if gcd(i,15) == 1:
...           print(i)
1
2
4
7
8
11
13
14记住缩进是如何决定if、for和while语句的块结构的:

sage: def legendre(a,p):
...       is_sqr_modp=-1
...       for i in range(p):
...           if a % p == i^2 % p:
...               is_sqr_modp=1
...       return is_sqr_modp
 
sage: legendre(2,7)
1
sage: legendre(3,7)
-1当然这不是实现勒让德符号的有效途径!它是用来说明Python/Sage编程的各个方面。Sage中的函数{kronecker}通过C语言库调用PARI来高效计算勒让德符号。

最后,我们注意到数字之间的比较运算符,例如==、!=、<=、>=、>、<,如果可能则会自动把两边的数字转换成同一种类型:

sage: 2 < 3.1; 3.1 <= 1
True
False
sage: 2/3 < 3/2;   3/2 < 3/1
True
True几乎任何两个对象都可以被比较;对象是全序的这一假设是不存在的。

sage: 2 < CC(3.1,1)
True
sage: 5 < VectorSpace(QQ,3)   # 输出是随机的
True对于符号不等式使用bool函数:

sage: x < x + 1
x < x + 1
sage: bool(x < x + 1)
True在Sage中比较不同类型对象时,在多数情况下Sage尝试将两个对象规范地强制转换到同一父类。如果成功,比较在两个强制转换后的对象间进行;如果不成功则对象被认为是不相等的。要检查两个变量是否引用同一对象,使用is。例如:

sage: 1 is 2/2
False
sage: 1 is 1
False
sage: 1 == 2/2
True在下面两行中,第一个等式是False,因为没有规范射[Math Processing Error],从而没有规范方法比较[Math Processing Error]中的1和[Math Processing Error]。相反,有规范映射[Math Processing Error],因此第二个比较是True。注意顺序并不重要。

sage: GF(5)(1) == QQ(1); QQ(1) == GF(5)(1)
False
False
sage: GF(5)(1) == ZZ(1); ZZ(1) == GF(5)(1)
True
True
sage: ZZ(1) == QQ(1)
True警告:Sage中的比较比Magma中更严格,Magma中[Math Processing Error]等于[Math Processing Error]。

sage: magma('GF(5)!1 eq Rationals()!1')            # 需要magma
true性能分析
本节作者:Martin Albrecht(malb@informatik.uni-bremen.de

“Premature optimization is the root of all evil.”——Donald Knuth

有时检查代码中的瓶颈对于明白哪个部分占用最多计算时间是很有用的;这样可以清楚了解要优化哪部分。Python和Sage提供几种性能分析(这个过程的称呼)选项。

最容易使用的是交互式Shell中的prun命令。它返回一个关于哪个函数使用多少计算时间的概要。例如,要分析8)有限域上的矩阵乘法,执行:

sage: k,a = GF(2**8, 'a').objgen()
sage: A = Matrix(k,10,10,[k.random_element() for _ in range(10*10)])sage: %prun B = A*A
       32893 function calls in 1.100 CPU seconds
 
Ordered by: internal time
 
ncalls tottime percall cumtime percall filename:lineno(function)
 12127  0.160   0.000   0.160  0.000 :0(isinstance)
  2000  0.150   0.000   0.280  0.000 matrix.py:2235(__getitem__)
  1000  0.120   0.000   0.370  0.000 finite_field_element.py:392(__mul__)
  1903  0.120   0.000   0.200  0.000 finite_field_element.py:47(__init__)
  1900  0.090   0.000   0.220  0.000 finite_field_element.py:376(__compat)
   900  0.080   0.000   0.260  0.000 finite_field_element.py:380(__add__)
     1  0.070   0.070   1.100  1.100 matrix.py:864(__mul__)
  2105  0.070   0.000   0.070  0.000 matrix.py:282(ncols)
  ...这里ncalls是调用的次数,tottime是给定函数用的总时间(不包括调用子函数的时间),percall是tottime除以ncalls的商。cumtime是这个函数和所有子函数所用的时间(即从开始到退出),percall是cumtime除以基础调用9)的商,并且filename:lineno(function)提供了每个函数的相关数据。这里的基本原则是10):函数在列表中越靠前,就越花时间。因此对优化更有吸引力。

通常,prun?提供如何使用分析器11)和理解输出结果的详细信息。

分析数据也可以写入一个对象中,以便更仔细地检查:

sage: %prun -r A*A
sage: stats = _
sage: stats?注意:输入stats = prun -r A\*A会显示语法错误信息,因为prun是IPython shell命令,不是正常的函数。

要得到分析数据漂亮的图像表示,您可以使用hotshot分析器,一个调用hotshot2cachetree和kcachegrind程序(仅Unix)的小脚本。用hotshot分析器分析同一个例子12):

sage: k,a = GF(2**8, 'a').objgen()
sage: A = Matrix(k,10,10,[k.random_element() for _ in range(10*10)])
sage: import hotshot
sage: filename = "pythongrind.prof"
sage: prof = hotshot.Profile(filename, lineevents=1)sage: prof.run("A*A")
<hotshot.Profile instance at 0x414c11ec>
sage: prof.close()这会在当前工作目录产生pythongrind.prof的文件。它现在可以转换为可视化的cachegrind格式。

在系统shell里输入

hotshot2calltree -o cachegrind.out.42 pythongrind.prof输出文件cachegrind.out.42现在可以用kcachegrind检查。请注意需要遵守命名习惯cachegrind.out.XX。

1) Enter键/回车键
2) 如果下面编译出现错误,请将“test.c”改为绝对路径,如“/home/username/test.c”
3) 如果返回“/usr/bin/env: sage -python: 没有该文件或目录 ”,则将第一句改为“#!/usr/bin/env sage”
4) 实际上只有对整数的因式分解可以运行,第二个例子会返回“NameError: name 'x' is not defined ”。此问题最近有人向Sage官方反映了,不过还没结论
5) 原文:list comprehensions。也有译为列表综合、列表包含、列表推导式的。我不确定哪个最合适
6) 原文:associative array
7) 原文:hashable
8) 原文:To profile
9) 原文:primitive calls
10) 原文:The rule of thumb here is
11) 原文:profiler
12) 其实我这里出错了……还没弄明白什么原因

posted on 2012-04-11 21:44  孤独的猫  阅读(1937)  评论(0编辑  收藏  举报