CoCalc Public Fileswww / talks / pyrex / pyrex.tex
Author: William A. Stein
1\documentclass[12pt, landscape]{article}
2%\documentclass[slides, landscape]{article}
3\usepackage{fullpage}
4\usepackage{fancybox}
5\author{William Stein\\
6{\tt http://modular.fas.harvard.edu}}
7\date{Massachussetts Python Meeting: November 11, 2004}
8\include{macros}
9\renewcommand{\dual}{\vee}
10\usepackage[hypertex]{hyperref}
11
12\setlength{\fboxsep}{1em}
13\setlength{\parindent}{0cm}
14\usepackage{pstricks}
15\usepackage{graphicx}
16\newrgbcolor{dblue}{0 0 0.8}
17\renewcommand{\hd}[1]{\begin{center}\LARGE\bf\dblue #1\vspace{-2.5ex}\end{center}}
18\newrgbcolor{dred}{0.7 0 0}
19\newrgbcolor{dgreen}{0 0.3 0}
20\newcommand{\rd}[1]{{\bf \dred #1}}
21\renewcommand{\defn}[1]{{\bf \dblue #1}}
22\newrgbcolor{purple}{0.3 0 0.3}
23\renewcommand{\magma}{{\purple\sc Magma}}
24
25\newcommand{\bd}[1]{{\bf\dred #1}}
26
27\newcommand{\mysection}[1]{\newpage\begin{center}{
28\section{\Huge{\dblue #1}}}\end{center}}
29
30\newcommand{\mysubsection}[1]{{\newpage\subsection{\Huge{\dblue #1}}}}
31
32
33\theoremstyle{plain}
34\newtheorem{listing}[theorem]{Listing}
35%\usepackage{python}
36\title{\Huge\blue A User's Perspective On Pyrex}
37
38\begin{document}
39\sf\Large
40\maketitle
41
42\begin{abstract}
43\vspace{1ex}
44\LARGE
45
46Pyrex is a language for writing Python extension classes that run at C
47speed.  In this 15 minute talk I will briefly describe Pyrex, why it
48useful for the software I'm currently writing, and how it compares to
49other systems for my application.
50
51The Pyrex manual describes Pyrex as
52\begin{center}
53\Huge\bf A smooth blend of the finest Python\\with the
54unsurpassed power of raw C.
55\end{center}
56
57\end{abstract}
58
59\mysection{Pyrex}
60
61\begin{itemize}
62\item The author of Pyrex is \rd{Greg Ewing}, of University of Canterbury, New Zealand.
63The Pyrex web page is:
64\begin{center}
65{\tt http://www.cosc.canterbury.ac.nz/\~{ }greg/python/Pyrex/}
66\end{center}
67\item There is an active mailing list:
68\begin{center}
69{\tt http://lists.copyleft.no/mailman/listinfo/pyrex}
70\end{center}
71
72\item Pyrex is a language very similar to Python for writing Python extension
73modules.  It is easy to use C functions from Pyrex, and much Python
74code also works in Pyrex.
75
76\item I'm just an enthusiastic Pyrex user.
77I have no connection with Pyrex development, and have not read
78  the source code, which appears to be about 13000 lines of pure
79  Python.
80
81\end{itemize}
82
83\mysection{MANIN: A Python Modular Forms package}
84
85\rd{MANIN:} I'm writing a Python package for number theorists; it's like
86the SciPy package, but for number theory (and cryptography?).
87The first version is aimed mainly at computing
88with modular forms'', which are holomorphic functions whose power
89series expansion encodes deep arithmetic information.
90I will say nothing further about modular forms here, except to
91remark that computing them involves sparse and dense linear algebra over
92exact base fields (e.g., finite fields, the rational numbers, and
93number fields).
94For my application, \rd{speed} and \rd{memory} {\bf efficiency} are
95crucial, and must meet or exceed what is already available in the
96competing modular forms programs that are partly written in C or C++.
97\vspace{2ex}
98
99\rd{MAGMA:} Currently the non-profit Australian computer algebra
100system MAGMA dominates in number theory computation.  It's expensive,
101and has several drawbacks (e.g., closed source, no Pickle'' of
102objects, no user defined classes, no graphics).  MAGMA is over 2
103million lines of C code and a few hundred thousand lines of code
104written in the MAGMA language.  Thus MAGMA is written in a hybrid of C
105and an interpreted language.  This hybrid approach seems to be the best way
106to write a fast sophisticated computer algebra system.
107
108\vspace{2ex}
109\rd{NZMATH:} A Japanese mathematician is trying to make a serious number
110theory package in \rd{pure} Python. See {\tt http://tnt.math.metro-u.ac.jp/nzmath/}.
111Pure Python seems like a  bad idea to me, since it is too slow.
112
113
115\begin{itemize}
116\item It's difficult to create a number theory package in pure Python,
117  because Python is too slow for \rd{certain} things.  For example,
118  multiplying two matrices over the integers modulo a small prime~$p$
119  is fast in C and slow in pure Python, as we will see below.
120
121\item I will now give a complete example of the simple
122  matrix-multiplication algorithm in pure Python and Pyrex.  (There
123  are better asymptotically fast divide and conquer'' matrix
124  multiplication algorithms that grew out of work of Strassen, which I
125  will not discuss.)
126
127\item Excellent support for exception handling in extension modules.
128\end{itemize}
129
130\mysection{Pure Python Implementation of Matrix Multiply}
131We multiply two square matrices modulo a small prime $p$.  I have
132not included error checking.
133
134{\small
135\begin{verbatim}
136# mult1.py
137class Matrix:
138    def __init__(self, p, n, entries=None):
139        """p -- prime
140           n -- positive integer
141           entries -- entries of the matrix (defaults to None, which means 0 matrix).
142        """
143        self.__n = n; self.__p = p
144        if entries != None:
145            self.__entries = [x % p for x in entries]
146        else:
147            self.__entries = [0 for _ in range(n*n)]
148
149    def __repr__(self):
150        s = ""; n = self.__n
151        for i in range(n):
152            for j in range(n):
153                s += "%s, "%self.__entries[n*i + j]
154            s += "\n"
155        return s
156
157    def __mul__(self, other):
158        ans = []; n = self.__n
159        for i in range(n):
160            for j in range(n):
161                v = [self.__entries[i*n+k] * other.__entries[k*n+j] for k in range(n)]
162                ans.append(sum(v) % self.__p)
163        return Matrix(self.__p, n, ans)
164\end{verbatim}
165}
166Example usage:
167\begin{verbatim}
168      >>> import mult1
169      >>> a = mult1.Matrix(101,3,range(9))
170      >>> a
171      0, 1, 2,
172      3, 4, 5,
173      6, 7, 8,
174      >>> a*a
175      15, 18, 21,
176      42, 54, 66,
177      69, 90, 10,
178\end{verbatim}
179
180
181\mysection{Pyrex Implementation of Matrix Multiply}
182I copied {\tt mult1.py} to {\tt mult2.pyx}, and made
183the appropriate modifications.
184{\small
185\begin{verbatim}
186
187# mult2.pyx
188cdef extern from "Python.h":    # we will use these C functions below
189    void* PyMem_Malloc(int)
190    void PyMem_Free(void *p)
191
192cdef class Matrix:
193    cdef int *entries
194    cdef int p, n
195
196    def __new__(self, int p, int n, entries=None):
197        self.p = p; self.n = n
198        self.entries = <int*> PyMem_Malloc(sizeof(int)*n*n)   # cast to int pointer
199
200    def __dealloc__(self):
201        PyMem_Free(self.entries)         # using a C function
202
203    def __init__(self, int p, int n, entries=None):
204        """ p -- prime
205            n -- positive integer
206            entries -- entries of the matrix (defaults to None, which means 0 matrix).
207        """
208        cdef int i, j, k, x
209        if entries != None:
210            for i from 0 <= i < self.n:         # Pyrex's version of for loop;
211                for j from 0 <= j < self.n:     # get converted to C with no Python calls
212                    k = i*self.n + j
213                    x = entries[k] % p
214                    if x<0: x = x + p
215                    self.entries[k] = x
216        else:
217            for i from 0 <= i < self.n:
218                for j from 0 <= j < self.n:
219                    self.entries[i*self.n + j] = 0
220
221    def __repr__(self):
222        cdef int i, j, n
223        s = ""
224        n = self.n
225        for i from 0 <= i < n:
226            for j from 0 <= j < n:
227                s = s + "%s, "%self.entries[n*i + j]
228            s = s + "\n"
229        return s
230
231    cdef Matrix __mul__(Matrix self, Matrix B):
232        cdef int s, i, j, k, n
233        cdef Matrix ans
234        ans = Matrix(self.p, self.n)
235        n = self.n
236        for i from 0 <= i < n:
237            for j from 0 <= j < n:
238                # The i,j entry of the product
239                s = 0
240                for k from 0 <= k < n:
241                    s = (s + (self.entries[i*n+k] * B.entries[k*n+j])) % self.p
242                    if s < 0: s = s + self.p
243                    ans.entries[i*n+j] = s
244        return ans
245\end{verbatim}
246}
247
248We first convert {\tt mult2.pyx} to a C program, then compile
249it to a shared object extension library:
250\begin{verbatim}
251# pyrex mult2.pyx         # convert pyx file to a 677 line pure C file.
252# ls -lh mult2.c
253-rw-r--r--  1 was was 23K Nov 11 15:16 mult2.c
254# gcc -shared -O3 -fPIC -I/home/was/local/include/python2.4 mult2.c -o mult2.so
255\end{verbatim}
256{\small
257\begin{verbatim}
258mult2.c: In function __pyx_f_5mult2_6Matrix___mul__':
259mult2.c:302: warning: use of cast expressions as lvalues is deprecated   # yikes!
260mult2.c:316: warning: use of cast expressions as lvalues is deprecated   # many warnings...
261mult2.c:362: warning: use of cast expressions as lvalues is deprecated
262mult2.c: At top level:
263mult2.c:427: warning: initialization from incompatible pointer type      # yikes!
264\end{verbatim}
265}
266\begin{verbatim}
267# ls -lh mult2.so
268-rwxr-xr-x  1 was was 17K Nov 11 15:16 mult2.so
269\end{verbatim}
270
271Now we can use the mult2 module, which has the same interface as mult.
272\begin{verbatim}
273>>> import mult2
274>>> a = mult2.Matrix(101,3,range(9))
275>>> a*a
27615, 18, 21,
27742, 54, 66,
27869, 90, 10,
279\end{verbatim}
280
281\mysection{Some Generated C Code}
282{\small
283\begin{verbatim}
284/* Generated by Pyrex 0.9.3 on Thu Nov 11 15:16:11 2004 */
285#include "Python.h"
286#include "structmember.h"
287typedef struct {PyObject **p; char *s;} __Pyx_InternTabEntry; /*proto*/
288...
289/* Declarations from mult2 */
290staticforward PyTypeObject __pyx_type_5mult2_Matrix;
291struct __pyx_obj_5mult2_Matrix {
293  struct __pyx_vtabstruct_5mult2_Matrix *__pyx_vtab;
294  int (*entries);
295  int p;
296  int n;
297};
298...
299static int __pyx_f_5mult2_6Matrix___new__(PyObject *__pyx_v_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
300static int __pyx_f_5mult2_6Matrix___new__(PyObject *__pyx_v_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
301  int __pyx_v_p;
302  int __pyx_v_n;
303  PyObject *__pyx_v_entries = 0;
304  int __pyx_r;
305  static char *__pyx_argnames[] = {"p","n","entries",0};
306  __pyx_v_entries = __pyx_k1;
307  if (!PyArg_ParseTupleAndKeywords(__pyx_args, __pyx_kwds, "ii|O", __pyx_argnames, &__pyx_v_p, &__pyx_v_n, &__pyx_v_entries)) return -1;
308  Py_INCREF(__pyx_v_self);
309  Py_INCREF(__pyx_v_entries);
310  /* "/home/was/talks/pyrex/mult2.pyx":22 */
311  ((struct __pyx_obj_5mult2_Matrix *)__pyx_v_self)->n = __pyx_v_n;
312...
313\end{verbatim}
314}
315
316
317\mysection{Comparing Timings}
318I also wrote a \rd{Psyco} version {\tt mult3.py}, which is {\tt mult1.py}, but
319with the line {\tt import pysco; psyco.full()} appended at the top.
320
321I used the following function to test the speed of each matrix multiply routine.
322\begin{verbatim}
323
324import timeit
325def matmul(tries=10, p=97, n=20):
326    T0 = timeit.Timer("a=mult1.Matrix(%s,%s,range(%s*%s)); a*a"%\
327                      (p,n,n,n), "import mult1")
328    T1 = timeit.Timer("a=mult2.Matrix(%s,%s,range(%s*%s)); a*a"%\
329                      (p,n,n,n), "import mult2")
330    T2 = timeit.Timer("a=mult3.Matrix(%s,%s,range(%s*%s)); a*a"%\
331                      (p,n,n,n), "import mult3")
332    t0 = T0.timeit(tries)
333    t1 = T1.timeit(tries)
334    t2 = T2.timeit(tries)
335    print "Pure Python: %s\n\nPyrex: %s (Speedup: %s)\n"%(t0, t1, t0/t1)
336    print "Psyco: %s (Speedup: %s)"%(t2, t0/t2)
337\end{verbatim}
338
339\newpage
340Here are some timing results (under Python 2.3, since my Pysco install
341only supports Python 2.3), on my Pentium-M 1.6Ghz laptop.
342\begin{verbatim}
343>>> prof.matmul(1000)
344Pure Python: 4.03225588799
345Pyrex: 0.380643129349 (Speedup: 10.5932711695)
346Psyco: 3.22930121422 (Speedup: 1.24864657104)
347>>> prof.matmul(200, n=30)
348Pure Python: 2.62041187286
349Pyrex: 0.219752788544 (Speedup: 11.9243623265)
350Psyco: 2.07065701485 (Speedup: 1.26549778842)
351>>> prof.matmul(200, n=50)
352Pure Python: 11.7379591465
353Pyrex: 0.95321393013 (Speedup: 12.3140868754)
354Psyco: 9.31631708145 (Speedup: 1.25993555649)
355\end{verbatim}
356We gain a full order of magnitude (base 10) improvement by using Pyrex!
357That's very significant for my application, where matrix multiply is
358at the core of many algorithms.  Also, Psyco doesn't seem to help much (this
359could be my ignorance about how best to use Pysco!).
360
361\mysection{Comparison with Psyco, SWIG, Boost, C}
362\begin{itemize}
363
364\item{}[{\bf Psyco}] Psyco is a just in time compiler''.  As the above timings
365  indicate, it is not as fast as Pyrex in our matrix multiplication
366  example.  Also, the Pysco web page says Psyco currently uses a lot
367  of memory. It only runs on Intel 386-compatible processors (under
368  any OS) right now.''  These are both major obstructions for my application.
369
370\item{}[{\bf SWIG}] SWIG wraps C++ classes as C-extension modules; the classes are
371  also wrapped by pure Python classes, which can be a performance hit.
372  For example, imagine implementing integers modulo $p$ in C++/SWIG
373  versus implementing them using Pyrex.  Every multiply would be
374  costly in SWIG.  It is also painful to access Python objects from
375  C++/SWIG (one uses typemaps'').  Finally, it is difficult \rd{for me}
376  to constantly switch between programming in C++ and programming in
377  Python.
378
379  The documentation for SWIG is {\em excellent}, it is a rapidly
380  evolving system, and SWIG supports exposing C++ code to many other
381  languages besides Python.
382
383\item{}[{\bf Boost}] Boost.Python seems really neat when I read about it, but I found
384  it way too complicated for my needs when I actually use it.  I just
385  couldn't get my head around it.
386
387%\item{}[{\bf PyInline}] This project seems dead (nothing new since 2001).
388
389\item{}[{\bf C}] Implementing extensions directly in C is painful because of all
390  the explicit reference counting.  Also, again it is confusing (to
391  me) to constantly switch between programming in C and programming
392  in Python.
393\end{itemize}
394
395
396\mysection{Drawbacks to Using Pyrex}
397
398\begin{itemize}
399\item Important to be familiar with the books that come with
400  Python called \rd{The Python/C API} and \rd{Extending and
401    Embedding}.
402%Otherwise when you try to read the C code produced by
403%  Pyrex in order to do further optimization or debuging, you might not
404%  understand it.
405
406\item Several basic Python constructs are not yet supported in Pyrex.
407None of the following work:
408\begin{verbatim}
409       a += 2
410       a = [i*i for i in range(10)]
411\end{verbatim}
412
413\item No direct support for C++ libraries or C++ code.
414
415\item Pyrex does not parse C header files.  That's why we had
416{\small
417\begin{verbatim}
418cdef extern from "Python.h":
419    void* PyMem_Malloc(int)
420    void PyMem_Free(void *p)
421\end{verbatim}
422}
423at the beginning of {\tt mult2.pyx}, instead of just {\tt cdef include "Python.h"}''
424
425\item Debuging can be confusing.
426
427\item It is sometimes tricky to divide Pyrex code into many files and
428  call Pyrex C functions from one file in another file.
429
430\end{itemize}
431
432
433
434\end{document}
435
`