The Sage Mathematical Software System
What is Sage?
Sage is a software distribution containing almost 100 software packages, and a huge Python library (the "Sage library") to tie them all together into a (somewhat) unified user interface.
Sage command language
{{{id=2|
', '.join(['Hello', "world!"])
///
}}}
{{{id=5|
%python
5/3
///
}}}
{{{id=7|
5/3
///
}}}
{{{id=8|
preparse('5/3')
///
}}}
{{{id=11|
%python
# This is bitwise XOR
2^7
///
}}}
{{{id=9|
# This is exponentiation
2^7
///
}}}
Cython demo
{{{id=12|
def py_gcd(a, b):
if a < b:
a, b = b, a
while b != 0:
a, b = b, a%b
return a
def py_totient(n):
count = 0
for i in range(1, n):
if py_gcd(i, n) == 1:
count += 1
return count
///
}}}
{{{id=16|
timeit('py_totient(123)')
///
}}}
{{{id=17|
%cython
cpdef unsigned long cy_gcd(unsigned long a, unsigned long b):
if a < b:
a, b = b, a
while b != 0:
a, b = b, a%b
return a
def cy_totient(unsigned long n):
cdef unsigned long count = 0
cdef unsigned long i
for i in range(1, n):
if cy_gcd(i, n) == 1:
count += 1
return count
///
}}}
{{{id=19|
timeit('cy_totient(123)')
///
}}}
{{{id=18|
timeit('cy_totient(123456)')
///
}}}
{{{id=14|
timeit('euler_phi(123456)')
///
}}}
{{{id=20|
%cython
from sage.libs.mpfr cimport *
def pi_str(int bits):
cdef mpfr_t pi
mpfr_init2(pi, bits)
mpfr_const_pi(pi, GMP_RNDN)
cdef mp_exp_t exp
cdef char *out = mpfr_get_str(NULL, &exp, 10, 0, pi, GMP_RNDN)
result = '%s * 10^%d' % (str(out), exp-int(strlen(out)))
mpfr_free_str(out)
return result
///
}}}
{{{id=25|
pi_str(30)
///
}}}
{{{id=26|
timeit('pi_str(1024)')
///
}}}
{{{id=21|
pi_str(1024)
///
}}}
3D graphics
{{{id=27|
var('x,y')
P = plot3d(x^2 + 15*sin(y), (x, -10, 10), (y, -10, 10))
P.show(viewer='tachyon')
///
}}}
{{{id=24|
P.show()
///
}}}
{{{id=29|
var('x,y,z')
T = RDF(golden_ratio)
p = 2 - (cos(x + T*y) + cos(x - T*y) + cos(y + T*z) + cos(y - T*z) + cos(z - T*x) + cos(z + T*x))
rad = 4.77
implicit_plot3d(p, (x, -rad, rad), (y, -rad, rad), (z, -rad, rad), plot_points=80).show()
///
}}}
Some of my contributions to Sage
Repeatable pseudo-random numbers
{{{id=30|
ZZ.random_element(10^6)
///
}}}
{{{id=32|
set_random_seed(1234)
ZZ.random_element(10^6)
///
}}}
{{{id=33|
set_random_seed(1234)
ZZ.random_element(10^6)
///
}}}
{{{id=34|
G = PermutationGroup([[(1,2,3),(4,5)], [(1,2)]])
set_random_seed(5678)
G.random_element()
///
}}}
{{{id=35|
set_random_seed(1234)
G.random_element()
///
}}}
{{{id=36|
set_random_seed(1234)
G.random_element()
///
}}}
set_random_seed() actually controls 7 random number generators: GMP, Python, NTL, libpari, pari/gp, GAP (2 generators); Singular is missing
RealIntervalField (wrapper for MPFI)
{{{id=38|
RIF128 = RealIntervalField(128)
RIF128.pi()
///
}}}
{{{id=37|
i3 = RIF128(3)
i4 = RIF128(4)
i3/i4*i4
///
}}}
{{{id=40|
i4/i3*i3
///
}}}
Real root isolation
{{{id=41|
x = polygen(QQ)
(x^2-x-1).roots(ring=RIF)
///
}}}
{{{id=44|
(x^495 * (x^2 - 9999)^2 - 1)
///
}}}
{{{id=45|
rts = (x^495 * (x^2 - 9999)^2 - 1).roots(ring=RIF)
///
}}}
{{{id=46|
rts
///
}}}
{{{id=47|
rts[2][0] - rts[1][0]
///
}}}
{{{id=48|
rts[2][0]^2
///
}}}
Field of algebraic numbers (real and complex)
{{{id=49|
rt2 = AA(sqrt(2))
rt2, rt2^2
///
}}}
{{{id=51|
rt2^2 == 2
///
}}}
{{{id=52|
QQbar.zeta(5)
///
}}}
{{{id=53|
QQbar.zeta(5).real()
///
}}}
{{{id=54|
foo = AA(sqrt(3) + sqrt(5))
///
}}}
{{{id=55|
number_field_elements_from_algebraics([foo])
///
}}}
{{{id=56|
foo = QQbar(sqrt(3) + I*sqrt(5))
///
}}}
{{{id=57|
number_field_elements_from_algebraics([foo])
///
}}}
Interface to QEPCAD B
{{{id=58|
var('x,y')
qf = qepcad_formula
ellipse = 3*x^2 + 2*x*y + y^2 - x + y - 7
ellplot = implicit_plot(ellipse, (x, -6, 6), (y, -6, 6))
ellplot.show(aspect_ratio=1)
///
}}}
{{{id=60|
F = qf.exists(y, ellipse == 0); F
///
}}}
{{{id=61|
qepcad(F)
///
}}}
{{{id=62|
xx = polygen(QQ)
rts = (8*xx^2 - 8*xx - 29).roots(ring=RDF, multiplicities=False)
rts
///
}}}
{{{id=65|
def vline(k):
return implicit_plot(x - k, (x, -6, 6), (y, -6, 6))
///
}}}
{{{id=63|
(ellplot + vline(rts[0]) + vline(rts[1])).show(aspect_ratio=1)
///
}}}
{{{id=64|
circle = x^2 + y^2 - 3
cirplot = implicit_plot(circle, (x, -6, 6), (y, -6, 6))
(cirplot + ellplot).show(aspect_ratio=1)
///
}}}
{{{id=66|
F = qf.exactly_k(3, y, circle * ellipse == 0); F
///
}}}
{{{id=67|
pts = qepcad(F, solution='all-points'); pts
///
}}}
{{{id=68|
(cirplot + ellplot + sum(map(lambda r: vline(RDF(r['x'])), pts))).show(aspect_ratio=1)
///
}}}
A tour of the Sage packages
Now I'll go through most of the packages included in Sage, and include at least one computation that makes use of each package. (I'll skip some of the non-mathematical packages.)
ATLAS (floating point matrix manipulation)
{{{id=69|
M = random_matrix(RDF, 3, 3)
///
}}}
{{{id=71|
M, M.inverse()
///
}}}
{{{id=72|
M * M.inverse()
///
}}}
cddlib
{{{id=73|
P = Polyhedron(ieqs = [[0,1,0],[0,0,1],[1,-1,-1]])
P.Hrepresentation()
///
}}}
{{{id=75|
P.Vrepresentation()
///
}}}
cliquer
{{{id=76|
G = Graph({0:[1,2,3], 1:[2], 3:[0,1]})
G.show(figsize=[4,3])
clique_number(G)
///
}}}
Conway polynomial database
{{{id=78|
conway_polynomial(37, 5)
///
}}}
cvxopt
We minimize $-4x_1 - 5x_2$ subject to $2x_1 + x_2 \leq 3$, $x_1 + 2x_2 \leq 3$, $x_1 \geq 0$, and $x_2 \geq 0$.
{{{id=80|
c=vector(RDF,[-4,-5])
G=matrix(RDF,[[2,1],[1,2],[-1,0],[0,-1]])
h=vector(RDF,[3,3,0,0])
sol=linear_program(c,G,h)
sol['x']
///
}}}
ecl
{{{id=83|
%lisp
(defun hello (x)
(concatenate 'string "Hello, " x))
///
}}}
{{{id=82|
%lisp
(hello "world")
///
}}}
eclib
{{{id=85|
E = EllipticCurve('11a1')
EE = E.mwrank_curve()
EE
///
}}}
{{{id=87|
EE.isogeny_class()
///
}}}
ecm
{{{id=88|
import sage.libs.libecm
from sage.libs.libecm import ecmfactor
ecmfactor(999, 0.00, verbose=True)
///
}}}
elliptic curves database
{{{id=90|
elliptic_curves.rank(n=5, rank=3, tors=2, labels=True)
///
}}}
FLINT
{{{id=92|
K = PolynomialRing(ZZ, 'x', implementation='FLINT')
p = K([-1, -1, 1]); (p, type(p))
///
}}}
flintqs
{{{id=98|
k = 19; n = next_prime(10^k)*next_prime(10^(k+1))
n
///
}}}
{{{id=97|
qsieve(n, time=True)
///
}}}
GAP
{{{id=99|
A4 = PermutationGroup([[(1,2,3)],[(2,3,4)]]); A4
///
}}}
{{{id=100|
A4.center()
///
}}}
{{{id=103|
n = gap(20062006); n
///
}}}
{{{id=104|
n.Factors()
///
}}}
{{{id=105|
%gap
Factors(2010)
///
}}}
genus2reduction
{{{id=106|
x = polygen(QQ)
R = genus2reduction(x^3 - 2*x^2 - 2*x + 1, -5*x^5)
R
///
}}}
gfan
{{{id=108|
R. = PolynomialRing(QQ, 2)
gf = R.ideal([x^3 - y, y^3 - x - 1]).groebner_fan()
print gf.gfan()
///
}}}
givaro
{{{id=110|
K. = GF(3^4)
///
}}}
{{{id=114|
(a+a^2)^7 / (a-1)
///
}}}
glpk
{{{id=115|
g = graphs.PetersenGraph()
p = MixedIntegerLinearProgram(maximization=True)
b = p.new_variable()
p.set_objective(sum([b[v] for v in g]))
for (u,v) in g.edges(labels=None):
p.add_constraint(b[u] + b[v], max=1)
p.set_binary(b)
p.solve(objective_only=True)
///
}}}
graphs database
{{{id=117|
Q = GraphQuery(display_cols=['graph6','num_vertices','degree_sequence'],num_edges=['<=',5],min_degree=1)
Q.number_of()
///
}}}
{{{id=119|
Q.show()
///
}}}
GSL
{{{id=120|
RDF.factorial(100)
///
}}}
IML
{{{id=122|
A = matrix(QQ, [[2,1,-5,-8],[-1,-1,4,6],[1,0,-1,-2]]); A
///
}}}
{{{id=124|
A.right_kernel(algorithm="iml")
///
}}}
lcalc
{{{id=125|
lcalc.zeros(4)
///
}}}
libfplll
{{{id=127|
A = Matrix(ZZ,3,3,range(1,10))
A.LLL()
///
}}}
libm4ri
{{{id=129|
M1 = random_matrix(GF(2), 1000, 1000)
M2 = random_matrix(GF(2), 1000, 1000)
M1 * M2
///
}}}
linbox
{{{id=131|
A = matrix(ZZ,6, range(36))
A.minpoly()
///
}}}
matplotlib
{{{id=133|
var('x')
plot(sin(x^2), (x, -10, 10))
///
}}}
maxima
{{{id=135|
var('x')
integral(sin(x^2), x)
///
}}}
MPFI
{{{id=137|
reset(['pi'])
RealIntervalField(1024)(pi)
///
}}}
MPFR
{{{id=139|
RealField(1024)(pi)
///
}}}
MPIR (fork of GMP)
{{{id=141|
factorial(1024)
///
}}}
mpmath
{{{id=143|
Ei(3+I).n()
///
}}}
networkx
{{{id=145|
graphs.BarbellGraph(5, 7).show()
///
}}}
NTL
{{{id=147|
K = PolynomialRing(ZZ, 'x', implementation='NTL')
p = K([-1, -1, 1]); (p, type(p))
///
}}}
numpy
{{{id=149|
V = vector(RDF, range(1, 7))
V.standard_deviation()
///
}}}
palp
{{{id=151|
o = lattice_polytope.octahedron(3)
lattice_polytope.all_nef_partitions([o])
o.nef_partitions()
///
}}}
pari
{{{id=153|
factorial(35).factor(algorithm="pari")
///
}}}
{{{id=155|
%gp
factor(binomial(50, 17))
///
}}}
polybori
{{{id=156|
sr = mq.SR(1, 1, 1, 4, gf2=True, polybori=True)
P = sr.vector([0, 0, 1, 0])
K = sr.vector([1, 0, 0, 1])
F, s = sr.polynomial_system(P, K)
F.groebner_basis()
///
}}}
polytopes database
{{{id=158|
ReflexivePolytopes(2)
///
}}}
pynac (fork of ginac)
{{{id=160|
var('x,y,z')
expand((x+y)^7 - x^7)
///
}}}
R
{{{id=162|
x = r([10.4,5.6,3.1,6.4,21.7]); x
///
}}}
{{{id=164|
x.var()
///
}}}
ratpoints
{{{id=165|
from sage.libs.ratpoints import ratpoints
ratpoints([1..6], 200)
///
}}}
rubik's cube solvers
{{{id=167|
from sage.interfaces.rubik import *
solver = DikSolver()
C = RubiksCube("R U F L B D")
solver.solve(C.facets())
///
}}}
scipy
{{{id=169|
bessel_I(1,1,"scipy")
///
}}}
singular
{{{id=171|
K. = QQ[]
(x+y)^7
///
}}}
symmetrica
{{{id=176|
SemistandardTableaux([3,2,1], [2, 2, 2]).cardinality()
///
}}}
sympow
{{{id=178|
E = EllipticCurve('5077a'); E
///
}}}
{{{id=180|
E.modular_degree()
///
}}}
sympy
{{{id=181|
var('x,y,z')
(x^y - z).integrate(y, algorithm="sympy")
///
}}}
zn_poly
{{{id=183|
P. = PolynomialRing(GF(next_prime(2^30)))
f = P.random_element(20)
g = P.random_element(20)
f._mul_zn_poly(g)
///
}}}
{{{id=185|
///
}}}