Cython, Part 1

In addition to the Sage-related examples and discussion below, to really learn Cython, I strongly recommend that you read straight as much of the Cython Users Guide (see http://docs.cython.org/src/userguide/) as you can trying out everything, and the other Cython documentation (see http://docs.cython.org/) as well.

Recall this example of Cython code, which illustrates speeding up a Python program using Cython by declaring types of variables (important reminder: %cython must be the first line of the input cell!)

{{{id=1| %cython def mysum(n): cdef long i, s s = 0 for i in range(1, n+1): s += i return s /// }}} {{{id=4| timeit('mysum(10^5)') /// 625 loops, best of 3: 68.8 µs per loop }}} {{{id=5| %python def pysum(n): s = 0 for i in range(1, n+1): s += i return s /// }}} {{{id=6| timeit('pysum(10^5)') /// 25 loops, best of 3: 14.3 ms per loop }}}

It's faster than it's Python counterpart, because of the explicit datatypes:

{{{id=7| 14.3/0.0688 /// 207.848837209302 }}}

Another important application of Cython is that it makes it possible to very efficiently make use of data types and functions defined in any C/C++ library.   Since there is an enormous amount of useful, fast, debugged C/C++ code out there, this is an extremely useful feature.  

Here's a first very simple example.   Type "man random" on the command line (or google it). 

{{{id=11| %cython cdef extern from "stdlib.h": # (1) long random() # (2) def random_nums(int n): # (3) cdef int i # (4) v = [random() for i in range(n)] # (5) return v /// }}} {{{id=25| random_nums(5) /// [1315705257, 1147455227, 1571270137, 1106977565, 1805149207] }}} {{{id=10| timeit('v = random_nums(10^5)') /// 125 loops, best of 3: 5.56 ms per loop }}} {{{id=8| %python import random k = 2**31-1 def py_random_nums(n): return [random.randint(0,k) for i in range(n)] /// }}} {{{id=26| py_random_nums(5) /// [317567506, 1289482476, 1766134327, 1216261810, 1427493671] }}} {{{id=13| timeit('v = py_random_nums(10^5)') /// 5 loops, best of 3: 251 ms per loop }}} {{{id=27| 251/5.56 /// 45.1438848920863 }}} {{{id=15| /// }}}

Now for a more mathematical example: arithmetic with arbitrary precision rational numbers.     The MPIR C library (which is in Sage) provides highly optimized arithmetic with arbitrary precision rationals.  (Note: MPIR and GMP are basically the same for our discussion; they are "forks" of each other, but export essentially the same functions.)  We could make use of MPIR by reading the documentation for MPIR and using "cdef extern" as above.  Fortunately, all of the necessary cdef extern declarations needed to use MPIR are already declared in Sage.   You can view all the declarations from the notebook by navigating to /src/libs/gmp [link opens in new window].

Let's use MPIR directly to create two rational numbers and add them together.  This code is complicated and illustrates all kinds of issues, so we will explain it in depth.  Once you understand all of this, you can probably deal with many issues that will come up with Cython.

{{{id=14| %cython from sage.libs.gmp.all cimport * # (1) def add_rationals(bytes a, bytes b): # (2) cdef mpq_t x, y, z # (3) mpq_init(x); mpq_init(y); mpq_init(z) # (4) mpq_set_str(x, a, 10) # base 10 string # (5) mpq_set_str(y, b, 10) mpq_add(z, x, y) # (6) cdef int n = (mpz_sizeinbase (mpq_numref(z), 10) # (7) + mpz_sizeinbase (mpq_denref(z), 10) + 3) cdef char* s = sage_malloc(sizeof(char)*n) # (8) if not s: raise MemoryError # (9) cdef bytes c = mpq_get_str(s, 10, z) # (10) mpq_clear(x); mpq_clear(y); mpq_clear(z) # (11) sage_free(s) # (12) return c /// }}} {{{id=23| add_rationals('2/3', '-5/21') /// '3/7' }}} {{{id=22| 2/3 - 5/21 /// 3/7 }}} {{{id=28| add_rationals('1/29048203984092834823049', '-394/29302938402384092834') /// '-11444963066794174536188472/851197732045760533660225724673976778930866' }}} {{{id=20| timeit("add_rationals('2/3', '-5/21')") /// 625 loops, best of 3: 1.29 µs per loop }}} {{{id=21| timeit('2/3 - 5/21') /// 625 loops, best of 3: 2.16 µs per loop }}}

A simplistic check that we probably didn't screw up and introduce any memory leaks.  (Go up to the code and comment out some frees to see how this changes.)

{{{id=16| print get_memory_usage() timeit("add_rationals('9012038409238411/13', '-4/9082309482309487')",number=10^6) get_memory_usage() /// 917.5625 1000000 loops, best of 3: 1.72 µs per loop 917.5625 }}} {{{id=17| /// }}} {{{id=18| /// }}}