Math 581d: Understanding How to Use Cython in Mathematics in Practice



<h1 style="text-align: center;">Math 581d: Understanding How to Use Cython in Mathematics in Practice</h1>
<p>The Problem: Implement arithmetic with "golden elements", i.e., elements of $R = \ZZ[\gamma]$ where $\gamma=\frac{1+\sqrt{5}}{2}$ is the golden ratio.&nbsp;</p>
<p>We will follow through this example at many, many different levels of detail and optimization, starting first with "playing around in the notebook", then making a standalone Cython/Python package, and culminating with adding new code to the Sage library. &nbsp; Your homework (number 7) will be to generalize these examples to the field $K=\QQ(\gamma)$.&nbsp;</p>
<p>Math:</p>
<ul>
<li>The elements of $R$ are all numbers of the form $a+b\gamma$ with $a,b\in\ZZ$. &nbsp;</li>
<li>Addition: $(a+b\gamma) + (c+d\gamma) = (a+c) + (b+d)\gamma$.</li>
<li>Multiplication? &nbsp;We have $\gamma^2 = \gamma+1$, so $$(a+b\gamma) \cdot (c+d\gamma) = (ac + bd) + (ad+bc+bd)\gamma.$$</li>
</ul>

<p>Figure out multiplication using Sage's commutative algebra (or you could just do it with pencil and paper...):</p>

{{{id=7|
R.<a,b,c,d,x> = QQ[]
S.<a,b,c,d,gamma> = R.quotient([x^2 - x - 1])
(a+b*gamma)*(c+d*gamma)
///
b*c*gamma + a*d*gamma + b*d*gamma + a*c + b*d
}}}

{{{id=12|

///
}}}

<h2>Python Class</h2>
<p>Before using Cython, let's make a first easy Python implementation.</p>

{{{id=11|
# Simple Python version
class Elt:    
    def __init__(self, a, b):
        self.a = a; self.b = b
    def __repr__(self): 
        return '%s + %s*gamma'%(self.a, self.b)
    def __add__(left, right):
        return Elt(left.a + right.a, left.b + right.b)
    def __mul__(left, right): 
        a, b, c, d = left.a, left.b, right.a, right.b
        return Elt(a*c + b*d, b*c + a*d + b*d)
///
}}}

<p>Let's try out our new class:</p>

{{{id=10|
z = Elt(3, 5); w = Elt(2, 8); z, w
///
(3 + 5*gamma, 2 + 8*gamma)
}}}

{{{id=42|
type(z)
///
<type 'instance'>
}}}

{{{id=9|
z + w
///
5 + 13*gamma
}}}

{{{id=16|
z * w
///
46 + 74*gamma
}}}

<p>Awesome! &nbsp; Back make sure you understand that the below should fail (and why). &nbsp;We will officially not worry about such things until probably next week (when we learn about "the coercion model").</p>

