MathLink and PARI

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.

Writing the code

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]

  qfbp[a,b,c,n]

  qfbnupow(a,b,c,n)

  nupow(f,nn)

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

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.

Further Reading

• The MathLink Reference Guide is about as helpful as most Wolfram documentation.
• Better is A MathLink Tutorial, by Todd Gayley.
• 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'

Sample Mathematica session

  (* The Mathematica session then looks similar to GP: *)
  In[1]: Directory[]
  Out[1]=/Users/stopple

  In[2]: SetDirectory["~/parishell/build"]
  Out[2]=/Users/stopple/parishell/build

  In[3]: FileNames[]
  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 *)

  In[5]: ?Kronecker
  Kronecker[x,y] gives the Kronecker symbol (x/y).

  In[6]: Kronecker[-23,13]
  Out[6]=1

  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]
  Out[15]=Qfb[13,9,2]

  In[16]: ?QfbComp
  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).

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

  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.

  In[21]: ?QfbPow
  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]
  Out[24]=./parishell
  (* You need to Uninstall the link when you are finished to kill the C program. *)

Implementation / Template (parishell.tm.c)

  /* 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, p],
             Message[QfbPrimeForm::notone, d]]
        ,
          If[JacobiSymbol[d, p] != -1, qfbpf[d, p]
          ,
            Message[QfbPrimeForm::notsquare, 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]]
  

  long kroneckersymbol(long x, long y);
  void qfbnupow(long a, long b, long c, long n );
  void qfbnucomp(long a, long b, long c, long r, long s, long t, double l);
  void qfbprimeform(long d, long p);
  void quadclassunit(long d);

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

  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;
  }

  void qfbprimeform(long d, long p)
  {
    long ltop = avma;
    GEN 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;
  }

  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;
  }

  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;
  }

  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: