4.3 The PARI C-library Interface

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.

Some 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 -- this is the preferred way as it is much faster. Thus this section focuses on using PARI.

So we will add a new method to the gen class: this is the abstract representation of all PARI library objects. That means that once we add a method to this class, every PARI object, whether it is a number, polynomial or matrix, will have our new method. So you can do pari(1).matfrobenius(), but since PARI wants to apply matfrobenius to matrices, not numbers, you will receive a PariError in this case.

The gen class is defined in SAGE_ROOT/devel/sage/sage/libs/pari/gen.pyx, and this is 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. Note that self.new_gen() calls a closing _sig_off macro. These two must always come in pairs, i.e. every _sig_on must be matched by a closing _sig_off. 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 about which function to call and how to call it can be retrieved from the PARI user's manual (note: Sage includes the development version of PARI, so check that version of the user's manual). Looking for matfrobenius you can find: "The library syntax is matfrobenius(M,flag)".

In case you are familiar with gp, please note that the PARI C function may have a different name than the corresponding gp function (for example, see mathnf), so always check the manual.

At this point we are done with the PARI interface and can add some more interfaces around it for convenience: first we add a functional representation of the method to SAGE_ROOT/devel/sage/sage/libs/pari/functional.py so we can do matfrobenius(<some.pari.gen>) as well as <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 to the matrix after doing some sanity checking. Then we convert output from PARI to Sage objects:

    def frobenius(self,flag=0):
        """
	If flag is 0 (the default value), return the Frobenius 
            form of this matrix.
        If flag is 1, return only the elementary divisors.
        If flag is 2, return a two-component 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