Sage Demo - before talk
system:sage


<h2>The Sage Mathematical Software System</h2>
<h3>What is Sage?</h3>
<p>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.</p>
<h3>Sage command language</h3>

{{{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
///
}}}

<h3>Cython demo</h3>

{{{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)
///
}}}

<h3>3D graphics</h3>

{{{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()
///
}}}

<h2>Some of my contributions to Sage</h2>
<h3>Repeatable pseudo-random numbers</h3>

{{{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()
///
}}}

<p>set_random_seed() actually controls 7 random number generators: GMP, Python, NTL, libpari, pari/gp, GAP (2 generators); Singular is missing</p>
<h3>RealIntervalField (wrapper for MPFI)</h3>

{{{id=38|
RIF128 = RealIntervalField(128)
RIF128.pi()
///
}}}

{{{id=37|
i3 = RIF128(3)
i4 = RIF128(4)
i3/i4*i4
///
}}}

{{{id=40|
i4/i3*i3
///
}}}

<h3>Real root isolation</h3>

{{{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
///
}}}

<h3>Field of algebraic numbers (real and complex)</h3>

{{{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])
///
}}}

<h3>Interface to QEPCAD B</h3>

{{{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)
///
}}}

<h2>A tour of the Sage packages</h2>
<p>Now I'll go through most of the packages included in Sage, and include at least one computation that makes use of each package.&nbsp; (I'll skip some of the non-mathematical packages.)</p>
<h3>ATLAS (floating point matrix manipulation)</h3>

{{{id=69|
M = random_matrix(RDF, 3, 3)
///
}}}

{{{id=71|
M, M.inverse()
///
}}}

{{{id=72|
M * M.inverse()
///
}}}

<h3>cddlib</h3>

{{{id=73|
P = Polyhedron(ieqs = [[0,1,0],[0,0,1],[1,-1,-1]])
P.Hrepresentation()
///
}}}

{{{id=75|
P.Vrepresentation()
///
}}}

<h3>cliquer</h3>

{{{id=76|
G = Graph({0:[1,2,3], 1:[2], 3:[0,1]})
G.show(figsize=[4,3])
clique_number(G)
///
}}}

<h3>Conway polynomial database</h3>

{{{id=78|
conway_polynomial(37, 5)
///
}}}

<h3>cvxopt</h3>
<p>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$.</p>

{{{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']
///
}}}

<h3>ecl</h3>

{{{id=83|
%lisp
(defun hello (x)
  (concatenate 'string "Hello, " x))
///
}}}

{{{id=82|
%lisp
(hello "world")
///
}}}

<h3>eclib</h3>

{{{id=85|
E = EllipticCurve('11a1')
EE = E.mwrank_curve()
EE
///
}}}

{{{id=87|
EE.isogeny_class()
///
}}}

<h3>ecm</h3>

{{{id=88|
import sage.libs.libecm
from sage.libs.libecm import ecmfactor
ecmfactor(999, 0.00, verbose=True)
///
}}}

<h3>elliptic curves database</h3>

{{{id=90|
elliptic_curves.rank(n=5, rank=3, tors=2, labels=True)
///
}}}

<h3>FLINT</h3>

{{{id=92|
K = PolynomialRing(ZZ, 'x', implementation='FLINT')
p = K([-1, -1, 1]); (p, type(p))
///
}}}

<h3>flintqs</h3>

{{{id=98|
k = 19; n = next_prime(10^k)*next_prime(10^(k+1))
n
///
}}}

{{{id=97|
qsieve(n, time=True)
///
}}}

<h3>GAP</h3>

{{{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)
///
}}}

<h3>genus2reduction</h3>

{{{id=106|
x = polygen(QQ)
R = genus2reduction(x^3 - 2*x^2 - 2*x + 1, -5*x^5)
R
///
}}}

<h3>gfan</h3>

{{{id=108|
R.<x,y> = PolynomialRing(QQ, 2)
gf = R.ideal([x^3 - y, y^3 - x - 1]).groebner_fan()
print gf.gfan()
///
}}}

<h3>givaro</h3>

{{{id=110|
K.<a> = GF(3^4)
///
}}}

{{{id=114|
(a+a^2)^7 / (a-1)
///
}}}

<h3>glpk</h3>

{{{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)
///
}}}

