% for book1\documentclass{report}2\voffset=-0.05\textheight \textheight=1.1\textheight3\hoffset=-0.05\textwidth \textwidth=1.1\textwidth45% for in-class presentation6%\documentclass[landscape,12pt]{report}7%\voffset=-0.3\textheight \textheight=1.5\textheight8%\hoffset=-0.35\textwidth \textwidth=1.7\textwidth91011\bibliographystyle{amsalpha}12\include{macros}13\DeclareMathOperator{\Li}{Li}1415\renewcommand{\todo}[1]{}161718\usepackage{pdfsync}19\usepackage[hypertex]{hyperref}2021\usepackage{graphicx}22\newcommand{\sagegraphic}[2]{\begin{center}\includegraphics[width=#1\textwidth]{graphics/#2}\end{center}}2324\title{Sage for Power Users}2526\author{William Stein}2728\begin{document}29%\sf\LARGE % for in-class presentation3031\maketitle32\tableofcontents3334%\part{Sage, Python and Cython}3536\chapter*{Preface}37This is a book about Sage \url{http://sagemath.org}, which is a large38free open source software project that I started in 2005, whose39``mission statement'' is to create a viable free open source40alternative to the commercial programs Magma, Maple, Mathematica, and41Matlab. I have given many talks, tutorials, and workshops on Sage,42and this book records what I have personally found to be the most43important key ideas that are needed to make effective use of Sage. My44intention is that you read the whole book cover-to-cover, and have45thus kept the book intentionally short.464748% Background49I assume that you have some previous computer programming experience,50but not necessarily in Python. Though I'll make reference to many51mathematical topics when illustrating how to use Sage, do not worry if52some are not familiar to you.5354% License55This book is licensed under the Creative Commons Attribution 3.056license\footnote{See \url{http://creativecommons.org/licenses/by/3.0/}.}, so57it will always be freely available in many formats.5859% Overview of book60\todo{Put an overview of the book -- chapter by chapter -- here.}616263\chapter{Introduction to Sage}64\section{Motivation}6566I started the Sage project in early 2005 in order to create a viable67free open source mathematical software package that I could use for my68research. I was frustrated with not being allowed to easily change or69understand the internals of closed source systems\footnote{For me,70this was a powerful niche program called ``Magma''.} and I had a71deep concern that my students and colleagues could not easily use the72commercially distributed software that I had spent many years73developing (and contributing). I started Sage as a new project74instead of switching to another system, since the capabilities of any75available software for number theory at the time were far behind76many key features of commercial systems.\footnote{For example, Magma's tools77for linear algebra over the rational numbers and finite fields were78vastly superior to anything available anywhere else.} Several79hundred people have since become involved in Sage development, and the80goals have broadened substantially.8182Sage uses a mainstream programming language, unlike all popular83mathematics software, including Maple,84Mathematica, and Matlab, which each use their own85special-purpose languages written just for mathematics. One works86with Sage using Python, which is one of the world's most popular87general purpose scripting languages. By using Python, one can use88almost anything ever written in Python directly in Sage. And89there is much useful Python code out there that addresses a wide90range of application areas.9192Instead of writing many of the core libraries from scratch like most93math software systems have done in Sage I assembled together the best94open source software out there, and built on it\footnote{Sage includes95over 500,000 lines of code that does not come from third party96projects.}. Also, the complete system is easily buildable from97source on a range of computers. There are challenges: some of the98upstream libraries can be difficult to understand, are written in a99range of languages, and have different conventions than Sage. Thus it100is important to strongly encouraging good relations with the projects101that create many of the components of Sage.102103A wide and vibrant community of developers and users have become104involved with Sage. Due to the broad interests of this large105community of developers, Sage has grown into a project with the106following specific goal:107\begin{quote}108{\bf Mission Statement:} Provide a viable free open source alternative109to Magma, Maple, Mathematica, and Matlab.110\end{quote}111%The Sage development model has matured. There are regular releases of112%Sage, and all code that goes into Sage is peer reviewed.113114\section{What is Sage?}115116Sage is a free open-source mathematics software system licensed under117the GNU Public License (GPL). It combines the power of about 100118open-source packages with a large amount of new code to provide a free119open source platform for mathematical computation. Sage has many120notable features.121122\begin{itemize}123\item Sage is free, due mainly to the volunteer effort of hundreds of124people and generous funding from the National Science Foundation,125private donations, and other organizations such as Google and126Microsoft. There are no license codes or copy protection. Sage is127also open source, so there are absolutely no secret or proprietary128algorithms anywhere in Sage. There is nothing that you are not129allowed to see or change.130131\item Sage uses the mainstream programming language Python. Learning132Sage will make you proficient in this popular, widely used, and well133supported free programming language, which you will likely also use134for other non-mathematics projects. Moreover, Sage features the135Cython compiler, which allows one to combine Python, C/C++/Fortran136libraries, and native machine types for potentially huge speedups.137138\item Sage is uniquely able to combine functionality from dozens of139other mathematical software programs and programming languages via140smart psuedoterminal interfaces. You can combine Lisp, Mathematica,141and C code to attack a single problem.142143\item Sage has both a sophisticated multiuser web-based graphical user144interface and a powerful command line interface. Sage can also be145made to work with any other Python interactive development environment (IDE).146147\item Sage may have the widest range of mathematical capabilities of148any single mathematical software system available. Sage and its149components are developed by an active and enthusiastic worldwide150community of people from many areas of mathematics, science, and151engineering.152153\item Modifications to Sage are publicly peer reviewed, and what goes154into Sage is decided via community discussions; no matter who you155are, if you have a brilliant idea, the energy, and can clearly argue156that something should go into Sage, it probably will. Known bugs in157Sage, and all discussions about them are available for all to see.158\end{itemize}159160Sage is nothing like Magma, Maple, Mathematica, and Matlab, in which161details of their implementations of algorithms is secret, their list162of bugs is concealed, how they decided what got included in each163release is under wraps, their custom programming language locks you164in, and you must fight with license codes, copy protection and165intentionally crippled web interfaces.166167\section{``This unique American idea of the168entrepreneurial company.''}169170The Mathematica documentation has an argument for why looking at the internals171of mathematical software is not necessary.172\begin{quote}173``Particularly in more advanced applications of Mathematica, it may174sometimes seem worthwhile to try to analyze internal algorithms in175order to predict which way of doing a given computation will be the176most efficient. And there are indeed occasionally major improvements177that you will be able to make in specific computations as a result178of such analyses.179180But most often the analyses will not be worthwhile. For the181internals of Mathematica are quite complicated, and even given a182basic description of the algorithm used for a particular purpose, it183is usually extremely difficult to reach a reliable conclusion about184how the detailed implementation of this algorithm will actually185behave in particular circumstances.''186187-- {\tiny \url{http://reference.wolfram.com/mathematica/tutorial/WhyYouDoNotUsuallyNeedToKnowAboutInternals.html}}188\end{quote}189190Wolfram, who founded the company that sells Mathematica, admits that191the mathematical community hates some of what he has done, arguing192that a closed source commercial model is the only approach that193can possibly work.194195\begin{quote}196``There's another thing, quite honestly, that that community has a197hard time with. They sort of hate one aspect of what I have done,198which is to take intellectual developments and make a company out of199them and sell things to people.200201My own view of that, which has hardened over the years, is, my god,202that's the right thing to do. If you look at what's happened with203TeX, for example, which went in the other direction... well,204Mathematica could not have been brought to where it is today if it205had not been done as a commercial effort. The amount of money that206has to be spent to do all the details of development, you just can't207support that in any other way than this unique American idea of the208entrepreneurial company.''209210-- Stephen Wolfram, 1993, Doctor Dobbs Journal Interview211\end{quote}212213214\noindent{}For the last 20 years, Matlab, Mathematica, and the other215commercial systems have dominanted with on the order of a hundred216million dollars a year in revenue. If the Sage project succeeds at217its goals (still a big if), it will have proved that Wolfram is wrong218and radically change the landscape of computational mathematics.219220\section{Getting Started}221The easiest way to get started with Sage right {\em now} is to visit222\url{http://480.sagenb.org} and login using OpenID by clicking one of223the buttons at the bottom right. This should work with nearly any224operating system and browser combination\footnote{I recommend that you225avoid using Internet Explorer if at all possible.}. Using Sage via226the above webpage is fine if you just want to use Sage via the227notebook, e.g., for learning Python (Chapter~\ref{ch:python}) and228Cython (Chapter~\ref{ch:cython}).229230231There are some situations where you will instead want to install Sage232on your own computer, or get an account on a command-line server on233which Sage is installed:234235\begin{enumerate}236\item You want to use the Sage command line interface.237\item You want to use the interactive command line profiler and238debugger, which haven't been properly ported to the notebook yet239(see Chapter~\ref{ch:profdeb}).240\item You want to modify Sage and contribute back new code (see241Chapter~\ref{ch:contrib}).242\item You want to interface nonfree software with Sage (see243Chapter~\ref{ch:interfaces}). It would be illegal for me to allow244just anybody to run Maple/Mathematica/etc. code at245\url{http://480.sagenb.org}.246\item You do not have access to the Internet.247248\end{enumerate}249250251\begin{remark}252Eliminating all but the last reason above are current goals of the253Sage project. A command line interface should be added to the254notebook, and it should support the profiler and debugger. It255should be possible to edit all files in the source code of Sage, use256revision control systems, etc., completely via the web. Even the257legal issue involving nonfree software could be resolved by hooking258into our University's authentication system, just as you259authenticate for off-campus access to library resources.260\end{remark}261262\section{A Tour}263264Sage uses the basic user-interface principle of ``question and265answer'' found in many other mathematical software systems. You enter266input written in a well-defined language and, after pressing the267\fbox{\sf return} key in the command line interface or pressing268\fbox{\sf shift+return} in the notebook interface, Sage evaluates your269input and returns the result.270271A traditional test that Sage is working is to compute 2+2:272\begin{lstlisting}273sage: 2 + 22744275\end{lstlisting}276277We factor a whole number.278\begin{lstlisting}279sage: factor(2012)2802^2 * 503281\end{lstlisting}282Thus $2012 = 2\times 2 \times 503$.283Sage can also factor negative numbers and rational numbers:284\begin{lstlisting}285sage: factor(-2012/2015)286-1 * 2^2 * 5^-1 * 13^-1 * 31^-1 * 503287\end{lstlisting}288289The language that Sage uses is almost the same as the290Python programming language.291One difference between Sage and Python is that \verb|^| means292exponentiation in Sage but exclusive or in Python.293Another difference is that integer division results in a rational294number in Sage, but is floor division in Python.295\begin{lstlisting}296sage: 2^32978298sage: 2012/62991006/3300\end{lstlisting}301302We can also factor symbolic expressions using Sage. To introduce a303symbolic variable, use the {\sf var} command.304\begin{lstlisting}305sage: var('x,y')306(x, y)307sage: F = factor(x^2 - 4*sin(y)^2); F308(x - 2*sin(y))*(x + 2*sin(y))309\end{lstlisting}%link310311If you want to put any result in a \LaTeX{} document\footnote{\LaTeX{}312is the dominant tool for producing {\em professional} quality mathematical313papers and books; it is free and open source and you should learn314it.}, use the {\sf latex} command:315%link316\begin{lstlisting}317sage: latex(F)318{\left(x - 2 \, \sin\left(y\right)\right)} {\left(x + 2 \,319\sin\left(y\right)\right)}320\end{lstlisting}321which looks like this:322$$323{\left(x - 2 \, \sin\left(y\right)\right)} {\left(x + 2 \,324\sin\left(y\right)\right)}325$$326327Sage knows Calculus:328\begin{lstlisting}329sage: integrate(e^x * sin(x), x)3301/2*(sin(x) - cos(x))*e^x331sage: derivative(1/2*(sin(x) - cos(x))*e^x).expand()332e^x*sin(x)333\end{lstlisting}334335Sage can plot functions:336\begin{lstlisting}337sage: g = plot(sin(x) + (1-x^2), (x, 0, 2)); g338\end{lstlisting}%link339\sagegraphic{.4}{plot1}340To include this plot in a document, save it as a PDF file:341%link342\begin{lstlisting}343sage: g.save('plot1.pdf')344\end{lstlisting}345346We numerically find a root of $\sin(x) + (1-x^2)$ between 0 and 2, as follows:347\begin{lstlisting}348sage: find_root(sin(x) + (1-x^2), 0, 2)3491.4096240040025754350\end{lstlisting}351352You can use some other programming languages directly from Sage, such as Lisp:353\begin{lstlisting}354sage: s = "(defun factorial(n)"355sage: s += " (if (= n 1) 1 (* n (factorial (- n 1)))))"356sage: lisp(s)357FACTORIAL358sage: lisp('(factorial 10)')3593628800360\end{lstlisting}361362Or use Mathematica (this won't work if you don't have Mathematica):363\begin{lstlisting}364sage: mathematica('Integrate[Sin[x^2],x]') # optional - Mathematica365Sqrt[Pi/2]*FresnelS[Sqrt[2/Pi]*x]366\end{lstlisting}367368Or use Magma, over the web (this should work as long as you have an369Internet connection, since it just uses370\url{http://magma.maths.usyd.edu.au/calc/}):371\begin{lstlisting}372sage: magma_free("Factorisation(2012)")373[ <2, 2>, <503, 1> ]374\end{lstlisting}375376% \chapter{Using Sage: The Notebook and Command Line}377378379% \section{The Notebook}380381382383% \subsection{Documentation}384385% \subsection{Creating interacts}386387% \section{The Command Line}\label{sec:line}388389\part{Programming Sage}390391\chapter{Python}\label{ch:python}392393Sage uses the Python programming language, which is relatively easy to394learn and fun to use. This chapter is a quick Sage-oriented395introduction to Python, which you should supplement with a book.396Fortunately, the two best books on Python are free: {\em The Python397Tutorial} (see \url{http://docs.python.org/}) and {\em Dive Into398Python} (see \url{http://www.diveintopython.net/}).399%The best non-free Python reference manual is {\em Python in a Nutshell} (see400%\url{http://oreilly.com/catalog/9780596001889}).401402\section{What is Python?}403Python is a popular free open source language with {\em no} particular404company pushing it. Many big companies such as Google use Python, and405support its development. From \url{http://python.org}:406\begin{quote}407``Python is a programming language that lets you work more quickly408and integrate your systems more effectively. You can learn to use409Python and see almost immediate gains in productivity and lower410maintenance costs.''411\end{quote}412\begin{itemize}413\item {\em Work more quickly}: you get stuff done instead of414fighting with the language and environment for silly reasons415\item {\em Integrate your systems}: Python is particular useful at416creating big systems out of possibly messy collections of software417tools.418\item {\em Maintenance costs}: Python code is more likely to be readable419and hackable.420\end{itemize}421422Sage is a big integrated system built out of {\em several million}423lines of possibly messy software, code written using Sage tends to424be readable and hackable, and people use Sage since it helps them425get stuff done immediately.426427\section{The Sage Preparser}428When you type commands into Sage, the computer programming language429you use is (almost) Python. Each line of code gets automatically run430through a preparser before it is sent to the Python interpreter. To431see exactly what changes occur, use the {\tt preparse} command:432\begin{lstlisting}433sage: preparse('a = 2.5^3')434"a = RealNumber('2.5')**Integer(3)"435\end{lstlisting}436As you can see, decimal literals get wrapped using the {\tt437RealNumber} command, so when you type 2.5, Python will see {\tt438RealNumber('2.5')}. Similarly, integer literals get wrapped using439{\tt Integer}. Finally, the caret symbol \verb|^| is replaced by {\tt440**}, which is Python's exponentiation operator. One motivation for441doing all this is that in Magma, Maple, Mathematica, Matlab and \LaTeX the442\verb|^| operator is exponentiation, and making Sage have the same443behavior as those systems helps minimize confusion (whereas in Python \verb|^| is444``exclusive or''). The preparse does a few other things, but not much445more.446447If you want to turn off the preparser, type {\tt preparser(False)}:448%skip449\begin{lstlisting}450sage: preparser(False)451sage: 2/3 + 2^34521453\end{lstlisting}454\begin{lstlisting}455sage: preparser(True)456sage: 2/3 + 2^345726/3458\end{lstlisting}459460Read more about the preparser at \url{http://www.sagemath.org/doc/reference/sage/misc/preparser.html}.461462463\section{Variables}464%http://wiki.wstein.org/10/480b/lectures/lec2465466In Python you create a variable by writing {\tt var = expression}; for example,467\begin{lstlisting}468sage: a = 2469sage: b = 3470sage: a + b4715472\end{lstlisting}473474You can also include several assignment statements on the same line if you separate them475with a semicolon:476\begin{lstlisting}477sage: c = 7; d = 15; e = 5478sage: c + d + e47927480\end{lstlisting}481You do {\em not} have to end lines with a semicolon.482You can also assign the same value to several variables at once:483\begin{lstlisting}484sage: c = d = 10485sage: c + d48620487\end{lstlisting}488489We have only used integers as the expression, but490Python supports many other types of objects, such as lists, which491we make using square brackets:492\begin{lstlisting}493sage: v = [7, 15, 5]; v494[7, 15, 5]495\end{lstlisting}496497The most important gotcha (and feature) of variables in Python is that498variables are a {\em reference} to a Python object, not a new copy of499that object. Thus in the example below, both {\tt v} and {\tt w}500``reference'' exactly the same Python list:501\begin{lstlisting}502sage: v = [1, 2, 3]503sage: w = v504sage: w[0] = 10505sage: v506[10, 2, 3]507\end{lstlisting}508Continuing the above example:509\begin{lstlisting}510sage: v[1] = 5511sage: w512[10, 5, 3]513\end{lstlisting}514515If you want a copy of an object, use the \verb|copy| command.516\begin{lstlisting}517sage: v = [1,2,3]518sage: w = v519sage: z = copy(w)520sage: v[0] = 10521sage: print w522[10, 2, 3]523sage: z524[1, 2, 3]525\end{lstlisting}526527The {\tt copy} function only copies the {\em references} in the list:528\begin{lstlisting}529sage: v = [[1,2], 3, 4]530sage: w = copy(v)531sage: w[1] = 10532sage: w[0][0] = 5533sage: v534[[5, 2], 3, 4]535sage: w536[[5, 2], 10, 4]537\end{lstlisting}538539To recursively make a new copy of everything (as much as possible), use540{\tt deepcopy}:541\begin{lstlisting}542sage: v = [[1,2], 3, 4]543sage: w = deepcopy(v)544sage: w[1] = 10; w[0][0] = 5545sage: v546[[1, 2], 3, 4]547sage: w548[[5, 2], 10, 4]549\end{lstlisting}550551You probably won't have to use {\tt deepcopy} often. In over 500,000552lines of code in the core Sage library, deepcopy is used around 177553times:554\begin{lstlisting}555sage -grep deepcopy |wc -l556177557\end{lstlisting}558559The main reason many people are very confused by variables being560references in Python is that most other mathematical software561(including both Mathematica and Matlab) works differently. For562example, in Matlab assignment to a variable creates a new copy. For563example, noting that arrays in Matlab are 1-based instead of 0-based,564\begin{lstlisting}565$ matlab566>> v = [1,2,3]567v =5681 2 3569>> w = v570w =5711 2 3572>> v(1) = 10573v =57410 2 3575>> w576w =5771 2 3578\end{lstlisting}579580And in Mathematica,581\begin{lstlisting}582In[27]:= v = {1,2,3};583In[28]:= w = v;584In[29]:= v[[1]] = 10;585In[30]:= v586Out[30]= {10, 2, 3}587In[31]:= w588Out[31]= {1, 2, 3}589\end{lstlisting}590591But of course in Sage:592\begin{lstlisting}593sage: v = [1,2,3]594sage: w = v595sage: v[0] = 10596sage: v597[10, 2, 3]598sage: w599[10, 2, 3]600\end{lstlisting}601602\begin{remark}603Another subtle difference in various computer604languages is that exponentiation may associate either605left to right or right to left. For example,606\begin{lstlisting}607sage: 3^3^36087625597484987609\end{lstlisting}610But in Matlab, we have611\begin{lstlisting}612$ matlab613>> 3^3^361419683615\end{lstlisting}616Finally, in Maple we have617\begin{lstlisting}618> 3^3^3619syntax error, ambiguous use of `^`, please use parentheses:620\end{lstlisting}621Thus watch out: of the two possible design choices about the622meaning of \verb|3^3^3|, we quickly find three design decisions623made in practice!624\end{remark}625626Like in Magma, Maple, Matlab, and Mathematica, you do not have to627explicitly declare the type of a variable, and it can have several628different types in a single snippet of code. This is629different to the situation with C, C++ and Java\footnote{Sage630is ``dynamically typed'', whereas C, C++ and Java are631``statically typed''.}.632Use the {\tt633type} function to determine the type of a variable,634and the {\tt id} function to find out the memory location635of the object that a variable references.636\begin{lstlisting}637sage: a = 10638sage: type(a)639<type 'sage.rings.integer.Integer'>640sage: id(a) # random; memory location a points at6414468006416642sage: a = "hello world"643sage: type(a)644<type 'str'>645sage: id(a) # random; new memory location a now points at6464507478816647\end{lstlisting}648649\section{Control Flow}650The basic control flow statements in Python are {\tt if}, {\tt while} and {\tt for}.651The {\bf if} statement lets you choose between alternative code at runtime.652Here is an example:653%skip654\begin{lstlisting}655a = 2; b = 3656if a > b:657print(1)658print("-----")659elif a == b:660print(2)661else:662print(3)663\end{lstlisting}664The Python interpreter evaluates the expression right after {\tt if}665and before the colon, and if it evaluates to {\tt True}, then all of666the code that is indented before the {\tt elif} or {\tt else} is667executed. Otherwise, the expression right after {\tt elif} is668evaluated, and if {\tt True}, then the indented code directly below it669is evaluated. Otherwise the code under the final {\tt else} is670evaluated. The {\tt elif} and {\tt else} are optional, and you can have671any number of {\tt elif} blocks.672673674Unlike nearly every other programming language, there are no explicit675begin and end markers around the block of code that will get evaluated676when a branch of the if statement is satisfied. Instead the code is677indented. There are at least two advantages to Python's choice:678(1) you will type and read less, and679(2) you will not be fooled by680misleading indentation in code like this C code:681\begin{lstlisting}682if (a>b)683printf("1");684printf("-----");685\end{lstlisting}686687Python's {\tt while} statement repeatedly executes the code indented688below it until the expression between the {\tt while} and the colon689evaluates to False, or until an explicit break statement is executed.690Here is an example:691%skip692\begin{lstlisting}693i = 5694while i > 0:695print(i)696i = i - 1697if i == 20:698break699\end{lstlisting}700When you evaluate this code, you'll see the following output:701%skip702\begin{lstlisting}70357044705370627071708\end{lstlisting}709Each time the indented block of code is executed,710the number $i$ is printed out, then the line {\tt i = i - 1} replaces711{\tt i} by an integer that is one smaller. Once $0$ is reached, the while712loop terminates.713714If instead, we set {\tt i = 25} at the top, and evaluate the code, we see:715%skip716\begin{lstlisting}7172571824719237202272121722\end{lstlisting}723This is because the {\tt if} statement evaluates to {\tt True} once724$i$ hits $20$, and the {\tt break} statement causes725the while loop to terminate.726727Use the Python {\tt for} loop to iterate over each element in a728list or any other ``iterable'' object.729For example,730%skip731\begin{lstlisting}732for i in [1, 2, 3, 4, 5]:733print(i, i*i)734\end{lstlisting}735will make a table of squares:736%skip737\begin{lstlisting}738(1, 1)739(2, 4)740(3, 9)741(4, 16)742(5, 25)743\end{lstlisting}744You can also use {\tt break} in a {\tt for} loop.745746There are many ways to make lists to iterate over747(see Section~\ref{sec:lists}), for example:748\begin{lstlisting}749sage: range(10)750[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]751sage: range(5,20)752[5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]753sage: [1..10]754[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]755sage: [n^2 for n in [1..10]]756[1, 4, 9, 16, 25, 36, 49, 64, 81, 100]757sage: [1,3,..,10]758[1, 3, 5, 7, 9]759sage: xrange(10^10,10^10+10^9) # a "lazy" list760xrange(10000000000, 11000000000)761\end{lstlisting}762763For example,764%skip765\begin{lstlisting}766for i in xrange(10^10, 10^10+10^9):767print(i)768if i > 10^10 + 5: break769\end{lstlisting}770results in771\begin{lstlisting}77210000000000773100000000017741000000000277510000000003776100000000047771000000000577810000000006779\end{lstlisting}780781\section{Functions}782%http://uw.sagenb.org/home/pub/28/783784Use \verb|def| to define a function in Python.785%skip786\begin{lstlisting}787def foo(a, bar, w=10):788if a:789print bar790# a block of code that is indented791print a, bar, w792\end{lstlisting}793The syntax is similar to the syntax of {\tt if}, {\tt for}, and {\tt794while}: a keyword, something, a colon, then an indented block of795code that gets executed under certain circumstances. More precisely,796define a function put {\tt def}, then the name of the function, then797in parenthesis the inputs to the function with possible default values798(e.g., {\tt w=10} above makes {\tt w} default to {\tt 10} if {\tt w}799is not specified). When the function is called, the input variables800to the function are set to reference the inputs, and the code in the801body of the function is executed.802803\begin{lstlisting}804sage: foo(1, 'abc', 5)805abc8061 abc 5807sage: foo(1, 'xyz')808xyz8091 xyz 10810\end{lstlisting}811812You can explicitly specify the input variables as follows,813which can make reading your code later easier:814\begin{lstlisting}815sage: foo(bar='gold', a=False, w=3)816False gold 3817\end{lstlisting}818819Unlike the situation with C/C++/Java, there is {\bf absolutely no way}820in the Python language to explicitly declare that the types of the821inputs and outputs of a function. You can put constraints on types822explicitly using the {\tt isinstance} function, or using decorators823(see Section~\ref{sec:decorators}).824\begin{lstlisting}825def g(a, b):826if not isinstance(a, Integer):827raise TypeError828return a+b829\end{lstlisting}830Then we have:831\begin{lstlisting}832sage: g(2, 3)8335834sage: g('sage', 'math')835Traceback (click to the left of this block for traceback)836...837TypeError838\end{lstlisting}839840Returning to the function {\tt foo} defined above, it will841just work with any inputs for which {\tt +} is defined.842For example, {\tt a}, {\tt b}, and {\tt c} could be strings or lists:843\begin{lstlisting}844sage: foo('a', 'b', 'c')845a= a b= b c= c846'abc'847sage: f([1,2], [3,4], [5,6])848a= [1, 2] b= [3, 4] c= [5, 6]849[1, 2, 3, 4, 5, 6]850\end{lstlisting}851Thus illustrates something in Python called ``duck typing''. So long852as an object quacks like a duck (in our case, something that supports853+), then we just treat it like a duck. In this sense, all Python functions854are extremely generic.855856Any variables that are created in the body of the function are857{\em local to the function}, unless you explicitly use the {\tt858global} keyword. For example,859860%skip861\begin{lstlisting}862c = 1; d = 1863def bar(a, b):864global d865c = a; d = b866print c, d867\end{lstlisting}868869When we call {\tt bar}, the global variable {\tt d} gets changed, but870the global {\tt c} does not change:871%skip872\begin{lstlisting}873sage: bar(5, 10)8745 10875sage: print c, d8761 10877\end{lstlisting}878879You can also have functions nested within functions (etc.), where the880nested function is completey hidden within the scode of the function881that contains it:882\begin{lstlisting}883c = 1; d = 1884def bar(a, b):885global d # this is a rare beast.886c = a; d = b887print c, d888def foo(x, y):889c = 'fun'; d = 'stuff'890print c, d891foo(c,d)892print c,d893\end{lstlisting}894Running this, note that the global c is not changed, and locally895within {\tt foo} we have yet another pair of variables also called896{\tt c} and {\tt d} that have nothing to do with the global {\tt c},897{\tt d} or the {\tt c}, {\tt d} defined at the top of {\tt bar}.898\begin{lstlisting}899sage: bar(5,10)9005 10901fun stuff9025 10903sage: c,d904(1, 10)905\end{lstlisting}906907908As illustrated above, a Python function can have side effects, and909behave differently depending on global variables. Thus Python910``functions'' are different than the functions $f:X\to Y$ in911mathematics. In mathematics, $f(x)$ depends only on $x$, not on the912state of some global variable, the time of day, phase of the moon,913etc, but in Python {\tt f(x)} can depend on more than just {\tt x}.914The following Python function evaluates to $x^2$ when the number of915seconds since the beginning is even, and $x^3$ when it is odd:916917%skip918\begin{lstlisting}919def f(x):920import time921if int(time.time()) % 2 == 0:922return x^2923else:924return x^3925\end{lstlisting}926927Here we imported a Python module called {\tt time} using the Python928{\tt import} command.929In case you are wondering ``what the heck is 'time''', you can type930\begin{lstlisting}931sage: import time932sage: help(time)933\end{lstlisting}934into Sage. In the notebook, you'll get a link to a webpage all about935the {\tt time} Python module. Python includes an {\em enormous}936standard library of modules, and you should read all about them at937\url{http://docs.python.org/library/}. I have more than once938reimplemented functionality that is already in one of the standard939modules because I didn't think to look at the above web page. Want to940use Python to parse a webpage? create JSON or XML? use regular941expressions? walk a directory tree? compress a file? use a SQLite942database? Then consult the Python standard library.943944Returning to our function {\tt f} defined above, when we run it, we945might get:946%skip947\begin{lstlisting}948sage: f(7)94949950sage: f(7)951343952\end{lstlisting}953954955Sage (but not Python) also has a notion of Calculus-style functions.956For example,957\begin{lstlisting}958sage: f(x,y) = sin(x) + e^cos(y)959sage: f(2,pi)960e^(-1) + sin(2)961sage: f.integrate(x)962(x, y) |--> x*e^cos(y) - cos(x)963sage: plot3d(f, (x,-pi,pi), (y,-2*pi,2*pi), viewer='tachyon')964\end{lstlisting}965\sagegraphic{.6}{3dplot1}966967968\subsection{Further remarks on passing variables to functions}969We mentioned above that Python uses call by reference semantics.970The following example helps clarify this point very explicitly. First971we create a list and note where it is stored in memory (at address97269421536 on my computer right now).973\begin{lstlisting}974sage: v = [1,2,3]975sage: id(v) # random - memory location of v97669421536977\end{lstlisting}%link978Next we define a function that979prints where in memory its input {\tt w} is stored, and modifies {\tt w}:980%link981\begin{lstlisting}982sage: def foo(w):983... print "location of w =", id(w)984... w.append('hello')985... print "w =", w986\end{lstlisting}%link987When we call foo with {\tt v}, note that the variable {\tt w}988points to the same memory location as {\tt v}:989%link990\begin{lstlisting}991sage: foo(v)992location of w = 69421536993w = [1, 2, 3, 'hello']994\end{lstlisting}%link995Moreover, and it's critical you understand this, the list {\tt v}996has now changed!997%link998\begin{lstlisting}999sage: v1000[1, 2, 3, 'hello']1001\end{lstlisting}%link1002If we want foo to modify a copy of {\tt v} instead, we have to1003explicitly use the {\tt copy} function:1004%link1005\begin{lstlisting}1006sage: foo(copy(v))1007location of w = 695359361008w = [1, 2, 3, 'hello', 'hello']1009\end{lstlisting}%link1010And this worked fine, as expected:1011%link1012\begin{lstlisting}1013sage: v1014[1, 2, 3, 'hello']1015\end{lstlisting}10161017This illustrates part of the ``Zen of Python'':1018\begin{center}1019\LARGE\bf Explicit is better than implicit.1020\end{center}1021To see the rest of the Zen of Python, type {\tt import this} into Sage.10221023\subsection{Gotcha: Default Arguments}1024Consider the following function \verb|my_append(a,L)| which appends {\tt a} to {\tt L}, but whose second argument is optional, so \verb|my_append(a)|1025just creates a new list {\tt [a]}:1026\begin{lstlisting}1027sage: def my_append(a, L=[]):1028... L.append(a)1029... print L1030sage: my_append(1)1031[1]1032sage: my_append(2) # what?1033[1, 2]1034sage: my_append(3) # what? what?1035[1, 2, 3]1036\end{lstlisting}10371038What happened? You might have expected to see output1039of {\tt [1]}, then {\tt [2]}, then {\tt [3]}.1040Let's modify the function \verb|my_append| to also print the memory1041location of {\tt L}.1042\begin{lstlisting}1043sage: def my_append(a, L=[]):1044... L.append(a)1045... print L, id(L)1046sage: my_append(1) # random memory location1047[1] 694384241048sage: my_append(2) # same random memory location1049[1, 2] 694384241050sage: my_append(3) # same random memory location1051[1, 2, 3] 694384241052\end{lstlisting}10531054When the function \verb|my_append| is first encountered by the Python1055interpreter, it evaluates each of the default arguments. When Python1056sees {\tt L=[]}, it creates a list in memory at location 69438424.1057Each time you call \verb|my_append| and don't specify the second argument, that1058same list---at address 69438424---is used, {\em and modified} in1059this case.10601061\subsection{Gotcha: Recursion}10621063Python supports recursive functions, but they are cripled in that1064there is by default a fairly small limit on the depth of recursion1065(you can increase this). This is because Python does not have ``tail1066recursion'' like a language such as lisp.10671068%skip1069\begin{lstlisting}1070def my_factorial(n):1071if n == 1: return n1072assert n > 01073return n * my_factorial(n-1)1074\end{lstlisting}10751076\noindent{}This works fine:1077%skip1078\begin{lstlisting}1079sage: my_factorial(20)108024329020081766400001081\end{lstlisting}1082But:1083%skip1084\begin{lstlisting}1085sage: my_factorial(1000)1086Traceback (click to the left of this block for traceback)1087...1088RuntimeError: maximum recursion depth exceeded in cmp1089\end{lstlisting}10901091So be careful when writing recursive functions. Often recursive1092functions will never ever be called with a big depth. However, if you1093need to write a recursive function that will be called with a big1094depth, you can simply increase the {\tt recursionlimit} as illustrated1095below.1096%skip1097\begin{lstlisting}1098sage: import sys1099sage: sys.getrecursionlimit()110010001101sage: sys.setrecursionlimit(1000000)1102sage: a = my_factorial(1000) # works fine!1103\end{lstlisting}11041105As an aside, you can in fact write little lisp programs using Sage if1106you want, since Sage includes an embedded lisp interpreter.1107For example,1108\begin{lstlisting}1109sage: lisp.eval('(defun factorial (n) (if (= n 1) 1 (* n (factorial (- n 1)))))')1110'FACTORIAL'1111sage: lisp('(factorial 10)')111236288001113sage: lisp(10).factorial()111436288001115\end{lstlisting}111611171118\subsection{Style}1119There is a standard coding style that almost everybody uses when1120writing Python code. Read about it in the Python tutorial:1121\begin{center}\large1122\url{http://docs.python.org/tutorial/controlflow.html#intermezzo-coding-style}1123\end{center}11241125Here is a stylish example:1126\begin{lstlisting}1127def good_function(a, b = 10):1128"""1129This is a good function.11301131This function has a docstring and is named using1132lower_case_with_underscores.11331134It takes as input integers a and b and outputs something computed1135using them. (Notice that the above line is <= 79 characters.)1136"""1137c = 011381139for i in range(a):1140# add i-th power of b to a and1141# put spaces around operators (comment on line of its own).1142c = b**i + a11431144# Another block, and a correctly named class (CamelCase).1145class UselessWrapper(int):1146pass11471148return UselessWrapper(c)1149\end{lstlisting}1150115111521153\section{Classes}11541155Python classes are typically used to define your own new data type1156(though they can be used in other ways as well). New classes are easy1157to define, and support standard object-oriented features such as1158``multiple inheritance'' and ``operator overloading''.11591160Here is a trivial example of a class:1161%skip1162\begin{lstlisting}1163class CoolThing(object):1164def foo(self, xyz):1165print self, xyz1166\end{lstlisting}1167Let's try it out:1168\begin{lstlisting}1169sage: z = CoolThing()1170sage: z.foo('abc')1171<__main__.CoolThing object at 0x...> abc1172sage: type(z)1173<class '__main__.CoolThing'>1174\end{lstlisting}11751176The line {\tt class CoolThing(object):} starts declaration of1177the class {\tt CoolThing}, which derives from the builtin1178class {\tt object}. Typing {\tt z = CoolThing()} creates1179a new instance of the class with the variable {\tt z}1180referencing it. The {\tt foo} method defined above is a function1181that can only be used with instances, which we call1182by writing {\tt z.foo('abc')}. Note that the first argument1183to {\tt def foo(self, xyz)} is {\tt self}, which refers to the1184particular instance of the class.11851186Next, we make a more complicated class, which also illustrates1187how to customize creation of new objects using the1188\verb|__init__| ``dunder method'', define how our1189objects print themselves using \verb|__repr__|,1190and define how {\tt +} and {\tt *} implement1191arithmetic using \verb|__add__| and \verb|__mul__|.11921193% skip1194\begin{lstlisting}1195class MyRational:1196def __init__(self, n, d):1197self._n = Integer(n); self._d = Integer(d)1198def __repr__(self):1199return '%s/%s'%(self._n, self._d)1200def __add__(self, right):1201return MyRational(self._n*right._d + self._d*right._n,1202self._d*right._d)1203def __mul__(self, right):1204return MyRational(self._n*right._n, self._d*right._d)1205def reduced_form(self):1206"""Return the reduced form of this rational number."""1207a = self._n / self._d1208return MyRational(a.numerator(), a.denominator())1209\end{lstlisting}12101211Once we define the above class, we have our own new version of1212``rational numbers''.1213\begin{lstlisting}1214sage: a = MyRational(2,6); b = MyRational(2, 3)1215sage: print a, b12162/6 2/31217sage: a.reduced_form()12181/31219sage: c = a + b; c122018/181221sage: c.reduced_form()12221/11223\end{lstlisting}1224However, notice that subtraction doesn't work:1225\begin{lstlisting}1226sage: a - b1227Traceback (most recent call last):1228...1229TypeError: unsupported operand type(s) for -: 'instance' and 'instance'1230\end{lstlisting}12311232This is because we didn't define a \verb|__sub__| method. You can add1233that method, which looks just like the \verb|__add__| method, except1234with the {\tt +} replaced by a {\tt -}, and subtraction will work.1235Alternatively, we can define a derived class that also defines a1236\verb|__sub__| method, as follows:1237\begin{lstlisting}1238class MyRational2(MyRational): # inheritence (multiple also fully supported)1239def __sub__(self, right):1240return MyRational2(self._n*right._d - self._d*right._n,1241self._d*right._d)1242\end{lstlisting}1243This has absolutely no impact on the original {\tt MyRational} class:1244\begin{lstlisting}1245sage: MyRational(2,6) - MyRational(2, 3)1246Traceback (most recent call last):1247...1248TypeError: unsupported operand type(s) for -: 'MyRational' and 'MyRational'1249\end{lstlisting}1250However, instances of {\tt MyRational2} support subtraction, in addition1251to the multiplication and addition defined above:1252\begin{lstlisting}1253sage: a = MyRational2(2,6); b = MyRational2(2, 3)1254sage: print a, b12552/6 2/31256sage: a + b125718/181258sage: a - b1259-6/181260\end{lstlisting}12611262Big caveat (!): If you do {\tt a+b}, then the resulting object1263is an instance of {\tt MyRational}, not of {\tt MyRational2}!12641265\begin{lstlisting}1266sage: type(a-b)1267<class '__main__.MyRational2'>1268sage: type(a+b)1269<class '__main__.MyRational'>1270\end{lstlisting}12711272This is because the \verb|__add__| method is execute, which explicitly1273refers to MyRational. You can make the code more robust regarding1274derivation by using \verb|self.__class__|, as illustrated below:1275\begin{lstlisting}1276class MyRational3(object):1277def __init__(self, n, d): # called to initialize object1278self._n = Integer(n); self._d = Integer(d)1279def __add__(self, right): # called to implement self + right1280return self.__class__(self._n*right._d + self._d*right._n,1281self._d*right._d)12821283class MyRational4(MyRational3):1284def __sub__(self, right): # called to implement self + right1285return self.__class__(self._n*right._d - self._d*right._n,1286self._d*right._d)1287\end{lstlisting}12881289Now things work better:1290\begin{lstlisting}1291sage: a = MyRational4(2,6); b = MyRational4(2, 3)1292sage: type(a-b), type(a+b)1293(<class '__main__.MyRational4'>, <class '__main__.MyRational4'>)1294\end{lstlisting}129512961297Here is another example that illustrates a default class attribute.1298%skip1299\begin{lstlisting}1300class MyClass:1301"""1302A simple example class.1303"""1304# a Python object attribute; this is basically a default1305# piece of data that is available to each instance of the1306# class, but can be changed in the instance without changing1307# it in the class. (See example below.)1308i = 1234513091310# A function attribute. Again, this is available to each1311# instance, and can be changed in the instance without1312# changing the class object itself.1313def f(self):1314return 'hello world'1315\end{lstlisting}1316First notice that {\tt MyClass} itself is just another Python object1317(we can have variables reference it, pass it into functions, etc.):1318%skip1319\begin{lstlisting}1320sage: MyClass1321<class __main__.MyClass at 0x...>1322sage: MyClass.i1323123451324sage: MyClass.f1325<unbound method MyClass.f>1326sage: MyClass.__doc__1327'A simple example class.'1328\end{lstlisting}1329We ``call'' {\tt MyClass} to create an instance \verb|x| of it:1330%skip1331\begin{lstlisting}1332sage: x = MyClass(); x1333<__main__.MyClass instance at 0x...>1334\end{lstlisting}1335We can then call methods of the instance \verb|x| and get access to its attributes.1336%skip1337\begin{lstlisting}1338sage: x.f()1339'hello world'1340sage: x.i1341123451342\end{lstlisting}1343We can also change the attributes and methods of {\tt x}.1344%skip1345\begin{lstlisting}1346sage: x.i = 501347sage: def g(): return "goodbye"1348sage: x.f = g1349sage: x.f()1350'goodbye'1351\end{lstlisting}13521353This does not change the attributes or methods of {\tt MyClass} or1354new instances of {\tt MyClass}.1355%skip1356\begin{lstlisting}1357sage: y = MyClass(); y.i1358123451359sage: y.f()1360'hello world'1361\end{lstlisting}1362We could change those if we wanted to though, as follows:1363%skip1364\begin{lstlisting}1365sage: def g(self): return "goodbye"1366sage: MyClass.f = g1367sage: y = MyClass()1368sage: y.f()1369'goodbye'1370\end{lstlisting}137113721373As you can see, Python is a {\em dynamic} language. The above is all1374happening at runtime. This is different than static languages such as1375C/C++/Java. It has pros and cons, with the main con being that Python1376can be slower. We will learn about Cython soon, which is similar to1377Python but gives you the option of surrending some of the dynamic1378features of Python in exchange for faster (but less dynamic) static1379semantics.138013811382\subsection{Creating a Number}1383The next example illustrates how to use self and some ``dunder''1384(=double underscore) methods:13851386%skip1387\begin{lstlisting}1388class Number:1389def __init__(self, x):1390# called when Number is instantiated1391self.x = x1392def __repr__(self):1393# defines how Number prints1394return "The Number %s"%self.x1395def __add__(self, right):1396# defines how "+" works1397return Number(self.x + right.x)1398\end{lstlisting}139914001401Now we create a number {\tt n}, print it, and add it (using +) to1402another number.1403%skip1404\begin{lstlisting}1405sage: n = Number(37)1406sage: n1407The Number 371408sage: n + Number(15)1409The Number 521410\end{lstlisting}14111412Try to add subtraction and multiplication to the class {\tt Number}1413right now. The names of the relevant dunder methods are \verb|__sub__|1414and \verb|__mul__|.1415141614171418See \url{http://docs.python.org/reference/datamodel.html} for long1419lists of dunder methods.1420142114221423\section{Data Types: list, tuple, set, str, and file}\label{sec:datatypes}14241425% http://wiki.wstein.org/10/480b/lectures/lec71426% http://wiki.wstein.org/10/480b/lectures/lec814271428\subsection{Lists}\label{sec:lists}1429A list in Python is a finite ordered ``list'' of any Python objects at1430all. Many useful operations are supported, along with a handy ``list1431comprehension'' notation that makes building lists easy.14321433First we create a list, whose entries are an integer, a string, a data1434type, and another list with a list in it. Note that {\tt v} has type1435{\tt list}.1436\begin{lstlisting}1437sage: v = [3, 'hello', Integer, ['a', [1,2]]]1438sage: type(v)1439<type 'list'>1440sage: v1441[3, 'hello', <type 'sage.rings.integer.Integer'>, ['a', [1, 2]]]1442\end{lstlisting}%link1443Lists in Python are 0-based, in that {\tt v[0]} is1444the first entry in the list. {\em Remember this!}1445%link1446\begin{lstlisting}1447sage: v[0]144831449sage: v[1]1450'hello'1451\end{lstlisting}%link1452You can also index into the list from the other side by using negative numbers:1453%link1454\begin{lstlisting}1455sage: v[-1]1456['a', [1, 2]]1457sage: v[-2]1458<type 'sage.rings.integer.Integer'>1459\end{lstlisting}%link1460You can slice lists. When slicing you specify a start and stop point,1461and take all the elements between. Keep in mind that it includes the1462starting point you specify, but excludes the endpoint.1463%link1464\begin{lstlisting}1465sage: v[1:]1466['hello', <type 'sage.rings.integer.Integer'>, ['a', [1, 2]]]1467sage: v[0:3]1468[3, 'hello', <type 'sage.rings.integer.Integer'>]1469sage: v[0:3:2] # just the even-indexed positions1470[3, <type 'sage.rings.integer.Integer'>]1471\end{lstlisting}1472Use {\tt len} to get the length of a list. New Sage/Python users often get very frustrated trying to figure out how to find the length of a list. Just memorize this right now!1473%link1474\begin{lstlisting}1475sage: len(v)147641477\end{lstlisting}%link1478You can also sort, append to, delete elements from, extend, etc., lists.1479See the Python documentation.1480%link1481\begin{lstlisting}1482sage: w = copy(v)1483sage: w.sort(); w1484[3, ['a', [1, 2]], 'hello', <type 'sage.rings.integer.Integer'>]1485sage: w.extend([1,2,3,4]); w1486[3, ['a', [1, 2]], 'hello', <type 'sage.rings.integer.Integer'>,14871, 2, 3, 4]1488\end{lstlisting}14891490You can build lists in place using list comprehension, which is a lot like "set building notation" in mathematics. For example:1491%link1492\begin{lstlisting}1493sage: [n*(n+1)/2 for n in range(1, 10) if n%2 == 1]1494[1, 6, 15, 28, 45]1495\end{lstlisting}149614971498The basic structure of a list comprehension is the following (there are1499more complicated forms):1500%skip1501\begin{lstlisting}1502[ <expression(i)> for i in <iterable> <optional if condition> ]1503\end{lstlisting}15041505Notice above that {\tt for n in range(1,10)} and {\tt if n\%2 == 1}1506are both valid snippets of Python code. Aside from possible scoping1507issues, list comprehensions are basically equivalent to combining a1508{\tt for} loop with an {\tt if} statement in them, where you append to1509a list. To illustrate this, note that you can literally almost1510rearrange the code of such a {\tt for} loop into a list comprehension, for1511example:1512%skip1513\begin{lstlisting}1514z = []1515for n in range(1, 10):1516if n % 2 == 1:1517z.append(n*(n+1)/2)1518\end{lstlisting}15191520If you evaluate the above code, then print {\tt z}, you'll see1521%skip1522\begin{lstlisting}1523sage: z1524[1, 6, 15, 28, 45]1525\end{lstlisting}15261527If you want to be effective with Sage/Python, you must master lists.15281529\subsection{Tuples}15301531Tuples are similar to lists, except you can't change which objects are1532stored in a tuple. Also, there is no tuple-comprehension; you have to1533make a list {\tt v}, then change it into a tuple by typing {\tt1534tuple(v)}. You can however, change the objects themselves if they1535are mutable.15361537\begin{lstlisting}1538sage: v = (3, 'hello', Integer, ['a', [1,2]]); type(v)1539<type 'tuple'>1540sage: v[0] = 5 # nope!1541Traceback (most recent call last):1542...1543TypeError: 'tuple' object does not support item assignment1544sage: v[3].append('change a mutable entry'); v1545(3, 'hello', <type 'sage.rings.integer.Integer'>,1546['a', [1, 2], 'change a mutable entry'])1547\end{lstlisting}15481549{\bf BIG FAT WARNING:} The following looks like a ``tuple1550comprehension'' (if there were such a thing), but it isn't one:1551\begin{lstlisting}1552sage: w = (n*(n+1)/2 for n in range(1, 10) if n%2 == 1); type(w)1553<type 'generator'>1554\end{lstlisting}%link15551556Notice that you can't index into {\tt w}:1557%link1558\begin{lstlisting}1559sage: w[0]1560Traceback (click to the left of this block for traceback)1561...1562TypeError: 'generator' object is unsubscriptable1563\end{lstlisting}%link1564You can iterate over {\tt w} though:1565%link1566\begin{lstlisting}1567sage: for n in w: print n,15681 6 15 28 451569\end{lstlisting}%link1570Here, we get no output since {\tt w} is ``used up''.1571%link1572\begin{lstlisting}1573sage: for n in w: print n,1574\end{lstlisting}157515761577Anyway, if you want to make a tuple using a list comprehension, be1578explicit, like so:1579\begin{lstlisting}1580sage: tuple( n*(n+1)/2 for n in range(1, 10) if n%2 == 1 )1581(1, 6, 15, 28, 45)1582\end{lstlisting}1583158415851586\subsection{Strings}15871588A string is a finite immutable (unchangeable) sequence of characters.1589Python supports a wonderful range of string processing functions. To1590make a string literal:1591\begin{itemize}1592\item Enclose it is in either single or double quotes (just be1593consistent) -- if you use single quotes you can use double quotes in1594your string without escaping them, and vice versa.1595\item For a multiline string use three single or double quotes in a1596row -- then you can include newlines directly in your string.1597\item There are many escape characters for including special1598characters in strings, e.g., \verb|'\n'| for ``newline''. If you put1599the letter {\tt r} right before the quotes you get a raw string, for1600which a backslash just stays a backslash and you can't escape anything;1601this is often useful for \LaTeX{} code.1602\end{itemize}16031604The following examples illustrates some of the above ways of creating strings.16051606\begin{lstlisting}1607sage: s = "this is a string's string using double quotes"; s1608"this is a string's string using double quotes"1609sage: print s1610this is a string's string using double quotes1611sage: s = 'this is a string"s using single quotes'; s1612'this is a string"s using single quotes'1613\end{lstlisting}16141615%skip1616\begin{lstlisting}1617s = """this is a1618multiline string."""16191620s = r"""Consider \sin(x) +1621\cos(y) and add \pi."""1622\end{lstlisting}16231624Strings in Python are extremely flexible and easy to manipulate. You1625can slice them exactly like lists, find substrings, concatenate, etc.16261627\begin{lstlisting}1628sage: s = "This is a string."; s[:10]1629'This is a '1630sage: s[10:]1631'string.'1632sage: s[::2] # get just the even indexed characters1633'Ti sasrn.'1634sage: s.find('a')163581636sage: s + " Yes, a string."1637'This is a string. Yes, a string.'1638sage: s.replace('a', 'b')1639'This is b string.'1640\end{lstlisting}16411642The join method is also amazingly useful. If {\tt s} is a string,1643then {\tt s.join([list of strings])} joins together the list of1644strings putting {\tt s} between each.16451646\begin{lstlisting}1647sage: ', '.join(['Stein', 'William', 'Arthur'])1648'Stein, William, Arthur'1649\end{lstlisting}16501651Other useful methods are upper and capitalize:16521653\begin{lstlisting}1654sage: s = 'this is lower case'; s.upper()1655'THIS IS LOWER CASE'1656sage: s.capitalize()1657'This is lower case'1658\end{lstlisting}16591660Finally, the string formating operator \% appears constantly in Python1661code and is extremely useful to know about. Basically, you just put1662\%s's in your string, and these get replaced by the string1663representations of a tuple of Python objects. Here's how you use it:1664\begin{lstlisting}1665sage: 'Hi %s. Meet %s.'%('Mom', 2/3)1666'Hi Mom. Meet 2/3.'1667\end{lstlisting}16681669Really what just happened was we created a string and a tuple, and1670used the mod operator on them, as illustrated below.1671\begin{lstlisting}1672sage: s = 'Hi %s. Meet %s.'1673sage: t = ('Mom', 2/3)1674sage: s % t1675'Hi Mom. Meet 2/3.'1676\end{lstlisting}16771678There are many other formating options besides just \%s. E.g., \%f is useful for numerical computations.16791680\begin{lstlisting}1681sage: '%.2f %.3f'%(.5, 7/11)1682'0.50 0.636'1683\end{lstlisting}16841685Above, \%.2f formats the string with 2 decimal digits after the point,1686and \%.3f with 3 decimal digits.16871688\subsection{Sets}16891690A set consists of unique elements with no ordering. You know what is1691or isn't in the set and can interate over it. The elements of a set1692must be immutable, since otherwise there would be no way to guarantee1693objects stay unique after putting them together in a set. Lists are1694{\em not immutable} so can't be put in a set, but strings can1695be.1696\begin{lstlisting}1697sage: v = ['this', 'is', 'what', 'this', 'is']; v1698['this', 'is', 'what', 'this', 'is']1699sage: X = set(v); X1700set(['this', 'is', 'what'])1701sage: type(X)1702<type 'set'>1703sage: X[0] # makes no sense1704Traceback (most recent call last):1705...1706TypeError: 'set' object does not support indexing1707sage: 'this' in X1708True1709sage: 'that' in X1710False1711sage: for a in X: print a,1712this is what1713\end{lstlisting}1714Here is how to use the set data structure to obtain1715the distinct types appearing in a list:1716\begin{lstlisting}1717sage: v = [1/2, 5/8, 2.5, 5/2, 3.8]1718sage: t = [type(a) for a in v]; t1719[<type 'sage.rings.rational.Rational'>,1720<type 'sage.rings.rational.Rational'>,1721<type 'sage.rings.real_mpfr.RealLiteral'>,1722<type 'sage.rings.rational.Rational'>,1723<type 'sage.rings.real_mpfr.RealLiteral'>]1724sage: list(set(t))1725[<type 'sage.rings.real_mpfr.RealLiteral'>,1726<type 'sage.rings.rational.Rational'>]1727\end{lstlisting}17281729If you create your own class, you can decide whether or not Python1730should consider it immutable by whether or not you define a1731\verb|__hash__| dunder method. If defined, then your object is1732considered immutable, and is allowed in sets. First, notice that1733sets can't contain lists.1734\begin{lstlisting}1735sage: v = [[1,2], [1,4]]1736sage: set(v)1737Traceback (most recent call last):1738...1739TypeError: unhashable type: 'list'1740\end{lstlisting}1741However, nothing stops us from making a class that derives from list and has a hash method:1742\begin{lstlisting}1743class NaughtyList(list):1744def __hash__(self): # a 32 or 64-bit int; equal objects should have the same hash.1745return hash(str(self))1746\end{lstlisting}1747\begin{lstlisting}1748sage: v = [NaughtyList([1,2]), NaughtyList([1,4])]; v1749[[1, 2], [1, 4]]1750sage: X = set(v); X1751set([[1, 2], [1, 4]])1752\end{lstlisting}1753Do something naughty:1754\begin{lstlisting}1755sage: v[1][1] = 21756sage: X1757set([[1, 2], [1, 2]])1758sage: v[0] == v[1]1759True1760\end{lstlisting}1761The set doesn't know:1762\begin{lstlisting}1763sage: X1764set([[1, 2], [1, 2]])1765\end{lstlisting}1766176717681769\subsection{Files}17701771It is straightforward to open, read, write, append to, and close files1772on disk. For example, below we create a file {\tt foo}, write to it,1773cose it, open it, then read it.17741775\begin{lstlisting}1776sage: F = open('foo','w')1777sage: F1778<open file 'foo', mode 'w' at 0x...>1779sage: F.write('hello there')1780sage: F.close()1781sage: print open('foo').read()1782hello there1783\end{lstlisting}17841785In the Sage notebook each input cell is executed in a different1786directory. Thus if you just create a file in one cell, you can't1787easily open and read it in another cell. The best workaround is to use1788the {\tt DATA} variable, which is a string that contains the name of a1789single directory that all cells have access to, and which you can1790upload/download files to and from using the Data menu.17911792\begin{lstlisting}1793sage notebook: open(DATA + 'foo','w').write('hi')1794sage notebook: print open(DATA + 'foo').read()1795hi1796sage notebook: os.system('ls -l %s'%DATA)1797total 41798-rw-r--r-- 1 sagenbflask sagenbflask 2 ... ... foo179901800sage notebook: print DATA1801/sagenb/flask/sage_notebook.sagenb/home/.../.../data/1802\end{lstlisting}18031804Another important topic involving files is how to read in interesting1805files, e.g., png image files, wav audio files, csv files, Excel1806spreadsheets, etc. There are various ways of loading a huge range of1807interesting files into Sage, but unfortunately there is still no1808single simple command that parses them all. This would be a good idea1809for a student project.181018111812\section{Exception Handling}18131814Python fully supports exception handling, which allows us to raise and1815handle error conditions eloquently. The syntax in Python for1816exception handling is as simple and straightforward as you can1817possibly imagine.18181819We would like to write a function \verb|formal_sum| that1820takes as input two arbitrary objects, and adds them (using +) if {\em1821possible}, and if not creates their sum in some formal sense.1822Our first attempt, of course, does not just magically just work:1823\begin{lstlisting}1824def formal_sum(a, b):1825return a + b1826\end{lstlisting}1827Then:1828\begin{lstlisting}1829sage: formal_sum(2, 3) # good183051831sage: formal_sum(5, [1,2,3]) # nope1832Traceback (most recent call last):1833...1834TypeError: unsupported operand parent(s) for '+':1835'Integer Ring' and '<type 'list'>'1836\end{lstlisting}1837How can we {\em know} whether or not \verb|a + b| will work?1838You could try to do something really complicated by attempting1839to predict (via looking at code?) whether \verb|__add__| will work,1840but that way insanity lies. Instead, use {\em exception handling}:18411842\begin{lstlisting}1843class FormalSum(object):1844def __init__(self, a, b):1845self.a = a; self.b = b1846def __repr__(self):1847return '%s + %s'%(self.a, self.b)18481849def formal_sum(a, b):1850try:1851return a + b1852except TypeError:1853return FormalSum(a,b)1854\end{lstlisting}1855The class {\tt FormalSum} block above defines a new Python class whose instances1856represent the formal sum of the two attributes {\tt a} and {\tt b}, and which1857print themselves nicely. The function \verb|formal_sum| tries to add1858{\tt a} and {\tt b}, and if this causes a {\tt TypeError}, it instead creates1859a {\tt FormalSum} instance.1860\begin{lstlisting}1861sage: formal_sum(3, 8)1862111863sage: formal_sum(5, [1,2,3])18645 + [1, 2, 3]1865sage: formal_sum(5, 'five')18665 + five1867\end{lstlisting}1868186918701871For our next example, instead of catching an error, we create a1872function {\tt divide} that raises an error if the second input {\tt d}1873equals {\tt 0}.18741875%skip1876\begin{lstlisting}1877def divide(n, d):1878if d == 0:1879raise ZeroDivisionError, "divide by 0!?!"1880return n/d1881\end{lstlisting}18821883Try it out:18841885%skip1886\begin{lstlisting}1887sage: divide(5, 7)18885/71889sage: divide(5, 0)1890Traceback (most recent call last):1891...1892ZeroDivisionError: divide by 0!?!1893\end{lstlisting}18941895Typically, if you try to divide numbers at the denominator is 0, Sage1896will raise a {\tt ZeroDivisionError}. Just as above, we can catch1897this case if we want, and return something else, as illustrated below:18981899%skip1900\begin{lstlisting}1901def divide2(n, d):1902try:1903return divide(n, d) # or just put "n/d"1904except ZeroDivisionError:1905return 'infinity'1906\end{lstlisting}19071908%skip1909\begin{lstlisting}1910sage: divide2(5, 3)19115/31912sage: divide2(5, 0)1913'infinity'1914\end{lstlisting}19151916This web page \url{http://docs.python.org/lib/module-exceptions.html}1917lists all the standard builtin exceptions along with what each means.1918Some common exceptions that often appear in the context of mathematics1919are: {\tt TypeError, ZeroDivisionError, ArithmeticError, ValueError,1920RuntimeError, NotImplementedError, OverflowError, IndexError}.1921We illustrate each of these below:19221923\begin{lstlisting}1924sage: ''.join([1,2])1925Traceback (most recent call last):1926...1927TypeError: sequence item 0: expected string, sage.rings.integer.Integer found1928sage: 1/01929Traceback (most recent call last):1930...1931ZeroDivisionError: Rational division by zero1932sage: factor(0)1933Traceback (most recent call last):1934...1935ArithmeticError: Prime factorization of 0 not defined.1936sage: CRT(2, 1, 3, 3)1937Traceback (most recent call last):1938...1939ValueError: No solution to crt problem since gcd(3,3) does not divide 2-11940sage: find_root(SR(1), 0, 5)1941Traceback (most recent call last):1942...1943RuntimeError: no zero in the interval, since constant expression is not 0.1944sage: RealField(50)(brun)1945Traceback (most recent call last):1946...1947NotImplementedError: brun is only available up to 41 bits1948sage: float(5)^float(902830982304982)1949Traceback (most recent call last):1950...1951OverflowError: (34, 'Numerical result out of range')1952sage: v = [1,2,3]1953sage: v[10]1954Traceback (most recent call last):1955...1956IndexError: list index out of range1957\end{lstlisting}19581959The key points to remember about exceptions are:1960\begin{enumerate}1961\item Three keywords: {\tt try}, {\tt except}, {\tt raise}1962\item How to catch multiple possible exceptions correctly (there is a gotcha here -- see below!).1963\item One more keyword: {\tt finally}1964\end{enumerate}19651966There is more to exceptions, but these are the key points. We1967illustrate the last two below in a contrived example.196819691970%skip1971\begin{lstlisting}1972def divide(n, d):1973try:1974return n/d1975except (ZeroDivisionError, ValueError), msg:1976print msg1977return '%s/%s'%(n,d)1978except TypeError, NotImplementedError:1979# the above line is PURE EVIL(!)1980print "NotImplementedError is now '%s'"%NotImplementedError1981print "What have I done?!"1982finally:1983print "The finally block is *always* executed."1984\end{lstlisting}19851986Now try it out:1987%skip1988\begin{lstlisting}1989sage: divide(2,3)1990The finally block is *always* executed.19912/31992sage: divide(2, 0)1993Rational division by zero1994The finally block is *always* executed.1995'2/0'1996sage: divide('hi', 'mom')1997NotImplementedError is now 'unsupported operand type(s) for /: 'str' and 'str''1998What have I done?!1999The finally block is *always* executed.2000\end{lstlisting}20012002The form of the except statement is:2003%skip2004\begin{lstlisting}2005except [single exception], message2006\end{lstlisting}2007%skip2008\begin{lstlisting}2009except (tuple,of,exceptions), message</p>2010\end{lstlisting}20112012For example,2013\begin{lstlisting}2014try:2015import foobar20161/02017except (ZeroDivisionError, ImportError), msg:2018print "oops --", msg2019\end{lstlisting}2020outputs {\tt oops -- No module named foobar}2021and2022\begin{lstlisting}2023try:20241/02025import foobar2026except (ZeroDivisionError, ImportError), msg:2027print "oops --", msg2028\end{lstlisting}2029outputs {\tt oops -- Rational division by zero}.20302031An extremely confusing error, which has cost me hours of frustration,2032is to write2033%skip2034\begin{lstlisting}2035except exception1, exception2:2036\end{lstlisting}20372038The result is that if exception1 occurs, then exception2 is set to the2039error message. Don't make the same mistake.20402041For example, if we evaluate2042\begin{lstlisting}2043try:20441/02045import foobar2046except ZeroDivisionError, ImportError:2047print "oops --"2048\end{lstlisting}2049then evaluate2050\begin{lstlisting}2051try:2052import foobar20531/02054except ImportError, ZeroDivisionError:2055print "oops --"2056\end{lstlisting}2057we see2058\begin{lstlisting}2059Traceback (click to the left of this block for traceback)2060...2061ImportError: No module named foobar2062\end{lstlisting}20632064Wait, what just happened above? We appear to have totally broken2065Python!? Actually, we have smashed the {\tt ImportError} variable,2066making it point at the {\tt ZeroDivisionError} message above!20672068\begin{lstlisting}2069sage: ImportError2070ZeroDivisionError('Rational division by zero',)2071\end{lstlisting}2072We can fix this for now using the {\tt reset} command, which resets a variable2073to its default state when Sage started up:2074\begin{lstlisting}2075sage: reset('ImportError')2076sage: ImportError2077<type 'exceptions.ImportError'>2078\end{lstlisting}207920802081Another {\em major mistake} I made once\footnote{Actually, I made it2082several hundred times in 2005--2006!}2083with exceptions is illustrated in the following example code:2084%skip2085\begin{lstlisting}2086def divide(n, d):2087if d == 0:2088raise ZeroDivisionError, "error dividing n(=%s) by d(=%s)"%(n,d)2089\end{lstlisting}20902091It's so friendly and nice having a helpful error message that explains2092what went wrong in the division, right? (No!):2093%skip2094\begin{lstlisting}2095sage: divide(3948,0)2096Traceback (click to the left of this block for traceback)2097...2098ZeroDivisionError: error dividing n(=3948) by d(=0)2099\end{lstlisting}21002101But if we put a large value of $n$ as input, then several seconds (or2102minutes!) will be spent just creating the error message. It's2103ridiculous that divide2 below takes over 3 seconds, given that all the2104time is spent creating an error message that we just ignore.2105%skip2106\begin{lstlisting}2107def divide2(n,d):2108try:2109divide(n, d)2110except ZeroDivisionError, msg:2111return 'infinity'2112\end{lstlisting}21132114%skip2115\begin{lstlisting}2116sage: n = 3^(10^7)2117sage: time divide2(n, 0)2118'infinity'2119Time: CPU 3.45 s, Wall: 3.46 s2120\end{lstlisting}21212122Once the Sage developer David Harvey spent a long time tracking down2123why certain power series arithmetic in Sage was so slow for his2124application. It turned out that deep in the code there was a2125try/except block in which the error message itself took over a minute2126to construct, and then it was immediately discarded. {\bf Moral:}2127{\em be very careful when constructing the error message that you2128include along with an exception!}212921302131\section{Decorators}\label{sec:decorators}21322133The definition of decorators is remarkably simple,2134but using them is subtle, powerful, and potentially dangerous.2135From PEP 318 (see \url{http://www.python.org/dev/peps/pep-0318}), we have2136the following new notation in Python (note the first line with the mysterious2137@ sign):2138%skip2139\begin{lstlisting}2140@dec12141def func(arg1, arg2, ...):2142pass2143\end{lstlisting}2144This is equivalent to:21452146%skip2147\begin{lstlisting}2148def func(arg1, arg2, ...):2149pass2150func = dec2(dec1(func))2151\end{lstlisting}2152That's it!21532154To motivate the point of decorators, let's make a function called {\tt2155echo} that takes as input a function {\tt f}, and returns a new2156function that acts just like {\tt f}, except that it prints all of its2157inputs. Here we use {\tt *args} and {\tt **kwds}, which is something2158that we have not discussed before. In Python, use {\tt *args} to2159refer to all of the positional inputs to a function, and {\tt **kwds}2160to refer to all of the keyword inputs. When you do this, {\tt args}2161is a Python tuple containing the positional inputs in order, and {\tt kwds}2162is a dictionary of the {\tt keyword=value} pairs. You can pass args and2163kwds on to another function (as illustrated below) by typing {\tt *args}2164and {\tt **kwds}.21652166%skip2167\begin{lstlisting}2168def echo(f):2169def g(*args, **kwds):2170print "args =", args2171print "kwds =", kwds2172return f(*args, **kwds)2173return g2174\end{lstlisting}21752176Now, let's try it out. Define a function:2177%skip2178\begin{lstlisting}2179def add_em_up(a, b, c):2180return a + b + c2181\end{lstlisting}2182Now use it:2183%skip2184\begin{lstlisting}2185sage: add_em_up(1, 2, 3)218662187\end{lstlisting}21882189The following works, but it sort of looks funny.2190%skip2191\begin{lstlisting}2192sage: add_em_up = echo(add_em_up)2193sage: add_em_up(1, 2, 3)2194args = (1, 2, 3)2195kwds = {}219662197\end{lstlisting}21982199Using a decorator right when we define \verb|add_em_up| is much, much cleaner:2200%skip2201\begin{lstlisting}2202@echo2203def add_em_up(a, b, c):2204return a + b + c2205\end{lstlisting}2206Now we have:2207%skip2208\begin{lstlisting}2209sage: add_em_up(1, 2, 3)2210args = (1, 2, 3)2211kwds = {}221262213\end{lstlisting}22142215Here's another example of a very handy decorator (only available in the Sage notebook):22162217%skip2218\begin{lstlisting}2219@interact2220def add_em_up(a=1, b=[1..10], c=(1..10)):2221return a + b + c2222\end{lstlisting}2223\begin{center}2224\includegraphics[width=.5\textwidth]{graphics/interact}2225\end{center}22262227A hope you can {\em sense the possibilities...}. Here we do type checking:22282229%skip2230\begin{lstlisting}2231class returns:2232def __init__(self, typ):2233self._typ = typ2234def __call__(self, f):2235return lambda *args, **kwds: self._typ(f(*args, **kwds))223622372238@returns(float)2239def f(n,m):2240"""Returns n + m."""2241return n + m2242\end{lstlisting}2243Let's try it out:2244%skip2245\begin{lstlisting}2246sage: f(2,3)22475.02248sage: type(f(5,6))2249<type 'float'>2250sage: f('4', '123')22514123.02252\end{lstlisting}22532254Here's another example I use all the time. If you put {\tt2255@parallel(ncpus)} before a function {\em and} you call the function2256using a list as input, then the function gets evaluated at each2257element of the list in parallel, and the results are returned as an2258iterator. If you call the function without giving a list as input, it2259just works as usual (not in parallel).226022612262%skip2263\begin{lstlisting}2264@parallel(10)2265def f(n):2266sleep(1) # make function seem slow2267return n*(n+1)/22268\end{lstlisting}2269First, try it not in parallel, which takes a long time.2270%skip2271\begin{lstlisting}2272%time2273sage: for n in [1..10]: print n, f(n)22741 122752 322763 622774 1022785 1522796 2122807 2822818 3622829 45228310 552284CPU time: 0.00 s, Wall time: 10.00 s2285\end{lstlisting}22862287Now try it in parallel:2288%skip2289\begin{lstlisting}2290%time2291sage: for X in f([1..10]): print X2292(((1,), {}), 1)2293(((2,), {}), 3)2294(((3,), {}), 6)2295(((4,), {}), 10)2296(((5,), {}), 15)2297(((6,), {}), 21)2298(((7,), {}), 28)2299(((8,), {}), 36)2300(((9,), {}), 45)2301(((10,), {}), 55)2302CPU time: 0.19 s, Wall time: 1.32 s2303\end{lstlisting}2304230523062307\section{The Ecosystem}23082309The Sage distributuion itself consists of about 100 open source2310programs and libraries, which (like Linux) are developed by a loosely2311knit international group of programmers. Many of these programs are2312written as Python libraries.23132314Any software engineer knows that a programming language is much more2315than just the formal language specification or even a particular2316implementation. It's also the user community, the general pace of2317development, and---most importantly---the collections of tools and2318libraries that are available in that language, especially the free2319ones. Python excels in available tools, as the following list of many2320of the Python-based components of Sage attests:23212322\begin{itemize}2323\item Pycrypto -- fast implementations of many cryptosystems.2324\item Cython -- a Python compiler and tool for efficient use of C/C++2325libraries from Python. We will have much more to say about Cython2326in Chapter~\ref{ch:cython}.2327\item IPython -- interactive interpreter shell2328\item Jinja2 -- HTML and other templating tools; popular for web2329applications.2330\item Moinmoin -- a standalone wiki, e.g., the one used by \url{http://wiki.sagemath.org}.2331\item PIL -- Python imaging library (a ``programmable Photoshop'')2332\item Pygments -- HTML source code highlighting2333\item SQLalchemy -- abstracts interface to most SQL databases and an2334object:relational mapper.2335\item Sphinx -- ReST documentation system for Python, which is used by2336many Python projects (including Sage).2337\item Twisted -- a networking framework; everything from web applications2338to email to ssh servers are implemented in Twisted.2339\item ZODB -- The Zope object-oriented database2340\item arpack -- A sparse numerical linear algebra library.2341\item CVXopt -- A library for solving convex (and other) optimization problems.2342\item Docutils -- related to Python documentation2343\item easy-install -- you can do \verb|easy_install foobar| to2344install any of the over 13,000 Python packages available at2345\url{http://pypi.python.org/}.23462347\item gd -- very quickly draw png images with lines, arcs, etc.23482349\item matplotlib -- the canonical Python 2d graphics library23502351\item mpmath -- arbitrary precision floating point mathematics2352special functions, numerical integration, matrices, etc.23532354\item NumPy -- an $n$-dimensional array library, which is the2355fundamental package needed for scientific computing with Python.23562357\item pexpect -- control command-line subprocesses23582359\item rpy2 -- fast compiled interface to the R statistics program,2360which is also included in Sage.23612362\item sage -- the Sage library; mainly implements mathematical2363algorithms, especially symbolic ones.23642365\item sagenb -- the Sage notebook web application (can be used2366standalone separate from Sage).23672368\item sagetex -- allows you to embed Sage in \LaTeX{} documents23692370\item SciPy -- a large library of numerical functions that are useful2371in mathematics, science, and engineering, including numerical2372integration, optimization, statistics, differential equations, etc.23732374\item setuptools -- package for distributing and working with standalone python packages23752376\item SymPy -- a lightweight Python library for symbolic mathematics.2377\end{itemize}23782379\section{Exercise: Build Python from Source}2380If your computer operating system is Linux or OS X (with XCode2381installed), it is an easy ``exercise'' to build the Python language2382from source. This is particularly relevant if you want to understand2383Python more deeply, since you can change anything you want in the2384interpreter itself, recompile, and try out the result!23852386First, go to \url{http://python.org/download/} and download some2387version of Python. I am using OS X (with XCode installed) and choose2388Python 3.2. In a few seconds I have the file2389{\tt Python-3.2.tar.bz2} in my {\sf Downloads} folder.2390Using the {\sf Terminal} application, I navigate to that folder, extract2391Python, configure and build it, which takes under 2 minutes (!).2392%skip2393\begin{lstlisting}2394deep:~ wstein$ cd Downloads2395deep:Downloads wstein$ tar xf Python-3.2.tar.bz22396deep:Downloads wstein$ cd Python-3.22397deep:Python-3.2 wstein$ ./configure; time make -j82398...2399real 1m18.284s2400user 1m59.552s2401sys 0m9.980s2402deep:Python-3.2 wstein$2403\end{lstlisting}24042405And now let's try it out:2406\begin{lstlisting}2407deep:Python-3.2 wstein$ ./python.exe2408Python 3.2 (r32:88445, Mar 30 2011, 10:20:45)2409[GCC 4.2.1 (Apple Inc. build 5666) (dot 3)] on darwin2410Type "help", "copyright", "credits" or "license" for more2411information.2412>>> 2 + 2241342414\end{lstlisting}24152416For fun, let's change something in the core of Python, recompile, and2417observe our change. On line 288 of {\tt2418Python-3.2/Objects/listobject.c}, I insert a line that calls the C2419{\sf printf} function to print out some graffiti:2420\begin{lstlisting}2421...2422PyObject *2423PyList_GetItem(PyObject *op, Py_ssize_t i)2424{2425printf("Hi Mom!\n"); /* I added this! */2426if (!PyList_Check(op)) {2427PyErr_BadInternalCall();2428...2429\end{lstlisting}24302431I then type {\tt make} again, wait a few seconds, and try out Python again:2432\begin{lstlisting}2433deep:Python-3.2 wstein$ ./python.exe2434Python 3.2 (r32:88445, Mar 30 2011, 10:25:56)2435[GCC 4.2.1 (Apple Inc. build 5666) (dot 3)] on darwin2436Type "help", "copyright", "credits" or "license" for more2437information.2438Hi Mom!2439...2440Hi Mom!2441>>> v = [1,2,3]2442>>> v[0]244312444>>> v['a']2445Traceback (most recent call last):2446File "<stdin>", line 1, in <module>2447Hi Mom!2448Hi Mom!2449Hi Mom!2450Hi Mom!2451Hi Mom!2452TypeError: list indices must be integers, not str2453\end{lstlisting}24542455Interestingly, the function \verb|PyList_GetItem| appears to not be2456called when we use an integer to access a list, but it is used when we2457try to access the list with anything else.24582459For more information about how the Python source code is laid out,2460see the README file, especially the section at the end called2461``Distribution structure''.24622463\chapter{Cython}\label{ch:cython}24642465Cython is a compiled variant of the Python language, that can be used2466to write very, very fast code for Sage, and makes it possible to2467efficiently use code from existing C/C++ libraries in Sage, or write2468part of a program in C and the other part in Sage.24692470Additional references for learning Cython include the Cython Users2471Guide (see \url{http://docs.cython.org/src/userguide/}) as you can,2472and the other Cython documentation (see \url{http://docs.cython.org/}).24732474\section{An Example: Speeding up a Simple Function}24752476Let's start with a first simple example. We write a brute force2477Python program that computes a double precision (53-bit) floating2478point approximation to2479$$2480f(n) = \sin(1) + \sin(2) + \sin(3) + \cdots + \sin(n-1) + \sin(n),2481$$2482for $n$ any positive integer.2483Let's start with the most naive and straightforward2484implementation of this:2485\begin{lstlisting}2486def python_sum_symbolic(n):2487return float( sum(sin(n) for n in xrange(1,n+1)) )2488\end{lstlisting}24892490As a sanity check for each of our implementations, we compute $f(1000)$:2491\begin{lstlisting}2492sage: python_sum_symbolic(1000)24930.8139696340731642494\end{lstlisting}24952496Let's benchmark it (these benchmarks were done under OS X on a 2.66GHz2497Intel Core i7 laptop).2498\begin{lstlisting}2499sage: timeit('python_sum_symbolic(1000)')25005 loops, best of 3: 193 ms per loop2501\end{lstlisting}25022503Next, try a bigger input. We use time instead, since that only runs2504the function once, and this is going to take awhile:25052506\begin{lstlisting}2507sage: time python_sum_symbolic(10^4)25081.63389102179244672509Time: CPU 15.88 s, Wall: 17.64 s2510\end{lstlisting}251125122513This is really, really, shockingly slow. By the end of this section,2514you'll see how to compute $f(n)$ over {\bf 17 million times} faster!251525162517One reason \verb|python_sum_symbolic| is so bad, is that the2518\verb|sin| function is a ``symbolic function'', in that it keeps2519track of exact values, etc. For example,2520\begin{lstlisting}2521sage: sum(sin(n) for n in xrange(1, 10+1))2522sin(1) + sin(2) + sin(3) + sin(4) + sin(5) + sin(6) +2523sin(7) + sin(8) + sin(9) + sin(10)2524\end{lstlisting}25252526So the \verb|python_sum_symbolic| function is computing a huge formal2527sum, then converting each summand to a float, and finally adding them2528up. This is wasteful. Our next implementation using a different sin2529function, one that is standard with Python, which immediately turns2530its input into a float (53-bit double precision number), computes2531\verb|sin|, and returns the result as a float.25322533\begin{lstlisting}2534def python_sum(n):2535from math import sin2536return sum(sin(i) for i in xrange(1, n+1))2537\end{lstlisting}253825392540How does this do?2541\begin{lstlisting}2542sage: python_sum(1000) # right answer25430.81396963407316592544sage: timeit('k = python_sum(1000)')2545625 loops, best of 3: 185 microseconds per loop2546sage: timeit('k = python_sum(10^4)')2547125 loops, best of 3: 2.07 ms per loop2548\end{lstlisting}2549Thus \verb|python_sum| is over {\em one thousand times faster}2550than \verb|python_sum_symbolic|!25512552\begin{lstlisting}2553sage: 193e-3 / 185e-625541043.243243243242555sage: 15.88 / 2.07e-325567671.497584541062557\end{lstlisting}255825592560Perhaps there is some win in not using \verb|sum|?2561\begin{lstlisting}2562def python_sum2(n):2563from math import sin2564s = float(0)2565for i in range(1, n+1):2566s += sin(i)2567return s2568\end{lstlisting}25692570Try it out: nope, no win at all (!):2571\begin{lstlisting}2572sage: timeit('k = python_sum2(10^4)')2573125 loops, best of 3: 2.11 ms per loop2574\end{lstlisting}257525762577Next we try using Cython to make our code faster.2578Note that our2579``rewritten'' program looks identical---the only difference so far is2580that we told Sage to compile the program using Cython by putting {\tt2581\%cython} at the beginning of the block (if you are using the2582command line instead of the notebook, put the code without2583\verb|%cython|2584in a file \verb|foo.pyx|, then type \verb|load foo.pyx|.)25852586%skip2587\begin{lstlisting}2588%cython2589def cython_sum(n):2590from math import sin2591return sum(sin(i) for i in xrange(1, n+1))2592\end{lstlisting}25932594If you evaluate the above code in the Sage notebook, you'll see that2595two linked files appear after the input cell:2596\begin{enumerate}2597\item A file that ends in .c: this is the C program that the above2598code got turned into. This is compiled and linked automatically2599into the running copy of Sage.2600\item A file that ends in .html: this is an annotated version of the2601above Cython program; double click on a line to see2602the corresponding C code.2603\end{enumerate}26042605Is the Cython program any faster?2606%skip2607\begin{lstlisting}2608sage: cython_sum(1000) # it works26090.81396963407316592610sage: timeit('cython_sum(1000)')2611625 loops, best of 3: 144 microseconds per loop2612sage: timeit('k = cython_sum(10^4)')2613625 loops, best of 3: 1.52 ms per loop2614\end{lstlisting}2615It's faster, but only about 30\% faster. This is very typical2616of what you should expect by simply putting \verb|%cython| above Python2617code.2618\begin{lstlisting}2619sage: 2.07/1.5226201.361842105263162621\end{lstlisting}262226232624The Cython program is not {\em that} much faster, because the computer2625is doing essentially the same thing in both functions.2626In the case of \verb|python_sum|,2627the Python interpreter is carrying out a sequence of operations2628(calling functions in the Python C library), and in the case of2629\verb|cython_sum|, a C program is running (the compiled Cython2630module), which is simply calling pretty much the same functions in the2631Python C library.26322633To get a further speedup, we declare certain variables to have C data2634types and use a C version of the sin function.26352636%skip2637\begin{lstlisting}2638%cython2639cdef extern from "math.h":2640double sin(double)26412642def cython_sum_typed_lib(long n):2643cdef long i2644return sum(sin(i) for i in range(1, n+1))2645\end{lstlisting}2646The differences are that we declared {\tt n} to be long and we added a2647new line \verb|cdef long i|, which declares {\tt i} to2648also be long. This tells Cython to treat {\tt n} and {\tt i}2649as being of data type {\tt long}, which is a 32 or 64-bit integer,2650depending on the computer you're using.2651This is the same as the long2652datatype in C.2653We also call the \verb|sin| function in the math C library.2654Let's see if this is faster.26552656\begin{lstlisting}2657sage: cython_sum_typed_lib(1000)26580.81396963407316592659sage: timeit('cython_sum_typed_lib(1000)')2660625 loops, best of 3: 85.2 Âµs per loop2661sage: timeit('cython_sum_typed_lib(10^4)')2662625 loops, best of 3: 951 Âµs per loop2663sage: 2.07e-3 / 778e-626642.660668380462722665\end{lstlisting}26662667So now our coding is beating the pure Python code by a factor of nearly 3.26682669In general, you have to be more careful, e.g., {\tt long} integers2670{\em overflow}. They may be either 32 or 64-bit2671depending on the computer you are using. The following example2672illustrates overflow:2673%skip2674\begin{lstlisting}2675%cython2676def longmul(long a, long b):2677return a*b2678\end{lstlisting}2679Now let's try it:26802681%skip2682\begin{lstlisting}2683sage: longmul(2^10, 2^20)268410737418242685sage: longmul(2^20, 2^50) # overflows!268602687sage: 2^20 * 2^50268811805916207174113034242689\end{lstlisting}269026912692For our next optimization, we use a for loop instead of the2693sum command:2694\begin{lstlisting}2695%cython2696cdef extern from "math.h":2697double sin(double)26982699def cython_sum_typed_lib_loop(long n):2700cdef long i2701cdef double s = 02702for i in range(1, n+1):2703s += sin(i)2704return s2705\end{lstlisting}2706In Cython, this is really worth it.2707\begin{lstlisting}2708sage: cython_sum_typed_lib_loop(1000)27090.81396963407316592710sage: timeit('cython_sum_typed_lib_loop(10^4)')2711625 loops, best of 3: 298 Âµs per loop2712sage: 2.07e-3 / 298e-627136.946308724832212714\end{lstlisting}2715We thus obtain a speedup of a factor of about 7 by switching2716from Python to Cython, and implementing exactly the same algorithm.2717Way under the hood, the same \verb|sin| function (providing by the2718math library on the operating system) is being called, but2719the Cython version of the function avoids a lot of overhead.2720272127222723Another approach to this particular numerical problem is to use the2724numpy library, which allows one to evaluate a function on all entries2725in an array, and sum the entries. This is called ``vectorization''.2726Here's what we get in this particular example:27272728\begin{lstlisting}2729def sum_numpy(n):2730import numpy2731return sum(numpy.sin(numpy.array(range(1, n+1))))2732\end{lstlisting}27332734Let's try it out:2735\begin{lstlisting}2736sage: sum_numpy(1000)27370.813969634073165922738sage: timeit('sum_numpy(10^4)')2739625 loops, best of 3: 1.33 ms per loop2740sage: 2.07e-3 / 1.33e-327411.556390977443612742\end{lstlisting}2743So for this benchmark, our Numpy implementation is slightly better than2744pure Python, but Cython is still much faster.27452746Finally, we try a different algorithm, both using Python and Cython.2747The symbolic capabilities of Sage can be used to find closed2748form expressions for certain formal sums.2749\begin{lstlisting}2750sage: var('i, n')2751sage: f = sum(sin(i), i, 1, n).full_simplify(); f27521/2*((cos(1) - 1)*sin(n*arctan(sin(1)/cos(1))) +2753sin(1)*cos(n*arctan(sin(1)/cos(1))) - sin(1))/(cos(1) - 1)2754\end{lstlisting}2755Thus (and I had no idea before trying this),2756$$2757f(n) = \frac{{\left(\cos\left(1\right) - 1\right)} \sin\left(n \arctan\left(\frac{\sin\left(1\right)}{\cos\left(1\right)}\right)\right) + \sin\left(1\right) \cos\left(n \arctan\left(\frac{\sin\left(1\right)}{\cos\left(1\right)}\right)\right) - \sin\left(1\right)}{2 \, {\left(\cos\left(1\right) - 1\right)}}2758$$27592760Using this, we can give an algorithm to compute this sum that is much2761faster than anything above, and scales much better as well2762to larger $n$.2763\begin{lstlisting}2764def sum_formula(n):2765from math import sin, cos, atan as arctan2766return (1/2*((cos(1) - 1)*sin(n*arctan(sin(1)/cos(1))) +2767sin(1)*cos(n*arctan(sin(1)/cos(1))) - sin(1))/(cos(1) - 1))2768\end{lstlisting}2769How does it do?2770\begin{lstlisting}2771sage: sum_formula(1000)27720.81396963407316642773sage: timeit('sum_formula(10^4)')2774625 loops, best of 3: 36.5 Âµs per loop2775sage: 2.07e-3 / 36.5e-6277656.71232876712332777\end{lstlisting}2778Finally, we implement the closed formula, but instead in Cython.27792780\begin{lstlisting}2781%cython27822783cdef extern from "math.h":2784double sin(double)2785double cos(double)2786double atan(double)27872788def sum_formula_cython(double n):2789return (.5*((cos(1) - 1)*sin(n*atan(sin(1)/cos(1))) +2790sin(1)*cos(n*atan(sin(1)/cos(1))) - sin(1))/(cos(1) - 1))2791\end{lstlisting}27922793This is 40 times faster than our Python implementation!2794\begin{lstlisting}2795sage: sum_formula_cython(1000)27960.81396963407316642797sage: timeit('sum_formula_cython(10^4)')2798625 loops, best of 3: 906 ns per loop2799sage: 36.5e-6 / 906e-9280040.28697571743932801\end{lstlisting}2802And it is over 2000 times faster than \verb|python_sum|.28032804\begin{lstlisting}2805sage: 2.07e-3 / 906e-928062284.768211920532807\end{lstlisting}2808It is a whopping {\bf 17 million times} faster than our first attempt!2809\begin{lstlisting}2810sagre: 15.88 / 906e-928111.75275938189845e72812\end{lstlisting}2813281428152816\section{Using External C/C++ Code}28172818Cython is absolutely critical to the design of Sage, and potentially2819very important to your own work, because it makes it possible to2820efficiently make use of data types and functions defined in any C/C++2821library. Since there is an enormous amount of useful, fast, debugged2822C/C++ code out there, Cython gives your Sage and Python programs2823access to vast amounts of useful capabilities. Also, when used2824correctly, there is no overhead in calling out to the C libraries,2825unlike the situation with SWIG, ctypes, and many other approaches2826to writing C library wrappers.28272828\subsection{Simple random example}2829Here's a first simple example. Type {\tt man random} on the command line2830(or Google it) to find out about the random C library function:2831\begin{lstlisting}2832RANDOM(3) BSD Library Functions Manual RANDOM(3)28332834NAME2835initstate, random, setstate, srandom, srandomdev -- better2836random number generator; routines for changing generators28372838LIBRARY2839Standard C Library (libc, -lc)28402841SYNOPSIS2842#include <stdlib.h>28432844char *2845initstate(unsigned seed, char *state, size_t size);28462847long2848random(void);2849...2850\end{lstlisting}28512852Despite {\tt random} being a function defined in the standard C2853library, we can still call it from Cython, as follows:2854%skip2855\begin{lstlisting}2856%cython2857cdef extern from "stdlib.h": # (1)2858long random() # (2)28592860def random_nums(int n): # (3)2861cdef int i # (4)2862v = [random() for i in range(n)] # (5)2863return v2864\end{lstlisting}28652866Let's try it out:2867%skip2868\begin{lstlisting}2869sage: random_nums(5)2870[1315705257, 1147455227, 1571270137, 1106977565, 1805149207]2871sage: timeit('v = random_nums(10^5)')2872125 loops, best of 3: 5.56 ms per loop2873\end{lstlisting}28742875It's interesting to see how this compares to pure Python. Here's the same program in Python:2876%skip2877\begin{lstlisting}2878%python2879import random2880k = 2**31-12881def py_random_nums(n):2882return [random.randint(0,k) for i in range(n)]2883\end{lstlisting}2884So the speedup is by a factor of nearly 50:2885%skip2886\begin{lstlisting}2887sage: py_random_nums(5)2888[317567506, 1289482476, 1766134327, 1216261810, 1427493671]2889sage: timeit('v = random_nums(10^5)')28905 loops, best of 3: 251 ms per loop2891sage: 251/5.56289245.14388489208632893\end{lstlisting}28942895Finally we explain the above code line by line. (TODO)2896289728982899\subsection{Adding rational numbers using MPIR}29002901We next consider a more mathematical example: arithmetic with2902arbitrary precision rational numbers. The MPIR C library (which is2903included with Sage, but can also be downloaded separately for free for2904any standard operating system from \url{http://mpir.org/}) provides2905highly optimized arithmetic with arbitrary precision integers and2906rational numbers.\footnote{MPIR and GMP \url{http://gmplib.org/} are2907basically the same for our discussion; technically they are2908``forks'' of each other, but export essentially the same functions.}2909We could make use of MPIR by reading the documentation for MPIR and2910using {\tt cdef extern} as above. Fortunately, all of the necessary2911{\tt cdef extern} declarations needed to use MPIR are already declared2912in Sage. You can view all the declarations from the notebook by2913navigating to {\tt <url of notebook server>/src/libs/gmp}.29142915Let's use MPIR directly to create two rational numbers and add them2916together. The code below is complicated and illustrates many issues2917and techniques, so we will explain it in great depth. Once you2918understand this, you can deal with many issues that will come up with2919Cython.29202921%skip2922\begin{lstlisting}2923%cython2924from sage.libs.gmp.all cimport * # (1)2925def add_rationals(bytes a, bytes b): # (2)2926cdef mpq_t x, y, z # (3)2927mpq_init(x); mpq_init(y); mpq_init(z) # (4)2928mpq_set_str(x, a, 10) # base 10 string # (5)2929mpq_set_str(y, b, 10)2930mpq_add(z, x, y) # (6)2931cdef int n = (mpz_sizeinbase (mpq_numref(z), 10) # (7)2932+ mpz_sizeinbase (mpq_denref(z), 10) + 3)2933cdef char* s = <char*>sage_malloc(sizeof(char)*n) # (8)2934if not s: raise MemoryError # (9)2935cdef bytes c = mpq_get_str(s, 10, z) # (10)2936mpq_clear(x); mpq_clear(y); mpq_clear(z) # (11)2937sage_free(s) # (12)2938return c2939\end{lstlisting}29402941Now let's try it out:2942%skip2943\begin{lstlisting}2944sage: add_rationals('2/3', '-5/21')2945'3/7'2946sage: 2/3 - 5/2129473/72948sage: add_rationals('1/29048203984092834823049',2949'-394/29302938402384092834')2950'-11444963066794174536188472/851197732045760533660225724673976778930866'2951\end{lstlisting}29522953Timings suggest we didn't mess up:2954%skip2955\begin{lstlisting}2956sage: timeit("add_rationals('2/3', '-5/21')")2957625 loops, best of 3: 1.29 Âµs per loop2958sage: timeit('2/3 - 5/21')2959625 loops, best of 3: 2.16 Âµs per loop2960\end{lstlisting}2961Here's a simplistic check that we probably didn't screw up and2962introduce any memory leaks. (Go up to the code and comment out some2963frees to see how this changes.)29642965%skip2966\begin{lstlisting}2967sage: print get_memory_usage()2968sage: timeit("add_rationals('9012038409238411/13',2969'-4/9082309482309487')",number=10^6)2970sage: get_memory_usage()2971917.562529721000000 loops, best of 3: 1.72 Âµs per loop2973917.56252974\end{lstlisting}29752976Finally, we will go line by line through the code and explain exactly2977what is going on and why. TODO297829792980\section{Important Cython Language Constructions}2981In this section we systematically go through {\em the most important}2982standard Cython language constructions. We will not talk about using2983numpy from Cython, dynamic memory allocation, or subtleties of the C2984language in this section. Instead we cover declaring and using cdef'd2985variables, explicit type casts, declaring external data types and2986functions, defining new Cython cdef'd functions, and declaring new2987Cython cdef'd classes that can have C attributes.29882989\subsection{Declaring Cython Variables Using {\tt cdef}}2990\begin{lstlisting}2991cdef type_name variable_name1, variable_name2, ...2992\end{lstlisting}2993The single most important statement that Cython adds to Python is2994\begin{lstlisting}2995cdef type_name2996\end{lstlisting}2997This allows you to declare a variable to have a type. The possibilities for the type include:2998\begin{itemize}2999\item C data type: int, float, double, char. Each can be modified by: short, long, signed, unsigned.3000\item Certain Python types, including: list, dict, str, object (=Python object), etc.3001\item Name of a known {\tt cdef class} (see below). You may have to cimport the class.3002\item More complicated C/C++ data types: struct, C++ class, typedef, etc., that have been declared using some other method described below.3003\end{itemize}30043005\begin{lstlisting}3006%cython3007def C_type_example():3008# ^ = exclusive or -- no preparsering in Cython!3009cdef int n=5/3, x=2^33010cdef long int m=9082309482394893943011cdef float y=4.59693012cdef double z=2.133013cdef char c='c'3014cdef char* s="a C string"3015print n, x, m, y, z, c, s3016\end{lstlisting}3017When we run the above function, we get the following. Note the lack3018of preparsing, and that the char variable {\tt c} is treated as a3019number.3020\begin{lstlisting}3021sage: C_type_example()30221 1 908230948239489394 4.59689998627 2.13 99 a C string3023\end{lstlisting}30243025\begin{lstlisting}3026%cython3027def type_example2(x, y):3028cdef list v3029cdef dict z3030v = x3031z = y3032\end{lstlisting}30333034\begin{lstlisting}3035sage: type_example2([1,2], {'a':5})3036sage: type_example2(17, {'a':5})3037Traceback (most recent call last):3038...3039TypeError: Expected list, got sage.rings.integer.Integer3040sage: type_example2([1,2], 17)3041Traceback (most recent call last):3042...3043TypeError: Expected dict, got sage.rings.integer.Integer3044\end{lstlisting}304530463047We can also define a new Cython cdef'd class and use that type in cdef:3048\begin{lstlisting}3049%cython3050cdef class MyNumber:3051cdef int n3052def __init__(self, k):3053self.n = k3054def __repr__(self):3055return repr(self.n)30563057def type_example3(MyNumber z):3058cdef MyNumber w = MyNumber(5)3059print z, w3060\end{lstlisting}3061Use it:30623063\begin{lstlisting}3064sage: type_example3(MyNumber(10))306510 53066sage: type_example3(10)3067Traceback (most recent call last):3068...3069TypeError: Argument 'z' has incorrect type...3070\end{lstlisting}30713072Terrifying caveat!: Any type can also be {\tt None}, which will cause3073horrible segfaults. Use \verb|not None| to deal with this.30743075\begin{lstlisting}3076sage: type_example3(None) # this could be very bad3077None 53078\end{lstlisting}30793080This is how to be safe:3081\begin{lstlisting}3082%cython3083cdef class MyNumber:3084cdef int n3085def __init__(self, k):3086self.n = k3087def __repr__(self):3088return repr(self.n)30893090def type_example4(MyNumber z not None):3091cdef MyNumber w = MyNumber(5)3092print z, w3093\end{lstlisting}30943095Now None is rejected:3096\begin{lstlisting}3097sage: type_example4(MyNumber(4))30984 53099sage: type_example4(None)3100Traceback (most recent call last):3101...3102TypeError: Argument 'z' has incorrect type ...3103\end{lstlisting}310431053106For the Cython source code of Sage integers, in the Sage library see3107\verb|rings/integer.pxd| and \verb|rings/integer.pyx|. Also, browse3108\verb|libs/gmp/| for the definition of functions such as3109\verb|mpz_set| below.31103111\begin{lstlisting}3112%cython3113from sage.rings.integer cimport Integer # note the cimport!3114def unsafe_mutate(Integer n, Integer m):3115mpz_set(n.value, m.value)3116\end{lstlisting}31173118\begin{lstlisting}3119sage: n = 153120sage: print n, id(n)312115 548527523122sage: unsafe_mutate(n, 2011)3123sage: print n, id(n)31242011 548527523125\end{lstlisting}312631273128\subsection{Explicit casts}3129\begin{lstlisting}3130<data_type> foo3131\end{lstlisting}31323133If you need to force the compiler to treat a variable of one data type3134as another, you have to use an explicit cast. In Java and C/C++ you3135would use parenthesis around a type name, as follows:31363137\begin{lstlisting}3138int i = 1;3139long j = 3;3140i = (int)j;3141\end{lstlisting}31423143In Cython, you use angle brackets (note: in Cython this particular3144cast isn't strictly necessary to get the code to compile, but in Java3145it is):31463147\begin{lstlisting}3148%cython3149cdef int i = 13150cdef long j = 33151i = <int> j3152print i3153\end{lstlisting}31543155Here's an example where we convert a Python string to a char* (i.e., a3156pointer to an array of characters), then change one of the characters,3157thus mutating an immutable string.31583159\begin{lstlisting}3160%cython3161def unsafe_mutate_str(bytes s, n, c):3162cdef char* t = <char*>s3163t[n] = ord(c)3164\end{lstlisting}3165Try it out:3166\begin{lstlisting}3167sage: s = 'This is an immutable string.'3168sage: print s, id(s), hash(s)3169This is an immutable string. 72268152 -56549257170928878183170sage: unsafe_mutate_str(s, 9, ' ')3171sage: unsafe_mutate_str(s, 11, ' ')3172sage: unsafe_mutate_str(s, 12, ' ')3173print s, id(s), hash(s)3174This is a mutable string. 72268152 -56549257170928878183175sage: hash('This is a mutable string.')3176-74761660604858060823177\end{lstlisting}31783179\subsection{Declaring External Data Types and Functions}31803181In order for Cython to make use of a function or data type defined in3182external C/C++ library, Cython has to {\em explicitly} be told what3183the input and output types are for that function and what the function3184should be called. Cython will then generate appropriate C/C++ code3185and conversions based on these assumptions. There are a large number3186of files in Sage and Cython itself that declare all the functions3187provided by various standard libraries, but sometimes you want to make3188use of a function defined elsewhere, e.g., in your own C/C++ library,3189so you have to declare things yourself. The purpose of the following3190examples is to illustrate how to do this. It is also extremely useful3191to look at the Sage library source code for thousands of additional3192nontrivial working examples.3193\begin{lstlisting}3194cdef extern from "filename.h":3195declarations ...3196\end{lstlisting}31973198The following examples illustrates several different possible3199declarations. We'll describe each line in detail. This first example3200declares a single type of round function on doubles -- it's as3201straightforward as it gets.3202\begin{lstlisting}3203%cython3204cdef extern from "math.h":3205double round(double)32063207def f(double n):3208return round(n)3209\end{lstlisting}3210Try it out:3211\begin{lstlisting}3212sage: f(10.53595)321311.03214\end{lstlisting}32153216Now suppose we want a version of round that returns a long. By3217consulting the man page for round, we find that there is a round3218function declared as follows:3219\begin{lstlisting}3220long int lround(double x);3221\end{lstlisting}3222We can declare it exactly like the above, or we can use a C ``name3223specifier'', which let's us tell Cython we want to call the function3224{\tt round} in our Cython code, but when Cython generates code it should3225actually emit {\tt lround}. This is what we do below.32263227\begin{lstlisting}3228%cython3229cdef extern from "math.h":3230long int round "lround"(double)32313232def f(double n):3233return round(n)3234\end{lstlisting}32353236\begin{lstlisting}3237sage: f(10.53595)3238113239\end{lstlisting}32403241Another case when using C name specifiers is useful if you want to be3242able to call both a C library version of a function and a builtin3243Python function with the same name.32443245\begin{lstlisting}3246%cython3247cdef extern from "stdlib.h":3248int c_abs "abs"(int i)32493250def myabs(n):3251print abs(n)3252print c_abs(n)3253\end{lstlisting}3254Now use it:3255\begin{lstlisting}3256sage: myabs(-10)3257103258103259\end{lstlisting}32603261We can also declare data types and variables using {\tt cdef extern}.3262To write the code below, I used the {\tt man} command on my computer3263several times on each referenced function. I knew the relevant3264functions because I read a book on the C programming language when I3265was a freshman; learning the basics of the C programming language and3266standard libraries is a very good idea if you want to be able to make3267effective use of Cython... or computers in general, since most systems3268programming is done in C.32693270Coming up with the declarations below is a little bit of an art form,3271in that they are not exactly what is given from the man pages, though3272they are close. Just realize that the declarations you give here do3273exactly one thing: they inform Cython about what C code it should3274generate, e.g., it will convert the3275