Next we illustrate how to load programs written in a separate file
into SAGE. Create a file called example.sage
with the following
content:
print "Hello World" print 2^3
You can read in and execute example.sage
file using the load
command.
sage: load "example.sage" Hello World 8
You can also attach a SAGE file to a running session using
the attach
command:
sage: attach "example.sage" Hello World 8
example.sage
and enter one blank line into
SAGE (i.e., hit ``return''), then the contents of example.sage
will be automatically reloaded into SAGE.
In particular, attach has the side effect of (auto-reload), very handy when debugging code, while load does not.
When SAGE loads example.sage
it converts it to Python,
which is then executed by the Python interpreter. This conversion
is minimal; it mainly involves wrapping integer literals in ZZ()
,
floating point literals in RR()
, replacing ^
's by **
's,
and replacing e.g., R.2
by R.gen(2)
.
The converted version of example.sage
is contained in the same directory as example.sage
and is called example.sage.py
. This file contains the
following code:
print "Hello World" print ZZ(2)**ZZ(3)
^
is
replaced by a **
. (In Python ^
means ``exclusive
or'' and **
means ``exponentiation''.)
Note: This preparsing is implemented in sage/misc/interpreter.py.
You can paste multi-line indented code into SAGE as long as
there are newlines to make new blocks (this is not necessary in
files). However, the best way to enter such code into SAGE is to save
it to a file and use attach
, as described above.
Speed is crucial in mathematical computations. Though Python is a convenient very high-level language, certain calculations can be several orders of magnitude faster than in Python if they are implemented using static types in a compiled language. It is virtually impossible to write serious competitive computer algebra software if one restricts oneself to interpreted Python code.
To deal with this problem, SAGE supports a compiled ``version'' of
Python called Pyrex (see [Pyr]). Pyrex is simultaneously similar
to both Python and C. List comprehensions are not allowed and
expressions like +=
are not allowed, but most other Python
constructions are allowed, such as importing code that you have
written in other Python modules. Morever, one can declare arbitrary C
variables and arbitrary C library calls can be made directly. The
resulting code is converted to C and compiled using a C compiler.
In order to make your own compiled SAGE code, give the file an .spyx
extension (instead of .sage
). You can attach and load compiled
code exactly like with interpreted code. The actual compilation is
done ``behind the scenes'' without your having to do anything
explicitly. See SAGE_ROOT/examples/pyrex/factorial.spyx
for an
example of a very fast compiled implementation of the factorial
function that directly uses the GMP C library. To try this out for
yourself, cd to SAGE_ROOT/examples/pyrex/
, then do the following:
sage: load "factorial.spyx" *************************************************** Recompiling factorial.spyx *************************************************** sage: factorial(50) 1520704660085668902180630408303238442218882078448025600000000000000L sage: time n = factorial(10000) CPU times: user 0.06 s, sys: 0.00 s, total: 0.06 s Wall time: 0.06
Note that SAGE will not recompile factorial.spyx
unless you
change it, even if you quit and restart SAGE. The compiled shared
object library is stored under $HOME/.sage/spyx
.
Full SAGE preparsing is applied to spyx files unless they
contain the string __no_preparse__
. If foo
is a function in the SAGE library, it is available in
an spyx file via sage.foo( ... )
.
The following standalone SAGE script factors integers, polynomials, etc:
#!/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
must
be in your PATH. If the above script is called factor
,
here is an example usage:
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: s = "sage"; type(s) <type 'str'> sage: s = 'sage'; type(s) # you can use either single or double quotes <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'>
To this SAGE adds many other types. E.g., vector spaces:
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'>
Only certain functions can be called on
. In other math software
systems, these would be called using the ``functional'' notation foo(V,...).
In SAGE, certain functions are attached to the type (or class) of
, and are called using an object-oriented syntax like in Java or
C++, e.g.,
V.foo(...)
. This helps keep the global namespace
from being polluted with tens of thousands of functions, and means
that many different functions with different behavior can be named
foo, without having to use type-checking of arguments (or case
statements) to decide which to call. Also, if you reuse the name
of a function, that function is still available (e.g., if you
call something zeta
, then want to compute the value of the
Riemann-Zeta function at 0.5, you can still type (0.5).zeta()).
sage: zeta = -1 sage: (0.5).zeta() -1.4603545088095868
In some very common cases the usual functional notation is also supported for convenience and because mathematical expressions might look confusing using object-oriented notation. Here are some examples.
sage: n = 2; n.sqrt() 1.4142135623730951 sage: sqrt(2) 1.4142135623730951 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^2 + 2*x + 5 sage: charpoly(A) x^2 + 2*x + 5
To list all member functions for
, use tab completion. Just type
A.
, then type the [tab]
key on your keyboard, as
explained in Section 3.6.
The list data type stores elements of arbitrary type. Like in C, C++, etc. (but unlike most standard computer algebra systems), the elements of the list are indexed starting from 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
Note: When indexing into a list, the index must be a Python int! A SAGE Integer will not work. The preparser doesn't change integer literals to Python ints when they are indexing lists. However, if you index a list with a variable, you must be careful.
sage: v = [1,2,3] sage: v[2] 3 sage: n = 2 # SAGE Integer sage: v[n] # bad Traceback (most recent call last): ... TypeError: list indices must be integers sage: v[int(n)] # good 3
The range
function creates a list of Python int's (not SAGE Integers):
sage: range(1, 15) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
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]
List slicing is a wonderful feature. If L
is a list,
then L[m:n]
returns the sublist of L
obtained
by starting at the
th element and stopping at the
st element,
as illustrated below.
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]
Tuples are similar to lists, except they are immutable, meaning once they are created they can't be changed.
sage: v = (1,2,3,4); v (1, 2, 3, 4) sage: type(v) <type 'tuple'> sage: v[1] = 5 Traceback (most recent call last): ... TypeError: object does not support item assignment
Sequences are a third list-oriented SAGE type. Unlike lists and tuples,
Sequence is not built-in Python type. By default, a sequence is
mutable, but using the Sequence
class method
set_immutable
, it can be set to be immutable, as the following
example illustrates. All elements of a sequence have a common parent,
called the sequences 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 'rational.Rational'> sage: v.universe() Rational Field sage: v.is_immutable() False sage: v.set_immutable() sage: v[0] = 3 Traceback (most recent call last): ... 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'>
As another example, basis for vector spaces are immutable sequences, since it's important that you don't change them.
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] Traceback (most recent call last): ... ValueError: object is immutable; please change a copy instead. sage: B.universe() Vector space of dimension 3 over Rational Field
A dictionary (also sometimes called an associative array) is a mapping from hashable objects to arbitrary objects.
sage: d = {1:5, 'sage':17, ZZ:GF(7)} sage: type(d) <type 'dict'> sage: d.keys() [Integer Ring, 1, 'sage'] sage: d['sage'] 17 sage: d[ZZ] Finite Field of size 7 sage: d[1] 5
The third key illustrates that the indexes of a dictionary can be complicated, e.g., the ring of integers.
You can turn the above dictionary into a list with the same data:
sage: d.items() [(Integer Ring, Finite Field of size 7), (1, 5), ('sage', 17)]
sage: d = {2:4, 3:9, 4:16} sage: [a*b for a, b in d.iteritems()] [8, 27, 64]
A dictionary is unordered, as the last output illustrates.
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 also has its own set type that is (in some cases) implemented using the built-in Python set type, but has a little bit of extra SAGE-related functionality. Create a SAGE set using Set(...). For example,
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
sage: v = (n^2 for n in xrange(10000000)) sage: v.next() 0 sage: v.next() 1 sage: v.next() 4
sage: w = (4*p + 1 for p in Primes() if is_prime(4*p+1)) sage: w # random 0x number <generator object at 0xb0853d6c> sage: w.next() 13 sage: w.next() 29 sage: w.next() 53
Certain rings, e.g., finite fields and the integers have iterators associated to them:
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)
We have seen a few examples already of some common uses of
for
loops. In Python, a for
loop has an indented structure,
such as
>>> for i in range(5): print(i) 0 1 2 3 4
print(i)
. This indentation is important.
In SAGE, the indentation is automatically put in for you
when you hit enter
after a ``:'', as illustrated below.
sage: for i in range(5): print(i) # now hit enter twice 0 1 2 3 4 sage:
The symbol =
is used for assignment.
The symbol ==
is used to check for equality:
sage: for i in range(15): if gcd(i,15) == 1: print(i) 1 2 4 7 8 11 13 14
if
, for
, and while
statements:
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
kronecker
, which comes with
SAGE, computes the Legendre symbol efficiently
via a C-library call to PARI.
Finally, we note that comparisons,
such as ==
, !=
,
<=
, >=
, >
, <
,
between numbers will automatically convert both
numbers into the same type if possible:
sage: 2 < 3.1; 3.1 <= 1 True False sage: 2/3 < 3/2; 3/2 < 3/1 True True
sage: 2 < 3.1+0.0*i; 3.1+2*i<4+3*i; 4+3*i < 3.1+2*i True True False sage: 5 < VectorSpace(QQ,3) True
is
. For example:
sage: 1 is 2/2 False sage: 1 is 1 False sage: 1 == 2/2 True
False
because
there is no canonical morphism
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
WARNING: Comparison in SAGE is more restrictive than in
Magma, which declares the
equal to
.
sage: magma('GF(5)!1 eq Rationals()!1') # optional magma required true
If you want to extend SAGE by adding your own methods to a SAGE class there is a convenient way built into Python. You should understand the basics of object-oriented programming in Python in order to understand how to override class methods (see, e.g., the excellent free Python tutorial [PyT]). In Python nearly everything is changeable, so for instance you can alter class functions on the fly. We will use this by adding a rather silly method to the Matrix class in the following example. We start by writing a MatrixMixIn class:
class MatrixMixIn: def ncols_plus_nrows(self): """ Return the number of rows plus the number of columns. """ return self.ncols() + self.nrows()
Now we load the file where we defined this class into SAGE and add
this class as a superclass to the Matrix class (you could also
just paste the above code into your SAGE session). The superclasses are
stored in an attribute called __bases__
for every class (except
builtin compiled extension classes). It is
a tuple of classes, so we just add our class to this tuple:
sage: sage.matrix.matrix.Matrix.__bases__ += (MatrixMixIn,)
You can also get the class of an object using the
__class__
attribute of any object of that class:
sage: A = Matrix(Integers(8),2,2,[1,2,3,4]) sage: A.__class__ <class 'sage.matrix.matrix.Matrix_generic_dense'>
But note that we cannot write Matrix.__bases__
above as
Matrix
is a function to construct matrices and not the Matrix
class. This is not always the case, but it's just the way Matrices are
treated in SAGE:
sage: type(Matrix) <type 'function'> sage: type(sage.matrix.matrix.Matrix) <type 'type'> sage: type(Integer) <type 'type'>
Now every instance of this class and every instance of any class inheriting from this class has our new method.
sage: A = Matrix(GF(2**8),10,10,) sage: A.ncols_plus_nrows() 20 sage: A.ncols_plus_nrows? Type: instancemethod Base Class: <type 'instancemethod'> String Form: <bound method Matrix_generic_dense_field.ncols_plus_nrows of [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 <...> 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0] [0 0 0 0 0 0 0 0 0 0]> Namespace: Interactive File: /home/martin/Uni-Bremen/BES/trunk/.test.sage.py Definition: A.ncols_plus_nrows(self) Docstring: Return the number of rows plus the number of columns.
When you restart SAGE the new method will be gone and you will have to mix it in again. So this is a non intrusive way to patch SAGE temporarily. (And you could make it permanent for your personal use by putting the above code in a file and loading it when you start SAGE.)
If you add any functionality which might be of general interest,
please consider sending it to William Stein or posting it to
sage-forum@lists.sourceforge.net
, so he can include it into the
main SAGE distribution (see Section 6.2 later in this
tutorial).
``Premature optimization is the root of all evil.'' - Donald Knuth
Sometimes it is useful to check for bottlenecks in code to understand which parts take the most computational time; this can give a good idea of which parts to optimize. Python and therefore SAGE offers several profiling--as this process is called--options.
The simplest to use is the prun
command in the interactive
shell. It returns a summary describing which functions took how much
computational time. To profile (the currently slow! - as of version 1.0) matrix
multiplication over finite fields, for example, do:
sage: k, a = GF(2**8).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) ...
Here ncalls
is the number of calls, tottime
is the total time
spent in the given function (and excluding time made in calls to sub-functions),
percall
is the quotient of tottime
divided by ncalls
.
cumtime
is the total time spent in this and all sub-functions (i.e.,
from invocation until exit), percall
is the quotient of cumtime
divided by primitive calls, and filename:lineno(function)
provides
the respective data of each function. The rule of thumb here is: The
higher the function in that listing the more expensive it is. Thus it
is more interesting for optimization.
As usual, prun?
provides details on how to use the profiler and
understand the output.
The profiling data may be written to an object as well to allow closer examination:
sage: prun -r A*A sage: stats = _ sage: ?stats
Please note that you cannot do a stats = prun -r A*A
for some
internal reason.
For a nice graphical representation of profiling data you can use the
hotshot profiler, a small script called hotshot2cachetree
and
the program kcachegrind
(Unix only). The same example with the
hotshot profiler:
sage: k, a = GF(2**8).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()
This results in a file pythongrind.prof
in the current working
directory. It can now be converted to the cachegrind format for
visualization.
On a system shell type
hotshot2calltree -o cachegrind.out.42 pythongrind.prof
The output file cachegrind.out.42
can now be examined with
kcachegrind
. Please note that the naming convention
cachegrind.out.XX
needs to be obeyed.
See About this document... for information on suggesting changes.