6.1 The PARI C-library Interface

(This chapter was written by Martin Albrecht.)

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.