\documentclass[runningheads,a4paper]{llncs}1\usepackage{ucs}2\usepackage[utf8x,utf8]{inputenc}34\usepackage{amsmath,amssymb}5\setcounter{tocdepth}{3}6\usepackage[pdftex]{graphicx}78\usepackage{url}9\usepackage{hyperref}1011\newcommand{\todo}[1]{[[\rule{0.5em}{1ex}\hspace{1em}#1]]}12\newcommand{\todobf}[1]{ { \rule{0.5em}{1ex}\textbf{ #1}}}1314\newcommand{\goth}{\mathfrak}1516\usepackage{fancybox}171819\def\C{\mathbb C}20\def\Q{\mathbb Q}21\def\Z{\mathbb Z}22\def\R{\mathbb R}232425\usepackage{color}26\usepackage{listings}27\usepackage{xspace} % to allow for text macros that don't eat space2829\lstdefinelanguage{Sage}[]{Python}30{morekeywords={True,False,sage,singular},31sensitive=true}32\lstset{frame=none,33showtabs=False,34showspaces=False,35showstringspaces=False,36commentstyle={\ttfamily\color{dredcolor}},37keywordstyle={\ttfamily\color{dbluecolor}\bfseries},38stringstyle ={\ttfamily\color{dgraycolor}\bfseries},39language = Sage,40basicstyle={\small \ttfamily},41aboveskip=.3em,42belowskip=.1em43}44\definecolor{dblackcolor}{rgb}{0.0,0.0,0.0}45\definecolor{dbluecolor}{rgb}{.01,.02,0.7}46\definecolor{dredcolor}{rgb}{0.8,0,0}47\definecolor{dgraycolor}{rgb}{0.30,0.3,0.30}48\newcommand{\dblue}{\color{dbluecolor}\bf}49\newcommand{\dred}{\color{dredcolor}\bf}50\newcommand{\dblack}{\color{dblackcolor}\bf}5152\renewcommand{\emph}[1]{{\dblack{#1}}}5354\newcommand{\Sage}{{\color{dbluecolor}\sf Sage}\xspace}555657\newcommand{\keywords}[1]{\par\addvspace\baselineskip58\noindent\keywordname\enspace\ignorespaces#1}5960\begin{document}6162\mainmatter % start of an individual contribution6364% first the title is needed65\title{The Sage Project:\\66\large Unifying Free Mathematical Software to Create a Viable Alternative to67Magma, Maple, Mathematica and MATLAB}6869% a short form should be given in case it is too long for the running head70\titlerunning{Sage: Unifying Free Mathematical Software}7172% the name(s) of the author(s) follow(s) next73%74\author{Burçin Eröcal\inst{1}%75%\thanks{Please note that the LNCS Editorial assumes that all authors have used76%the western naming convention, with given names preceding surnames. This determines77%the structure of the names in the running heads and the author index.}%78\and William Stein\inst{2}}79%80%\authorrunning{}8182% the affiliations are given next; don't give your e-mail address83% unless you accept that it will be published84\institute{Research Institute for Symbolic Computation\\85Johannes Kepler University,\\86Linz, Austria\\87Supported by FWF grants P20347 and DK W1214.\\88\email{burcin@erocal.org}89\and90Department of Mathematics\\91University of Washington\\92\email{wstein@uw.edu}\\93Supported by NSF grant DMS-0757627 and DMS-0555776.}9495%96% NB: a more complex sample for affiliations and the mapping to the97% corresponding authors can be found in the file "llncs.dem"98% (search for the string "\mainmatter" where a contribution starts).99% "llncs.dem" accompanies the document class "llncs.cls".100%101102%\toctitle{Lecture Notes in Computer Science}103%\tocauthor{Authors' Instructions}104\maketitle105106107%The abstract should summarize the contents of the paper and should108%contain at least 70 and at most 150 words. It should be written using the109%\emph{abstract} environment.110\begin{abstract}111Sage is a free, open source,112self-contained distribution of mathematical software, including a113large library that provides a unified interface to the components of114this distribution. This library also builds on the components of115Sage to implement novel algorithms covering a broad range of116mathematical functionality from algebraic combinatorics to number117theory and arithmetic geometry.118119\keywords{Python, Cython, Sage, Open Source, Interfaces}120\end{abstract}121122%\todo{Building the Car - internals, design decisions etc. this section can be skipped before looking at the examples in the next section, but the examples might refer to it}123124%\todo{Interesting problems - (Test drive) - showcasing the capabilities (of the car)}125126127\section{Introduction}128129In order to use mathematical software for exploration, we often push130the boundaries of available computing resources and continuously try to131improve our implementations and algorithms. Most mathematical132algorithms require basic building blocks, such as multiprecision133numbers, fast polynomial arithmetic, exact or numeric linear algebra,134or more advanced algorithms such as Gr\"obner basis computation or135integer factorization. Though implementing some of these basic136foundations from scratch can be a good exercise, the resulting code137may be slow and buggy. Instead, one can build on existing optimized138implementations of these basic components, either by using a general139computer algebra system, such as Magma, Maple, Mathematica or MATLAB,140or by making use of the many high quality open source libraries that141provide the desired functionality. These two approaches both have142significant drawbacks. This paper is about143Sage,\footnote{\url{http://www.sagemath.org}} which provides an144alternative approach to this problem.145146Having to rely on a closed propriety system can be frustrating, since147it is difficult to gain access to the source code of the software,148either to correct a bug or include a simple optimization in an149algorithm. Sometimes this is by design:150\begin{quote}\small ``Indeed, in almost all practical uses of Mathematica,151issues about how Mathematica works inside turn out to be largely152irrelevant. You might think that knowing how Mathematica153works inside would be necessary [...]'' (See \cite{mathematica:internals}.)154\end{quote}155Even if we manage to contact the developers, and they find time to156make the changes we request, it might still take months or years157before these changes are made available in a new release.158159Fundamental questions of correctness, reproducibility and scientific160value arise when building a mathematical research program on top of161proprietary software (see, e.g., \cite{joyner-stein:notices}). There162are many published refereed papers containing results that163rely on computations performed in Magma, Maple, or164Mathematica.\footnote{Including by the second author of this paper,165e.g., \cite{conrad-edixhoven-stein:j1p}!} In some cases, a specific166version of Magma is the only software that can carry out the167computation. This is not the infrastructure on which we want to build168the future of mathematical research.169170In sharp contrast, open source libraries provide a great deal of171flexibility, since anyone can see and modify the source code as they172wish. However, functionality is often segmented into different173specialized libraries and advanced algorithms are hidden behind custom174interpreted languages. One often runs into trouble trying to install175dependencies before being able use an open source software package.176Also, converting the output of one package to the input format of177another package can present numerous difficulties and178introduce subtle errors.179180Sage, which started in 2005 (see \cite{joyner-stein:sigsam}),181attacks this problem by providing:182183\begin{enumerate}184\item a {\em self-contained distribution} of mathematical software that185installs from source easily, with the only dependency being compiler186tools,187\item {\em unified interfaces} to other mathematical software to make it188easier to use all these programs together, and189\item a {\em new library} that builds on the included software packages and190implements a broad range of mathematical functionality.191\end{enumerate}192193%\begin{center}194%\shadowbox{%195%\begin{minipage}{0.85\textwidth}196%{\bf The Sage Mission Statement:} Create a viable free open source %alternative197%to Magma, Maple, Mathematica and Matlab.198%\end{minipage}}199%\end{center}200201202203204%William Stein started the Sage project in early 2005, with financial205%support from the NSF, Google, Microsoft, University of Washington, and206%many other organizations and individuals.207208%Sage\footnote{\url{http://www.sagemath.org}} is a free, open source, self-contained209%distribution of mathematical software, including a large library210%that provides a unified interface to the components of this distribution. This library also builds upon211%the components of Sage to implement novel algorithms covering a broad212%range of mathematical functionality from algebraic combinatorics to213%number theory and arithmetic geometry.214215216The rest of this paper goes into more detail about Sage. In217Section~\ref{sec:notebook}, we describe the Sage graphical user218interface. Section~\ref{sec:dev} is about the Sage development219process, Sage days workshops, mailing lists,220and documentation. The subject of Section~\ref{sec:car} is the221sophisticated way in which Sage is built out of a wide range of open222source libraries and software. In Section~\ref{sec:interfaces} we223explain how we use Python and Cython as the glue that224binds the compendium of software included in Sage into a unified225whole. We then delve deeper into Python, Cython and the Sage226preparser in Section~\ref{sec:python}, and illustrate some227applications to mathematics in Section~\ref{sec:algebraic}. Sage is228actively used for research, and in229Section~\ref{sec:interesting} we describe some230capabilities of Sage in advanced areas of mathematics.231232\subsection{The Notebook}\label{sec:notebook}233\begin{figure}234\begin{center}235\shadowbox{\includegraphics[width=.8\textwidth]{fig/sagenb}}236\caption{The Sage Notebook\label{fig:sagenb}}237\end{center}238\end{figure}239As illustrated in Figure~\ref{fig:sagenb}, the graphical user240interface for Sage is a web application, inspired by Google Documents241\cite{googledocs}, which provides convenient access to all242capabilities of Sage, including 3D graphics. In single user mode,243Sage works like a regular application whose main window happens to be244your web browser. In multiuser mode, this architecture allows users to245easily set up servers for accessing their work over the Internet as246well as sharing and collaborating with colleagues. One can try the247Sage notebook by visiting \url{www.sagenb.org}, where there are over24830,000 user accounts and over 2,000 published worksheets.249250Users also download Sage to run it directly on their computers. We251track all downloads from \url{www.sagemath.org}, though there252are several other high-profile sites that provide mirrors of our253binaries. Recently, people download about 6,000 copies of Sage per254month directly from the Sage website.255256257\subsection{The Sage Development Process}\label{sec:dev}258259There are over 200 developers from across the world who have260contributed to the Sage project. People often contribute because they261write code using Sage as part of a research project, and in this262process find and fix bugs, speed up parts of Sage, or want the code263portion of their research to be peer reviewed. Each contribution to264Sage is first posted to the Sage Trac server265\url{trac.sagemath.org}; it is then peer reviewed, and finally266added to Sage after all issues have been sorted out and all267requirements are met. Nothing about this process is anonymous; every268step of the peer review process is recorded indefinitely for all to269see.270271The Sage Developer's Guide begins with an easy-to-follow tutorial that272guides developers through each step involved in contributing code to Sage.273Swift feedback is available through the \verb|sage-devel| mailing274list, and the \verb|#sage-devel| IRC chat room on275\verb|irc.freenode.net| (see276\url{www.sagemath.org/development.html}).277278%(#sage-devel on irc.freenode.net).279280281Much development of Sage has taken place at the Sage Days282workshops. There have been two dozen Sage Days283\cite{sagedays} and many more are planned. These are essential to284sustaining the momentum of the Sage project and also help ensure that285developers work together toward a common goal, rather than competing286with each other and fragmenting our limited community resources.287288A major goal is ensuring that there will be many Sage Days workshops289for the next couple of years. The topics will depend on funding, but290will likely include numerical computation, large-scale bug fixing,291$L$-functions and modular forms, function fields, symbolic292computation, topology, and combinatorics. The combination of293experienced developers with a group of enthusiastic mathematicians at294each of these workshops has rapidly increased the developer community,295and we hope that it will continue to do so.296297%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%298299300\section{Building the Car\dots}\label{sec:car}301\begin{center}302\includegraphics[height=3cm]{fig/sagecar_crop}303\end{center}304305306With the motto ``building the car instead of reinventing the wheel,''307Sage brings together numerous open source software308packages (see Table~\ref{tab:standard} and \cite{components}).309310\begin{table}311\begin{center}312\caption{Packages Included With Every Copy of Sage-4.4.2\label{tab:standard}}313\begin{tabular}{|l|l|l|l|l|}\hline314atlas& gap& libgcrypt& palp& scipy\_sandbox\\\hline315blas& gd& libgpg\_error& pari& scons\\\hline316boehm\_gc& gdmodule& libm4ri& pexpect& setuptools\\\hline317boost& genus2reduction& libpng& pil& singular\\\hline318cddlib& gfan& linbox& polybori& sphinx\\\hline319cliquer& ghmm& matplotlib& pycrypto& sqlalchemy\\\hline320cvxopt& givaro& maxima& pygments& sqlite\\\hline321cython& gnutls& mercurial& pynac& symmetrica\\\hline322docutils& gsl& moin& python& sympow\\\hline323ecl& iconv& mpfi& python\_gnutls& sympy\\\hline324eclib& iml& mpfr& r& tachyon\\\hline325ecm& ipython& mpir& ratpoints& termcap\\\hline326f2c& jinja& mpmath& readline& twisted\\\hline327flint& jinja2& networkx& rubiks& weave\\\hline328flintqs& lapack& ntl& sagenb& zlib\\\hline329fortran& lcalc& numpy& sagetex& zn\_poly\\\hline330freetype& libfplll& opencdk& scipy& zodb3\\\hline331\end{tabular}332\end{center}333\end{table}334335Many applications of Sage require using these336libraries together. Sage handles the conversion of data behind the337scenes, automatically using the best tool for the job, and allows338the user to concentrate on the problem at hand.339340In the following example, which we explain in detail below, Sage341uses the FLINT library \cite{flint} for univariate342polynomials over the ring $\Z$ of integers, whereas Singular343\cite{singular} is used for multivariate344polynomials. The option to use the NTL library \cite{ntl} for345univariate polynomials is still available, if the user so chooses.346347\begin{lstlisting}[numbers=left,numberstyle=\tiny,numbersep=0pt]348sage: R.<x> = ZZ[]349sage: type(R.an_element())350<type 'sage.rings...Polynomial_integer_dense_flint'>351sage: R.<x,y> = ZZ[]352sage: type(R.an_element())353<type 'sage.rings...MPolynomial_libsingular'>354sage: R = PolynomialRing(ZZ, 'x', implementation='NTL')355sage: type(R.an_element())356<type 'sage.rings...Polynomial_integer_dense_ntl'>357\end{lstlisting}358359The first line in the example above constructs the univariate polynomial ring360${\tt R} = \Z[x]$, and assigns the variable \verb|x| to be the generator of361this ring. Note that $\Z$ is represented by \verb|ZZ| in Sage. The362expression \verb|R.<x> = ZZ[]| is not valid Python, but can be used363in Sage code as a shorthand as explained in Section \ref{sec:preparser}. The364next line asks the ring \verb|R| for an element, using the365\verb|an_element| function, then uses the builtin Python function366\verb|type| to query its type. We learn that it is an instance of the class367\verb|Polynomial_integer_dense_flint|. Similarly line 4 constructs ${\tt R}368= \Z[x,y]$ and line 7 defines ${\tt R} = \Z[x]$, but this time using the369\verb|PolynomialRing| constructor explicitly and specifying that we want370the underlying implementation to use the NTL library.371372Often these interfaces are used under the hood, without the user373having to know anything about the corresponding systems. Nonetheless, there are easy374ways to find out what is used by inspecting the source code, and users375are strongly encouraged to cite components they use in published376papers. The following example illustrates another377way to get a list of378components used when a specific command is run.379380\begin{lstlisting}381sage: from sage.misc.citation import get_systems382sage: get_systems('integrate(x^2, x)')383['ginac', 'Maxima']384sage: R.<x,y,z> = QQ[]385sage: I = R.ideal(x^2+y^2, z^2+y)386sage: get_systems('I.primary_decomposition()')387['Singular']388\end{lstlisting}389390\subsection{Interfaces}\label{sec:interfaces}391392Sage makes it possible to use a wide range of mathematical software393packages together by providing a unified interface that394handles data conversion automatically. The complexity and395functionality of these interfaces varies greatly, from simple396text-based interfaces that call external software for an individual397computation, to using a library as the basis for an398arithmetic type. The interfaces can also run code from libraries399written in the interpreted language of another program.400Table~\ref{fig:interfaces} lists401the interfaces provided by Sage.402403\begin{table}404\begin{center}405\caption{Sage Interfaces to the above Mathematical Software\label{fig:interfaces}}406407\begin{tabular}{|l|l|}\hline408{\bf Pexpect} & axiom, ecm, fricas, frobby, gap, g2red, gfan,409gnuplot, gp, \\410& kash, lie, lisp, macaulay2, magma, maple, mathematica,411\\412& matlab, maxima, mupad, mwrank, octave, phc, polymake, \\413&povray, qepcad, qsieve, r,414rubik, scilab, singular, tachyon\\\hline415{\bf C Library} & eclib, fplll, gap (in progress), iml, linbox, maxima,\\416& ratpoints, r (via rpy2), singular, symmetrica\\\hline417{\bf C Library arithmetic} & flint, mpir, ntl, pari, polybori, pynac, singular\\\hline418\end{tabular}419\end{center}420\end{table}421422423The above interfaces are the result of many years writing424Python and Cython \cite{cython} code to adapt425Singular \cite{singular}, GAP \cite{gap}, Maxima \cite{maxima}, Pari426\cite{pari}, GiNaC/Pynac \cite{ginac}, NTL \cite{ntl}, FLINT427\cite{flint}, and many other libraries, so that they can be used428smoothly and efficiently in a unified way from Python \cite{python}.429Some of these programs were originally designed to be used only430through their own interpreter and made into a library by Sage431developers. For example libSingular was created by Martin Albrecht in432order to use the fast multivariate polynomial arithmetic in Singular433from Sage. The libSingular interface is now used by other434projects, including Macaulay2 \cite{M2} and GFan \cite{gfan}.435436There are other approaches to linking mathematical software together.437The recent paper \cite{science_openmath} reports on the state of the438art using OpenMath. Sage takes a439dramatically different approach to this problem. Instead of using a440general string-based XML protocol to communicate with other mathematical441software, Sage interfaces are tailor made to the specific software and442problem at hand. This results in far more efficient and flexible443interfaces. The main disadvantage compared to OpenMath is that the444interfaces all go through Sage.445446%Sage also allows easy use of commercial packages, by encapsulating the447%knowledge\todo{sounds funny} to transfer data to a format these packages can448%understand.449450Having access to many programs which can perform the same computation,451without having to worry about data conversion, also makes it easier to452double check results. For example, below we first use Maxima, an open453source symbolic computation package distributed with Sage, to454integrate a function, then perform the same computation using Maple455and Mathematica.456457\begin{lstlisting}458sage: var('x')459sage: integrate(sin(x^2), x)4601/8*((I - 1)*sqrt(2)*erf((1/2*I - 1/2)*sqrt(2)*x) + \461(I + 1)*sqrt(2)*erf((1/2*I + 1/2)*sqrt(2)*x))*sqrt(pi)462sage: maple(sin(x^2)).integrate(x)4631/2*2^(1/2)*Pi^(1/2)*FresnelS(2^(1/2)/Pi^(1/2)*x)464sage: mathematica(sin(x^2)).Integrate(x)465Sqrt[Pi/2]*FresnelS[Sqrt[2/Pi]*x]466\end{lstlisting}467468469470471472The most common type of interface, called a \verb|pexpect|473interface, communicates with another command line program by reading474and writing strings to a text console, as if another user was in front475of the terminal. Even though these are relatively simple to develop,476the overhead of having to print and parse strings to represent the477data makes this process potentially cumbersome and inefficient. This478is the default method of communication with most high level479mathematics software, including commercial and open source programs,480such as Maple, Mathematica, Magma, KASH or GAP.481482Sage provides a framework to represent elements over these interfaces, perform483arithmetic with them or apply functions to the given object, as well as using a484file to pass the data if the string representation is too big. The following485demonstrates arithmetic with GAP elements.486487\begin{lstlisting}488sage: a = gap('22')489sage: a*a490484491\end{lstlisting}492493It is also possible to use \verb|pexpect| interfaces over remote consoles. In494the following code, we connect to the \verb|localhost| as a different user495and call Mathematica functions. Note that the interface can handle indexing496vectors as well.497498\begin{lstlisting}499sage: mma = Mathematica(server="rmma60@localhost")500sage: mma("2+2")5014502sage: t = mma("Cos[x]")503sage: t.Integrate('x')504Sin[x]505sage: t = mma('{0,1,2,3}')506sage: t[2]5071508\end{lstlisting}509510Sage also includes specialized libraries that are linked directly from compiled511code written in Cython. These are used to handle specific problems, such as the512characteristic polynomial computation in the example below.513514\begin{lstlisting}515sage: M = Matrix(GF(5), 10, 10)516sage: M.randomize()517sage: M.charpoly(algorithm='linbox')518x^10 + 4*x^9 + 4*x^7 + 3*x^4 + 3*x^3 + 3*x^2 + 4*x + 3519\end{lstlisting}520521Many basic arithmetic types also use Cython to directly utilize data522structures from efficient arithmetic libraries, such as MPIR or FLINT. An523example of this can be seen at the beginning of this section, where elements524of the ring $\Z[x]$ are represented by the class525\verb|Polynomial_integer_dense_flint|.526527The Singular interface is one of the most advanced included in528Sage. Singular has a large library of code written in its own language.529Previously the only way to access these functions, which include530algorithms for Gr\"obner basis and primary decomposition, was to call531Singular through a \verb|pexpect| interface, passing data back and forth532using strings. Recently, due to work of Michael Brickenstein and Martin533Albrecht, Sage acquired the ability to call these functions directly.534535In the example below, we import the function \verb|primdecSY| from536\verb|primdec.lib|, and call it the same way we would call a Python537function. The interface handles the conversion of the data to Singular's format538and back. Since Sage already uses Singular data structures directly to represent539multivariate polynomials and ideals over multivariate polynomial rings, there540are no conversion costs. It is only a matter of passing the right pointer.541542\begin{lstlisting}543sage: pr = sage.libs.singular.ff.primdec__lib.primdecSY544sage: R.<x,y,z> = QQ[]545sage: p = z^2+1; q = z^3+2546sage: I = R.ideal([p*q^2,y-z^2])547sage: pr(I)548[[[z^2 - y, y^3 + 4*y*z + 4], \549[z^2 - y, y*z + 2, y^2 + 2*z]], \550[[y + 1, z^2 + 1], [y + 1, z^2 + 1]]]551\end{lstlisting}552553Efforts are under way to extend these capabilities to other programs,554for example to GAP which provides Sage's underlying group theory555functionality. Up to now, GAP was only available through its interpreter,556through a \verb|pexpect| interface that was written by Steve Linton. As the557following example demonstrates, the performance of this interface is558far from ideal.\footnote{All timings in this paper were performed on559an 2.66GHz Intel Xeon X7460 based computer.}560561\begin{lstlisting}562sage: b = gap('10')563sage: b*b564100565sage: timeit('b*b')566625 loops, best of 3: 289 microseconds per loop567\end{lstlisting}568569The code snippet above constructs the element \verb|b| in GAP using the570\verb|pexpect| interface, and measures the time it takes to square571\verb|b|. Compare these numbers to the following example, which uses the572library interface to GAP, recently developed by the second author573(but {\em not} included in Sage yet).574575\begin{lstlisting}576sage: import sage.libs.gap.gap as g577sage: a = g.libgap('10'); a57810579sage: type(a)580<type 'sage.libs.gap.gap.GapElement'>581sage: a*a582100583sage: timeit('a*a')584625 loops, best of 3: 229 nanoseconds per loop585\end{lstlisting}586The library interface is about 1,000 times faster than the pexpect interface.587588589\subsection{Python - a mainstream language}\label{sec:python}590591In line with the principle of not reinventing the wheel, Sage is built on the592mainstream programming language Python, both as the main development language593and the user language. This frees the Sage developers, who are mainly594mathematicians, from the troubles of language design, and gives access to595an immense array of general purpose Python libraries and tools.596597Python is an interpreted language with a clear, easy to read and learn598syntax. Since it is dynamically typed, it is ideal for rapid prototyping,599providing an environment to easily test new ideas and algorithms.600601%Note that many mathematical packages, such as GINV, PolyBoRi, or CVXOPT,602%already have a Python interface.603604\subsubsection{A fast interpreter}605606In the following Singular session, we first declare the ring ${\tt607r}=\Q[x,y,z]$ and the polynomial ${\tt f} \in {\tt r}$, then608measures the time to square ${\tt f}$ repeatedly, 10,000 times.609610\begin{lstlisting}611singular: int t = timer; ring r = 0,(x,y,z), dp;612singular: def f = y^2*z^2-x^2*y^3-x*z^3+x^3*y*z;613singular: int j; def g = f;614singular: for (j = 1; j <= 10^5; j++) { g = f*f; }615singular: (timer-t), system("--ticks-per-sec");616990 1000617\end{lstlisting}618619The elapsed time is 990 milliseconds. Next we use Sage to do the same620computation, using the same Singular data structures directly, but621without going through the interpreter.622623\begin{lstlisting}624sage: R.<x,y,z> = QQ[]625sage: f = y^2*z^2 - x^2*y^3 - x*z^3 + x^3*y*z; type(f)626<type 'sage.rings.polynomial...MPolynomial_libsingular'>627sage: timeit('for j in xrange(10^5): g = f*f')6285 loops, best of 3: 91.8 ms per loop629\end{lstlisting}630631Sage takes only 91.8 milliseconds for the same operation. This difference is632because the Python interpreter is more efficient at performing \verb|for| loops.633634635636\subsubsection{Cython - compiled extensions}637638Python alone is too slow to implement a serious mathematical software639system. Fortunately, Cython \cite{cython}640makes it easy to optimize parts of your641program or access existing C/C++ libraries. It can translate Python642code with annotations containing static type information to C/C++643code, which is then compiled as a Python extension module.644645Many of the basic arithmetic types in Sage are provided by Cython wrappers of646C libraries, such as FLINT for univariate polynomials over $\Z$, Singular for647multivariate polynomials, and Pynac for symbolic expressions.648649The code segment below defines a Python function to add integers from 0 to650$N$ and times the execution of this function with the argument \verb|10^7|.651652\begin{lstlisting}653sage: def mysum(N):654....: s = int(0)655....: for k in xrange(1,N): s += k656....: return s657....:658sage: time mysum(10^7)659CPU times: user 0.52 s, sys: 0.00 s, total: 0.52 s66049999995000000661\end{lstlisting}662663Here is the same function, but the loop index {\tt k} is declared to be664a C integer and the accumulator {\tt s} is a C {\tt long long}.665666\begin{lstlisting}667sage: cython("""668....: def mysum_cython(N):669....: cdef int k670....: cdef long long s = 0671....: for k in xrange(N): s += k672....: return s673....: """)674sage: time mysum_cython(10^7)675CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s67649999995000000L677\end{lstlisting}678679The code is compiled and linked to the interpreter on the fly, and the680function \verb|mysum_cython| is available immediately. Note that the run681time for the Cython function is 60 times faster than the Python equivalent.682683Cython also handles the conversion of Python types to C types automatically.684In the following example, we call the C function \verb|sinl| using Cython685to wrap it in a Python function named \verb|sin_c_wrap|.686687\begin{lstlisting}688sage: cython("""689....: cdef extern from "math.h":690....: long double sinl(long double)691....: def sin_c_wrap(a):692....: return sinl(a)693....: """)694sage: sin_c_wrap(3.14)6950.0015926529164868282696sage: sin_c_wrap(1)6970.8414709848078965698sage: sin_c_wrap(1r)6990.8414709848078965700\end{lstlisting}701702Note that the conversion of Sage types in the first two calls to703\verb|sin_c_wrap| or the Python type integer in the last call is performed704transparently by Cython.705706\subsubsection{The Preparser}\label{sec:preparser}707708While Python has many advantages as a programming and glue language, it also709has some undesirable features. Sage hides these problems by using a preparser710to change the commands passed to Python in an interactive session (or when running711a script with the \verb|.sage| extension). In order to712maintain compatibility with Python, changes performed by the preparser are713kept to a minimum. Moreover, the Sage library code is not preparsed, and is714written in Cython or Python directly.715716Python, like C and many other programming languages, performs integer floor717division. This means typing \verb|1/2| results in 0, not the rational number718$1/2$. Sage wraps all numeric literals entered in the command line or the719notebook with its own type declarations, which behave as expected with respect720to arithmetic and have the advantage that they are backed by efficient721multiprecision arithmetic libraries such as MPIR \cite{mpir} and MPFR \cite{mpfr}, which are thousands of times faster than Python for large722integer arithmetic.723724To call the preparser directly on a given string, use the725\verb|preparse| function.726727\begin{lstlisting}728sage: preparse("1/2")729'Integer(1)/Integer(2)'730sage: preparse("1.5")731"RealNumber('1.5')"732\end{lstlisting}733734Adding a trailing {\tt r} after a number indicates that the preparser should leave735that as the ``raw'' literal. The following illustrates division with Python integers.736737\begin{lstlisting}738sage: preparse("1r/2r")739'1/2'740sage: 1r/2r7410742\end{lstlisting}743Here is the result of performing the same division in Sage.744\begin{lstlisting}745sage: 1/27461/2747sage: type(1/2)748<type 'sage.rings.rational.Rational'>749sage: (1/2).parent()750Rational Field751\end{lstlisting}752753The preparser also changes the \verb|^| sign to the exponentiation operator754\verb|**| and provides a shorthand to create new mathematical domains and755name their generator in one command.756757\begin{lstlisting}758sage: preparse("2^3")759'Integer(2)**Integer(3)'760sage: preparse("R.<x,y> = ZZ[]")761"R = ZZ['x, y']; (x, y,) = R._first_ngens(2)"762\end{lstlisting}763764765\subsection{Algebraic, Symbolic and Numerical Tools}~\label{sec:algebraic}766767Sage combines algebraic, symbolic and numerical computation tools under one768roof, enabling users to choose the tool that best suits the problem. This769combination also makes Sage more accessible to a wide audience---scientists,770engineers, pure mathematicians and mathematics teachers771can all use the same platform for scientific computation.772773While not concentrating on only one of these domains might seem to divide774development resources unnecessarily, it actually results in a better overall775experience for everyone, since users do not have to come up with makeshift776solutions to compensate for the lack of functionality from a different field.777Moreover, because Sage is a distributed mostly-volunteer open source778project, widening our focus results in substantially more developer resources.779780\subsubsection{Algebraic Tools: The Coercion System}781782An algebraic framework, similar to that of Magma or Axiom, provides access to783efficient data structures and specialized algorithms associated to particular784mathematical domains. The Python language allows785classes to define how arithmetic786operations like \verb|+| and \verb|*| will be handled, in a787similar way to how C++ allows overloading of operators. However, the built-in support788for overloading in Python is too simple789to support operations with a range of objects in a mathematical type hierarchy.790791Sage abstracts the process of deciding what an arithmetic operation792means, or equivalently, in which domain the operation should be793performed, in a framework called the {\em coercion system}, which was794developed and implemented by Robert Bradshaw, David Roe,795and many others. Implementations of new mathematical objects only need796to define which other domains have a natural embedding to their797domain. When performing arithmetic with objects, the798coercion system will find a common domain where both arguments can be799canonically mapped, perform the necessary type conversions automatically,800thus allowing the implementation to only handle the case where both801objects have the same parent.802803In the following example, the variable \verb|t| is an element of $\Z$ whereas804\verb|u| is in $\Q$. In order to perform the addition, the coercion system805first deduces that the result should be in $\Q$ from the fact that \verb|t|806can be converted to the domain of \verb|u|, namely $\Q$, but canonical conversion807in the other direction is not possible. Then the addition is performed with808both operands having the same domain $\Q$.809810\begin{lstlisting}811sage: t = 1812sage: t.parent()813Integer Ring814sage: u = 1/2815sage: u.parent()816Rational Field817sage: v = t + u; v8183/2819sage: v.parent()820Rational Field821\end{lstlisting}822823Similarly, in the following example, the common domain $\Q[x]$ is found for824arguments from $\Z[x]$ and $\Q$. Note that in this case, the result is not in825the domain of either of the operands.826827\begin{lstlisting}828sage: R.<x> = ZZ[]829sage: r = x + 1/2830sage: r.parent()831Univariate Polynomial Ring in x over Rational Field832sage: 5*r8335*x + 5/2834\end{lstlisting}835836%Sage's coercion system is perhaps more systematic than Magma's; for example,837%Magma V2.16-8 does not handle the above example:838%\begin{lstlisting}839% > R<x> := PolynomialRing(Integers());840% > x + 1/2;841% ^842% Runtime error in '+': Bad argument types843% Argument types given: RngUPolElt[RngInt], FldRatElt844%\end{lstlisting}845846%\begin{lstlisting}847% sage: R.<x> = QQ[]848% sage: P.<x,y> = ZZ[]849% sage: r = R.random_element(); r850% -1/2*x^2 + 12*x + 4/3851% sage: res = r + y; res852% -1/2*x^2 + 12*x + y + 4/3853% sage: res.parent()854% Multivariate Polynomial Ring in x, y over Rational Field855%\end{lstlisting}856857%For more details on the coercion system, and the full range of supported858%features.see \todo{coercion section in the reference manual?}859860\subsubsection{Algebraic Tools: The Category Framework}861862Another abstraction to make implementing mathematical863structures easier is the {\em category framework}, whose development864was spearheaded by Nicolas Thi\'ery and Florent Hivert. Similar in spirit to the865mathematical programming facilities developed in Axiom and encapsulated866in Aldor, the category framework uses Python's dynamic class867creation capabilities to combine functions relevant for a mathematical object,868inherited through a mathematical hierarchy, into a class at run time.869870This process greatly simplifies the troubles of having to combine object-oriented programming concepts with mathematical structural concerns, while871keeping efficiency in mind. Efficient implementations can keep the872inheritance hierarchy imposed by the data structures, while generic methods873to compute basic properties are implemented in the {\em category} and874automatically attached to the element classes when they are needed.875876%\todo{refer to877%http://www.sagemath.org/doc/reference/sage/categories/primer.html?}878879\subsubsection{Symbolic Tools}880881The symbolic subsystem of Sage provides an environment similar to Maple or882Mathematica, where the input is treated only as an expression883without any concern about the underlying mathematical structure.884885886%Pynac (http://pynac.sagemath.org/), which was started in August 2008887%by William Stein and Burcin Erocal, is a hybrid C++ and Cython library888%built on top of GiNaC. Describing itself as "an open framework for889%symbolic computation with the C++ programming language," GiNaC is890%indeed a powerful C++ library for symbolic manipulation that was first891%released in 1999. All of GiNaC's basic arithmetic is done by the CLN892%library. For Sage, we replaced all dependence on CLN by calls to893%Cython functions via the Cython public C API, thus simultaneously894%replacing the depence on CLN by a dependence on Python, and making it895%possible for GiNaC to directly manipulate symbolic expressions896%involving arbitrary Python objects.897898Sage uses Pynac \cite{pynac}, a hybrid C++ and Cython library built on top of GiNaC \cite{ginac},899to work with symbolic expressions. High level symbolic calculus problems900including symbolic integration, solution of differential equations and Laplace901transforms are solved using Maxima behind the scenes.902903Here is an example of how to use the symbolic computation facilities in Sage.904Note that in contrast to other symbolic software such as Maple, variables must905be declared before they are used.906\begin{lstlisting}907sage: x,y,z = var('x,y,z')908sage: sin(x).diff(x)909cos(x)910sage: psi(x).series(x,4)911(-1)*x^(-1) + (-euler_gamma) + (1/6*pi^2)*x + \912(-zeta(3))*x^2 + (1/90*pi^4)*x^3 + Order(x^4)913sage: w = SR.wild() # wildcard for symbolic substitutions914sage: ((x^2+y^2+z^2)*zeta(x)).subs({w^2:5})91515*zeta(x)916\end{lstlisting}917918%In order to develop new symbolics functionality in Sage, in most cases919%you do not need to modify the Pynac library directly. Everything can,920%and should, be done from the Python (and Cython) layer in Sage. Pynac921%only provides fast arithmetic and pattern matching support in the922%backend. Since its first release as part of Sage one year ago, few923%changes in Pynac were needed, which shows the API between Sage and924%Pynac is very stable.925926927%Maxima: a Lisp CAS linked to Sage928%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~929930%Maxima is the GPL-licensed free software descendant of the Macsyma931%project, which was developed for the DOE at the Massachusetts932%Institute of Technology starting in the late 1960s. Maxima has a wide933%range of capabilities. Many algorithms for symbolic computation, such934%as Gosper's algorithm for indefinite summation of hypergeometric935%functions, were first implemented in Maxima.936937%Sage uses Maxima as the basis for certain symbolic calculus938%facilities, Maxima is designed as an939%interactive command line program, where the user types in input at the940%terminal and answers questions that come up during computation as941%necessary. The Sage interface to Maxima hides this interaction and942%string passing behind the scenes and provides a unified library943%interface as a part of the other symbolics functions in Sage.944%945%\begin{lstlisting}946% sage: var('x,z,s')947% (x, z, s)948% sage: (z + exp(x)).laplace(x, s)949% z/s + 1/(s - 1)950%\end{lstlisting}951952%Recently, a library interface to Maxima, using the foreign function953%interface of the ECL lisp environment, was developed by Nils Bruin and954%included in Sage. While still not used by default in Sage, this955%interface provides a promising approach to further enhance the956%efficiency and robustness of using Maxima from within Sage. For957%example, we illustrate the relative speed of computing 2+2 using the958%library interface to ECL, using a pexpect interface to ECL, and using959%our pexpect-based interface to Maxima below.960%961%\begin{lstlisting}962% sage: import sage.libs.ecl963% sage: timeit("sage.libs.ecl.ecl_eval('(+ 2 2)')")964% 625 loops, best of 3: 8.71 microseconds per loop965% sage: timeit("lisp('(+ 2 2)')")966% 625 loops, best of 3: 480 microseconds per loop967% sage: timeit("maxima('2 + 2')")968% 125 loops, best of 3: 1330 microseconds per loop969%\end{lstlisting}970971\subsubsection{Numerical Tools}972In addition to code for symbolic computation, the standard numerical973Python packages NumPy, SciPy, and Matplotlib are included in Sage,974along with the numerical libraries cvxopt, GSL, Mpmath, and R.975976For numerical applications, Robert Bradshaw and Carl Witty developed977a compiler for Sage that converts978symbolic expressions into an internal format suitable for blazingly fast979floating point evaluation.980981\begin{lstlisting}982sage: f(x,y) = sqrt(x^2 + y^2)983sage: a = float(2)984sage: timeit('float(f(a,a))')985625 loops, best of 3: 216 microseconds per loop986sage: g = fast_float(f)987sage: timeit('float(g(a,a))')988625 loops, best of 3: 0.406 microseconds per loop989\end{lstlisting}990The \verb|fast_float| feature is automatically used by the \verb|minimize| command.991\begin{lstlisting}992sage: minimize(f, (a,a))993(-5.65756135618e-05, -5.65756135618e-05)994\end{lstlisting}995Performance is typically within a factor of996two from what one gets using a direct implementation in C or Fortran.997998999\section{Afterword}\label{sec:interesting}10001001In this article, we have showed that1002Sage is a powerful platform for developing sophisticated mathematical1003software. Sage is actively used in research mathematics,1004and people use Sage to develop state-of-the-art algorithms.1005Sage is particularly strong in number theory, algebraic combinatorics,1006and graph theory. For further examples, see the100753 published articles, 11 Ph.D. theses, 10 books, and 30 preprints at1008\url{www.sagemath.org/library-publications.html}.10091010For example, Sage has extensive functionality for computations related to the1011Birch and Swinnerton-Dyer conjecture. In addition to Mordell-Weil1012group computations using \cite{mwrank} and point counting over large1013finite fields using the SEA package in \cite{pari}, there is much1014novel elliptic curve code written directly for Sage. This includes the1015fastest known algorithm for computation of $p$-adic heights1016\cite{harvey:padicheights,mazur-stein-tate:padic}, and code for computing $p$-adic1017$L$-series of elliptic curves at ordinary, supersingular, and split1018multiplicative primes. Sage combines these capabilities to1019compute explicitly bounds on Shafarevich-Tate groups of elliptic1020curves \cite{stein-wuthrich}. Sage also has code for computation1021with modular forms, modular abelian varieties, and ideal class groups in1022quaternion algebras.10231024The MuPAD-combinat project, which was started by Florent Hivert and1025Nicolas M. Thi\'ery in 2000, built the world's preeminent system for1026algebraic combinatorics on top of MuPAD (see \cite{icms:mupad2006} and1027\cite{mupad-hivert-thiery}). Page 54 of \cite{mupad-hivert-thiery}:1028``They [MuPAD] also have promised to release the code source of the1029library under a well known open-source license, some day.'' In 2008,1030MuPAD was instead purchased by MathWorks (makers of MATLAB), so MuPAD1031is no longer available as a separate product, and will probably never1032be open source. Instead it now suddenly costs \$30001033(commercial) or \$700 (academic).1034%3000=2000+10001035%700=500+2001036%\begin{quote}1037%{\small1038%``{\bf 2. Can I purchase MuPAD Pro licenses?}\\1039%MuPAD Pro is no longer sold as a standalone product. However, The1040%MathWorks has incorporated MuPAD technology into Symbolic Math1041%Toolbox, which can be purchased for use with MATLAB.''1042%-- See \cite{matlabmupad}. }1043%\end{quote}10441045As a result, the MuPAD-combinat group has spent several years1046reimplementing everything in Sage (see \cite{sagecombinatroadmap} for1047the current status). The MuPAD-combinat1048group was not taken by surprise by the failure of MuPAD, but instead1049were concerned from the beginning by the inherent risk in building1050their research program on top of MuPAD. In fact, they decided to1051switch to Sage two months before the bad news hit, and have made1052tremendous progress porting:1053\begin{quote}1054{1055``It has been such a relief during the last two years not to have1056this Damocles sword on our head!''10571058\hspace{5em} -- Nicolas Thi\'ery}1059\end{quote}106010611062% \begin{verbatim}1063% talk about the main areas of current devel focus for sage, since current1064% momentum best predictor for future:10651066% - number theory and arith geom1067% - algebraic combinatorics1068% - symbolic comp1069% - graph theory1070% - algebraic topology (?)10711072% KILLER FEATURES:1073% E.g., unique algorithms, like p-adic l-series.1074% flint.1075% interfaces (we already discuss that).1076% p-adic cohomology to get zeta functions.1077% modular forms1078% fast computation of det's.1079% - multiply polys over ZZ1080% - graph isomorphism1081% - algebraic topology -- mention Steenrod algebras10821083% People love learning about "killer features", since they are useful to know about when you attack a problem.1084% You can think "aha, sage can do that".1085% and algebraic combinatorics.1086% we can mention our graph isomorphism algorithm.1087% and algebraic topology -- mention Steenrod algebras.1088% convex polytopes1089% \end{verbatim}109010911092109310941095\bibliographystyle{amsalpha}1096\bibliography{sage}1097\end{document}109810991100