Here is the step-by-step guide to adding a new PARI functions to Sage. We use the Frobenius form of a matrix as an example.
The heavy lifting for matrices over integers is implemented using the
PARI library. To compute the frobenius form in PARI the
matfrobenius
function is used.
There are two ways to interact with the PARI library from Sage: The gp interface uses the gp interpreter and the PARI interface uses direct calls to the PARI C functions and is the preferred way as it is much faster.
So we need to add a new method to the gen class which is the abstract
representation of all PARI library objects. That means that if we add a
method to this class every PARI object regardless if it is a number,
polynomial or matrix will have our new method. So you can do
pari(1).matfrobenius()
but you will receive a PariError in this case.
The gen class is defined in sage/libs/pari/gen.pyx where we add the method
matfrobenius
:
def matfrobenius(self, flag=0): """ matfrobenius(M,{flag}): Return the Frobenius form of the square matrix M. If flag is 1, return only the elementary divisors. If flag is 2, return a two-components vector [F,B] where F is the Frobenius form and B is the basis change so that M=B^-1*F*B. """ _sig_on return self.new_gen(matfrobenius(self.g, flag))
The _sig_on
statement is some magic to prevent SIGSEGVs from the pari C
library to crash the Sage interpreter by catching these signals. The
self.new_gen()
call constructs a new SAGE-python-gen object from a given
pari-C-gen where the pari-C-gen is stored as the SAGE-python-gen.g attribute.
The matfrobenius call is just a call to the pari C library function
"matfrobenius" with the appropriate parameters.
The information which function to call and how to call it can be retrieved
from the PARI users manual (note: Sage includes the development version of
PARI so check the development version's user manual). Looking for
matfrobenius
you can find:
"The library syntax is matfrobenius(M,flag)"
Please not that the pari C function may have a different name than gp function
(see e.g. mathnf
), so always check with the manual.
At this point we are done with the PARI interface and may add some more interfaces around it for convenience:
First we add a functional representation of the method to
sage/libs/pari/functional.py so we can do:
matfrobenius(<some.pari.gen>)
additionally to
<some.pari.gen>.matfrobenius()
.
def matfrobenius(self): return pari(self).matfrobenius()
Then we also add a frobenius(flag) method to the matrix_integer class where we
call the matfrobenius()
method on the _pari_()
object associated with the
matrix after doing some sanity checking. Then we need to convert the pari gen
to some meaningful Sage objects depending on the return value as described in
the PARI user's manual.
def frobenius(self,flag=0): """ Return the Frobenius form of this matrix. If flag is 1, return only the elementary divisors. If flag is 2, return a two-components vector [F,B] where F is the Frobenius form and B is the basis change so that M=B^-1*F*B. INPUT: flag -- 0,1 or 2 as described above ALGORITHM: uses pari's matfrobenius() EXAMPLE: sage: A = MatrixSpace(IntegerRing(), 3)(range(9)) sage: A.frobenius(0) [ 0 0 0] [ 1 0 18] [ 0 1 12] sage: A.frobenius(1) [x3 - 12*x2 - 18*x] sage: A.frobenius(2) ([ 0 0 0] [ 1 0 18] [ 0 1 12], [ -1 2 -1] [ 0 23/15 -14/15] [ 0 -2/15 1/15]) """ if self.nrows()!=self.ncols(): raise ArithmeticError, "frobenius matrix of non-square matrix not defined." v = self._pari_().matfrobenius(flag) if flag==0: return self.matrix_space()(v.python()) elif flag==1: r = polynomial_ring.PolynomialRing(self.base_ring()) #BUG: this should be handled in PolynomialRing not here return [eval(str(x).replace("^","**"),{},r.gens_dict()) for x in v.python_list()] elif flag==2: F = matrix_space.MatrixSpace(rational_field.RationalField(),self.nrows())(v[0].python()) B = matrix_space.MatrixSpace(rational_field.RationalField(),self.nrows())(v[1].python()) return F,B
Then we add some examples, wait for make test to complete without errors and commit.
See About this document... for information on suggesting changes.