Sharedwww / talks / pyrex / pyrex.texOpen in CoCalc
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
46
Pyrex is a language for writing Python extension classes that run at C
47
speed. In this 15 minute talk I will briefly describe Pyrex, why it
48
useful for the software I'm currently writing, and how it compares to
49
other systems for my application.
50
51
The Pyrex manual describes Pyrex as
52
\begin{center}
53
\Huge\bf A smooth blend of the finest Python\\with the
54
unsurpassed 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.
63
The 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
73
modules. It is easy to use C functions from Pyrex, and much Python
74
code also works in Pyrex.
75
76
\item I'm just an enthusiastic Pyrex user.
77
I 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
86
the SciPy package, but for number theory (and cryptography?).
87
The first version is aimed mainly at computing
88
with ``modular forms'', which are holomorphic functions whose power
89
series expansion encodes deep arithmetic information.
90
I will say nothing further about modular forms here, except to
91
remark that computing them involves sparse and dense linear algebra over
92
exact base fields (e.g., finite fields, the rational numbers, and
93
number fields).
94
For my application, \rd{speed} and \rd{memory} {\bf efficiency} are
95
crucial, and must meet or exceed what is already available in the
96
competing 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
100
system MAGMA dominates in number theory computation. It's expensive,
101
and has several drawbacks (e.g., closed source, no ``Pickle'' of
102
objects, no user defined classes, no graphics). MAGMA is over 2
103
million lines of C code and a few hundred thousand lines of code
104
written in the MAGMA language. Thus MAGMA is written in a hybrid of C
105
and an interpreted language. This hybrid approach seems to be the best way
106
to write a fast sophisticated computer algebra system.
107
108
\vspace{2ex}
109
\rd{NZMATH:} A Japanese mathematician is trying to make a serious number
110
theory package in \rd{pure} Python. See {\tt http://tnt.math.metro-u.ac.jp/nzmath/}.
111
Pure Python seems like a bad idea to me, since it is too slow.
112
113
114
\mysection{How Pyrex is Helpful}
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}
131
We multiply two square matrices modulo a small prime $p$. I have
132
not included error checking.
133
134
{\small
135
\begin{verbatim}
136
# mult1.py
137
class 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
}
166
Example 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}
182
I copied {\tt mult1.py} to {\tt mult2.pyx}, and made
183
the appropriate modifications.
184
{\small
185
\begin{verbatim}
186
187
# mult2.pyx
188
cdef extern from "Python.h": # we will use these C functions below
189
void* PyMem_Malloc(int)
190
void PyMem_Free(void *p)
191
192
cdef 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
248
We first convert {\tt mult2.pyx} to a C program, then compile
249
it 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}
258
mult2.c: In function `__pyx_f_5mult2_6Matrix___mul__':
259
mult2.c:302: warning: use of cast expressions as lvalues is deprecated # yikes!
260
mult2.c:316: warning: use of cast expressions as lvalues is deprecated # many warnings...
261
mult2.c:362: warning: use of cast expressions as lvalues is deprecated
262
mult2.c: At top level:
263
mult2.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
271
Now 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
276
15, 18, 21,
277
42, 54, 66,
278
69, 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"
287
typedef struct {PyObject **p; char *s;} __Pyx_InternTabEntry; /*proto*/
288
...
289
/* Declarations from mult2 */
290
staticforward PyTypeObject __pyx_type_5mult2_Matrix;
291
struct __pyx_obj_5mult2_Matrix {
292
PyObject_HEAD
293
struct __pyx_vtabstruct_5mult2_Matrix *__pyx_vtab;
294
int (*entries);
295
int p;
296
int n;
297
};
298
...
299
static int __pyx_f_5mult2_6Matrix___new__(PyObject *__pyx_v_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
300
static 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}
318
I also wrote a \rd{Psyco} version {\tt mult3.py}, which is {\tt mult1.py}, but
319
with the line {\tt import pysco; psyco.full()} appended at the top.
320
321
I used the following function to test the speed of each matrix multiply routine.
322
\begin{verbatim}
323
324
import timeit
325
def 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
340
Here are some timing results (under Python 2.3, since my Pysco install
341
only supports Python 2.3), on my Pentium-M 1.6Ghz laptop.
342
\begin{verbatim}
343
>>> prof.matmul(1000)
344
Pure Python: 4.03225588799
345
Pyrex: 0.380643129349 (Speedup: 10.5932711695)
346
Psyco: 3.22930121422 (Speedup: 1.24864657104)
347
>>> prof.matmul(200, n=30)
348
Pure Python: 2.62041187286
349
Pyrex: 0.219752788544 (Speedup: 11.9243623265)
350
Psyco: 2.07065701485 (Speedup: 1.26549778842)
351
>>> prof.matmul(200, n=50)
352
Pure Python: 11.7379591465
353
Pyrex: 0.95321393013 (Speedup: 12.3140868754)
354
Psyco: 9.31631708145 (Speedup: 1.25993555649)
355
\end{verbatim}
356
We gain a full order of magnitude (base 10) improvement by using Pyrex!
357
That's very significant for my application, where matrix multiply is
358
at the core of many algorithms. Also, Psyco doesn't seem to help much (this
359
could 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.
407
None 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}
418
cdef extern from "Python.h":
419
void* PyMem_Malloc(int)
420
void PyMem_Free(void *p)
421
\end{verbatim}
422
}
423
at 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