CoCalc Shared Fileswww / papers / icms / icms_2010.tex
Author: William A. Stein
1\documentclass[runningheads,a4paper]{llncs}
2\usepackage{ucs}
3\usepackage[utf8x,utf8]{inputenc}
4
5\usepackage{amsmath,amssymb}
6\setcounter{tocdepth}{3}
7\usepackage[pdftex]{graphicx}
8
9\usepackage{url}
10\usepackage{hyperref}
11
12\newcommand{\todo}[1]{[[\rule{0.5em}{1ex}\hspace{1em}#1]]}
13\newcommand{\todobf}[1]{ { \rule{0.5em}{1ex}\textbf{ #1}}}
14
15\newcommand{\goth}{\mathfrak}
16
17\usepackage{fancybox}
18
19
20\def\C{\mathbb C}
21\def\Q{\mathbb Q}
22\def\Z{\mathbb Z}
23\def\R{\mathbb R}
24
25
26\usepackage{color}
27\usepackage{listings}
28\usepackage{xspace}  % to allow for text macros that don't eat space
29
30\lstdefinelanguage{Sage}[]{Python}
31{morekeywords={True,False,sage,singular},
32sensitive=true}
33\lstset{frame=none,
34          showtabs=False,
35          showspaces=False,
36          showstringspaces=False,
38          keywordstyle={\ttfamily\color{dbluecolor}\bfseries},
39          stringstyle ={\ttfamily\color{dgraycolor}\bfseries},
40          language = Sage,
41	  basicstyle={\small \ttfamily},
42	  aboveskip=.3em,
43	  belowskip=.1em
44          }
45\definecolor{dblackcolor}{rgb}{0.0,0.0,0.0}
46\definecolor{dbluecolor}{rgb}{.01,.02,0.7}
47\definecolor{dredcolor}{rgb}{0.8,0,0}
48\definecolor{dgraycolor}{rgb}{0.30,0.3,0.30}
49\newcommand{\dblue}{\color{dbluecolor}\bf}
50\newcommand{\dred}{\color{dredcolor}\bf}
51\newcommand{\dblack}{\color{dblackcolor}\bf}
52
53\renewcommand{\emph}[1]{{\dblack{#1}}}
54
55\newcommand{\Sage}{{\color{dbluecolor}\sf Sage}\xspace}
56
57
59\noindent\keywordname\enspace\ignorespaces#1}
60
61\begin{document}
62
63\mainmatter  % start of an individual contribution
64
65% first the title is needed
66\title{The Sage Project:\\
67\large Unifying Free Mathematical Software to Create a Viable Alternative to
68Magma, Maple, Mathematica and MATLAB}
69
70% a short form should be given in case it is too long for the running head
71\titlerunning{Sage: Unifying Free Mathematical Software}
72
73% the name(s) of the author(s) follow(s) next
74%
75\author{Burçin Eröcal\inst{1}%
76%\thanks{Please note that the LNCS Editorial assumes that all authors have used
77%the western naming convention, with given names preceding surnames. This determines
78%the structure of the names in the running heads and the author index.}%
79\and William Stein\inst{2}}
80%
81%\authorrunning{}
82
84% unless you accept that it will be published
85\institute{Research Institute for Symbolic Computation\\
86Johannes Kepler University,\\
87Linz, Austria\\
88Supported by FWF grants P20347 and DK W1214.\\
89\email{burcin@erocal.org}
90\and
91Department of Mathematics\\
92University of Washington\\
93\email{wstein@uw.edu}\\
94Supported by NSF grant DMS-0757627 and DMS-0555776.}
95
96%
97% NB: a more complex sample for affiliations and the mapping to the
98% corresponding authors can be found in the file "llncs.dem"
99% (search for the string "\mainmatter" where a contribution starts).
100% "llncs.dem" accompanies the document class "llncs.cls".
101%
102
103%\toctitle{Lecture Notes in Computer Science}
104%\tocauthor{Authors' Instructions}
105\maketitle
106
107
108%The abstract should summarize the contents of the paper and should
109%contain at least 70 and at most 150 words. It should be written using the
110%\emph{abstract} environment.
111\begin{abstract}
112  Sage is a free, open source,
113  self-contained distribution of mathematical software, including a
114  large library that provides a unified interface to the components of
115  this distribution.  This library also builds on the components of
116  Sage to implement novel algorithms covering a broad range of
117  mathematical functionality from algebraic combinatorics to number
118  theory and arithmetic geometry.
119
120  \keywords{Python, Cython, Sage, Open Source, Interfaces}
121\end{abstract}
122
123%\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}
124
125%\todo{Interesting problems - (Test drive) - showcasing the capabilities (of the car)}
126
127
128\section{Introduction}
129
130In order to use mathematical software for exploration, we often push
131the boundaries of available computing resources and continuously try to
132improve our implementations and algorithms.  Most mathematical
133algorithms require basic building blocks, such as multiprecision
134numbers, fast polynomial arithmetic, exact or numeric linear algebra,
135or more advanced algorithms such as Gr\"obner basis computation or
136integer factorization.  Though implementing some of these basic
137foundations from scratch can be a good exercise, the resulting code
138may be slow and buggy.  Instead, one can build on existing optimized
139implementations of these basic components, either by using a general
140computer algebra system, such as Magma, Maple, Mathematica or MATLAB,
141or by making use of the many high quality open source libraries that
142provide the desired functionality.  These two approaches both have
143significant drawbacks.  This paper is about
144Sage,\footnote{\url{http://www.sagemath.org}} which provides an
145alternative approach to this problem.
146
147Having to rely on a closed propriety system can be frustrating, since
148it is difficult to gain access to the source code of the software,
149either to correct a bug or include a simple optimization in an
150algorithm.  Sometimes this is by design:
151\begin{quote}\small Indeed, in almost all practical uses of Mathematica,
152issues about how Mathematica works inside turn out to be largely
153irrelevant.   You might think that knowing how Mathematica
154works inside would be necessary [...]'' (See \cite{mathematica:internals}.)
155\end{quote}
156Even if we manage to contact the developers, and they find time to
157make the changes we request, it might still take months or years
158before these changes are made available in a new release.
159
160Fundamental questions of correctness, reproducibility and scientific
161value arise when building a mathematical research program on top of
162proprietary software (see, e.g., \cite{joyner-stein:notices}).  There
163are many published refereed papers containing results that
164rely on computations performed in Magma, Maple, or
165Mathematica.\footnote{Including by the second author of this paper,
166  e.g., \cite{conrad-edixhoven-stein:j1p}!} In some cases, a specific
167version of Magma is the only software that can carry out the
168computation.  This is not the infrastructure on which we want to build
169the future of mathematical research.
170
171In sharp contrast, open source libraries provide a great deal of
172flexibility, since anyone can see and modify the source code as they
173wish.  However, functionality is often segmented into different
174specialized libraries and advanced algorithms are hidden behind custom
175interpreted languages. One often runs into trouble trying to install
176dependencies before being able use an open source software package.
177Also, converting the output of one package to the input format of
178another package can present numerous difficulties and
179introduce subtle errors.
180
181Sage, which started in 2005 (see \cite{joyner-stein:sigsam}),
182attacks this problem by providing:
183
184\begin{enumerate}
185\item a {\em self-contained distribution} of mathematical software that
186  installs from source easily, with the only dependency being compiler
187  tools,
188\item {\em unified interfaces} to other mathematical software to make it
189  easier to use all these programs together, and
190\item a {\em new library} that builds on the included software packages and
191  implements a broad range of mathematical functionality.
192\end{enumerate}
193
194%\begin{center}
196%\begin{minipage}{0.85\textwidth}
197%{\bf The Sage Mission Statement:} Create a viable free open source %alternative
198%to Magma, Maple, Mathematica and Matlab.
199%\end{minipage}}
200%\end{center}
201
202
203
204
205%William Stein started the Sage project in early 2005, with financial
206%support from the NSF, Google, Microsoft, University of Washington, and
207%many other organizations and individuals.
208
209%Sage\footnote{\url{http://www.sagemath.org}} is a free, open source, self-contained
210%distribution of mathematical software, including a large library
211%that provides a unified interface to the components of this distribution.  This library also builds upon
212%the components of Sage to implement novel algorithms covering a broad
213%range of mathematical functionality from algebraic combinatorics to
214%number theory and arithmetic geometry.
215
216
217The rest of this paper goes into more detail about Sage.  In
218Section~\ref{sec:notebook}, we describe the Sage graphical user
219interface.  Section~\ref{sec:dev} is about the Sage development
220process, Sage days workshops, mailing lists,
221and documentation. The subject of Section~\ref{sec:car} is the
222sophisticated way in which Sage is built out of a wide range of open
223source libraries and software.  In Section~\ref{sec:interfaces} we
224explain how we use Python and Cython as the glue that
225binds the compendium of software included in Sage into a unified
226whole.  We then delve deeper into Python, Cython and the Sage
227preparser in Section~\ref{sec:python}, and illustrate some
228applications to mathematics in Section~\ref{sec:algebraic}.  Sage is
229actively used for research, and in
230Section~\ref{sec:interesting} we describe some
231capabilities of Sage in advanced areas of mathematics.
232
233\subsection{The Notebook}\label{sec:notebook}
234\begin{figure}
235\begin{center}
237\caption{The Sage Notebook\label{fig:sagenb}}
238\end{center}
239\end{figure}
240As illustrated in Figure~\ref{fig:sagenb}, the graphical user
241interface for Sage is a web application, inspired by Google Documents
243capabilities of Sage, including 3D graphics.  In single user mode,
244Sage works like a regular application whose main window happens to be
246easily set up servers for accessing their work over the Internet as
247well as sharing and collaborating with colleagues.  One can try the
248Sage notebook by visiting \url{www.sagenb.org}, where there are over
24930,000 user accounts and over 2,000 published worksheets.
250
251Users also download Sage to run it directly on their computers.  We
253are several other high-profile sites that provide mirrors of our
255month directly from the Sage website.
256
257
258\subsection{The Sage Development Process}\label{sec:dev}
259
260There are over 200 developers from across the world who have
261contributed to the Sage project.  People often contribute because they
262write code using Sage as part of a research project, and in this
263process find and fix bugs, speed up parts of Sage, or want the code
264portion of their research to be peer reviewed.  Each contribution to
265Sage is first posted to the Sage Trac server
266\url{trac.sagemath.org}; it is then peer reviewed, and finally
267added to Sage after all issues have been sorted out and all
269step of the peer review process is recorded indefinitely for all to
270see.
271
272The Sage Developer's Guide begins with an easy-to-follow tutorial that
273guides developers through each step involved in contributing code to Sage.
274Swift feedback is available through the \verb|sage-devel| mailing
275list, and the \verb|#sage-devel| IRC chat room on
276\verb|irc.freenode.net| (see
277\url{www.sagemath.org/development.html}).
278
279%(#sage-devel on irc.freenode.net).
280
281
282Much development of Sage has taken place at the Sage Days
283workshops.  There have been two dozen Sage Days
284\cite{sagedays} and many more are planned.  These are essential to
285sustaining the momentum of the Sage project and also help ensure that
286developers work together toward a common goal, rather than competing
287with each other and fragmenting our limited community resources.
288
289A major goal is ensuring that there will be many Sage Days workshops
290for the next couple of years.  The topics will depend on funding, but
291will likely include numerical computation, large-scale bug fixing,
292$L$-functions and modular forms, function fields, symbolic
293computation, topology, and combinatorics. The combination of
294experienced developers with a group of enthusiastic mathematicians at
295each of these workshops has rapidly increased the developer community,
296and we hope that it will continue to do so.
297
298%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
299
300
301\section{Building the Car\dots}\label{sec:car}
302\begin{center}
303\includegraphics[height=3cm]{fig/sagecar_crop}
304\end{center}
305
306
307With the motto building the car instead of reinventing the wheel,''
308Sage  brings together numerous open source software
309packages (see Table~\ref{tab:standard} and \cite{components}).
310
311\begin{table}
312\begin{center}
313\caption{Packages Included With Every Copy of Sage-4.4.2\label{tab:standard}}
314\begin{tabular}{|l|l|l|l|l|}\hline
315atlas& gap& libgcrypt& palp& scipy\_sandbox\\\hline
316blas& gd& libgpg\_error& pari& scons\\\hline
317boehm\_gc& gdmodule& libm4ri& pexpect& setuptools\\\hline
318boost& genus2reduction& libpng& pil& singular\\\hline
319cddlib& gfan& linbox& polybori& sphinx\\\hline
320cliquer& ghmm& matplotlib& pycrypto& sqlalchemy\\\hline
321cvxopt& givaro& maxima& pygments& sqlite\\\hline
322cython& gnutls& mercurial& pynac& symmetrica\\\hline
323docutils& gsl& moin& python& sympow\\\hline
324ecl& iconv& mpfi& python\_gnutls& sympy\\\hline
325eclib& iml& mpfr& r& tachyon\\\hline
326ecm& ipython& mpir& ratpoints& termcap\\\hline
328flint& jinja2& networkx& rubiks& weave\\\hline
329flintqs& lapack& ntl& sagenb& zlib\\\hline
330fortran& lcalc& numpy& sagetex& zn\_poly\\\hline
331freetype& libfplll& opencdk& scipy& zodb3\\\hline
332\end{tabular}
333\end{center}
334\end{table}
335
336Many applications of Sage require using these
337libraries together. Sage handles the conversion of data behind the
338scenes, automatically using the best tool for the job, and allows
339the user to concentrate on the problem at hand.
340
341In the following example, which we explain in detail below, Sage
342uses the FLINT library \cite{flint} for univariate
343polynomials over the ring $\Z$ of integers, whereas Singular
344\cite{singular} is used for multivariate
345polynomials. The option to use the NTL library \cite{ntl} for
346univariate polynomials is still available, if the user so chooses.
347
348\begin{lstlisting}[numbers=left,numberstyle=\tiny,numbersep=0pt]
349    sage: R.<x> = ZZ[]
350    sage: type(R.an_element())
351    <type 'sage.rings...Polynomial_integer_dense_flint'>
352    sage: R.<x,y> = ZZ[]
353    sage: type(R.an_element())
354    <type 'sage.rings...MPolynomial_libsingular'>
355    sage: R = PolynomialRing(ZZ, 'x', implementation='NTL')
356    sage: type(R.an_element())
357    <type 'sage.rings...Polynomial_integer_dense_ntl'>
358\end{lstlisting}
359
360The first line in the example above constructs the univariate polynomial ring
361${\tt R} = \Z[x]$, and assigns the variable \verb|x| to be the generator of
362this ring. Note that $\Z$ is represented by \verb|ZZ| in Sage. The
363expression \verb|R.<x> = ZZ[]| is not valid Python, but can be used
364in Sage code as a shorthand as explained in Section \ref{sec:preparser}. The
365next line asks the ring \verb|R| for an element, using the
366\verb|an_element| function, then uses the builtin Python function
367\verb|type| to query its type. We learn that it is an instance of the class
368\verb|Polynomial_integer_dense_flint|. Similarly line 4 constructs ${\tt R} 369= \Z[x,y]$ and line 7 defines ${\tt R} = \Z[x]$, but this time using the
370\verb|PolynomialRing| constructor explicitly and specifying that we want
371the underlying implementation to use the NTL library.
372
373Often these interfaces are used under the hood, without the user
374having to know anything about the corresponding systems.  Nonetheless, there are easy
375ways to find out what is used by inspecting the source code, and users
376are strongly encouraged to cite components they use in published
377papers.  The following example illustrates another
378way to get a list of
379components used when a specific command is run.
380
381\begin{lstlisting}
382    sage: from sage.misc.citation import get_systems
383    sage: get_systems('integrate(x^2, x)')
384    ['ginac', 'Maxima']
385    sage: R.<x,y,z> = QQ[]
386    sage: I = R.ideal(x^2+y^2, z^2+y)
387    sage: get_systems('I.primary_decomposition()')
388    ['Singular']
389\end{lstlisting}
390
391\subsection{Interfaces}\label{sec:interfaces}
392
393Sage makes it possible to use a wide range of mathematical software
394packages together by providing a unified interface that
395handles data conversion automatically. The complexity and
396functionality of these interfaces varies greatly, from simple
397text-based interfaces that call external software for an individual
398computation, to using a library as the basis for an
399arithmetic type.  The interfaces can also run code from libraries
400written in the interpreted language of another program.
401Table~\ref{fig:interfaces} lists
402the interfaces provided by Sage.
403
404\begin{table}
405\begin{center}
406\caption{Sage Interfaces to the above Mathematical Software\label{fig:interfaces}}
407
408\begin{tabular}{|l|l|}\hline
409{\bf Pexpect} & axiom, ecm, fricas, frobby, gap, g2red, gfan,
410gnuplot, gp, \\
411& kash, lie, lisp, macaulay2, magma, maple, mathematica,
412\\
413& matlab, maxima, mupad, mwrank, octave, phc, polymake, \\
415rubik, scilab, singular, tachyon\\\hline
416{\bf C Library} & eclib, fplll, gap (in progress), iml, linbox, maxima,\\
417   & ratpoints, r (via rpy2), singular, symmetrica\\\hline
418{\bf C Library arithmetic} & flint, mpir, ntl, pari, polybori, pynac, singular\\\hline
419\end{tabular}
420\end{center}
421\end{table}
422
423
424The above interfaces are the result of many years writing
425Python and  Cython \cite{cython} code to adapt
426Singular \cite{singular}, GAP \cite{gap}, Maxima \cite{maxima}, Pari
427\cite{pari}, GiNaC/Pynac \cite{ginac}, NTL \cite{ntl}, FLINT
428\cite{flint}, and many other libraries, so that they can be used
429smoothly and efficiently in a unified way from Python \cite{python}.
430Some of these programs were originally designed to be used only
431through their own interpreter and made into a library by Sage
432developers. For example libSingular was created by Martin Albrecht in
433order to use the fast multivariate polynomial arithmetic in Singular
434from Sage. The libSingular interface is now used by other
435projects, including Macaulay2 \cite{M2} and GFan \cite{gfan}.
436
437There are other approaches to linking mathematical software together.
438The recent paper \cite{science_openmath} reports on the state of the
439art  using OpenMath.  Sage takes a
440dramatically different approach to this problem.  Instead of using a
441general string-based XML protocol to communicate with other mathematical
442software, Sage interfaces are tailor made to the specific software and
443problem at hand. This results in far more efficient and flexible
444interfaces.  The main disadvantage compared to OpenMath is that the
445interfaces all go through Sage.
446
447%Sage also allows easy use of commercial packages, by encapsulating the
448%knowledge\todo{sounds funny} to transfer data to a format these packages can
449%understand.
450
452without having to worry about data conversion, also makes it easier to
453double check results.  For example, below we first use Maxima, an open
454source symbolic computation package distributed with Sage, to
455integrate a function, then perform the same computation using Maple
456and Mathematica.
457
458\begin{lstlisting}
459    sage: var('x')
460    sage: integrate(sin(x^2), x)
461    1/8*((I - 1)*sqrt(2)*erf((1/2*I - 1/2)*sqrt(2)*x) + \
462      (I + 1)*sqrt(2)*erf((1/2*I + 1/2)*sqrt(2)*x))*sqrt(pi)
463    sage: maple(sin(x^2)).integrate(x)
464    1/2*2^(1/2)*Pi^(1/2)*FresnelS(2^(1/2)/Pi^(1/2)*x)
465    sage: mathematica(sin(x^2)).Integrate(x)
466    Sqrt[Pi/2]*FresnelS[Sqrt[2/Pi]*x]
467\end{lstlisting}
468
469
470
471
472
473The most common type of interface, called a \verb|pexpect|
474interface, communicates with another command line program by reading
475and writing strings to a text console, as if another user was in front
476of the terminal. Even though these are relatively simple to develop,
477the overhead of having to print and parse strings to represent the
478data makes this process potentially cumbersome and inefficient.  This
479is the default method of communication with most high level
480mathematics software, including commercial and open source programs,
481such as Maple, Mathematica, Magma, KASH or GAP.
482
483Sage provides a framework to represent elements over these interfaces, perform
484arithmetic with them or apply functions to the given object, as well as using a
485file to pass the data if the string representation is too big. The following
486demonstrates arithmetic with GAP elements.
487
488\begin{lstlisting}
489    sage: a = gap('22')
490    sage: a*a
491    484
492\end{lstlisting}
493
494It is also possible to use \verb|pexpect| interfaces over remote consoles. In
495the following code, we connect to the \verb|localhost| as a different user
496and call Mathematica functions. Note that the interface can handle indexing
497vectors as well.
498
499\begin{lstlisting}
500    sage: mma = Mathematica(server="rmma60@localhost")
501    sage: mma("2+2")
502    4
503    sage: t = mma("Cos[x]")
504    sage: t.Integrate('x')
505    Sin[x]
506    sage: t = mma('{0,1,2,3}')
507    sage: t[2]
508    1
509\end{lstlisting}
510
511Sage also includes specialized libraries that are linked directly from compiled
512code written in Cython. These are used to handle specific problems, such as the
513characteristic polynomial computation in the example below.
514
515\begin{lstlisting}
516    sage: M = Matrix(GF(5), 10, 10)
517    sage: M.randomize()
518    sage: M.charpoly(algorithm='linbox')
519    x^10 + 4*x^9 + 4*x^7 + 3*x^4 + 3*x^3 + 3*x^2 + 4*x + 3
520\end{lstlisting}
521
522Many basic arithmetic types also use Cython to directly utilize data
523structures from efficient arithmetic libraries, such as MPIR or FLINT. An
524example of this can be seen at the beginning of this section, where elements
525of the ring $\Z[x]$ are represented by the class
526\verb|Polynomial_integer_dense_flint|.
527
528The Singular interface is one of the most advanced included in
529Sage. Singular has a large library of code written in its own language.
530Previously the only way to access these functions, which include
531algorithms for  Gr\"obner basis and primary decomposition, was to call
532Singular through a \verb|pexpect| interface, passing data back and forth
533using strings. Recently, due to work of Michael Brickenstein and Martin
534Albrecht, Sage acquired the ability to call these functions directly.
535
536In the example below, we import the function \verb|primdecSY| from
537\verb|primdec.lib|, and call it the same way we would call a Python
538function. The interface handles the conversion of the data to Singular's format
539and back. Since Sage already uses Singular data structures directly to represent
540multivariate polynomials and ideals over multivariate polynomial rings, there
541are no conversion costs. It is only a matter of passing the right pointer.
542
543\begin{lstlisting}
544    sage: pr = sage.libs.singular.ff.primdec__lib.primdecSY
545    sage: R.<x,y,z> = QQ[]
546    sage: p = z^2+1; q = z^3+2
547    sage: I = R.ideal([p*q^2,y-z^2])
548    sage: pr(I)
549    [[[z^2 - y, y^3 + 4*y*z + 4], \
550	[z^2 - y, y*z + 2, y^2 + 2*z]], \
551	[[y + 1, z^2 + 1], [y + 1, z^2 + 1]]]
552\end{lstlisting}
553
554Efforts are under way to extend these capabilities to other programs,
555for example to GAP which provides Sage's underlying group theory
556functionality. Up to now, GAP was only available through its interpreter,
557through a \verb|pexpect| interface that was written by Steve Linton.  As the
558following example demonstrates, the performance of this interface is
559far from ideal.\footnote{All timings in this paper were performed on
560  an 2.66GHz Intel Xeon X7460 based computer.}
561
562\begin{lstlisting}
563    sage: b = gap('10')
564    sage: b*b
565    100
566    sage: timeit('b*b')
567    625 loops, best of 3: 289 microseconds per loop
568\end{lstlisting}
569
570The code snippet above constructs the element \verb|b| in GAP using the
571\verb|pexpect| interface, and measures the time it takes to square
572\verb|b|. Compare these numbers to the following example, which uses the
573library interface to GAP, recently developed by the second author
574(but {\em not} included in Sage yet).
575
576\begin{lstlisting}
577    sage: import sage.libs.gap.gap as g
578    sage: a = g.libgap('10'); a
579    10
580    sage: type(a)
581    <type 'sage.libs.gap.gap.GapElement'>
582    sage: a*a
583    100
584    sage: timeit('a*a')
585    625 loops, best of 3: 229 nanoseconds per loop
586\end{lstlisting}
587The library interface is about 1,000 times faster than the pexpect interface.
588
589
590\subsection{Python - a mainstream language}\label{sec:python}
591
592In line with the principle of not reinventing the wheel, Sage is built on the
593mainstream programming language Python, both as the main development language
594and the user language. This frees the Sage developers, who are mainly
596an immense array of general purpose Python libraries and tools.
597
598Python is an interpreted language with a clear, easy to read and learn
599syntax. Since it is dynamically typed, it is ideal for rapid prototyping,
600providing an environment to easily test new ideas and algorithms.
601
602%Note that many mathematical packages, such as GINV, PolyBoRi, or CVXOPT,
604
605\subsubsection{A fast interpreter}
606
607In the following Singular session, we first declare the ring ${\tt 608r}=\Q[x,y,z]$ and the polynomial ${\tt f} \in {\tt r}$, then
609measures the time to square ${\tt f}$ repeatedly, 10,000 times.
610
611\begin{lstlisting}
612    singular: int t = timer; ring r = 0,(x,y,z), dp;
613    singular: def f = y^2*z^2-x^2*y^3-x*z^3+x^3*y*z;
614    singular: int j; def g = f;
615    singular: for (j = 1; j <= 10^5; j++) { g = f*f; }
616    singular: (timer-t), system("--ticks-per-sec");
617    990 1000
618\end{lstlisting}
619
620The elapsed time is 990 milliseconds. Next we use Sage to do the same
621computation, using the same Singular data structures directly, but
622without going through the interpreter.
623
624\begin{lstlisting}
625    sage: R.<x,y,z> = QQ[]
626    sage: f = y^2*z^2 - x^2*y^3 - x*z^3 + x^3*y*z;  type(f)
627    <type 'sage.rings.polynomial...MPolynomial_libsingular'>
628    sage: timeit('for j in xrange(10^5): g = f*f')
629    5 loops, best of 3: 91.8 ms per loop
630\end{lstlisting}
631
632Sage takes only 91.8 milliseconds for the same operation. This difference is
633because the Python interpreter is more efficient at performing \verb|for| loops.
634
635
636
637\subsubsection{Cython - compiled extensions}
638
639Python alone is too slow to implement a serious mathematical software
640system.  Fortunately,  Cython \cite{cython}
641makes it easy to optimize parts of your
642program or access existing C/C++ libraries. It can translate Python
643code with annotations containing static type information to C/C++
644code, which is then compiled as a Python extension module.
645
646Many of the basic arithmetic types in Sage are provided by Cython wrappers of
647C libraries, such as FLINT for univariate polynomials over $\Z$, Singular for
648multivariate polynomials, and Pynac for symbolic expressions.
649
650The code segment below defines a Python function to add integers from 0 to
651$N$ and times the execution of this function with the argument \verb|10^7|.
652
653\begin{lstlisting}
654    sage: def mysum(N):
655    ....:     s = int(0)
656    ....:     for k in xrange(1,N): s += k
657    ....:     return s
658    ....:
659    sage: time mysum(10^7)
660    CPU times: user 0.52 s, sys: 0.00 s, total: 0.52 s
661    49999995000000
662\end{lstlisting}
663
664Here is the same function, but the loop index {\tt k} is declared to be
665a C integer and the accumulator {\tt s} is a C {\tt long long}.
666
667\begin{lstlisting}
668    sage: cython("""
669    ....: def mysum_cython(N):
670    ....:     cdef int k
671    ....:     cdef long long s = 0
672    ....:     for k in xrange(N): s += k
673    ....:     return s
674    ....: """)
675    sage: time mysum_cython(10^7)
676    CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
677    49999995000000L
678\end{lstlisting}
679
680The code is compiled and linked to the interpreter on the fly, and the
681function \verb|mysum_cython| is available immediately. Note that the run
682time for the Cython function is 60 times faster than the Python equivalent.
683
684Cython also handles the conversion of Python types to C types automatically.
685In the following example, we call the C function \verb|sinl| using Cython
686to wrap it in a Python function named \verb|sin_c_wrap|.
687
688\begin{lstlisting}
689    sage: cython("""
690    ....: cdef extern from "math.h":
691    ....:         long double sinl(long double)
692    ....: def sin_c_wrap(a):
693    ....:         return sinl(a)
694    ....: """)
695    sage: sin_c_wrap(3.14)
696    0.0015926529164868282
697    sage: sin_c_wrap(1)
698    0.8414709848078965
699    sage: sin_c_wrap(1r)
700    0.8414709848078965
701\end{lstlisting}
702
703Note that the conversion of Sage types in the first two calls to
704\verb|sin_c_wrap| or the Python type integer in the last call is performed
705transparently by Cython.
706
707\subsubsection{The Preparser}\label{sec:preparser}
708
709While Python has many advantages as a programming and glue language, it also
710has some undesirable features. Sage hides these problems by using a preparser
711to change the commands passed to Python in an interactive session (or when running
712a script with the \verb|.sage| extension). In order to
713maintain compatibility with Python, changes performed by the preparser are
714kept to a minimum.  Moreover, the Sage library code is not preparsed, and is
715written in Cython or Python directly.
716
717Python, like C and many other programming languages, performs integer floor
718division. This means typing \verb|1/2| results in 0, not the rational number
719$1/2$. Sage wraps all numeric literals entered in the command line or the
720notebook with its own type declarations, which behave as expected with respect
721to arithmetic and have the advantage that they are backed by efficient
722multiprecision arithmetic libraries such as MPIR \cite{mpir} and MPFR \cite{mpfr}, which are thousands of times faster than Python for large
723integer arithmetic.
724
725To call the preparser directly on a given string, use the
726\verb|preparse| function.
727
728\begin{lstlisting}
729    sage: preparse("1/2")
730    'Integer(1)/Integer(2)'
731    sage: preparse("1.5")
732    "RealNumber('1.5')"
733\end{lstlisting}
734
735Adding a trailing {\tt r} after a number indicates that the preparser should leave
736that as the raw'' literal.  The following illustrates division with Python integers.
737
738\begin{lstlisting}
739    sage: preparse("1r/2r")
740    '1/2'
741    sage: 1r/2r
742    0
743\end{lstlisting}
744Here is the result of performing the same division in Sage.
745\begin{lstlisting}
746    sage: 1/2
747    1/2
748    sage: type(1/2)
749    <type 'sage.rings.rational.Rational'>
750    sage: (1/2).parent()
751    Rational Field
752\end{lstlisting}
753
754The preparser also changes the \verb|^| sign to the exponentiation operator
755\verb|**| and provides a shorthand to create new mathematical domains and
756name their generator in one command.
757
758\begin{lstlisting}
759    sage: preparse("2^3")
760    'Integer(2)**Integer(3)'
761    sage: preparse("R.<x,y> = ZZ[]")
762    "R = ZZ['x, y']; (x, y,) = R._first_ngens(2)"
763\end{lstlisting}
764
765
766\subsection{Algebraic, Symbolic and Numerical Tools}~\label{sec:algebraic}
767
768Sage combines algebraic, symbolic and numerical computation tools under one
769roof, enabling users to choose the tool that best suits the problem. This
770combination also makes Sage more accessible to a wide audience---scientists,
771engineers, pure mathematicians and mathematics teachers
772can all use the same platform for scientific computation.
773
774While not concentrating on only one of these domains might seem to divide
775development resources unnecessarily, it actually results in a better overall
776experience for everyone, since users do not have to come up with makeshift
777solutions to compensate for the lack of functionality from a different field.
778Moreover, because Sage is a distributed mostly-volunteer open source
779project, widening our focus results in substantially more developer resources.
780
781\subsubsection{Algebraic Tools: The Coercion System}
782
783An algebraic framework, similar to that of Magma or Axiom, provides access to
784efficient data structures and specialized algorithms associated to particular
785mathematical domains.   The Python language allows
786classes to define how arithmetic
787operations like \verb|+| and \verb|*| will be handled, in a
788similar way to how C++ allows overloading of operators. However, the built-in support
790to support operations with a range of objects in a mathematical type hierarchy.
791
792Sage abstracts the process of deciding what an arithmetic operation
793means, or equivalently, in which domain the operation should be
794performed, in a framework called the {\em coercion system}, which was
795developed and implemented by Robert Bradshaw, David Roe,
796and many others. Implementations of new mathematical objects only need
797to define which other domains have a natural embedding to their
798domain. When performing arithmetic with objects, the
799coercion system will find a common domain where both arguments can be
800canonically mapped, perform the necessary type conversions automatically,
801thus allowing the implementation to only handle the case where both
802objects have the same parent.
803
804In the following example, the variable \verb|t| is an element of $\Z$ whereas
805\verb|u| is in $\Q$. In order to perform the addition, the coercion system
806first deduces that the result should be in $\Q$ from the fact that \verb|t|
807can be converted to the domain of \verb|u|, namely $\Q$, but canonical conversion
808in the other direction is not possible. Then the addition is performed with
809both operands having the same domain $\Q$.
810
811\begin{lstlisting}
812    sage: t = 1
813    sage: t.parent()
814    Integer Ring
815    sage: u = 1/2
816    sage: u.parent()
817    Rational Field
818    sage: v = t + u; v
819    3/2
820    sage: v.parent()
821    Rational Field
822\end{lstlisting}
823
824Similarly, in the following example, the common domain $\Q[x]$ is found for
825arguments from $\Z[x]$ and $\Q$. Note that in this case, the result is not in
826the domain of either of the operands.
827
828\begin{lstlisting}
829    sage: R.<x> = ZZ[]
830    sage: r = x + 1/2
831    sage: r.parent()
832    Univariate Polynomial Ring in x over Rational Field
833    sage: 5*r
834    5*x + 5/2
835\end{lstlisting}
836
837%Sage's coercion system is perhaps more systematic than Magma's; for example,
838%Magma V2.16-8 does not handle the above example:
839%\begin{lstlisting}
840%    > R<x> := PolynomialRing(Integers());
841%    > x + 1/2;
842%        ^
843%    Runtime error in '+': Bad argument types
844%    Argument types given: RngUPolElt[RngInt], FldRatElt
845%\end{lstlisting}
846
847%\begin{lstlisting}
848%    sage: R.<x> = QQ[]
849%    sage: P.<x,y> = ZZ[]
850%    sage: r = R.random_element(); r
851%    -1/2*x^2 + 12*x + 4/3
852%    sage: res = r + y; res
853%    -1/2*x^2 + 12*x + y + 4/3
854%    sage: res.parent()
855%    Multivariate Polynomial Ring in x, y over Rational Field
856%\end{lstlisting}
857
858%For more details on the coercion system, and the full range of supported
859%features.see \todo{coercion section in the reference manual?}
860
861\subsubsection{Algebraic Tools: The Category Framework}
862
863Another abstraction  to make implementing mathematical
864structures easier is the {\em category framework}, whose development
865was spearheaded by Nicolas Thi\'ery and Florent Hivert. Similar in spirit to the
866mathematical programming facilities developed in Axiom and encapsulated
867in Aldor, the category framework uses Python's dynamic class
868creation capabilities to combine functions relevant for a mathematical object,
869inherited through a mathematical hierarchy, into a class at run time.
870
871This process greatly simplifies the troubles of having to combine object-oriented programming concepts with mathematical structural concerns, while
872keeping efficiency in mind. Efficient implementations can keep the
873inheritance hierarchy imposed by the data structures, while generic methods
874to compute basic properties are implemented in the {\em category} and
875automatically attached to the element classes when they are needed.
876
877%\todo{refer to
878%http://www.sagemath.org/doc/reference/sage/categories/primer.html?}
879
880\subsubsection{Symbolic Tools}
881
882The symbolic subsystem of Sage provides an environment similar to Maple or
883Mathematica, where the input is treated only as an expression
884without any concern about the underlying mathematical structure.
885
886
887%Pynac (http://pynac.sagemath.org/), which was started in August 2008
888%by William Stein and Burcin Erocal, is a hybrid C++ and Cython library
889%built on top of GiNaC.  Describing itself as "an open framework for
890%symbolic computation with the C++ programming language," GiNaC is
891%indeed a powerful C++ library for symbolic manipulation that was first
892%released in 1999.  All of GiNaC's basic arithmetic is done by the CLN
893%library. For Sage, we replaced all dependence on CLN by calls to
894%Cython functions via the Cython public C API, thus simultaneously
895%replacing the depence on CLN by a dependence on Python, and making it
896%possible for GiNaC to directly manipulate symbolic expressions
897%involving arbitrary Python objects.
898
899Sage uses Pynac \cite{pynac}, a hybrid C++ and Cython library built on top of GiNaC \cite{ginac},
900to work with symbolic expressions. High level symbolic calculus problems
901including symbolic integration, solution of differential equations and Laplace
902transforms are solved using Maxima behind the scenes.
903
904Here is an example of how to use the symbolic computation facilities in Sage.
905Note that in contrast to other symbolic software such as Maple, variables must
906be declared before they are used.
907\begin{lstlisting}
908    sage: x,y,z = var('x,y,z')
909    sage: sin(x).diff(x)
910    cos(x)
911    sage: psi(x).series(x,4)
912    (-1)*x^(-1) + (-euler_gamma) + (1/6*pi^2)*x + \
913	(-zeta(3))*x^2 + (1/90*pi^4)*x^3 + Order(x^4)
914    sage: w = SR.wild()  # wildcard for symbolic substitutions
915    sage: ((x^2+y^2+z^2)*zeta(x)).subs({w^2:5})
916    15*zeta(x)
917\end{lstlisting}
918
919%In order to develop new symbolics functionality in Sage, in most cases
920%you do not need to modify the Pynac library directly. Everything can,
921%and should, be done from the Python (and Cython) layer in Sage. Pynac
922%only provides fast arithmetic and pattern matching support in the
923%backend.  Since its first release as part of Sage one year ago, few
924%changes in Pynac were needed, which shows the API between Sage and
925%Pynac is very stable.
926
927
928%Maxima: a Lisp CAS linked to Sage
929%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
930
931%Maxima is the GPL-licensed free software descendant of the Macsyma
932%project, which was developed for the DOE at the Massachusetts
933%Institute of Technology starting in the late 1960s. Maxima has a wide
934%range of capabilities. Many algorithms for symbolic computation, such
935%as Gosper's algorithm for indefinite summation of hypergeometric
936%functions, were first implemented in Maxima.
937
938%Sage uses Maxima as the basis for certain symbolic calculus
939%facilities,   Maxima is designed as an
940%interactive command line program, where the user types in input at the
941%terminal and answers questions that come up during computation as
942%necessary. The Sage interface to Maxima hides this interaction and
943%string passing behind the scenes and provides a unified library
944%interface as a part of the other symbolics functions in Sage.
945%
946%\begin{lstlisting}
947%    sage: var('x,z,s')
948%    (x, z, s)
949%    sage: (z + exp(x)).laplace(x, s)
950%    z/s + 1/(s - 1)
951%\end{lstlisting}
952
953%Recently, a library interface to Maxima, using the foreign function
954%interface of the ECL lisp environment, was developed by Nils Bruin and
955%included in Sage. While still not used by default in Sage, this
956%interface provides a promising approach to further enhance the
957%efficiency and robustness of using Maxima from within Sage.  For
958%example, we illustrate the relative speed of computing 2+2 using the
959%library interface to ECL, using a pexpect interface to ECL, and using
960%our pexpect-based interface to Maxima below.
961%
962%\begin{lstlisting}
963%    sage: import sage.libs.ecl
964%    sage: timeit("sage.libs.ecl.ecl_eval('(+ 2 2)')")
965%    625 loops, best of 3: 8.71 microseconds per loop
966%    sage: timeit("lisp('(+ 2 2)')")
967%    625 loops, best of 3: 480 microseconds per loop
968%    sage: timeit("maxima('2 + 2')")
969%    125 loops, best of 3: 1330 microseconds per loop
970%\end{lstlisting}
971
972\subsubsection{Numerical Tools}
973In addition to code for symbolic computation, the standard numerical
974Python packages NumPy, SciPy, and Matplotlib are included in Sage,
975along with the numerical libraries cvxopt, GSL, Mpmath, and R.
976
977For numerical applications, Robert Bradshaw and Carl Witty developed
978a compiler for Sage that converts
979symbolic expressions into an internal format suitable for blazingly fast
980floating point evaluation.
981
982\begin{lstlisting}
983    sage: f(x,y) = sqrt(x^2 + y^2)
984    sage: a = float(2)
985    sage: timeit('float(f(a,a))')
986    625 loops, best of 3: 216 microseconds per loop
987    sage: g = fast_float(f)
988    sage: timeit('float(g(a,a))')
989    625 loops, best of 3: 0.406 microseconds per loop
990\end{lstlisting}
991The \verb|fast_float| feature is automatically used by the \verb|minimize| command.
992\begin{lstlisting}
993    sage: minimize(f, (a,a))
994    (-5.65756135618e-05, -5.65756135618e-05)
995\end{lstlisting}
996Performance is typically within a factor of
997two from what one gets using a direct implementation in C or Fortran.
998
999
1000\section{Afterword}\label{sec:interesting}
1001
1003Sage is a powerful platform for developing sophisticated mathematical
1004software.  Sage is actively used in research mathematics,
1005and  people use Sage to develop state-of-the-art algorithms.
1006Sage is particularly strong in number theory, algebraic combinatorics,
1007and graph theory.  For further examples, see the
100853 published articles, 11 Ph.D. theses, 10 books, and 30 preprints at
1009\url{www.sagemath.org/library-publications.html}.
1010
1011For example, Sage has extensive functionality for computations related to the
1012Birch and Swinnerton-Dyer conjecture.  In addition to Mordell-Weil
1013group computations using \cite{mwrank} and point counting over large
1014finite fields using the SEA package in \cite{pari}, there is much
1015novel elliptic curve code written directly for Sage.  This includes the
1016fastest known algorithm for computation of $p$-adic heights
1017\cite{harvey:padicheights,mazur-stein-tate:padic}, and code for computing $p$-adic
1018$L$-series of elliptic curves at ordinary, supersingular, and split
1019multiplicative primes.  Sage combines these capabilities to
1020compute explicitly bounds on Shafarevich-Tate groups of elliptic
1021curves \cite{stein-wuthrich}.  Sage also has code for computation
1022with modular forms, modular abelian varieties, and ideal class groups in
1023quaternion algebras.
1024
1025The MuPAD-combinat project, which was started by Florent Hivert and
1026Nicolas M. Thi\'ery in 2000, built the world's preeminent system for
1029They [MuPAD] also have promised to release the code source of the
1030library under a well known open-source license, some day.''  In 2008,
1032is no longer available as a separate product, and will probably never
1033be open source.  Instead it now suddenly costs \$3000 1034(commercial) or \$700 (academic).
1035%3000=2000+1000
1036%700=500+200
1037%\begin{quote}
1038%{\small
1039%{\bf 2. Can I purchase MuPAD Pro licenses?}\\
1040%MuPAD Pro is no longer sold as a standalone product.  However, The
1041%MathWorks has incorporated MuPAD technology into Symbolic Math
1042%Toolbox, which can be purchased for use with MATLAB.''
1044%\end{quote}
1045
1046As a result, the MuPAD-combinat group has spent several years
1047reimplementing everything in Sage (see \cite{sagecombinatroadmap} for
1049group was not taken by surprise by the failure of MuPAD, but instead
1050were concerned from the beginning by the inherent risk in building
1051their research program on top of MuPAD.  In fact, they decided to
1052switch to Sage two months before the bad news hit, and have made
1053tremendous progress porting:
1054\begin{quote}
1055{
1056  It has been such a relief during the last two years not to have
1057  this Damocles sword on our head!''
1058
1059\hspace{5em} -- Nicolas Thi\'ery}
1060\end{quote}
1061
1062
1063% \begin{verbatim}
1064% talk about the main areas of current devel focus for sage, since current
1065% momentum best predictor for future:
1066
1067%     - number theory and arith geom
1068%     - algebraic combinatorics
1069%     - symbolic comp
1070%     - graph theory
1071%     - algebraic topology (?)
1072
1073% KILLER FEATURES:
1074% E.g., unique algorithms, like p-adic l-series.
1075% flint.
1076% interfaces (we already discuss that).
1077% p-adic cohomology to get zeta functions.
1078% modular forms
1079% fast computation of det's.
1080%      - multiply polys over ZZ
1081%      - graph isomorphism
1082%      - algebraic topology -- mention Steenrod algebras
1083
1084% People love learning about "killer features", since they are useful to know about when you attack a problem.
1085% You can think "aha, sage can do that".
1086% and algebraic combinatorics.
1087% we can mention our graph isomorphism algorithm.
1088% and algebraic topology -- mention Steenrod algebras.
1089% convex polytopes
1090% \end{verbatim}
1091
1092
1093
1094
1095
1096\bibliographystyle{amsalpha}
1097\bibliography{sage}
1098\end{document}
1099
1100