Using External Libraries and Interfaces

When writing code for Sage, use Python for the basic structure and interface. For speed, efficiency, or convenience, you can implement parts of the code using any of the following languages: Cython, C/C++, Fortran 95, GAP, Common Lisp, Singular, and PARI/GP. You can also use all C/C++ libraries included with Sage, see Packages and Features. And if you are okay with your code depending on external programs, you can use Octave, or even Magma, Mathematica, or Maple.

In this chapter, we discuss interfaces between Sage and PARI, GAP and Singular.

The PARI C library interface

Here is a step-by-step guide to adding 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. 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.

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/src/sage/libs/cypari2/gen.pyx, and this is where we add the method matfrobenius:

def matfrobenius(self, flag=0):
    r"""
    M.matfrobenius(flag=0): Return the Frobenius form of the square
    matrix M. If flag is 1, return only the elementary divisors (a list
    of polynomials). 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`.

    EXAMPLES::

        sage: a = pari('[1,2;3,4]')
        sage: a.matfrobenius()
        [0, 2; 1, 5]
        sage: a.matfrobenius(flag=1)
        [x^2 - 5*x - 2]
        sage: a.matfrobenius(2)
        [[0, 2; 1, 5], [1, -1/3; 0, 1/3]]
    """
    sig_on()
    return self.new_gen(matfrobenius(self.g, flag, 0))

Note the use of the sig_on() statement.

The matfrobenius call is just a call to the PARI C library function matfrobenius with the appropriate parameters.

The self.new_gen(GEN x) call constructs a new Sage gen object from a given PARI GEN where the PARI GEN is stored as the .g attribute. Apart from this, self.new_gen() calls a closing sig_off() macro and also clears the PARI stack so it is very convenient to use in a return statement as illustrated above. So after self.new_gen(), all PARI GEN’s which are not converted to Sage gen’s are gone. There is also self.new_gen_noclear(GEN x) which does the same as self.new_gen(GEN x) except that it does not call sig_off() nor clear the PARI stack.

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 GEN matfrobenius(GEN M, long flag, long v = -1), where v is a variable number.

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

