Main
Packages
Funding
SEARCH

Help / Community
FAQ
Documentation
Tutorials
Mailing Lists
Bugs
Timeline
Ateliers PARI/GP

Library
Publications
Contributed GP scripts
Fun!

Development
Latest Changes
Version Control
Coding Guidelines
PariDroid
Logo

Tests & benchmarks
Buildlogs
Coverage Report
Doc Coverage
Refcards test
Benchmarks

WWW Stats

Written by Jeffrey Stopple on pari-users (2002)

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]
```
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.

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

(* 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

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.

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

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

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.

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)."

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."

If[(d==-3)||(d==-4),{1, {}, {}, 1.},
If[(Mod[d, 4] == 0) || (Mod[d, 4] == 1), qcu[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);

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

{
GEN v;
long ltop = avma;
int i, lenny;

lenny=lg(v[2])-1;
for (i=1; i<=lenny;i++) {
}
for (i=1; i<=lenny;i++) {
}
avma = ltop;
}

void qfbprimeform(long d, long p)
{
long ltop = avma;
GEN q = primeform(stoi(d),stoi(p),1);
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));
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);
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:
: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:
```

PARI/GP Development