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:
Figure out multiplication using Sage's commutative algebra (or you could just do it with pencil and paper...):
{{{id=7| R. = QQ[] S. = 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| /// }}}Before using Cython, let's make a first easy Python implementation.
{{{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) /// }}}Let's try out our new class:
{{{id=10| z = Elt(3, 5); w = Elt(2, 8); z, w /// (3 + 5*gamma, 2 + 8*gamma) }}} {{{id=42| type(z) ///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").
{{{id=18| z + 1 /// Traceback (most recent call last): File "What about speed? It's not so bad, really:
{{{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| /// }}}Everything else we will do is about different versions of the above fairly simple class.
{{{id=27| /// }}}The only change is that I changed the name of the class and put %cython at the top of the block.
{{{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) /// }}}Let's try it out:
{{{id=25| z = Elt2(3, 5); w = Elt2(2, 8); z, w /// (3 + 5*gamma, 2 + 8*gamma) }}} {{{id=41| type(z) ///The timings are identical to the Python version.
{{{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| /// }}}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.
{{{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 "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:
Note that the type of an instance of a cdef'd class is "type" rather than "instance".
{{{id=35| type(z) ///Dang! Thwarted again! This time the problem is:
Let's try again...
{{{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) ///Finally! And is it any faster?
{{{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 }}}Daaaaamn! It's way faster. Excellent.
{{{id=63| z + 2 /// Traceback (most recent call last): File "The above is as expected. But...
{{{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). ------------------------------------------------------------ }}}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:
{{{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 "How does the speed do?
{{{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| /// }}}Nice and fast.
{{{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 }}}... and that is as good as you're going to get without creating an "object pool" and/or doing something with stack based allocation...
{{{id=90| /// }}} {{{id=89| /// }}} {{{id=88| /// }}}