MathLink and PARI, part II

With advice from David Withoff at Wolfram and Bill Allombert, I've figured out how to use Mathematica's MathLink capability to install PARI functions in a Mathematica session. This is a summary with some examples which you can modify.

MathLink includes a lot of low level functions for communication back and forth, but this can be avoided by the use of MathLink 'templates'. The template file consists of C functions which in this case are mainly calls to the PARI library functions, and the Mathematica defined functions which will correspond to them. The mprep application takes the template file and generates C code to set up the MathLink and handle the function calls.

The template file can also include Mathematica expressions preceded by :Evaluate: which mprep passes to the C code as character string. The C code passes them to Mathematica when you Install, and and tells Mathematica to evaluate them. This is useful to wrap Mathematica code around your installed function to check for errors, and to convert complicated Mathematica structures into things C will appreciate, such as reals or integers.

Each function thus typically has four incarnations. For example, below is a Mathematica function QfbPow[f_Qfb, n_Integer] which does some error checking and converts an f of the form Qfb[a,b,c] into a list of integers {a,b,c} which is passed to the Installed function qfbp[a,b,c,n]. The MathLink connection matches this up to the C function qfbnupow(a,b,c,n) which creates PARI objects f = Qfb(a,b,c) of type t_QFI and nn of type t_INT and calls the PARI library function nupow(f,nn). The output is converted back to C longs and then MathLinked manually back to Mathematica as a Qfb[a',b',c'].

Although error checking is done, the big drawback is that while Mathematica and PARI deal in arbitrary precision reals and arbitrarily long integers, C does not. The link will crash if some parameter does not fit in a C long. This could be dealt with by more preprocessing, at the cost of speed.

The C code is easily produced with mprep:

./mprep ./parishell.tm -o parishell.tm.c

Notice that unlike the examples in the MathLink documentation, there is no need to have your C code in a separate file.

The C code is then compiled along with the headers pari.h and mathlink.h and the libraries libpari.a and libML.a. Both libpari.a and libML.a caused my compiler to give 'out of date' error messages, which the ranlib utility in Unix fixed.

(* The Mathematica session then looks similar to GP: *)

Out[1]=/Users/stopple

In[2]: SetDirectory["~/parishell/build"]

Out[2]=/Users/stopple/parishell/build

Out[3]={.DS_Store,parishell,parishell.build}

In[4]: link=Install["./parishell",LinkMode->Launch]

Out[4]=LinkObject[./parishell,2,2]

(* the C program parishell is now running in the background, waiting for calls from the Mathematica kernel *)

Kronecker[x,y] gives the Kronecker symbol (x/y).

In[6]: Kronecker[-23,13]

In[7]: ?QuadClassUnit

QuadClassUnit[d] gives the class number, a vector giving the structure of the

class group as a product of cyclic groups, a vector of binary quadratic forms

which generate those cyclic groups, and a real which measures the

correctness of the result. If it is close to 1., the result is correct

(under GRH). If it is close to a larger integer, the class number is off by

a factor equal to this integer.

In[8]: QuadClassUnit[-100023]

Out[8]={168,{42,2,2},{Qfb[2,1,12503],Qfb[77,77,344],Qfb[33,33,766]},0.997724}

In[9]:=QuadClassUnit[-13]

From In[9]:=QuadClassUnit::notdisc: -13 is not a discriminant.

In[10]: QuadClassUnit[-3]

Out[10]={1,{},{},1.}

In[11]: ?QfbPrimeForm

QfbPrimeForm[d,p] returns a prime binary quadratic form of discriminant d

whose first coefficient is the prime number p. Error if d is not a

discriminant, or is not a quadratic residue mod p.

In[12]: QfbPrimeForm[-21,5]

From In[12]:=QfbPrimeForm::notdisc: -21 is not a discriminant.

In[13]: QfbPrimeForm[-23,5]

From In[13]:=QfbPrimeForm::notsquare: -23 is not congruent to a square modulo 5.

In[14]: QfbPrimeForm[-23,6]

From In[14]:=QfbPrimeForm::notprime: 6 is not prime.

In[15]: QfbPrimeForm[-23,13]

QfbComp[f,g] is the composition of the primitive positive definite binary

quadratic forms f and g, using the NUCOMP and NUDUPL algorithms of Shanks (a

In[17]: QfbComp[Qfb[2,1,3],Qfb[2,-1,3]]

In[18]: QfbComp[Qfb[2,1,12503],Qfb[77,77,344]]

Out[18]=Qfb[154,77,172]

In[19]: QfbComp[Qfb[2,1,3],Qfb[2,1,7]]

From In[19]:=QfbComp::notequal: Qfb[2,1,3] and Qfb[2,1,7] have different discriminants.

In[20]: QfbComp[Qfb[2,1,-3],Qfb[2,-1,-3]]

From In[20]:=QfbComp::notneg: The discriminant 25 is non-negative.

QfbPow[f,n] is the n-th power of the primitive positive definite binary

quadratic form f using the NUCOMP and NUDUPL algorithms.

In[22]: QfbPow[Qfb[2,1,12503],168]

Out[22]=Qfb[1,1,25006]

In[23]: QfbPow[Qfb[2,1,-3],7]

From In[23]:=QfbPow::notneg: The discriminant Qfb[2,1,-3] is non-negative.

In[24]: Uninstall[link]

(* You need to Uninstall the link when you are finished to kill the C program. *)

/* Here is the template file parishell.tm. The :Evaluate: commands are all first, then the C code, then the templates. */

:Evaluate:disc[f_Qfb] := Apply[((#2)^2 - 4(#1)(#3)) &, f]

:Evaluate: Kronecker::usage = "Kronecker[x,y] gives the Kronecker symbol (x/y)."

:Evaluate: QuadClassUnit::usage = "QuadClassUnit[d] gives the class number, a vector giving the structure of the class group as a product of cyclic groups, a vector of binary quadratic forms which generate those cyclic groups, and a real which measures the correctness of the result. If it is close to 1., the result is correct (under GRH). If it is close to a larger integer, the class number is off by a factor equal to this integer."

:Evaluate:QuadClassUnit::notdisc = "`1` is not a discriminant."

:Evaluate:QuadClassUnit[d_Integer] :=
If[(d==-3)||(d==-4),{1, {}, {}, 1.},
If[(Mod[d, 4] == 0) || (Mod[d, 4] == 1), qcu[d],

        Message[QuadClassUnit::notdisc, d]]]

:Evaluate: QfbPrimeForm::usage = "QfbPrimeForm[d,p] returns a prime binary quadratic form of discriminant d whose first coefficient is the prime number p. Error if d is not a discriminant, or is not a quadratic residue mod p."

:Evaluate:QfbPrimeForm::notdisc = "`1` is not a discriminant."

:Evaluate:QfbPrimeForm::notprime = "`1` is not prime."

:Evaluate:QfbPrimeForm::notone = "`1` is not congruent to 1 modulo 8 or 0 modulo 4."

:Evaluate:QfbPrimeForm::notsquare = "`1` is not congruent to a square modulo `2`."

:Evaluate:QfbPrimeForm[d_Integer, p_Integer] :=
        If[(Mod[d, 4] == 2) || (Mod[d, 4] ==3),Message[QfbPrimeForm::notdisc, d],
       If[!PrimeQ[p], Message[QfbPrimeForm::notprime, p],
      If[p == 2, If[(Mod[d, 8] == 1) || (Mod[d, 4] == 0), qfbpf[d, 2],

Message[QfbPrimeForm::notone, d]],
If[JacobiSymbol[d, p] == -1, Message[QfbPrimeForm::notsquare, d, p],

qfbpf[d, p]]]]];

:Evaluate:QfbComp::usage = "QfbComp[f,g] is the composition of the primitive positive definite binary quadratic forms f and g, using the NUCOMP and NUDUPL algorithms of Shanks (a la Atkin)."

:Evaluate:QfbComp::notequal = "`1` and `2` have different discriminants."

:Evaluate:QfbComp::notneg = "The discriminant `1` is non-negative."

:Evaluate:QfbComp[f_Qfb, g_Qfb] :=
If[disc[f] == disc[g],

If[disc[f] < 0,

Apply[qfbc, Append[Join[Apply[List, f],Apply[List, g]], N[Abs[disc[f]]^(1/4)]]],

Message[QfbComp::notneg, disc[f]]],

Message[QfbComp::notequal, f, g]]

:Evaluate:QfbPow::usage = "QfbPow[f,n] is the n-th power of the primitive positive definite binary quadratic form f using the NUCOMP and NUDUPL algorithms."

:Evaluate:QfbPow::notneg = "The discriminant `1` is non-negative."

:Evaluate:QfbPow[f_Qfb, n_Integer] :=
If[disc[f] < 0,

Apply[qfbp, Append[Apply[List, f], n]],

Message[QfbPow::notneg, f]]

#include <pari/pari.h>

long kroneckersymbol(long x, long y);

long kroneckersymbol(long x, long y)

{
{
    return itos(gkronecker(stoi(x),stoi(y)));
}
return 0;
}

void quadclassunit(long d);

void quadclassunit(long d)
{
GEN v;
long ltop=avma;
int i,lenny;

v = quadclassunit0(stoi(d),0,NULL,2);
MLPutFunction(stdlink,"List",4);
MLPutLongInteger(stdlink,itos((GEN) v[1]));
lenny=lg(v[2])-1;
MLPutFunction(stdlink,"List",lenny);
for (i=1; i<=lenny;i++) {
     MLPutLongInteger(stdlink,itos((GEN) ((GEN) v[2])[i]));
}
MLPutFunction(stdlink,"List",lenny);
for (i=1; i<=lenny;i++) {
     MLPutFunction(stdlink,"Qfb",3);
     MLPutLongInteger(stdlink,itos((GEN)((GEN) ((GEN) v[3])[i])[1]));
     MLPutLongInteger(stdlink,itos((GEN)((GEN) ((GEN) v[3])[i])[2]));
     MLPutLongInteger(stdlink,itos((GEN)((GEN) ((GEN) v[3])[i])[3]));
}
MLPutDouble(stdlink, rtodbl((GEN) v[5]));
avma=ltop;
return;
}

void qfbprimeform(long d, long p);

void qfbprimeform(long d, long p)
{

GEN q;
long ltop = avma;

q=primeform(stoi(d),stoi(p),1);
MLPutFunction(stdlink,"Qfb",3);
MLPutLongInteger(stdlink,itos((GEN) q[1]));
MLPutLongInteger(stdlink,itos((GEN) q[2]));
MLPutLongInteger(stdlink,itos((GEN) q[3]));
avma=ltop;
return;
}

void qfbnucomp(long a, long b, long c, long r, long s, long t, double l );

void qfbnucomp(long a, long b, long c, long r, long s, long t, double l )
{

GEN f, g, h;
long ltop = avma;

f=qfi(stoi(a),stoi(b),stoi(c));
g=qfi(stoi(r),stoi(s),stoi(t));
h=nucomp(f,g,dbltor(l));
MLPutFunction(stdlink,"Qfb",3);
MLPutLongInteger(stdlink,itos((GEN) h[1]));
MLPutLongInteger(stdlink,itos((GEN) h[2]));
MLPutLongInteger(stdlink,itos((GEN) h[3]));
avma=ltop;
return;
}

void qfbnupow(long a, long b, long c, long n );

void qfbnupow(long a, long b, long c, long n )
{

GEN f, nn, h;
long ltop = avma;

f=qfi(stoi(a),stoi(b),stoi(c));
nn=stoi(n);
h=nupow(f,nn);
MLPutFunction(stdlink,"Qfb",3);
MLPutLongInteger(stdlink,itos((GEN) h[1]));
MLPutLongInteger(stdlink,itos((GEN) h[2]));
MLPutLongInteger(stdlink,itos((GEN) h[3]));
avma=ltop;

return;
}

int main(int argc, char* argv[])
{
pari_init(4000000, 500000);

        return MLMain(argc, argv);
}

:Begin:
:Function:       kroneckersymbol
:Pattern:        Kronecker[i_Integer, j_Integer]
:Arguments:      { i, j }
:ArgumentTypes: { LongInteger, LongInteger }
:ReturnType:     LongInteger
:End:

:Begin:
:Function:       quadclassunit
:Pattern:        qcu[d_Integer]
:Arguments:      { d }
:ArgumentTypes: { LongInteger }
:ReturnType:     Manual
:End:

:Begin:
:Function:       qfbprimeform
:Pattern:        qfbpf[d_Integer,p_Integer]
:Arguments:      { d, p }
:ArgumentTypes: { LongInteger, LongInteger }
:ReturnType:     Manual
:End:

:Begin:
:Function:       qfbnucomp
:Pattern:        qfbc[a_Integer,b_Integer,c_Integer,r_Integer,s_Integer,t_Integer, l_Real]
:Arguments:      { a, b, c, r, s, t, l }
:ArgumentTypes: { LongInteger, LongInteger, LongInteger, LongInteger, LongInteger, LongInteger, Real }
:ReturnType:     Manual
:End:

:Begin:
:Function:       qfbnupow
:Pattern:        qfbp[a_Integer,b_Integer,c_Integer,n_Integer]

:Arguments: { a, b, c, n }
:ArgumentTypes: { LongInteger, LongInteger, LongInteger, LongInteger }
:ReturnType: Manual
:End:

"The MathLink Reference Guide", at

http://www.mathsource.com/Content22/Enhancements/MathLink/0204-398

is about as helpful as most Wolfram documentation. Better is "A MathLink Tutorial" by Todd Gayley, which can be found at

http://www.mathsource.com/Content/Enhancements/MathLink/0206-693

The Mathematica notebook DeveloperGuide.nb in the MathLink distribution guides you through the examples. Some of what it says, at least for Mac OS X, is incorrect.

The book "MathLink: Network Programming with Mathematica" by Miyaji and Abbott, is the most detailed source.

I found it helpful to compile and run both the MathLink example 'addtwo' and also the PARI example 'matexp.c'