CoCalc Shared Fileswww / papers / icms / icms_2010.texOpen in CoCalc with one click!
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},
32
sensitive=true}
33
\lstset{frame=none,
34
showtabs=False,
35
showspaces=False,
36
showstringspaces=False,
37
commentstyle={\ttfamily\color{dredcolor}},
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
58
\newcommand{\keywords}[1]{\par\addvspace\baselineskip
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
68
Magma, 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
83
% the affiliations are given next; don't give your e-mail address
84
% unless you accept that it will be published
85
\institute{Research Institute for Symbolic Computation\\
86
Johannes Kepler University,\\
87
Linz, Austria\\
88
Supported by FWF grants P20347 and DK W1214.\\
89
\email{burcin@erocal.org}
90
\and
91
Department of Mathematics\\
92
University of Washington\\
93
\email{wstein@uw.edu}\\
94
Supported 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
130
In order to use mathematical software for exploration, we often push
131
the boundaries of available computing resources and continuously try to
132
improve our implementations and algorithms. Most mathematical
133
algorithms require basic building blocks, such as multiprecision
134
numbers, fast polynomial arithmetic, exact or numeric linear algebra,
135
or more advanced algorithms such as Gr\"obner basis computation or
136
integer factorization. Though implementing some of these basic
137
foundations from scratch can be a good exercise, the resulting code
138
may be slow and buggy. Instead, one can build on existing optimized
139
implementations of these basic components, either by using a general
140
computer algebra system, such as Magma, Maple, Mathematica or MATLAB,
141
or by making use of the many high quality open source libraries that
142
provide the desired functionality. These two approaches both have
143
significant drawbacks. This paper is about
144
Sage,\footnote{\url{http://www.sagemath.org}} which provides an
145
alternative approach to this problem.
146
147
Having to rely on a closed propriety system can be frustrating, since
148
it is difficult to gain access to the source code of the software,
149
either to correct a bug or include a simple optimization in an
150
algorithm. Sometimes this is by design:
151
\begin{quote}\small ``Indeed, in almost all practical uses of Mathematica,
152
issues about how Mathematica works inside turn out to be largely
153
irrelevant. You might think that knowing how Mathematica
154
works inside would be necessary [...]'' (See \cite{mathematica:internals}.)
155
\end{quote}
156
Even if we manage to contact the developers, and they find time to
157
make the changes we request, it might still take months or years
158
before these changes are made available in a new release.
159
160
Fundamental questions of correctness, reproducibility and scientific
161
value arise when building a mathematical research program on top of
162
proprietary software (see, e.g., \cite{joyner-stein:notices}). There
163
are many published refereed papers containing results that
164
rely on computations performed in Magma, Maple, or
165
Mathematica.\footnote{Including by the second author of this paper,
166
e.g., \cite{conrad-edixhoven-stein:j1p}!} In some cases, a specific
167
version of Magma is the only software that can carry out the
168
computation. This is not the infrastructure on which we want to build
169
the future of mathematical research.
170
171
In sharp contrast, open source libraries provide a great deal of
172
flexibility, since anyone can see and modify the source code as they
173
wish. However, functionality is often segmented into different
174
specialized libraries and advanced algorithms are hidden behind custom
175
interpreted languages. One often runs into trouble trying to install
176
dependencies before being able use an open source software package.
177
Also, converting the output of one package to the input format of
178
another package can present numerous difficulties and
179
introduce subtle errors.
180
181
Sage, which started in 2005 (see \cite{joyner-stein:sigsam}),
182
attacks 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}
195
%\shadowbox{%
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
217
The rest of this paper goes into more detail about Sage. In
218
Section~\ref{sec:notebook}, we describe the Sage graphical user
219
interface. Section~\ref{sec:dev} is about the Sage development
220
process, Sage days workshops, mailing lists,
221
and documentation. The subject of Section~\ref{sec:car} is the
222
sophisticated way in which Sage is built out of a wide range of open
223
source libraries and software. In Section~\ref{sec:interfaces} we
224
explain how we use Python and Cython as the glue that
225
binds the compendium of software included in Sage into a unified
226
whole. We then delve deeper into Python, Cython and the Sage
227
preparser in Section~\ref{sec:python}, and illustrate some
228
applications to mathematics in Section~\ref{sec:algebraic}. Sage is
229
actively used for research, and in
230
Section~\ref{sec:interesting} we describe some
231
capabilities of Sage in advanced areas of mathematics.
232
233
\subsection{The Notebook}\label{sec:notebook}
234
\begin{figure}
235
\begin{center}
236
\shadowbox{\includegraphics[width=.8\textwidth]{fig/sagenb}}
237
\caption{The Sage Notebook\label{fig:sagenb}}
238
\end{center}
239
\end{figure}
240
As illustrated in Figure~\ref{fig:sagenb}, the graphical user
241
interface for Sage is a web application, inspired by Google Documents
242
\cite{googledocs}, which provides convenient access to all
243
capabilities of Sage, including 3D graphics. In single user mode,
244
Sage works like a regular application whose main window happens to be
245
your web browser. In multiuser mode, this architecture allows users to
246
easily set up servers for accessing their work over the Internet as
247
well as sharing and collaborating with colleagues. One can try the
248
Sage notebook by visiting \url{www.sagenb.org}, where there are over
249
30,000 user accounts and over 2,000 published worksheets.
250
251
Users also download Sage to run it directly on their computers. We
252
track all downloads from \url{www.sagemath.org}, though there
253
are several other high-profile sites that provide mirrors of our
254
binaries. Recently, people download about 6,000 copies of Sage per
255
month directly from the Sage website.
256
257
258
\subsection{The Sage Development Process}\label{sec:dev}
259
260
There are over 200 developers from across the world who have
261
contributed to the Sage project. People often contribute because they
262
write code using Sage as part of a research project, and in this
263
process find and fix bugs, speed up parts of Sage, or want the code
264
portion of their research to be peer reviewed. Each contribution to
265
Sage is first posted to the Sage Trac server
266
\url{trac.sagemath.org}; it is then peer reviewed, and finally
267
added to Sage after all issues have been sorted out and all
268
requirements are met. Nothing about this process is anonymous; every
269
step of the peer review process is recorded indefinitely for all to
270
see.
271
272
The Sage Developer's Guide begins with an easy-to-follow tutorial that
273
guides developers through each step involved in contributing code to Sage.
274
Swift feedback is available through the \verb|sage-devel| mailing
275
list, 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
282
Much development of Sage has taken place at the Sage Days
283
workshops. There have been two dozen Sage Days
284
\cite{sagedays} and many more are planned. These are essential to
285
sustaining the momentum of the Sage project and also help ensure that
286
developers work together toward a common goal, rather than competing
287
with each other and fragmenting our limited community resources.
288
289
A major goal is ensuring that there will be many Sage Days workshops
290
for the next couple of years. The topics will depend on funding, but
291
will likely include numerical computation, large-scale bug fixing,
292
$L$-functions and modular forms, function fields, symbolic
293
computation, topology, and combinatorics. The combination of
294
experienced developers with a group of enthusiastic mathematicians at
295
each of these workshops has rapidly increased the developer community,
296
and 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
307
With the motto ``building the car instead of reinventing the wheel,''
308
Sage brings together numerous open source software
309
packages (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
315
atlas& gap& libgcrypt& palp& scipy\_sandbox\\\hline
316
blas& gd& libgpg\_error& pari& scons\\\hline
317
boehm\_gc& gdmodule& libm4ri& pexpect& setuptools\\\hline
318
boost& genus2reduction& libpng& pil& singular\\\hline
319
cddlib& gfan& linbox& polybori& sphinx\\\hline
320
cliquer& ghmm& matplotlib& pycrypto& sqlalchemy\\\hline
321
cvxopt& givaro& maxima& pygments& sqlite\\\hline
322
cython& gnutls& mercurial& pynac& symmetrica\\\hline
323
docutils& gsl& moin& python& sympow\\\hline
324
ecl& iconv& mpfi& python\_gnutls& sympy\\\hline
325
eclib& iml& mpfr& r& tachyon\\\hline
326
ecm& ipython& mpir& ratpoints& termcap\\\hline
327
f2c& jinja& mpmath& readline& twisted\\\hline
328
flint& jinja2& networkx& rubiks& weave\\\hline
329
flintqs& lapack& ntl& sagenb& zlib\\\hline
330
fortran& lcalc& numpy& sagetex& zn\_poly\\\hline
331
freetype& libfplll& opencdk& scipy& zodb3\\\hline
332
\end{tabular}
333
\end{center}
334
\end{table}
335
336
Many applications of Sage require using these
337
libraries together. Sage handles the conversion of data behind the
338
scenes, automatically using the best tool for the job, and allows
339
the user to concentrate on the problem at hand.
340
341
In the following example, which we explain in detail below, Sage
342
uses the FLINT library \cite{flint} for univariate
343
polynomials over the ring $\Z$ of integers, whereas Singular
344
\cite{singular} is used for multivariate
345
polynomials. The option to use the NTL library \cite{ntl} for
346
univariate 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
360
The 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
362
this ring. Note that $\Z$ is represented by \verb|ZZ| in Sage. The
363
expression \verb|R.<x> = ZZ[]| is not valid Python, but can be used
364
in Sage code as a shorthand as explained in Section \ref{sec:preparser}. The
365
next 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
371
the underlying implementation to use the NTL library.
372
373
Often these interfaces are used under the hood, without the user
374
having to know anything about the corresponding systems. Nonetheless, there are easy
375
ways to find out what is used by inspecting the source code, and users
376
are strongly encouraged to cite components they use in published
377
papers. The following example illustrates another
378
way to get a list of
379
components 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
393
Sage makes it possible to use a wide range of mathematical software
394
packages together by providing a unified interface that
395
handles data conversion automatically. The complexity and
396
functionality of these interfaces varies greatly, from simple
397
text-based interfaces that call external software for an individual
398
computation, to using a library as the basis for an
399
arithmetic type. The interfaces can also run code from libraries
400
written in the interpreted language of another program.
401
Table~\ref{fig:interfaces} lists
402
the 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,
410
gnuplot, gp, \\
411
& kash, lie, lisp, macaulay2, magma, maple, mathematica,
412
\\
413
& matlab, maxima, mupad, mwrank, octave, phc, polymake, \\
414
&povray, qepcad, qsieve, r,
415
rubik, 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
424
The above interfaces are the result of many years writing
425
Python and Cython \cite{cython} code to adapt
426
Singular \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
429
smoothly and efficiently in a unified way from Python \cite{python}.
430
Some of these programs were originally designed to be used only
431
through their own interpreter and made into a library by Sage
432
developers. For example libSingular was created by Martin Albrecht in
433
order to use the fast multivariate polynomial arithmetic in Singular
434
from Sage. The libSingular interface is now used by other
435
projects, including Macaulay2 \cite{M2} and GFan \cite{gfan}.
436
437
There are other approaches to linking mathematical software together.
438
The recent paper \cite{science_openmath} reports on the state of the
439
art using OpenMath. Sage takes a
440
dramatically different approach to this problem. Instead of using a
441
general string-based XML protocol to communicate with other mathematical
442
software, Sage interfaces are tailor made to the specific software and
443
problem at hand. This results in far more efficient and flexible
444
interfaces. The main disadvantage compared to OpenMath is that the
445
interfaces 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
451
Having access to many programs which can perform the same computation,
452
without having to worry about data conversion, also makes it easier to
453
double check results. For example, below we first use Maxima, an open
454
source symbolic computation package distributed with Sage, to
455
integrate a function, then perform the same computation using Maple
456
and 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
473
The most common type of interface, called a \verb|pexpect|
474
interface, communicates with another command line program by reading
475
and writing strings to a text console, as if another user was in front
476
of the terminal. Even though these are relatively simple to develop,
477
the overhead of having to print and parse strings to represent the
478
data makes this process potentially cumbersome and inefficient. This
479
is the default method of communication with most high level
480
mathematics software, including commercial and open source programs,
481
such as Maple, Mathematica, Magma, KASH or GAP.
482
483
Sage provides a framework to represent elements over these interfaces, perform
484
arithmetic with them or apply functions to the given object, as well as using a
485
file to pass the data if the string representation is too big. The following
486
demonstrates arithmetic with GAP elements.
487
488
\begin{lstlisting}
489
sage: a = gap('22')
490
sage: a*a
491
484
492
\end{lstlisting}
493
494
It is also possible to use \verb|pexpect| interfaces over remote consoles. In
495
the following code, we connect to the \verb|localhost| as a different user
496
and call Mathematica functions. Note that the interface can handle indexing
497
vectors 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
511
Sage also includes specialized libraries that are linked directly from compiled
512
code written in Cython. These are used to handle specific problems, such as the
513
characteristic 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
522
Many basic arithmetic types also use Cython to directly utilize data
523
structures from efficient arithmetic libraries, such as MPIR or FLINT. An
524
example of this can be seen at the beginning of this section, where elements
525
of the ring $\Z[x]$ are represented by the class
526
\verb|Polynomial_integer_dense_flint|.
527
528
The Singular interface is one of the most advanced included in
529
Sage. Singular has a large library of code written in its own language.
530
Previously the only way to access these functions, which include
531
algorithms for Gr\"obner basis and primary decomposition, was to call
532
Singular through a \verb|pexpect| interface, passing data back and forth
533
using strings. Recently, due to work of Michael Brickenstein and Martin
534
Albrecht, Sage acquired the ability to call these functions directly.
535
536
In 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
538
function. The interface handles the conversion of the data to Singular's format
539
and back. Since Sage already uses Singular data structures directly to represent
540
multivariate polynomials and ideals over multivariate polynomial rings, there
541
are 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
554
Efforts are under way to extend these capabilities to other programs,
555
for example to GAP which provides Sage's underlying group theory
556
functionality. Up to now, GAP was only available through its interpreter,
557
through a \verb|pexpect| interface that was written by Steve Linton. As the
558
following example demonstrates, the performance of this interface is
559
far 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
570
The 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
573
library 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}
587
The library interface is about 1,000 times faster than the pexpect interface.
588
589
590
\subsection{Python - a mainstream language}\label{sec:python}
591
592
In line with the principle of not reinventing the wheel, Sage is built on the
593
mainstream programming language Python, both as the main development language
594
and the user language. This frees the Sage developers, who are mainly
595
mathematicians, from the troubles of language design, and gives access to
596
an immense array of general purpose Python libraries and tools.
597
598
Python is an interpreted language with a clear, easy to read and learn
599
syntax. Since it is dynamically typed, it is ideal for rapid prototyping,
600
providing an environment to easily test new ideas and algorithms.
601
602
%Note that many mathematical packages, such as GINV, PolyBoRi, or CVXOPT,
603
%already have a Python interface.
604
605
\subsubsection{A fast interpreter}
606
607
In the following Singular session, we first declare the ring ${\tt
608
r}=\Q[x,y,z]$ and the polynomial ${\tt f} \in {\tt r}$, then
609
measures 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
620
The elapsed time is 990 milliseconds. Next we use Sage to do the same
621
computation, using the same Singular data structures directly, but
622
without 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
632
Sage takes only 91.8 milliseconds for the same operation. This difference is
633
because the Python interpreter is more efficient at performing \verb|for| loops.
634
635
636
637
\subsubsection{Cython - compiled extensions}
638
639
Python alone is too slow to implement a serious mathematical software
640
system. Fortunately, Cython \cite{cython}
641
makes it easy to optimize parts of your
642
program or access existing C/C++ libraries. It can translate Python
643
code with annotations containing static type information to C/C++
644
code, which is then compiled as a Python extension module.
645
646
Many of the basic arithmetic types in Sage are provided by Cython wrappers of
647
C libraries, such as FLINT for univariate polynomials over $\Z$, Singular for
648
multivariate polynomials, and Pynac for symbolic expressions.
649
650
The 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
664
Here is the same function, but the loop index {\tt k} is declared to be
665
a 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
680
The code is compiled and linked to the interpreter on the fly, and the
681
function \verb|mysum_cython| is available immediately. Note that the run
682
time for the Cython function is 60 times faster than the Python equivalent.
683
684
Cython also handles the conversion of Python types to C types automatically.
685
In the following example, we call the C function \verb|sinl| using Cython
686
to 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
703
Note 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
705
transparently by Cython.
706
707
\subsubsection{The Preparser}\label{sec:preparser}
708
709
While Python has many advantages as a programming and glue language, it also
710
has some undesirable features. Sage hides these problems by using a preparser
711
to change the commands passed to Python in an interactive session (or when running
712
a script with the \verb|.sage| extension). In order to
713
maintain compatibility with Python, changes performed by the preparser are
714
kept to a minimum. Moreover, the Sage library code is not preparsed, and is
715
written in Cython or Python directly.
716
717
Python, like C and many other programming languages, performs integer floor
718
division. 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
720
notebook with its own type declarations, which behave as expected with respect
721
to arithmetic and have the advantage that they are backed by efficient
722
multiprecision arithmetic libraries such as MPIR \cite{mpir} and MPFR \cite{mpfr}, which are thousands of times faster than Python for large
723
integer arithmetic.
724
725
To 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
735
Adding a trailing {\tt r} after a number indicates that the preparser should leave
736
that 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}
744
Here 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
754
The preparser also changes the \verb|^| sign to the exponentiation operator
755
\verb|**| and provides a shorthand to create new mathematical domains and
756
name 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
768
Sage combines algebraic, symbolic and numerical computation tools under one
769
roof, enabling users to choose the tool that best suits the problem. This
770
combination also makes Sage more accessible to a wide audience---scientists,
771
engineers, pure mathematicians and mathematics teachers
772
can all use the same platform for scientific computation.
773
774
While not concentrating on only one of these domains might seem to divide
775
development resources unnecessarily, it actually results in a better overall
776
experience for everyone, since users do not have to come up with makeshift
777
solutions to compensate for the lack of functionality from a different field.
778
Moreover, because Sage is a distributed mostly-volunteer open source
779
project, widening our focus results in substantially more developer resources.
780
781
\subsubsection{Algebraic Tools: The Coercion System}
782
783
An algebraic framework, similar to that of Magma or Axiom, provides access to
784
efficient data structures and specialized algorithms associated to particular
785
mathematical domains. The Python language allows
786
classes to define how arithmetic
787
operations like \verb|+| and \verb|*| will be handled, in a
788
similar way to how C++ allows overloading of operators. However, the built-in support
789
for overloading in Python is too simple
790
to support operations with a range of objects in a mathematical type hierarchy.
791
792
Sage abstracts the process of deciding what an arithmetic operation
793
means, or equivalently, in which domain the operation should be
794
performed, in a framework called the {\em coercion system}, which was
795
developed and implemented by Robert Bradshaw, David Roe,
796
and many others. Implementations of new mathematical objects only need
797
to define which other domains have a natural embedding to their
798
domain. When performing arithmetic with objects, the
799
coercion system will find a common domain where both arguments can be
800
canonically mapped, perform the necessary type conversions automatically,
801
thus allowing the implementation to only handle the case where both
802
objects have the same parent.
803
804
In 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
806
first deduces that the result should be in $\Q$ from the fact that \verb|t|
807
can be converted to the domain of \verb|u|, namely $\Q$, but canonical conversion
808
in the other direction is not possible. Then the addition is performed with
809
both 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
824
Similarly, in the following example, the common domain $\Q[x]$ is found for
825
arguments from $\Z[x]$ and $\Q$. Note that in this case, the result is not in
826
the 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
863
Another abstraction to make implementing mathematical
864
structures easier is the {\em category framework}, whose development
865
was spearheaded by Nicolas Thi\'ery and Florent Hivert. Similar in spirit to the
866
mathematical programming facilities developed in Axiom and encapsulated
867
in Aldor, the category framework uses Python's dynamic class
868
creation capabilities to combine functions relevant for a mathematical object,
869
inherited through a mathematical hierarchy, into a class at run time.
870
871
This process greatly simplifies the troubles of having to combine object-oriented programming concepts with mathematical structural concerns, while
872
keeping efficiency in mind. Efficient implementations can keep the
873
inheritance hierarchy imposed by the data structures, while generic methods
874
to compute basic properties are implemented in the {\em category} and
875
automatically 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
882
The symbolic subsystem of Sage provides an environment similar to Maple or
883
Mathematica, where the input is treated only as an expression
884
without 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
899
Sage uses Pynac \cite{pynac}, a hybrid C++ and Cython library built on top of GiNaC \cite{ginac},
900
to work with symbolic expressions. High level symbolic calculus problems
901
including symbolic integration, solution of differential equations and Laplace
902
transforms are solved using Maxima behind the scenes.
903
904
Here is an example of how to use the symbolic computation facilities in Sage.
905
Note that in contrast to other symbolic software such as Maple, variables must
906
be 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}
973
In addition to code for symbolic computation, the standard numerical
974
Python packages NumPy, SciPy, and Matplotlib are included in Sage,
975
along with the numerical libraries cvxopt, GSL, Mpmath, and R.
976
977
For numerical applications, Robert Bradshaw and Carl Witty developed
978
a compiler for Sage that converts
979
symbolic expressions into an internal format suitable for blazingly fast
980
floating 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}
991
The \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}
996
Performance is typically within a factor of
997
two from what one gets using a direct implementation in C or Fortran.
998
999
1000
\section{Afterword}\label{sec:interesting}
1001
1002
In this article, we have showed that
1003
Sage is a powerful platform for developing sophisticated mathematical
1004
software. Sage is actively used in research mathematics,
1005
and people use Sage to develop state-of-the-art algorithms.
1006
Sage is particularly strong in number theory, algebraic combinatorics,
1007
and graph theory. For further examples, see the
1008
53 published articles, 11 Ph.D. theses, 10 books, and 30 preprints at
1009
\url{www.sagemath.org/library-publications.html}.
1010
1011
For example, Sage has extensive functionality for computations related to the
1012
Birch and Swinnerton-Dyer conjecture. In addition to Mordell-Weil
1013
group computations using \cite{mwrank} and point counting over large
1014
finite fields using the SEA package in \cite{pari}, there is much
1015
novel elliptic curve code written directly for Sage. This includes the
1016
fastest 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
1019
multiplicative primes. Sage combines these capabilities to
1020
compute explicitly bounds on Shafarevich-Tate groups of elliptic
1021
curves \cite{stein-wuthrich}. Sage also has code for computation
1022
with modular forms, modular abelian varieties, and ideal class groups in
1023
quaternion algebras.
1024
1025
The MuPAD-combinat project, which was started by Florent Hivert and
1026
Nicolas M. Thi\'ery in 2000, built the world's preeminent system for
1027
algebraic combinatorics on top of MuPAD (see \cite{icms:mupad2006} and
1028
\cite{mupad-hivert-thiery}). Page 54 of \cite{mupad-hivert-thiery}:
1029
``They [MuPAD] also have promised to release the code source of the
1030
library under a well known open-source license, some day.'' In 2008,
1031
MuPAD was instead purchased by MathWorks (makers of MATLAB), so MuPAD
1032
is no longer available as a separate product, and will probably never
1033
be 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.''
1043
%-- See \cite{matlabmupad}. }
1044
%\end{quote}
1045
1046
As a result, the MuPAD-combinat group has spent several years
1047
reimplementing everything in Sage (see \cite{sagecombinatroadmap} for
1048
the current status). The MuPAD-combinat
1049
group was not taken by surprise by the failure of MuPAD, but instead
1050
were concerned from the beginning by the inherent risk in building
1051
their research program on top of MuPAD. In fact, they decided to
1052
switch to Sage two months before the bad news hit, and have made
1053
tremendous 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