{{{id=18|
z + 1
///
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_62.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("eiArIDE="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "/private/var/folders/7y/7y-O1iZOGTmMUMnLq7otq++++TI/-Tmp-/tmp8Zx42s/___code___.py", line 3, in <module>
    exec compile(u'z + _sage_const_1 ' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "", line 7, in __add__
    
  File "element.pyx", line 306, in sage.structure.element.Element.__getattr__ (sage/structure/element.c:2666)
  File "parent.pyx", line 272, in sage.structure.parent.getattr_from_other_class (sage/structure/parent.c:2840)
  File "parent.pyx", line 170, in sage.structure.parent.raise_attribute_error (sage/structure/parent.c:2611)
AttributeError: 'sage.rings.integer.Integer' object has no attribute 'a'
}}}

<p>What about speed? &nbsp;It's not so bad, really:</p>

{{{id=17|
timeit('z+w', number=10^5)
///
100000 loops, best of 3: 1.89 µs per loop
}}}

{{{id=22|
timeit('z*w', number=10^5)
///
100000 loops, best of 3: 2.47 µs per loop
}}}

{{{id=28|

///
}}}

<p>Everything else we will do is about different versions of the above fairly simple class.</p>

{{{id=27|

///
}}}

<h2>Next, our First Naive <span style="text-decoration: underline;">Cython Version</span></h2>
<p>The <em>only</em>&nbsp;change is that I changed the name of the class and put %cython at the top of the block.&nbsp;</p>

{{{id=26|
%cython

# First naive (and pointless) use of Cython.
class Elt2:    
    def __init__(self, a, b):
        self.a = a; self.b = b
    def __repr__(self): 
        return '%s + %s*gamma'%(self.a, self.b)
    def __add__(left, right):
        return Elt2(left.a + right.a, left.b + right.b)
    def __mul__(left, right): 
        a, b, c, d = left.a, left.b, right.a, right.b
        return Elt2(a*c + b*d, b*c + a*d + b*d)
///
}}}

<p>Let's try it out:</p>

{{{id=25|
z = Elt2(3, 5); w = Elt2(2, 8); z, w
///
(3 + 5*gamma, 2 + 8*gamma)
}}}

{{{id=41|
type(z)
///
<type 'instance'>
}}}

{{{id=24|
z + w
///
5 + 13*gamma
}}}

{{{id=23|
z*w
///
46 + 74*gamma
}}}

<p>The timings are identical to the Python version.</p>

{{{id=33|
timeit('z+w', number=10^5)
///
100000 loops, best of 3: 1.85 µs per loop
}}}

{{{id=34|
timeit('z*w', number=10^5)
///
100000 loops, best of 3: 2.29 µs per loop
}}}

{{{id=39|

///
}}}

<h2>A Less Naive Cython Version: Make a Cdef'd Extension Class</h2>
<p>This looks for now just like what we did above, but we put "cdef class" instead of just "class". &nbsp;There are some major changes behind the scenes though.</p>

{{{id=38|
%cython

# A cdef'd class
cdef class Elt3:    
    def __init__(self, a, b):
        self.a = a; self.b = b
    def __repr__(self): 
        return '%s + %s*gamma'%(self.a, self.b)
    def __add__(left, right):
        return Elt3(left.a + right.a, left.b + right.b)
    def __mul__(left, right): 
        a, b, c, d = left.a, left.b, right.a, right.b
        return Elt3(a*c + b*d, b*c + a*d + b*d)
///
}}}

{{{id=37|
z = Elt3(3, 5); w = Elt3(2, 8); z, w
///
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_73.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("eiA9IEVsdDMoMywgNSk7IHcgPSBFbHQzKDIsIDgpOyB6LCB3"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "/private/var/folders/7y/7y-O1iZOGTmMUMnLq7otq++++TI/-Tmp-/tmp3JJIKH/___code___.py", line 3, in <module>
    exec compile(u'z = Elt3(_sage_const_3 , _sage_const_5 ); w = Elt3(_sage_const_2 , _sage_const_8 ); z, w' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage80_spyx_0.pyx", line 9, in _Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage80_spyx_0.Elt3.__init__ (_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage80_spyx_0.c:480)
    self.a = a; self.b = b
AttributeError: '_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage80_spyx_0.Elt3' object has no attribute 'a'
}}}

<p><strong>Dang: </strong>&nbsp;The first difference is that it doesn't work at all! &nbsp;</p>
<p>cdef'd classes are different than normal Python classes in a few ways:</p>
<ol>
<li>All their attributes <em>must</em>&nbsp;be declared at compile time. &nbsp;The plus is you can declare C data types (like int, long, etc.), which is impossible with Python classes.</li>
<li>Only single inheritance is supported.</li>
<li>You can declare "cdef'd methods", which are methods of the cdef'd class that can take C data types as input, return C data types as output, and for which calling them is extremely efficient. &nbsp; (Calling a Python function usually involves creating a tuple and/or dict, putting the inputs in the tuple, unpacking the tuple, etc. &nbsp;Lots of overhead and red tape!)</li>
</ol>

{{{id=44|
%cython

cdef class Elt4:    
    cdef object a, b
    def __init__(self, a, b):
        self.a = a; self.b = b
    def __repr__(self): 
        return '%s + %s*gamma'%(self.a, self.b)
    def __add__(left, right):
        return Elt4(left.a + right.a, left.b + right.b)
    def __mul__(left, right): 
        a, b, c, d = left.a, left.b, right.a, right.b
        return Elt4(a*c + b*d, b*c + a*d + b*d)
///
}}}

{{{id=43|
z = Elt4(3, 5); w = Elt4(2, 8); z, w
///
(3 + 5*gamma, 2 + 8*gamma)
}}}

<p>Note that the type of an instance of a cdef'd class is "type" rather than "instance". &nbsp;</p>

{{{id=35|
type(z)
///
<type '_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage83_spyx_0.Elt4'>
}}}

{{{id=46|
z+w
///
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_77.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("eit3"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "/private/var/folders/7y/7y-O1iZOGTmMUMnLq7otq++++TI/-Tmp-/tmpBX248K/___code___.py", line 2, in <module>
    exec compile(u'z+w' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage83_spyx_0.pyx", line 13, in _Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage83_spyx_0.Elt4.__add__ (_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage83_spyx_0.c:578)
    return Elt4(left.a + right.a, left.b + right.b)
AttributeError: '_Users_wstein__sage_psage_notebook_sagenb_home_adm' object has no attribute 'a'
}}}

<p><strong>Dang! &nbsp; </strong>Thwarted again! &nbsp;This time the problem is:</p>
<ul>
<li>Attributes of cdef'd classes are private by default. &nbsp;To access them you have to explicitly declare the type of the object to Cython... or make the attribute public (usually a bad idea, since it can result in very slow code).&nbsp;</li>
</ul>
<p>Let's try again...</p>

{{{id=49|
%cython

cdef class Elt5:    
    cdef object a, b
    def __init__(self, a, b):
        self.a = a; self.b = b
    def __repr__(self): 
        return '%s + %s*gamma'%(self.a, self.b)
    def __add__(Elt5 left, Elt5 right):
        return Elt5(left.a + right.a, left.b + right.b)
    def __mul__(Elt5 left, Elt5 right): 
        a, b, c, d = left.a, left.b, right.a, right.b
        return Elt5(a*c + b*d, b*c + a*d + b*d)
///
}}}

{{{id=51|
z = Elt5(3, 5); w = Elt5(2, 8); z, w
///
(3 + 5*gamma, 2 + 8*gamma)
}}}

{{{id=52|
type(z)
///
<type '_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage90_spyx_0.Elt5'>
}}}

{{{id=53|
z + w
///
5 + 13*gamma
}}}

{{{id=54|
z * w
///
46 + 74*gamma
}}}

<p>Finally! &nbsp; And is it any faster?</p>

{{{id=55|
timeit('z+w', number=10^5)
///
100000 loops, best of 3: 247 ns per loop
}}}

{{{id=58|
timeit('z*w', number=10^5)
///
100000 loops, best of 3: 520 ns per loop
}}}

<p><strong>Daaaaamn! &nbsp; </strong>It's way faster. &nbsp;Excellent.</p>

{{{id=63|
z + 2 
///
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_4.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("eiArIDI="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "/private/var/folders/7y/7y-O1iZOGTmMUMnLq7otq++++TI/-Tmp-/tmpseviB1/___code___.py", line 3, in <module>
    exec compile(u'z + _sage_const_2 ' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage2_spyx_0.pyx", line 12, in _Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage2_spyx_0.Elt5.__add__ (_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage2_spyx_0.c:572)
    def __add__(Elt5 left, Elt5 right):
TypeError: Argument 'right' has incorrect type (expected _Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage2_spyx_0.Elt5, got sage.rings.integer.Integer)
}}}

<p>The above is as expected. &nbsp; &nbsp;But...</p>

{{{id=60|
z + None
///
------------------------------------------------------------
Unhandled SIGSEGV: A segmentation fault occurred in Sage.
This probably occurred because a *compiled* component
of Sage has a bug in it (typically accessing invalid memory)
or is not properly wrapped with _sig_on, _sig_off.
You might want to run Sage under gdb with 'sage -gdb' to debug this.
Sage will now terminate (sorry).
------------------------------------------------------------
}}}

<h3><strong>Python totally crashed. &nbsp;Ouch. &nbsp;Sorry to burst your bubble, but -- snap! &nbsp; &nbsp;What just happened??!?</strong></h3>
<ul>
<li>In Cython when you declare an input variable to have a cdef class type, e.g., <strong>Elt5 right, </strong>then Cython converts if possible the input (right) to the type Elt5. &nbsp;It either gives a type error or succeeds. &nbsp;The funny thing is that by default None can be converted to any type, but if you try to access any methods or attributes of the resulting object, you get a segfault. &nbsp; &nbsp; The trick is to either explicitly check for None in your code (ugly) or do what is illustrated below.</li>
</ul>

{{{id=59|
%cython

cdef class Elt6:    
    cdef object a, b
    def __init__(self, a, b):
        self.a = a; self.b = b
    def __repr__(self): 
        return '%s + %s*gamma'%(self.a, self.b)
    def __add__(Elt6 left, Elt6 right not None):
        return Elt6(left.a + right.a, left.b + right.b)
    def __mul__(Elt6 left, Elt6 right not None): 
        a, b, c, d = left.a, left.b, right.a, right.b
        return Elt6(a*c + b*d, b*c + a*d + b*d)
///
}}}

{{{id=57|
z = Elt6(3, 5); w = Elt6(2, 8); z, w
///
(3 + 5*gamma, 2 + 8*gamma)
}}}

{{{id=65|
z + w
///
5 + 13*gamma
}}}

{{{id=66|
z + None
///
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_8.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("eiArIE5vbmU="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "/private/var/folders/7y/7y-O1iZOGTmMUMnLq7otq++++TI/-Tmp-/tmppzmbSx/___code___.py", line 2, in <module>
    exec compile(u'z + None' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage6_spyx_0.pyx", line 12, in _Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage6_spyx_0.Elt6.__add__ (_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage6_spyx_0.c:572)
    def __add__(Elt6 left, Elt6 right not None):
TypeError: Argument 'right' has incorrect type (expected _Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage6_spyx_0.Elt6, got NoneType)
}}}

{{{id=80|
timeit('z+w', number=10^5)
///
100000 loops, best of 3: 244 ns per loop
}}}

{{{id=81|
timeit('z*w', number=10^5)
///
100000 loops, best of 3: 543 ns per loop
}}}

{{{id=70|

///
}}}

<h2>Using cdef to declare C data structures for speed purposes</h2>
<p>Suppose we want our golden ring to be superfast, but we know we will only be dealing with small elements. &nbsp;E.g., suppose we are a table of something up to some bound, and only small elements are involve. &nbsp;Instead of "cdef object", we can replace object by any C datatype. &nbsp;Let's start with long:</p>

{{{id=69|
%cython

cdef class Elt7:    
    cdef long a, b
    def __init__(self, a, b):
        self.a = a; self.b = b
    def __repr__(self): 
        return '%s + %s*gamma'%(self.a, self.b)
    def __add__(Elt7 left, Elt7 right not None):
        return Elt7(left.a + right.a, left.b + right.b)
    def __mul__(Elt7 left, Elt7 right not None): 
        a, b, c, d = left.a, left.b, right.a, right.b
        return Elt7(a*c + b*d, b*c + a*d + b*d)
///
}}}

{{{id=68|
z = Elt7(3, 5); w = Elt7(2, 8); z, w
///
(3 + 5*gamma, 2 + 8*gamma)
}}}

{{{id=67|
z + w
///
5 + 13*gamma
}}}

{{{id=73|
z * w
///
46 + 74*gamma
}}}

{{{id=75|
a = Elt7(993098220348,2903842)
a*a
///
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_34.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("YSA9IEVsdDcoOTkzMDk4MjIwMzQ4LDI5MDM4NDIpCmEqYQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "/private/var/folders/7y/7y-O1iZOGTmMUMnLq7otq++++TI/-Tmp-/tmpTvf0YW/___code___.py", line 4, in <module>
    exec compile(u'a*a' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage13_spyx_0.pyx", line 16, in _Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage13_spyx_0.Elt7.__mul__ (_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage13_spyx_0.c:713)
    return Elt7(a*c + b*d, b*c + a*d + b*d)
  File "_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage13_spyx_0.pyx", line 9, in _Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage13_spyx_0.Elt7.__init__ (_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage13_spyx_0.c:486)
    self.a = a; self.b = b
OverflowError: long int too large to convert to int
}}}

<p>How does the speed do?</p>

{{{id=74|
timeit('z+w', number=10^5)
///
100000 loops, best of 3: 143 ns per loop
}}}

{{{id=79|
timeit('z*w', number=10^5)
///
100000 loops, best of 3: 309 ns per loop
}}}

{{{id=83|

///
}}}

<p>Nice and fast.</p>

{{{id=82|

///
}}}

{{{id=78|
%cython

cdef class Elt8:    
    cdef long a, b
    def __init__(self, a, b):
        self.a = a; self.b = b
    def __repr__(self): 
        return '%s + %s*gamma'%(self.a, self.b)
    def __add__(Elt8 left, Elt8 right not None):
        cdef Elt8 z = PY_NEW(Elt8)
        z.a = left.a + right.a
        z.b = left.b + right.b
        return z
    def __mul__(Elt8 left, Elt8 right not None): 
        cdef long a=left.a, b=left.b, c=right.a, d=right.b
        cdef Elt8 z = PY_NEW(Elt8)
        z.a = a*c + b*d
        z.b =b*c + a*d + b*d
        return z
///
}}}

{{{id=87|
z = Elt8(3, 5); w = Elt8(2, 8); z, w
///
(3 + 5*gamma, 2 + 8*gamma)
}}}

{{{id=86|
timeit('z+w', number=10^5)
///
100000 loops, best of 3: 104 ns per loop
}}}

{{{id=85|
timeit('z*w', number=10^5)
///
100000 loops, best of 3: 106 ns per loop
}}}

<p>... and that is as good as you're going to get without creating an "object pool" and/or doing something with stack based allocation...</p>

{{{id=90|

///
}}}

{{{id=89|

///
}}}

{{{id=88|

///
}}}