MainDownload Packages Timeline Search SupportFAQ Documentation Tutorials Ateliers PARI/GP Mailing Lists GP scripts libraryContributed scripts DevelopmentBugs Latest Changes Version Control Coding Guidelines Tests & benchmarksBuildlogs Coverage report Benchmarks MiscellaneousWWW Stats Logo Fun! Links |
## MathLink and PARIWritten by Jeffrey Stopple (pari-users-555)
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 codeMathLink 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.cNotice 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. ## 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] := 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: |

PARI/GP Development

Last Modified: 2013-02-02 19:53:39

Copyleft © 2003-2013
the PARI group.