<h3>graphs database</h3>

{{{id=117|
Q = GraphQuery(display_cols=['graph6','num_vertices','degree_sequence'],num_edges=['<=',5],min_degree=1)
Q.number_of()
///
}}}

{{{id=119|
Q.show()
///
}}}

<h3>GSL</h3>

{{{id=120|
RDF.factorial(100)
///
}}}

<h3>IML</h3>

{{{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")
///
}}}

<h3>lcalc</h3>

{{{id=125|
lcalc.zeros(4)
///
}}}

<h3>libfplll</h3>

{{{id=127|
A = Matrix(ZZ,3,3,range(1,10))
A.LLL()
///
}}}

<h3>libm4ri</h3>

{{{id=129|
M1 = random_matrix(GF(2), 1000, 1000)
M2 = random_matrix(GF(2), 1000, 1000)
M1 * M2
///
}}}

<h3>linbox</h3>

{{{id=131|
A = matrix(ZZ,6, range(36))
A.minpoly()
///
}}}

<h3>matplotlib</h3>

{{{id=133|
var('x')
plot(sin(x^2), (x, -10, 10))
///
}}}

<h3>maxima</h3>

{{{id=135|
var('x')
integral(sin(x^2), x)
///
}}}

<h3>MPFI</h3>

{{{id=137|
reset(['pi'])
RealIntervalField(1024)(pi)
///
}}}

<h3>MPFR</h3>

{{{id=139|
RealField(1024)(pi)
///
}}}

<h3>MPIR (fork of GMP)</h3>

{{{id=141|
factorial(1024)
///
}}}

<h3>mpmath</h3>

{{{id=143|
Ei(3+I).n()
///
}}}

<h3>networkx</h3>

{{{id=145|
graphs.BarbellGraph(5, 7).show()
///
}}}

<h3>NTL</h3>

{{{id=147|
K = PolynomialRing(ZZ, 'x', implementation='NTL')
p = K([-1, -1, 1]); (p, type(p))
///
}}}

<h3>numpy</h3>

{{{id=149|
V = vector(RDF, range(1, 7))
V.standard_deviation()
///
}}}

<h3>palp</h3>

{{{id=151|
o = lattice_polytope.octahedron(3)
lattice_polytope.all_nef_partitions([o])
o.nef_partitions()
///
}}}

<h3>pari</h3>

{{{id=153|
factorial(35).factor(algorithm="pari")
///
}}}

{{{id=155|
%gp
factor(binomial(50, 17))
///
}}}

<h3>polybori</h3>

{{{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()
///
}}}

<h3>polytopes database</h3>

{{{id=158|
ReflexivePolytopes(2)
///
}}}

<h3>pynac (fork of ginac)</h3>

{{{id=160|
var('x,y,z')
expand((x+y)^7 - x^7)
///
}}}

<h3>R</h3>

{{{id=162|
x = r([10.4,5.6,3.1,6.4,21.7]); x
///
}}}

{{{id=164|
x.var()
///
}}}

<h3>ratpoints</h3>

{{{id=165|
from sage.libs.ratpoints import ratpoints
ratpoints([1..6], 200)
///
}}}

<h3>rubik's cube solvers</h3>

{{{id=167|
from sage.interfaces.rubik import *
solver = DikSolver()
C = RubiksCube("R U F L B D")
solver.solve(C.facets())
///
}}}

<h3>scipy</h3>

{{{id=169|
bessel_I(1,1,"scipy")
///
}}}

<h3>singular</h3>

{{{id=171|
K.<x,y> = QQ[]
(x+y)^7
///
}}}

<h3>symmetrica</h3>

{{{id=176|
SemistandardTableaux([3,2,1], [2, 2, 2]).cardinality()
///
}}}

<h3>sympow</h3>

{{{id=178|
E = EllipticCurve('5077a'); E
///
}}}

{{{id=180|
E.modular_degree()
///
}}}

<h3>sympy</h3>

{{{id=181|
var('x,y,z')
(x^y - z).integrate(y, algorithm="sympy")
///
}}}

<h3>zn_poly</h3>

{{{id=183|
P.<x> = PolynomialRing(GF(next_prime(2^30)))
f = P.random_element(20)
g = P.random_element(20)
f._mul_zn_poly(g)
///
}}}

{{{id=185|

///
}}}