581d-2010-11-10-sqrt5

561 days ago by wstein

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

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. 

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.   Your homework (number 7) will be to generalize these examples to the field K=\QQ(\gamma). 

Math:

  • The elements of R are all numbers of the form a+b\gamma with a,b\in\ZZ.  
  • Addition: (a+b\gamma) + (c+d\gamma) = (a+c) + (b+d)\gamma.
  • Multiplication?  We have \gamma^2 = \gamma+1, so
    (a+b\gamma) \cdot (c+d\gamma) = (ac + bd) + (ad+bc+bd)\gamma.

Figure out multiplication using Sage's commutative algebra (or you could just do it with pencil and paper...):
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
b*c*gamma + a*d*gamma + b*d*gamma + a*c + b*d
 
       

Python Class

Before using Cython, let's make a first easy Python implementation.
# 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) 
       

Let's try out our new class:
z = Elt(3, 5); w = Elt(2, 8); z, w 
        
(3 + 5*gamma, 2 + 8*gamma)
(3 + 5*gamma, 2 + 8*gamma)
type(z) 
        
<type 'instance'>
<type 'instance'>
z + w 
        
5 + 13*gamma
5 + 13*gamma
z * w 
        
46 + 74*gamma
46 + 74*gamma

Awesome!   Back make sure you understand that the below should fail (and why).  We will officially not worry about such things until probably next week (when we learn about "the coercion model").
z + 1 
        
Traceback (click to the left of this block for traceback)
...
AttributeError: 'sage.rings.integer.Integer' object has no attribute 'a'
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"))
  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 
  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'

What about speed?  It's not so bad, really:
timeit('z+w', number=10^5) 
        
100000 loops, best of 3: 1.89 µs per loop
100000 loops, best of 3: 1.89 µs per loop
timeit('z*w', number=10^5) 
        
100000 loops, best of 3: 2.47 µs per loop
100000 loops, best of 3: 2.47 µs per loop
 
       

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

Next, our First Naive Cython Version

The only change is that I changed the name of the class and put %cython at the top of the block. 
%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) 
       

Let's try it out:
z = Elt2(3, 5); w = Elt2(2, 8); z, w 
        
(3 + 5*gamma, 2 + 8*gamma)
(3 + 5*gamma, 2 + 8*gamma)
type(z) 
        
<type 'instance'>
<type 'instance'>
z + w 
        
5 + 13*gamma
5 + 13*gamma
z*w 
        
46 + 74*gamma
46 + 74*gamma

The timings are identical to the Python version.
timeit('z+w', number=10^5) 
        
100000 loops, best of 3: 1.85 µs per loop
100000 loops, best of 3: 1.85 µs per loop
timeit('z*w', number=10^5) 
        
100000 loops, best of 3: 2.29 µs per loop
100000 loops, best of 3: 2.29 µs per loop
 
       

A Less Naive Cython Version: Make a Cdef'd Extension Class

This looks for now just like what we did above, but we put "cdef class" instead of just "class".  There are some major changes behind the scenes though.
%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) 
       
z = Elt3(3, 5); w = Elt3(2, 8); z, w 
        
Traceback (click to the left of this block for traceback)
...
AttributeError:
'_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage80_spy\
x_0.Elt3' object has no attribute 'a'
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"))
  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
  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'

Dang:  The first difference is that it doesn't work at all!  

cdef'd classes are different than normal Python classes in a few ways:

  1. All their attributes must be declared at compile time.  The plus is you can declare C data types (like int, long, etc.), which is impossible with Python classes.
  2. Only single inheritance is supported.
  3. 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.   (Calling a Python function usually involves creating a tuple and/or dict, putting the inputs in the tuple, unpacking the tuple, etc.  Lots of overhead and red tape!)
%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) 
       
z = Elt4(3, 5); w = Elt4(2, 8); z, w 
        
(3 + 5*gamma, 2 + 8*gamma)
(3 + 5*gamma, 2 + 8*gamma)

Note that the type of an instance of a cdef'd class is "type" rather than "instance".  
type(z) 
        
<type
'_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage83_spy\
x_0.Elt4'>
<type '_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage83_spyx_0.Elt4'>
z+w 
        
