\documentclass[12pt, landscape]{article}
%\documentclass[slides, landscape]{article}
\usepackage{fullpage}
\usepackage{fancybox}
\author{William Stein\\
{\tt http://modular.fas.harvard.edu}}
\date{Massachussetts Python Meeting: November 11, 2004}
\include{macros}
\renewcommand{\dual}{\vee}
\usepackage[hypertex]{hyperref}
\setlength{\fboxsep}{1em}
\setlength{\parindent}{0cm}
\usepackage{pstricks}
\usepackage{graphicx}
\newrgbcolor{dblue}{0 0 0.8}
\renewcommand{\hd}[1]{\begin{center}\LARGE\bf\dblue #1\vspace{-2.5ex}\end{center}}
\newrgbcolor{dred}{0.7 0 0}
\newrgbcolor{dgreen}{0 0.3 0}
\newcommand{\rd}[1]{{\bf \dred #1}}
\renewcommand{\defn}[1]{{\bf \dblue #1}}
\newrgbcolor{purple}{0.3 0 0.3}
\renewcommand{\magma}{{\purple\sc Magma}}
\newcommand{\bd}[1]{{\bf\dred #1}}
\newcommand{\mysection}[1]{\newpage\begin{center}{
\section{\Huge{\dblue #1}}}\end{center}}
\newcommand{\mysubsection}[1]{{\newpage\subsection{\Huge{\dblue #1}}}}
\theoremstyle{plain}
\newtheorem{listing}[theorem]{Listing}
%\usepackage{python}
\title{\Huge\blue A User's Perspective On Pyrex}
\begin{document}
\sf\Large
\maketitle
\begin{abstract}
\vspace{1ex}
\LARGE
Pyrex is a language for writing Python extension classes that run at C
speed. In this 15 minute talk I will briefly describe Pyrex, why it
useful for the software I'm currently writing, and how it compares to
other systems for my application.
The Pyrex manual describes Pyrex as
\begin{center}
\Huge\bf A smooth blend of the finest Python\\with the
unsurpassed power of raw C.
\end{center}
\end{abstract}
\mysection{Pyrex}
\begin{itemize}
\item The author of Pyrex is \rd{Greg Ewing}, of University of Canterbury, New Zealand.
The Pyrex web page is:
\begin{center}
{\tt http://www.cosc.canterbury.ac.nz/\~{ }greg/python/Pyrex/}
\end{center}
\item There is an active mailing list:
\begin{center}
{\tt http://lists.copyleft.no/mailman/listinfo/pyrex}
\end{center}
\item Pyrex is a language very similar to Python for writing Python extension
modules. It is easy to use C functions from Pyrex, and much Python
code also works in Pyrex.
\item I'm just an enthusiastic Pyrex user.
I have no connection with Pyrex development, and have not read
the source code, which appears to be about 13000 lines of pure
Python.
\end{itemize}
\mysection{MANIN: A Python Modular Forms package}
\rd{MANIN:} I'm writing a Python package for number theorists; it's like
the SciPy package, but for number theory (and cryptography?).
The first version is aimed mainly at computing
with ``modular forms'', which are holomorphic functions whose power
series expansion encodes deep arithmetic information.
I will say nothing further about modular forms here, except to
remark that computing them involves sparse and dense linear algebra over
exact base fields (e.g., finite fields, the rational numbers, and
number fields).
For my application, \rd{speed} and \rd{memory} {\bf efficiency} are
crucial, and must meet or exceed what is already available in the
competing modular forms programs that are partly written in C or C++.
\vspace{2ex}
\rd{MAGMA:} Currently the non-profit Australian computer algebra
system MAGMA dominates in number theory computation. It's expensive,
and has several drawbacks (e.g., closed source, no ``Pickle'' of
objects, no user defined classes, no graphics). MAGMA is over 2
million lines of C code and a few hundred thousand lines of code
written in the MAGMA language. Thus MAGMA is written in a hybrid of C
and an interpreted language. This hybrid approach seems to be the best way
to write a fast sophisticated computer algebra system.
\vspace{2ex}
\rd{NZMATH:} A Japanese mathematician is trying to make a serious number
theory package in \rd{pure} Python. See {\tt http://tnt.math.metro-u.ac.jp/nzmath/}.
Pure Python seems like a bad idea to me, since it is too slow.
\mysection{How Pyrex is Helpful}
\begin{itemize}
\item It's difficult to create a number theory package in pure Python,
because Python is too slow for \rd{certain} things. For example,
multiplying two matrices over the integers modulo a small prime~$p$
is fast in C and slow in pure Python, as we will see below.
\item I will now give a complete example of the simple
matrix-multiplication algorithm in pure Python and Pyrex. (There
are better asymptotically fast ``divide and conquer'' matrix
multiplication algorithms that grew out of work of Strassen, which I
will not discuss.)
\item Excellent support for exception handling in extension modules.
\end{itemize}
\mysection{Pure Python Implementation of Matrix Multiply}
We multiply two square matrices modulo a small prime $p$. I have
not included error checking.
{\small
\begin{verbatim}
# mult1.py
class Matrix:
def __init__(self, p, n, entries=None):
"""p -- prime
n -- positive integer
entries -- entries of the matrix (defaults to None, which means 0 matrix).
"""
self.__n = n; self.__p = p
if entries != None:
self.__entries = [x % p for x in entries]
else:
self.__entries = [0 for _ in range(n*n)]
def __repr__(self):
s = ""; n = self.__n
for i in range(n):
for j in range(n):
s += "%s, "%self.__entries[n*i + j]
s += "\n"
return s
def __mul__(self, other):
ans = []; n = self.__n
for i in range(n):
for j in range(n):
v = [self.__entries[i*n+k] * other.__entries[k*n+j] for k in range(n)]
ans.append(sum(v) % self.__p)
return Matrix(self.__p, n, ans)
\end{verbatim}
}
Example usage:
\begin{verbatim}
>>> import mult1
>>> a = mult1.Matrix(101,3,range(9))
>>> a
0, 1, 2,
3, 4, 5,
6, 7, 8,
>>> a*a
15, 18, 21,
42, 54, 66,
69, 90, 10,
\end{verbatim}
\mysection{Pyrex Implementation of Matrix Multiply}
I copied {\tt mult1.py} to {\tt mult2.pyx}, and made
the appropriate modifications.
{\small
\begin{verbatim}
# mult2.pyx
cdef extern from "Python.h": # we will use these C functions below
void* PyMem_Malloc(int)
void PyMem_Free(void *p)
cdef class Matrix:
cdef int *entries
cdef int p, n
def __new__(self, int p, int n, entries=None):
self.p = p; self.n = n
self.entries = PyMem_Malloc(sizeof(int)*n*n) # cast to int pointer
def __dealloc__(self):
PyMem_Free(self.entries) # using a C function
def __init__(self, int p, int n, entries=None):
""" p -- prime
n -- positive integer
entries -- entries of the matrix (defaults to None, which means 0 matrix).
"""
cdef int i, j, k, x
if entries != None:
for i from 0 <= i < self.n: # Pyrex's version of for loop;
for j from 0 <= j < self.n: # get converted to C with no Python calls
k = i*self.n + j
x = entries[k] % p
if x<0: x = x + p
self.entries[k] = x
else:
for i from 0 <= i < self.n:
for j from 0 <= j < self.n:
self.entries[i*self.n + j] = 0
def __repr__(self):
cdef int i, j, n
s = ""
n = self.n
for i from 0 <= i < n:
for j from 0 <= j < n:
s = s + "%s, "%self.entries[n*i + j]
s = s + "\n"
return s
cdef Matrix __mul__(Matrix self, Matrix B):
cdef int s, i, j, k, n
cdef Matrix ans
ans = Matrix(self.p, self.n)
n = self.n
for i from 0 <= i < n:
for j from 0 <= j < n:
# The i,j entry of the product
s = 0
for k from 0 <= k < n:
s = (s + (self.entries[i*n+k] * B.entries[k*n+j])) % self.p
if s < 0: s = s + self.p
ans.entries[i*n+j] = s
return ans
\end{verbatim}
}
We first convert {\tt mult2.pyx} to a C program, then compile
it to a shared object extension library:
\begin{verbatim}
# pyrex mult2.pyx # convert pyx file to a 677 line pure C file.
# ls -lh mult2.c
-rw-r--r-- 1 was was 23K Nov 11 15:16 mult2.c
# gcc -shared -O3 -fPIC -I/home/was/local/include/python2.4 mult2.c -o mult2.so
\end{verbatim}
{\small
\begin{verbatim}
mult2.c: In function `__pyx_f_5mult2_6Matrix___mul__':
mult2.c:302: warning: use of cast expressions as lvalues is deprecated # yikes!
mult2.c:316: warning: use of cast expressions as lvalues is deprecated # many warnings...
mult2.c:362: warning: use of cast expressions as lvalues is deprecated
mult2.c: At top level:
mult2.c:427: warning: initialization from incompatible pointer type # yikes!
\end{verbatim}
}
\begin{verbatim}
# ls -lh mult2.so
-rwxr-xr-x 1 was was 17K Nov 11 15:16 mult2.so
\end{verbatim}
Now we can use the mult2 module, which has the same interface as mult.
\begin{verbatim}
>>> import mult2
>>> a = mult2.Matrix(101,3,range(9))
>>> a*a
15, 18, 21,
42, 54, 66,
69, 90, 10,
\end{verbatim}
\mysection{Some Generated C Code}
{\small
\begin{verbatim}
/* Generated by Pyrex 0.9.3 on Thu Nov 11 15:16:11 2004 */
#include "Python.h"
#include "structmember.h"
typedef struct {PyObject **p; char *s;} __Pyx_InternTabEntry; /*proto*/
...
/* Declarations from mult2 */
staticforward PyTypeObject __pyx_type_5mult2_Matrix;
struct __pyx_obj_5mult2_Matrix {
PyObject_HEAD
struct __pyx_vtabstruct_5mult2_Matrix *__pyx_vtab;
int (*entries);
int p;
int n;
};
...
static int __pyx_f_5mult2_6Matrix___new__(PyObject *__pyx_v_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static int __pyx_f_5mult2_6Matrix___new__(PyObject *__pyx_v_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
int __pyx_v_p;
int __pyx_v_n;
PyObject *__pyx_v_entries = 0;
int __pyx_r;
static char *__pyx_argnames[] = {"p","n","entries",0};
__pyx_v_entries = __pyx_k1;
if (!PyArg_ParseTupleAndKeywords(__pyx_args, __pyx_kwds, "ii|O", __pyx_argnames, &__pyx_v_p, &__pyx_v_n, &__pyx_v_entries)) return -1;
Py_INCREF(__pyx_v_self);
Py_INCREF(__pyx_v_entries);
/* "/home/was/talks/pyrex/mult2.pyx":22 */
((struct __pyx_obj_5mult2_Matrix *)__pyx_v_self)->n = __pyx_v_n;
...
\end{verbatim}
}
\mysection{Comparing Timings}
I also wrote a \rd{Psyco} version {\tt mult3.py}, which is {\tt mult1.py}, but
with the line {\tt import pysco; psyco.full()} appended at the top.
I used the following function to test the speed of each matrix multiply routine.
\begin{verbatim}
import timeit
def matmul(tries=10, p=97, n=20):
T0 = timeit.Timer("a=mult1.Matrix(%s,%s,range(%s*%s)); a*a"%\
(p,n,n,n), "import mult1")
T1 = timeit.Timer("a=mult2.Matrix(%s,%s,range(%s*%s)); a*a"%\
(p,n,n,n), "import mult2")
T2 = timeit.Timer("a=mult3.Matrix(%s,%s,range(%s*%s)); a*a"%\
(p,n,n,n), "import mult3")
t0 = T0.timeit(tries)
t1 = T1.timeit(tries)
t2 = T2.timeit(tries)
print "Pure Python: %s\n\nPyrex: %s (Speedup: %s)\n"%(t0, t1, t0/t1)
print "Psyco: %s (Speedup: %s)"%(t2, t0/t2)
\end{verbatim}
\newpage
Here are some timing results (under Python 2.3, since my Pysco install
only supports Python 2.3), on my Pentium-M 1.6Ghz laptop.
\begin{verbatim}
>>> prof.matmul(1000)
Pure Python: 4.03225588799
Pyrex: 0.380643129349 (Speedup: 10.5932711695)
Psyco: 3.22930121422 (Speedup: 1.24864657104)
>>> prof.matmul(200, n=30)
Pure Python: 2.62041187286
Pyrex: 0.219752788544 (Speedup: 11.9243623265)
Psyco: 2.07065701485 (Speedup: 1.26549778842)
>>> prof.matmul(200, n=50)
Pure Python: 11.7379591465
Pyrex: 0.95321393013 (Speedup: 12.3140868754)
Psyco: 9.31631708145 (Speedup: 1.25993555649)
\end{verbatim}
We gain a full order of magnitude (base 10) improvement by using Pyrex!
That's very significant for my application, where matrix multiply is
at the core of many algorithms. Also, Psyco doesn't seem to help much (this
could be my ignorance about how best to use Pysco!).
\mysection{Comparison with Psyco, SWIG, Boost, C}
\begin{itemize}
\item{}[{\bf Psyco}] Psyco is a ``just in time compiler''. As the above timings
indicate, it is not as fast as Pyrex in our matrix multiplication
example. Also, the Pysco web page says ``Psyco currently uses a lot
of memory. It only runs on Intel 386-compatible processors (under
any OS) right now.'' These are both major obstructions for my application.
\item{}[{\bf SWIG}] SWIG wraps C++ classes as C-extension modules; the classes are
also wrapped by pure Python classes, which can be a performance hit.
For example, imagine implementing integers modulo $p$ in C++/SWIG
versus implementing them using Pyrex. Every multiply would be
costly in SWIG. It is also painful to access Python objects from
C++/SWIG (one uses ``typemaps''). Finally, it is difficult \rd{for me}
to constantly switch between programming in C++ and programming in
Python.
The documentation for SWIG is {\em excellent}, it is a rapidly
evolving system, and SWIG supports exposing C++ code to many other
languages besides Python.
\item{}[{\bf Boost}] Boost.Python seems really neat when I read about it, but I found
it way too complicated for my needs when I actually use it. I just
couldn't get my head around it.
%\item{}[{\bf PyInline}] This project seems dead (nothing new since 2001).
\item{}[{\bf C}] Implementing extensions directly in C is painful because of all
the explicit reference counting. Also, again it is confusing (to
me) to constantly switch between programming in C and programming
in Python.
\end{itemize}
\mysection{Drawbacks to Using Pyrex}
\begin{itemize}
\item Important to be familiar with the books that come with
Python called \rd{The Python/C API} and \rd{Extending and
Embedding}.
%Otherwise when you try to read the C code produced by
% Pyrex in order to do further optimization or debuging, you might not
% understand it.
\item Several basic Python constructs are not yet supported in Pyrex.
None of the following work:
\begin{verbatim}
a += 2
a = [i*i for i in range(10)]
\end{verbatim}
\item No direct support for C++ libraries or C++ code.
\item Pyrex does not parse C header files. That's why we had
{\small
\begin{verbatim}
cdef extern from "Python.h":
void* PyMem_Malloc(int)
void PyMem_Free(void *p)
\end{verbatim}
}
at the beginning of {\tt mult2.pyx}, instead of just ``{\tt cdef include "Python.h"}''
\item Debuging can be confusing.
\item It is sometimes tricky to divide Pyrex code into many files and
call Pyrex C functions from one file in another file.
\end{itemize}
\end{document}