Grand Tour - part 2 system:sage

The Sage Library -- a Grand Tour -- Part 2

June 2009

William Stein

{{{id=242| def ls(module): #auto os.system('ls $SAGE_ROOT/devel/sage/sage/%s'%module) ls('') /// __init__.py coding geometry logic parallel stats __init__.pyc combinat graphs matrix plot structure algebras crypto groups media probability symbolic all.py databases gsl misc quadratic_forms tests all_cmdline.py ext homology modular rings version.py all_notebook.py finance interfaces modules schemes calculus functions lfunctions monoids server categories games libs numerical sets }}} {{{id=4| /// }}}

The Sage library consists of well over 350,000 lines of new code written by Sage developers.  It provides a unified interface to most existing mathematical software and fills in the many gaps in functionality of open source mathematical software.

The Sage library is divided up into the following 39 major modules, and in this talk we will give a tour of the following subset of them:

Part 2

{{{id=382| /// }}} {{{id=381| /// }}}

 

Lfunctions

The lfunctions module provides a Sage interface to the following included Sage components:

There is other L-functions capabilities in Sage, but it is in other parts of the library, e.g., Christian Wuthrich and William Stein have implemented the world's most comprehensive package for $p$-adic $L$-functions of elliptic curves.

PROJECT: Dokchitser's algorithm is being completely coded from scratch in Cython by Robert Bradshaw.  This is the first non-Dokchitser implementation of the algorithm, and mistakes have been corrected in the paper.  He is also working on a provably correct version.

PROJECT: Work under way to wrap Lcalc more directly as a C++ library in Cython for enhanced speed and flexibility.

{{{id=262| ls('lfunctions') /// __init__.py all.py dokchitser.py lcalc.py sympow.py }}}

We compute on the fly imaginary parts of the nontrivial zeros of the Riemann zeta function $\zeta(s)=\sum \frac{1}{n^s}$:

{{{id=220| time lcalc.zeros(10) /// [14.1347251, 21.0220396, 25.0108576, 30.4248761, 32.9350616, 37.5861782, 40.9187190, 43.3270733, 48.0051509, 49.7738325] Time: CPU 0.01 s, Wall: 0.04 s }}}

We compute an elliptic curve $L$-series at $i=\sqrt{-1}$ using Dokchitser (implicitly):

{{{id=384| L = EllipticCurve('37a').lseries() /// }}} {{{id=385| L(I) /// -0.596082884215534 - 0.544769058656713*I }}} {{{id=43| /// }}}

Modular

"The are five fundamental operations of arithmetic: addition, subtraction, multiplication, division, and modular forms." -- Eichler

A modular form (of weight 2) and level $N$ is a holomorphic function $f(z)$ on the extended complex upper half plane such that $f(z)dz$ is invariant under the linear fractional transformation action of the group $\Gamma_0(N)$ of $2\times 2$ integer matrices with determinant $1$ that are upper triangle modulo $N$.  Modular forms played a key role in Wiles' proof of Fermat's Last Theorem, and are at the center of much of modern number theory.

The modular module:

I was teaching an introductory lecture course on modular forms, and I decided I would teach Sage how to do my problem sheets. So I systematically went through the basic things a student might want to calculate, and added the necessary functionality if it wasn't already there and/or fixed it where it was broken. I think most of what I had in mind there is now complete." -- David Loeffler


PROJECTS:

  1. Create a special class for level 1 modular forms, and directly compute all relevant information for these spaces via $q$-expansions, instead of modular symbols. (Craig Citro has various code sitting around to make this yet faster...). 
  2. Compute the ring of modular forms of a given level, and generate $q$-expansions.
  3. Write some code to compute with overconvergent modular symbols (similar to Rob Pollack's Magma code).
  4. Optimize and extend the modular abelian varieties package (not all known algorithms are implemented).
  5. Implement quaternion algebras over number fields, then Dembele's algorithm for Hilbert modular forms.
  6. Period lattices of modular abelian varieties and higher weight modular forms.
  7. In the long term, David Loeffler is interested in
  8. Compute automorphic forms for reductive groups other than $\text{GL}_2(\mathbf{Q})$, e.g. by merging in the code David Loeffler has written for definite unitary groups.
{{{id=268| ls('modular') /// __init__.py congroup_element.py hecke abvar curves modform all.py cusps.py modsym arithgroup dims.py overconvergent buzzard.py dirichlet.py quatalg congroup.py etaproducts.py ssmod }}} {{{id=269| ls('modular/modsym') /// __init__.py g1list.py p1list.c all.py ghlist.py p1list.pxd ambient.py hecke_operator.py p1list.pyx apply.c heilbronn.c relation_matrix.py apply.pxd heilbronn.pyx solaris_fix.h apply.pyx manin_symbols.py space.py boundary.py modsym.py subspace.py element.py modular_symbols.py tests.py }}} {{{id=388| ls('modular/modform') /// __init__.py defaults.py numerical.py all.py eis_series.py periods.py ambient.py eisenstein_submodule.py space.py ambient_R.py element.py submodule.py ambient_eps.py find_generators.py test.py ambient_g0.py half_integral.py tests.py ambient_g1.py hecke_operator_on_qexp.py theta.py constructor.py j_invariant.py vm_basis.py cuspidal_submodule.py notes.py }}}

Basis of modular symbols (written as Manin symbols):

{{{id=116| ModularSymbols(37).basis() /// ((1,0), (1,23), (1,32), (1,34), (1,35)) }}} {{{id=390| /// }}}

Basis of modular forms:

{{{id=115| ModularForms(37).basis() /// [ q + q^3 - 2*q^4 + O(q^6), q^2 + 2*q^3 - 2*q^4 + q^5 + O(q^6), 1 + 2/3*q + 2*q^2 + 8/3*q^3 + 14/3*q^4 + 4*q^5 + O(q^6) ] }}}

Picture of a modular form

{{{id=128| f = CuspForms(37,2).basis()[0].qexp(200).truncate() I = CDF.0; pi = CDF(pi) complex_plot(lambda z: f((2*pi*I*z).exp()), (0,2), (0, .5), plot_points=100).show( frame=True, axes=False) /// }}}

Computing with a modular abelian variety

{{{id=394| J = J0(37); J /// Abelian variety J0(37) of dimension 2 }}} {{{id=117| J.decomposition() /// [ Simple abelian subvariety 37a(1,37) of dimension 1 of J0(37), Simple abelian subvariety 37b(1,37) of dimension 1 of J0(37) ] }}} {{{id=396| /// }}}

Computing a Brandt Module

{{{id=118| B = BrandtModule(97,5); B /// Brandt module of dimension 48 of level 97*5 of weight 2 over Rational Field }}} {{{id=119| B.hecke_matrix(2).fcp() /// (x - 3) * (x - 1)^2 * x^2 * (x^2 - 5) * (x^3 + 2*x^2 - 5*x - 8) * (x^3 + 4*x^2 + 3*x - 1)^2 * (x^4 + x^3 - 4*x^2 - 2*x + 3) * (x^4 - 3*x^3 - x^2 + 6*x - 1)^2 * (x^6 + x^5 - 9*x^4 - 9*x^3 + 17*x^2 + 14*x + 1) * (x^7 - 2*x^6 - 10*x^5 + 18*x^4 + 26*x^3 - 35*x^2 - 21*x + 7) * (x^7 + x^6 - 9*x^5 - 7*x^4 + 23*x^3 + 12*x^2 - 15*x - 8) }}}

Supersingular Module

{{{id=123| X = SupersingularModule(97); X /// Module of supersingular points on X_0(1)/F_97 over Integer Ring }}} {{{id=122| X.supersingular_points() /// ([1, 20, 81*a + 84, 16*a + 68, 20*a + 35, 77*a + 55, 12*a + 75, 85*a + 87], {1: 0, 77*a + 55: 5, 81*a + 84: 2, 12*a + 75: 6, 85*a + 87: 7, 16*a + 68: 3, 20*a + 35: 4, 20: 1}) }}} {{{id=400| /// }}}

An Eta Product

{{{id=121| e = EtaProduct(3, {3:12, 1:-12}); e /// Eta product of level 3 : (eta_1)^-12 (eta_3)^12 }}} {{{id=125| e.qexp(10) /// q + 12*q^2 + 90*q^3 + 508*q^4 + 2391*q^5 + 9828*q^6 + 36428*q^7 + 124188*q^8 + 395199*q^9 + 1186344*q^10 + O(q^11) }}} {{{id=55| /// }}} {{{id=3| /// }}}

Quadratic forms

The quadratic forms code was written mostly by Jon Hanke (with contributions from Gonzalo Tornaria). It:

PROJECT: Expand this code to include QuadraticLattice and QuadraticSpace classes which should allow computations with basic lattice operations, and local-global operations for rational quadratic spaces over number fields.

{{{id=377| ls('quadratic_forms') /// __init__.py all.py binary_qf.py constructions.py count_local_2.c count_local_2.pyx extras.py genera quadratic_form.py quadratic_form__automorphisms.py quadratic_form__count_local_2.py quadratic_form__equivalence_testing.py quadratic_form__evaluate.c quadratic_form__evaluate.pyx quadratic_form__genus.py quadratic_form__local_density_congruence.py quadratic_form__local_density_interfaces.py quadratic_form__local_field_invariants.py quadratic_form__local_normal_form.py quadratic_form__local_representation_conditions.py quadratic_form__mass.py quadratic_form__mass__Conway_Sloane_masses.py quadratic_form__mass__Siegel_densities.py quadratic_form__neighbors.py quadratic_form__reduction_theory.py quadratic_form__siegel_product.py quadratic_form__split_local_covering.py quadratic_form__ternary_Tornaria.py quadratic_form__theta.py quadratic_form__variable_substitutions.py random_quadraticform.py special_values.py }}} {{{id=376| s = QuadraticForm(ZZ, 4, [1..10]); s /// Quadratic form in 4 variables over Integer Ring with coefficients: [ 1 2 3 4 ] [ * 5 6 7 ] [ * * 8 9 ] [ * * * 10 ] }}} {{{id=378| factor(s.disc()) /// 3 * 11 * 53 }}} {{{id=380| s.signature() /// 4 }}} {{{id=379| /// }}} {{{id=371| /// }}}

Matrix

The matrix module provides optimized sparse and dense matrices over exact and numerical base rings.  It contains much code by Robert Bradshaw and I, which was later greatly improved by Jason Grout, Josh Kantor, Rob Beezer, and others.

The matrix module is fairly mature.  Some new features in the last two years:

PROJECTS:

  1. Add new functionality in a wider range of cases, e.g., LU decomposition.
  2. Arbitrary precision LAPACK: Optimized high-precision numerical linear algebra -- we have high precision numerical linear algebra, but it uses completely generic algorithms, so is potentially very slow. 
  3. Optimize sparse linear algebra.
{{{id=265| ls('matrix') /// __init__.py matrix_integer_dense.pyx action.c matrix_integer_dense_hnf.py action.pxd matrix_integer_dense_saturation.py action.pyx matrix_integer_sparse.c all.py matrix_integer_sparse.pxd benchmark.py matrix_integer_sparse.pyx berlekamp_massey.py matrix_misc.py change_ring.c matrix_mod2_dense.c change_ring.pyx matrix_mod2_dense.pxd constructor.py matrix_mod2_dense.pyx docs.py matrix_modn_dense.c matrix.c matrix_modn_dense.pxd matrix.pxd matrix_modn_dense.pyx matrix.pyx matrix_modn_sparse.c matrix0.c matrix_modn_sparse.pxd matrix0.pxd matrix_modn_sparse.pyx matrix0.pyx matrix_mpolynomial_dense.cpp matrix1.c matrix_mpolynomial_dense.pxd matrix1.pxd matrix_mpolynomial_dense.pyx matrix1.pyx matrix_rational_dense.c matrix2.c matrix_rational_dense.pxd matrix2.pxd matrix_rational_dense.pyx matrix2.pyx matrix_rational_dense.pyx.orig matrix_complex_double_dense.c matrix_rational_sparse.c matrix_complex_double_dense.pxd matrix_rational_sparse.pxd matrix_complex_double_dense.pyx matrix_rational_sparse.pyx matrix_cyclo_dense.cpp matrix_real_double_dense.c matrix_cyclo_dense.pxd matrix_real_double_dense.pxd matrix_cyclo_dense.pyx matrix_real_double_dense.pyx matrix_dense.c matrix_space.py matrix_dense.pxd matrix_sparse.c matrix_dense.pyx matrix_sparse.pxd matrix_domain_dense.pxd matrix_sparse.pyx matrix_domain_dense.pyx matrix_symbolic_dense.c matrix_domain_sparse.pxd matrix_symbolic_dense.pxd matrix_domain_sparse.pyx matrix_symbolic_dense.pyx matrix_double_dense.c matrix_window.c matrix_double_dense.pxd matrix_window.pxd matrix_double_dense.pyx matrix_window.pyx matrix_generic_dense.c matrix_window_modn_dense.c matrix_generic_dense.pxd matrix_window_modn_dense.pxd matrix_generic_dense.pyx matrix_window_modn_dense.pyx matrix_generic_sparse.c misc.c matrix_generic_sparse.pxd misc.pyx matrix_generic_sparse.pyx padics matrix_integer_2x2.c strassen.c matrix_integer_2x2.pxd strassen.pyx matrix_integer_2x2.pyx symplectic_basis.py matrix_integer_dense.c template.pxd matrix_integer_dense.pxd tests.py }}} {{{id=345| A = random_matrix(ZZ,4); A /// [-1 -1 5 8] [-4 2 0 4] [ 0 0 0 -1] [ 2 1 6 0] }}} {{{id=223| A.hermite_form() /// [ 1 0 11 0] [ 0 1 60 0] [ 0 0 76 0] [ 0 0 0 1] }}} {{{id=346| A.det() /// -76 }}} {{{id=347| A = random_matrix(ZZ,500) time A.det() /// -448443024043381226917541943805423596583910939202237039654563385356224926721064839131254043268871001051590216113926485961496449908943073211017899233945179059993592738263782803938034257230613839005118770516758031996463030017000250861144331392077975596505551629794019734301908054013201907719854147738311382789573167119088253929605320555898618452194115301568218416789142691203687887679466662624480800908710720994067838169671700629380272751942663883724413150746360002976072950851528518118313998787194581792360365716135806985530238301835159253505987777828066846714951224785904107728561841821790157068451259646590328509126638514531426420741167765245890406517146754422307981135560845822581919422536042428691845630837196677237112704884356023878932575395554408423269144955308566235745524691016306064352372529341455423915447638405388211972862716372459557393587162670875028746749965820557123191285310795050077515228894351174256684234941106265979408502311840497226316909420947118244723240284530070426818553475963412591321529636507311167494971999225174284283605671199436780136439737896135815318505321511907312194446121916930891509888889668652451509953933168678835419075692841071934245534876255209046827647319948628267065822893761021963463777399888304709710113944941134995315299832273837948855988907611711456708521527004451737482586813671100556322290683626802078021996754891508 Time: CPU 3.09 s, Wall: 2.96 s }}} {{{id=348| A = random_matrix(RDF,1000) time B = A*A /// Time: CPU 0.31 s, Wall: 0.18 s }}} {{{id=49| /// }}}

Games

The games module:
Sage also has extensive code for working with Rubik's cube (4 different solvers, cool graphics, etc.), but it is in the groups module.

PROJECT: Rob Beezer just posted a much more sophisticated Sudoku solver, which will replace mine.

PROJECT: The games module should grow to include massive amounts of code for solving a wide range of puzzles.   Kiran Kedlaya volunteered a bunch of Python code that he used for the MIT Mystery Hunt.  When Sage gets a SAT solver (say minisat), then it can easily be used to solve lots of games and puzzles.
{{{id=254| ls('games') /// __init__.py all.py sudoku.py }}} {{{id=210| A = matrix(ZZ,9,[5,0,0, 0,8,0, 0,4,9, 0,0,0, 5,0,0, 0,3,0, 0,6,7, 3,0,0, 0,0,1, 1,5,0, 0,0,0, 0,0,0, 0,0,0, 2,0,8, 0,0,0, 0,0,0, 0,0,0, 0,1,8, 7,0,0, 0,0,4, 1,5,0, 0,3,0, 0,0,2, 0,0,0, 4,9,0, 0,5,0, 0,0,3]) /// }}} {{{id=211| A /// [5 0 0 0 8 0 0 4 9] [0 0 0 5 0 0 0 3 0] [0 6 7 3 0 0 0 0 1] [1 5 0 0 0 0 0 0 0] [0 0 0 2 0 8 0 0 0] [0 0 0 0 0 0 0 1 8] [7 0 0 0 0 4 1 5 0] [0 3 0 0 0 2 0 0 0] [4 9 0 0 5 0 0 0 3] }}} {{{id=161| sudoku(A) /// [5 1 3 6 8 7 2 4 9] [8 4 9 5 2 1 6 3 7] [2 6 7 3 4 9 5 8 1] [1 5 8 4 6 3 9 7 2] [9 7 4 2 1 8 3 6 5] [3 2 6 7 9 5 4 1 8] [7 8 2 9 3 4 1 5 6] [6 3 5 1 7 2 8 9 4] [4 9 1 8 5 6 7 2 3] }}} {{{id=212| matrix_plot(sudoku(A),cmap='Pastel1') /// }}} {{{id=214| /// }}} {{{id=29| /// }}}

Logic

The formal logic module implements formal logic (e.g., truth tables, logical formulas) mainly for educational applications.  (Written by Chris Gorecki.)

PROJECT: Extend so usable for research applications by wrapping MINISAT.

{{{id=264| ls('logic') /// __init__.py booleval.py logic.py logictable.py all.py boolformula.py logicparser.py propcalc.py }}} {{{id=222| import sage.logic.propcalc as propcalc s = propcalc.formula("a&b|~(c|a)") s.truthtable() /// a b c value False False False True False False True False False True False True False True True False True False False False True False True False True True False True True True True True }}} {{{id=300| /// }}} {{{id=47| /// }}}

Homology

The new homology module by algebraic topologist John Palmieri computes homology of simplicial complexes.

PROJECT: Several topologists attended Sage Days 15, and this module is now being extremely actively developed.  Sage will soon include a very wide range of powerful code for algebraic topology computations.

{{{id=260| ls('homology') /// __init__.py chain_complex.py simplicial_complex.py all.py examples.py }}} {{{id=218| X = SimplicialComplex(3, [[0,1], [1,2], [2,3], [3,0]]); X /// Simplicial complex with vertex set (0, 1, 2, 3) and facets {(1, 2), (2, 3), (0, 3), (0, 1)} }}} {{{id=353| X.stanley_reisner_ring() /// Quotient of Multivariate Polynomial Ring in x0, x1, x2, x3 over Integer Ring by the ideal (x1*x3, x0*x2) }}} {{{id=355| X.cohomology() /// {0: 0, 1: Z} }}} {{{id=354| S = SimplicialComplex(3, [[0,1], [1,2], [0,2]]) # circle T = S.product(S) # torus T /// Simplicial complex with 16 vertices and 18 facets }}} {{{id=352| T.cohomology() /// {0: 0, 1: Z x Z, 2: Z} }}} {{{id=402| T.graph().plot(layout='circular') /// }}} {{{id=368| /// }}} {{{id=39| /// }}}

Geometry

The geometry module currently contains code for working mainly with polytopes.  It's especially strong for working with polytopes with integers vertices (lattice polytopes).

PROJECT: Marshall Hampton is mainly working on this module right now.

{{{id=255| ls('geometry') /// __init__.py lattice_polytope.py polytope.py all.py polyhedra.py }}} {{{id=215| m = matrix(ZZ, [[1, 1, 2, -1, 3, 0], [0, 5, 0, 0, -1, 0], [0, 0, 1, 0, 0, -7]]) p = LatticePolytope(m); p /// A lattice polytope: 3-dimensional, 5 vertices. }}} {{{id=216| p.plot3d() /// }}} {{{id=31| /// }}}

Schemes

The schemes directory is divided up into subdirectories:

 

History:

  1. David Kohel and I wrote the very first version of the schemes code in 2005. 
  2. Vast amount of functionality added to the elliptic and hyperelliptic curve directories by David Harvey, John Cremona, Robert Bradshaw, David Roe, Christian Wuthrich, and many other people. 
  3. Alex Ghitza is currently working on code not in the elliptic or hyperelliptic curve directories. 
  4. David Harvey designed and implemented fast algorithm for computation of Frobenius on hyperelliptic curves, which is relevant to cryptography and number theory.

PROJECT: The goal for general schemes is of course to be able to tick off as many items as possible on the list at http://wiki.sagemath.org/AlgebraicGeometrySEP . Alex Ghitza plan:

  1. bring doctest coverage to 100 percent
  2. use the combinatorics code to implement some enumerative geometry, which people now do in Maple
  3. use the new algebraic topology code and your code for finitely generated modules to implement sheaf cohomology
  4. use the polytopes code to implement toric varieties, which will give good examples on which to doctest more advanced functionality

PROJECT: John Cremona is working on $S$-integral points over number fields.

PROJECT: Robert Bradshaw and William Stein working on a general Heegner points package.

PROJECT: Dan Shumow currently implementing everything he can related to isogenies of elliptic curves.

 

{{{id=284| ls('schemes') /// __init__.py generic plane_curves all.py hyperelliptic_curves plane_quartics elliptic_curves jacobians readme.py }}} {{{id=286| E = EllipticCurve('37a'); show(E) ///
y^2 + y = x^3 - x
}}} {{{id=310| plot(E,thickness=3).show(figsize=[6,3]) /// }}} {{{id=404| S = E.sha(); S /// Shafarevich-Tate group for the Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field }}} {{{id=406| S.an() /// 1 }}} {{{id=407| S.bound_kolyvagin() /// ([2], 1) }}} {{{id=408| S.two_selmer_bound() /// 0 }}} {{{id=309| /// }}}

Some general schemes

{{{id=287| Spec(ZZ) /// Spectrum of Integer Ring }}} {{{id=288| PP. = ProjectiveSpace(QQ,2); PP /// Projective Space of dimension 2 over Rational Field }}} {{{id=289| Z = PP.subscheme([y^2*z + y*z^2 - (x^3 - x*z^2), x-y]); Z /// Closed subscheme of Projective Space of dimension 2 over Rational Field defined by: -x^3 + y^2*z + x*z^2 + y*z^2 x - y }}} {{{id=283| Z.irreducible_components() /// [ Closed subscheme of Projective Space of dimension 2 over Rational Field defined by: y x, Closed subscheme of Projective Space of dimension 2 over Rational Field defined by: y - 2*z x - 2*z, Closed subscheme of Projective Space of dimension 2 over Rational Field defined by: y + z x + z ] }}} {{{id=290| /// }}}

David Harvey's Hyperelliptic Jacobian Point Counting

{{{id=414| from sage.schemes.hyperelliptic_curves.hypellfrob import hypellfrob R. = PolynomialRing(ZZ) f = x^5 + 2*x^2 + x + 1; p = 2003 time M = hypellfrob(p, 4, f); M /// Time: CPU 0.36 s, Wall: 0.37 s [ 4169095036140 + O(2003^4) 9084480306994 + O(2003^4) 6847436440991 + O(2003^4) 13920997340731 + O(2003^4)] [ 7328714194382 + O(2003^4) 2587287412441 + O(2003^4) 4450185817450 + O(2003^4) 3312129622624 + O(2003^4)] [ 4058918492851 + O(2003^4) 12645738016730 + O(2003^4) 7269250999683 + O(2003^4) 3307639750279 + O(2003^4)] [ 2811015576052 + O(2003^4) 4803825776245 + O(2003^4) 10902661795104 + O(2003^4) 2070582767902 + O(2003^4)] }}} {{{id=415| -M.trace() /// 16096216215996 + O(2003^4) }}} {{{id=416| ZZ(M.det()) /// 4012009 }}} {{{id=412| /// }}} {{{id=411| /// }}} {{{id=410| /// }}} {{{id=409| /// }}} {{{id=73| /// }}}