Traceback (click to the left of this block for traceback)
...
AttributeError: '_Users_wstein__sage_psage_notebook_sagenb_home_adm'
object has no attribute 'a'
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"))
  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
  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'

Dang!   Thwarted again!  This time the problem is:

  • Attributes of cdef'd classes are private by default.  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). 

Let's try again...
%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) 
       
z = Elt5(3, 5); w = Elt5(2, 8); z, w 
        
(3 + 5*gamma, 2 + 8*gamma)
(3 + 5*gamma, 2 + 8*gamma)
type(z) 
        
<type
'_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage90_spy\
x_0.Elt5'>
<type '_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage90_spyx_0.Elt5'>
z + w 
        
5 + 13*gamma
5 + 13*gamma
z * w 
        
46 + 74*gamma
46 + 74*gamma

Finally!   And is it any faster?
timeit('z+w', number=10^5) 
        
100000 loops, best of 3: 247 ns per loop
100000 loops, best of 3: 247 ns per loop
timeit('z*w', number=10^5) 
        
100000 loops, best of 3: 520 ns per loop
100000 loops, best of 3: 520 ns per loop

Daaaaamn!   It's way faster.  Excellent.
z + 2 
        
Traceback (click to the left of this block for traceback)
...
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)
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"))
  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 
  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)

The above is as expected.    But...
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).
------------------------------------------------------------
------------------------------------------------------------
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).
------------------------------------------------------------

Python totally crashed.  Ouch.  Sorry to burst your bubble, but -- snap!    What just happened??!?

  • In Cython when you declare an input variable to have a cdef class type, e.g., Elt5 right, then Cython converts if possible the input (right) to the type Elt5.  It either gives a type error or succeeds.  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.     The trick is to either explicitly check for None in your code (ugly) or do what is illustrated below.
%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) 
       
z = Elt6(3, 5); w = Elt6(2, 8); z, w 
        
(3 + 5*gamma, 2 + 8*gamma)
(3 + 5*gamma, 2 + 8*gamma)
z + w 
        
5 + 13*gamma
5 + 13*gamma
z + None 
        
Traceback (click to the left of this block for traceback)
...
TypeError: Argument 'right' has incorrect type (expected
_Users_wstein__sage_psage_notebook_sagenb_home_admin_22_code_sage6_spyx_\
0.Elt6, got NoneType)
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"))
  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
  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)
timeit('z+w', number=10^5) 
        
100000 loops, best of 3: 244 ns per loop
100000 loops, best of 3: 244 ns per loop
timeit('z*w', number=10^5) 
        
100000 loops, best of 3: 543 ns per loop
100000 loops, best of 3: 543 ns per loop
 
       

Using cdef to declare C data structures for speed purposes

Suppose we want our golden ring to be superfast, but we know we will only be dealing with small elements.  E.g., suppose we are a table of something up to some bound, and only small elements are involve.  Instead of "cdef object", we can replace object by any C datatype.  Let's start with long:
%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) 
       
z = Elt7(3, 5); w = Elt7(2, 8); z, w 
        
(3 + 5*gamma, 2 + 8*gamma)
(3 + 5*gamma, 2 + 8*gamma)
z + w 
        
5 + 13*gamma
5 + 13*gamma
z * w 
        
46 + 74*gamma
46 + 74*gamma
a = Elt7(993098220348,2903842) a*a 
        
Traceback (click to the left of this block for traceback)
...
OverflowError: long int too large to convert to int
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"))
  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
  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

How does the speed do?
timeit('z+w', number=10^5) 
        
100000 loops, best of 3: 143 ns per loop
100000 loops, best of 3: 143 ns per loop
timeit('z*w', number=10^5) 
        
100000 loops, best of 3: 309 ns per loop
100000 loops, best of 3: 309 ns per loop
 
       

Nice and fast.
 
       
%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 
       
z = Elt8(3, 5); w = Elt8(2, 8); z, w 
        
(3 + 5*gamma, 2 + 8*gamma)
(3 + 5*gamma, 2 + 8*gamma)
timeit('z+w', number=10^5) 
        
100000 loops, best of 3: 104 ns per loop
100000 loops, best of 3: 104 ns per loop
timeit('z*w', number=10^5) 
        
100000 loops, best of 3: 106 ns per loop
100000 loops, best of 3: 106 ns per loop

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