We can also add a frobenius_form(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_form(self, flag=0, var='x'):
    """
    Return the Frobenius form (rational canonical form) of this matrix.

    INPUT:

    -  ``flag`` -- 0 (default), 1 or 2 as follows:

        -  ``0`` -- (default) return the Frobenius form of this
           matrix.

        -  ``1`` -- return only the elementary divisor
           polynomials, as polynomials in var.

        -  ``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}FB`.

    -  ``var`` -- a string (default: 'x')

    ALGORITHM: uses PARI's :pari:`matfrobenius`

    EXAMPLES::

        sage: A = MatrixSpace(ZZ, 3)(range(9))
        sage: A.frobenius_form(0)
        [ 0  0  0]
        [ 1  0 18]
        [ 0  1 12]
        sage: A.frobenius_form(1)
        [x^3 - 12*x^2 - 18*x]
        sage: A.frobenius_form(1, var='y')
        [y^3 - 12*y^2 - 18*y]
    """
    if not self.is_square():
        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 = PolynomialRing(self.base_ring(), names=var)
        retr = []
        for f in v:
            retr.append(eval(str(f).replace("^","**"), {'x':r.gen()}, r.gens_dict()))
        return retr
    elif flag == 2:
        F = matrix_space.MatrixSpace(QQ, self.nrows())(v[0].python())
        B = matrix_space.MatrixSpace(QQ, self.nrows())(v[1].python())
        return F, B

GAP

Wrapping a GAP function in Sage can be done in two different ways. The first one uses the pexpect interface to pipe various commands to GAP and read back the input into Sage. The second way is to access GAP through its library interface libgap. We recommend the second one for the sake of efficiency.

For example, suppose we want to make a wrapper for the computation of the Cartan matrix of a simple Lie algebra. The Cartan matrix of \(G_2\) is available in GAP using the commands:

gap> L:= SimpleLieAlgebra( "G", 2, Rationals );
<Lie algebra of dimension 14 over Rationals>
gap> R:= RootSystem( L );
<root system of rank 2>
gap> CartanMatrix( R );

In Sage, one can access these commands by typing:

Sage

sage: L = libgap.SimpleLieAlgebra("G", 2, QQ); L
<Lie algebra of dimension 14 over Rationals>
sage: R = L.RootSystem(); R
<root system of rank 2>
sage: R.CartanMatrix()
[ [ 2, -1 ], [ -3, 2 ] ]

Python

>>> from sage.all import *
>>> L = libgap.SimpleLieAlgebra("G", Integer(2), QQ); L
<Lie algebra of dimension 14 over Rationals>
>>> R = L.RootSystem(); R
<root system of rank 2>
>>> R.CartanMatrix()
[ [ 2, -1 ], [ -3, 2 ] ]

The purpose of this section is to use this example to show how one might write a Python/Sage program whose input is, say, ('G',2) and whose output is the matrix above (but as a Sage Matrix—see the code in the directory SAGE_ROOT/src/sage/matrix/ and the corresponding parts of the Sage reference manual).

First, the input must be converted into a libgap object, either by applying libgap.eval on a string, or by using a GAP command CMD as a libgap.CMD method. Then one can work with these object using other GAP commands as libgap methods. At the end, one can convert back to sage using the method sage if it works.

def cartan_matrix(typ, rank):
    """
    Return the Cartan matrix of given Chevalley type and rank.

    INPUT:

    - type -- a Chevalley letter name, as a string, for
      a family type of simple Lie algebras
    - rank -- an integer (legal for that type)

    EXAMPLES::

        sage: cartan_matrix("A",5)
        [ 2 -1  0  0  0]
        [-1  2 -1  0  0]
        [ 0 -1  2 -1  0]
        [ 0  0 -1  2 -1]
        [ 0  0  0 -1  2]
        sage: cartan_matrix("G",2)
        [ 2 -1]
        [-3  2]
    """
    L = libgap.SimpleLieAlgebra(typ, rank, libgap.Rationals)
    R = L.RootSystem()
    sM = R.CartanMatrix()
    ans = sM.sage()
    MS = MatrixSpace(QQ, rank)
    return MS(ans)

The output ans is a Python list. The last two lines convert that list to an instance of the Sage class Matrix.

Defining “easy” and “hard” is subjective, but here is one definition. Wrapping a GAP function is “easy” if there is already a corresponding class in Python or Sage for the output data type of the GAP function you are trying to wrap. For example, wrapping any GUAVA (GAP’s error-correcting codes package) function is “easy” since error-correcting codes are vector spaces over finite fields and GUAVA functions return one of the following data types:

• vectors over finite fields,

• polynomials over finite fields,

• matrices over finite fields,

• permutation groups or their elements,

• integers.

Sage already has classes for each of these.

A “hard” example is left as an exercise! Here are a few ideas.

• Write a wrapper for GAP’s FreeLieAlgebra function (or, more generally, all the finitely presented Lie algebra functions in GAP). This would require creating new Python objects.

• Write a wrapper for GAP’s FreeGroup function (or, more generally, all the finitely presented groups functions in GAP). This would require writing some new Python objects.

• Write a wrapper for GAP’s character tables. Though this could be done without creating new Python objects, to make the most use of these tables, it probably would be best to have new Python objects for this.

LibGAP

The disadvantage of using other programs through interfaces is that there is a certain unavoidable latency (of the order of 10ms) involved in sending input and receiving the result. If you have to call functions in a tight loop this can be unacceptably slow. Calling into a shared library has much lower latency and furthermore avoids having to convert everything into a string in-between. This is why Sage includes a shared library version of the GAP kernel, available as \(libgap\) in Sage. The libgap analogue of the first example in GAP is:

Sage

sage: SimpleLieAlgebra = libgap.function_factory('SimpleLieAlgebra')
sage: L = SimpleLieAlgebra('G', 2, QQ)
sage: R = L.RootSystem();  R
<root system of rank 2>
sage: R.CartanMatrix()    # output is a GAP matrix
[ [ 2, -1 ], [ -3, 2 ] ]
sage: matrix(R.CartanMatrix())   # convert to Sage matrix
[ 2 -1]
[-3  2]

Python

>>> from sage.all import *
>>> SimpleLieAlgebra = libgap.function_factory('SimpleLieAlgebra')
>>> L = SimpleLieAlgebra('G', Integer(2), QQ)
>>> R = L.RootSystem();  R
<root system of rank 2>
>>> R.CartanMatrix()    # output is a GAP matrix
[ [ 2, -1 ], [ -3, 2 ] ]
>>> matrix(R.CartanMatrix())   # convert to Sage matrix
[ 2 -1]
[-3  2]

Singular

Using Singular functions from Sage is not much different conceptually from using GAP functions from Sage. As with GAP, this can range from easy to hard, depending on how much of the data structure of the output of the Singular function is already present in Sage.

First, some terminology. For us, a curve \(X\) over a finite field \(F\) is an equation of the form \(f(x,y) = 0\), where \(f \in F[x,y]\) is a polynomial. It may or may not be singular. A place of degree \(d\) is a Galois orbit of \(d\) points in \(X(E)\), where \(E/F\) is of degree \(d\). For example, a place of degree \(1\) is also a place of degree \(3\), but a place of degree \(2\) is not since no degree \(3\) extension of \(F\) contains a degree \(2\) extension. Places of degree \(1\) are also called \(F\)-rational points.

As an example of the Sage/Singular interface, we will explain how to wrap Singular’s NSplaces, which computes places on a curve over a finite field. (The command closed_points also does this in some cases.) This is “easy” since no new Python classes are needed in Sage to carry this out.

Here is an example on how to use this command in Singular:

 A Computer Algebra System for Polynomial Computations   /   version 3-0-0
                                                       0<
     by: G.-M. Greuel, G. Pfister, H. Schoenemann        \   May 2005
FB Mathematik der Universitaet, D-67653 Kaiserslautern    \
> LIB "brnoeth.lib";
[...]
> ring s=5,(x,y),lp;
> poly f=y^2-x^9-x;
> list X1=Adj_div(f);
Computing affine singular points ...
Computing all points at infinity ...
Computing affine singular places ...
Computing singular places at infinity ...
Computing non-singular places at infinity ...
Adjunction divisor computed successfully

The genus of the curve is 4
> list X2=NSplaces(1,X1);
Computing non-singular affine places of degree 1 ...
> list X3=extcurve(1,X2);

Total number of rational places : 6

> def R=X3[1][5];
> setring R;
> POINTS;
[1]:
   [1]:
      0
   [2]:
      1
   [3]:
      0
[2]:
   [1]:
      -2
   [2]:
      1
   [3]:
      1
[3]:
   [1]:
      -2
   [2]:
      1
   [3]:
      1
[4]:
   [1]:
      -2
   [2]:
      -1
   [3]:
      1
[5]:
   [1]:
      2
   [2]:
      -2
   [3]:
      1
[6]:
   [1]:
      0
   [2]:
      0
   [3]:
      1

Here is another way of doing this same calculation in the Sage interface to Singular:

Sage

sage: singular.LIB("brnoeth.lib")
sage: singular.ring(5,'(x,y)','lp')
    polynomial ring, over a field, global ordering
    //   coefficients: ZZ/5
    //   number of vars : 2
    //        block   1 : ordering lp
    //                  : names    x y
    //        block   2 : ordering C
sage: f = singular('y^2-x^9-x')
sage: print(singular.eval("list X1=Adj_div(%s);"%f.name()))
Computing affine singular points ...
Computing all points at infinity ...
Computing affine singular places ...
Computing singular places at infinity ...
Computing non-singular places at infinity ...
Adjunction divisor computed successfully

The genus of the curve is 4
sage: print(singular.eval("list X2=NSplaces(1,X1);"))
Computing non-singular affine places of degree 1 ...
sage: print(singular.eval("list X3=extcurve(1,X2);"))

Total number of rational places : 6

sage: singular.eval("def R=X3[1][5];")
''
sage: singular.eval("setring R;")
''
sage: L = singular.eval("POINTS;")

sage: print(L) # random
[1]:
   [1]:
      0
   [2]:
      1
   [3]:
      0
...

Python

>>> from sage.all import *
>>> singular.LIB("brnoeth.lib")
>>> singular.ring(Integer(5),'(x,y)','lp')
    polynomial ring, over a field, global ordering
    //   coefficients: ZZ/5
    //   number of vars : 2
    //        block   1 : ordering lp
    //                  : names    x y
    //        block   2 : ordering C
>>> f = singular('y^2-x^9-x')
>>> print(singular.eval("list X1=Adj_div(%s);"%f.name()))
Computing affine singular points ...
Computing all points at infinity ...
Computing affine singular places ...
Computing singular places at infinity ...
Computing non-singular places at infinity ...
Adjunction divisor computed successfully
<BLANKLINE>
The genus of the curve is 4
>>> print(singular.eval("list X2=NSplaces(1,X1);"))
Computing non-singular affine places of degree 1 ...
>>> print(singular.eval("list X3=extcurve(1,X2);"))
<BLANKLINE>
Total number of rational places : 6
<BLANKLINE>
>>> singular.eval("def R=X3[1][5];")
''
>>> singular.eval("setring R;")
''
>>> L = singular.eval("POINTS;")

>>> print(L) # random
[1]:
   [1]:
      0
   [2]:
      1
   [3]:
      0
...

From looking at the output, notice that our wrapper function will need to parse the string represented by \(L\) above, so let us write a separate function to do just that. This requires figuring out how to determine where the coordinates of the points are placed in the string \(L\). Python has some very useful string manipulation commands to do just that.

def points_parser(string_points, F):
    """
    This function will parse a string of points
    of X over a finite field F returned by Singular's NSplaces
    command into a Python list of points with entries from F.

    EXAMPLES::

        sage: F = GF(5)
        sage: points_parser(L,F)
        ((0, 1, 0), (3, 4, 1), (0, 0, 1), (2, 3, 1), (3, 1, 1), (2, 2, 1))
    """
    Pts = []
    n = len(L)
    # start block to compute a pt
    L1 = L
    while len(L1) > 32:
        idx =L1.index("     ")
        pt = []
        # start block1 for compute pt
        idx = L1.index("     ")
        idx2 = L1[idx:].index("\n")
        L2 = L1[idx:idx+idx2]
        pt.append(F(eval(L2)))
        # end block1 to compute pt
        L1 = L1[idx+8:] # repeat block 2 more times
        # start block2 for compute pt
        idx = L1.index("     ")
        idx2 = L1[idx:].index("\n")
        L2 = L1[idx:idx+idx2]
        pt.append(F(eval(L2)))
        # end block2 to compute pt
        L1=L1[idx+8:] # repeat block 1 more time
        # start block3 for compute pt
        idx=L1.index("     ")
        if "\n" in L1[idx:]:
            idx2 = L1[idx:].index("\n")
        else:
            idx2 = len(L1[idx:])
        L2 = L1[idx:idx+idx2]
        pt.append(F(eval(L2)))
        # end block3 to compute pt
        # end block to compute a pt
        Pts.append(tuple(pt))  # repeat until no more pts
        L1 = L1[idx+8:] # repeat block 2 more times
    return tuple(Pts)

Now it is an easy matter to put these ingredients together into a Sage function which takes as input a triple \((f,F,d)\): a polynomial \(f\) in \(F[x,y]\) defining \(X:\ f(x,y)=0\) (note that the variables \(x,y\) must be used), a finite field \(F\) of prime order, and the degree \(d\). The output is the number of places in \(X\) of degree \(d=1\) over \(F\). At the moment, there is no “translation” between elements of \(GF(p^d)\) in Singular and Sage unless \(d=1\). So, for this reason, we restrict ourselves to points of degree one.

def places_on_curve(f, F):
    """
    INPUT:

    - f -- element of F[x,y], defining X: f(x,y)=0
    - F -- a finite field of *prime order*

    OUTPUT:

    integer -- the number of places in X of degree d=1 over F

    EXAMPLES::

        sage: F = GF(5)
        sage: R = PolynomialRing(F,2,names=["x","y"])
        sage: x,y = R.gens()
        sage: f = y^2-x^9-x
        sage: places_on_curve(f,F)
        ((0, 1, 0), (3, 4, 1), (0, 0, 1), (2, 3, 1), (3, 1, 1), (2, 2, 1))
    """
    d = 1
    p = F.characteristic()
    singular.eval('LIB "brnoeth.lib";')
    singular.eval("ring s="+str(p)+",(x,y),lp;")
    singular.eval("poly f="+str(f))
    singular.eval("list X1=Adj_div(f);")
    singular.eval("list X2=NSplaces("+str(d)+",X1);")
    singular.eval("list X3=extcurve("+str(d)+",X2);")
    singular.eval("def R=X3[1][5];")
    singular.eval("setring R;")
    L = singular.eval("POINTS;")
    return points_parser(L,F)

Note that the ordering returned by this Sage function is exactly the same as the ordering in the Singular variable POINTS.

One more example (in addition to the one in the docstring):

Sage

sage: F = GF(2)
sage: R = MPolynomialRing(F,2,names = ["x","y"])
sage: x,y = R.gens()
sage: f = x^3*y+y^3+x
sage: places_on_curve(f,F)
((0, 1, 0), (1, 0, 0), (0, 0, 1))

Python

>>> from sage.all import *
>>> F = GF(Integer(2))
>>> R = MPolynomialRing(F,Integer(2),names = ["x","y"])
>>> x,y = R.gens()
>>> f = x**Integer(3)*y+y**Integer(3)+x
>>> places_on_curve(f,F)
((0, 1, 0), (1, 0, 0), (0, 0, 1))

Singular: another approach

There is also a more Python-like interface to Singular. Using this, the code is much simpler, as illustrated below. First, we demonstrate computing the places on a curve in a particular case:

Sage

sage: singular.lib('brnoeth.lib')
sage: R = singular.ring(5, '(x,y)', 'lp')
sage: f = singular.new('y^2 - x^9 - x')
sage: X1 = f.Adj_div()
sage: X2 = singular.NSplaces(1, X1)
sage: X3 = singular.extcurve(1, X2)
sage: R = X3[1][5]
sage: singular.set_ring(R)
sage: L = singular.new('POINTS')

Python

>>> from sage.all import *
>>> singular.lib('brnoeth.lib')
>>> R = singular.ring(Integer(5), '(x,y)', 'lp')
>>> f = singular.new('y^2 - x^9 - x')
>>> X1 = f.Adj_div()
>>> X2 = singular.NSplaces(Integer(1), X1)
>>> X3 = singular.extcurve(Integer(1), X2)
>>> R = X3[Integer(1)][Integer(5)]
>>> singular.set_ring(R)
>>> L = singular.new('POINTS')

Note that these elements of L are defined modulo 5 in Singular, and they compare differently than you would expect from their print representation:

Sage

sage: sorted([(L[i][1], L[i][2], L[i][3]) for i in range(1,7)])
[(0, 0, 1), (0, 1, 0), (2, 2, 1), (2, -2, 1), (-2, 1, 1), (-2, -1, 1)]

Python

>>> from sage.all import *
>>> sorted([(L[i][Integer(1)], L[i][Integer(2)], L[i][Integer(3)]) for i in range(Integer(1),Integer(7))])
[(0, 0, 1), (0, 1, 0), (2, 2, 1), (2, -2, 1), (-2, 1, 1), (-2, -1, 1)]

Next, we implement the general function (for brevity we omit the docstring, which is the same as above). Note that the point_parser function is not required:

def places_on_curve(f, F):
    p = F.characteristic()
    if F.degree() > 1:
        raise NotImplementedError
    singular.lib('brnoeth.lib')
    R = singular.ring(5, '(x,y)', 'lp')
    f = singular.new('y^2 - x^9 - x')
    X1 = f.Adj_div()
    X2 = singular.NSplaces(1, X1)
    X3 = singular.extcurve(1, X2)
    R = X3[1][5]
    singular.setring(R)
    L = singular.new('POINTS')
    return [(int(L[i][1]), int(L[i][2]), int(L[i][3])) \
             for i in range(1,int(L.size())+1)]

This code is much shorter, nice, and more readable. However, it depends on certain functions, e.g. singular.setring having been implemented in the Sage/Singular interface, whereas the code in the previous section used only the barest minimum of that interface.

Creating a new pseudo-TTY interface

You can create Sage pseudo-tty interfaces that allow Sage to work with almost any command line program, and which do not require any modification or extensions to that program. They are also surprisingly fast and flexible (given how they work!), because all I/O is buffered, and because interaction between Sage and the command line program can be non-blocking (asynchronous). A pseudo-tty Sage interface is asynchronous because it derives from the Sage class Expect, which handles the communication between Sage and the external process.

For example, here is part of the file SAGE_ROOT/src/sage/interfaces/octave.py, which defines an interface between Sage and Octave, an open source program for doing numerical computations, among other things:

import os
from expect import Expect, ExpectElement

class Octave(Expect):
    ...

The first two lines import the library os, which contains operating system routines, and also the class Expect, which is the basic class for interfaces. The third line defines the class Octave; it derives from Expect as well. After this comes a docstring, which we omit here (see the file for details). Next comes:

def __init__(self, script_subdirectory="", logfile=None,
             server=None, server_tmpdir=None):
    Expect.__init__(self,
                    name = 'octave',
                    prompt = '>',
                    command = "octave --no-line-editing --silent",
                    server = server,
                    server_tmpdir = server_tmpdir,
                    script_subdirectory = script_subdirectory,
                    restart_on_ctrlc = False,
                    verbose_start = False,
                    logfile = logfile,
                    eval_using_file_cutoff=100)

This uses the class Expect to set up the Octave interface:

def set(self, var, value):
    """
    Set the variable var to the given value.
    """
    cmd = '%s=%s;' % (var,value)
    out = self.eval(cmd)
    if out.find("error") != -1:
        raise TypeError("Error executing code in Octave\nCODE:\n\t%s\nOctave ERROR:\n\t%s"%(cmd, out))

def get(self, var):
    """
    Get the value of the variable var.
    """
    s = self.eval('%s' % var)
    i = s.find('=')
    return s[i+1:]

def console(self):
    octave_console()

These let users type octave.set('x', 3), after which octave.get('x') returns ' 3'. Running octave.console() dumps the user into an Octave interactive shell:

def solve_linear_system(self, A, b):
    """
    Use octave to compute a solution x to A*x = b, as a list.

    INPUT:

    - A -- mxn matrix A with entries in QQ or RR
    - b -- m-vector b entries in QQ or RR (resp)

    OUTPUT:

    A list x (if it exists) which solves M*x = b

    EXAMPLES::

        sage: M33 = MatrixSpace(QQ,3,3)
        sage: A   = M33([1,2,3,4,5,6,7,8,0])
        sage: V3  = VectorSpace(QQ,3)
        sage: b   = V3([1,2,3])
        sage: octave.solve_linear_system(A,b)    # optional - octave
        [-0.333333, 0.666667, 0]

    AUTHOR: David Joyner and William Stein
    """
    m = A.nrows()
    n = A.ncols()
    if m != len(b):
        raise ValueError("dimensions of A and b must be compatible")
    from sage.matrix.all import MatrixSpace
    from sage.rings.all import QQ
    MS = MatrixSpace(QQ, m, 1)
    b  = MS(list(b))  # converted b to a "column vector"
    sA = self.sage2octave_matrix_string(A)
    sb = self.sage2octave_matrix_string(b)
    self.eval("a = " + sA )
    self.eval("b = " + sb )
    soln = octave.eval("c = a \\ b")
    soln = soln.replace("\n\n ", "[")
    soln = soln.replace("\n\n", "]")
    soln = soln.replace("\n", ",")
    sol  = soln[3:]
    return eval(sol)

This code defines the method solve_linear_system, which works as documented.

These are only excerpts from octave.py; check that file for more definitions and examples. Look at other files in the directory SAGE_ROOT/src/sage/interfaces/ for examples of interfaces to other software packages.