\chapter{Programming PARI in Library Mode} \section{Introduction: initializations, universal objects, input and output.} To be able to use PARI in \ref{library mode}, you must write a C program and link it to the PARI library and the PARI include files. See the installation guide (given in Appendix A) on how to create and install the PARI library and include files. A sample Makefile is also given in Appendix B. Probably the best way to understand how programming is done is to work through a complete example. We will write such a program in section 4.3. Before doing this, a few explanations are in order. \subsec{Initializations and universal objects.} First, one must explain to the outside world what kind of objects and programs we are going to use. This is done simply with the statement {\tt \#include } Note that this is usually a link, created when you make the library, to a file called either {\tt genpari68k.h} (for 680x0-based machines with x at least 2) or {\tt genpariother.h} (for all other machines). This file {\tt \ref{genpari.h}}, first exports all the necessary constants, variables and functions, defines some important macros, and also defines the fundamental type for all PARI objects: the type {\bf \ref{GEN}}, which is simply a pointer to {\tt long}. {\bf Technical note}: we would have liked to define a type GEN to be a pointer on itself. This unfortunately is not possible in C, except by using structures, but then the names become unwieldy. On the other hand, by a simple trick, it can be done in Pascal for example. The result of this is that when we will use a component of a PARI object, it will be a {\sl long}, hence will need to be typecasted to a GEN again if we want to avoid warnings from the compiler. This will sometimes be quite tedious, but of course is trivially done. To take an example, a polynomial $P$ of degree 2 will be represented by a chunk of memory pointed to by the GEN $P$. $P[0]$ and $P[1]$ contain code information, in particular the type of the object, the degree of $P$, etc\dots. $P[2]$, $P[3]$, $P[4]$ contain pointers to the coefficients of degree 0, 1, and 2 of $P$ respectively (note the ascending order). This is where typecasting will be necessary: in principle $P[i]$ (for $i=2,3,4$) is a long, but we will want to use it as a GEN. The coefficients $P[i]$ themselves are in chunks of memory whose complexity depends on the types of the coefficients, and so on. Now we must state the most important law about programming in PARI, which must be respected if one wants to avoid disasters: {\bf Apart from universal objects {\rm (see below)} the chunks of memory used by a given PARI object must be in consecutive memory locations}. Don't panic: let's see the reason and the meaning of this, and how it can be achieved. When doing large computations, unwanted intermediate results clutter up memory very fast so some kind of garbage collecting is needed. Most large systems do garbage collecting when the cluttering gets heavy, and this slows down the performance. In PARI we have taken a different approach: you must do your own cleaning up as soon as the intermediate results are not needed. Special purpose routines have been written to do this, but the primary requirement is exactly as stated above: a PARI object must be (essentially) connected. As a consequence of this explanation, one also sees that there is an evident exception to the above law: if your computation is small enough so that you don't need to do any garbage collecting, then just go ahead, PARI won't mind disconnected objects in most cases. However, since PARI routines do their own garbage collecting, watch carefully what you are doing. The notion of \ref{universal object} alluded to above is quite simple: during the execution of your program, a number of objects will have been defined (by the system or by yourself) with the idea that they stay permanently with the same values. Examples are the integer 1, the fraction $\dfrac 1 2$, the polynomial $X$, or a prime $p$ which is used as a base modulus for integermods or $p$-adics. These universal objects are of course allowed to be disconnected from the other PARI objects. After declaring the use of the file genpari.h, the first executable statement of a main program should be to initialize the PARI system, and in particular the PARI stack which will contain all the computations. This is simply done with a call to the two variable function {\bf \ref{init}}, like {\tt init(4000000,100000)}. The first argument (here 4 million) is the number of bytes given to PARI to work with (it should not reasonably drop under 500000), and the second is the upper limit on a precomputed prime number table. If you don't want prime numbers, just put 2, but put an argument anyway because {\tt init()} expects one. We have now at our disposal: $\bullet$ the following universal objects: the integer 0 ({\bf \ref{gzero}} as a GEN, {\tt \ref{zero}} as a long), the integer 1 ({\bf \ref{gun}} as a GEN, {\tt \ref{un}} as a long), the integer 2 ({\bf \ref{gdeux}} as a GEN, {\tt \ref{deux}} as a long), the fraction $\dfrac 1 2$ ({\bf \ref{ghalf}} as a GEN, {\tt \ref{lhalf}} as a long), the complex number $i$, {\bf \ref{gi}} as a GEN. In addition, space is reserved for the polynomials 1 ({\bf \ref{polun[]}} as a GEN, {\tt \ref{lpolun}} as a long), and the polynomials $x_v$, ({\bf \ref{polx[]}} as GENs, {\tt \ref{lpolx}} as longs), where $x_v$ is the name of variable number $v$, where $0\le v\le 255$. However, they are not created upon initialization, and it is the programmer's responsibility to fill them before use. Since this is not very easy, we advise the user to use the function {\tt \ref{lisseq}} which has essentially the same effect as {\tt \ref{lisexpr}} except that it can execute a sequence of expressions and not only a single expression. For example, to prepare for use the variables $a$,$b$,$c$,$x$, write \centerline{\tt lisseq("x;a;b;c")} Note that {\tt polun} and {\tt polx} are arrays, the index being the polynomial variable number. $\bullet$ a large PARI \ref{stack} containing nothing but (in the present version) the 167 long words (668 bytes) of the predefined universal objects. $\bullet$ a \ref{heap} which is dealt with in a different way from the stack, and will contain other permanent universal objects. $\bullet$ a table of primes. $\bullet$ access to all the built-in functions of the PARI library. We have already described many of these functions in the preceding chapters. However some of them are specific to library mode and thus will be explained in this chapter. \subsec{Input and output.} Two important aspects have not yet been explained since they are specific to library mode: input and output of PARI objects. For \ref{input}, PARI provides you with two powerful high level functions which enables you to input your objects as if you were under GP. In fact, the second one {\sl is \/} essentially the GP syntactical parser, hence you can use it not only for input but for any computation that you can do under GP. These functions are called {\bf \ref{lisexpr}} and {\bf \ref{lisseq}}. The first one has the following syntax: {\tt GEN lisexpr(char* s);} Its effect is to analyze the input string s and to compute the result as in GP. However it is limited to one expression. If you want to read and evaluate a sequence of expressions, use {\tt GEN lisseq(char* s);} Warning: there is a slight difference between these functions and the GP syntactical parser: the expressions and sequences which you use must not contain any spaces. Once in a while, it may be useful to have the evaluation of the string involving a call to a function you have defined in C. The function {\bf \ref{install}} allows you to give a name to a function taking 0, 1, 2 or 3 GEN arguments and returning a single GEN. The syntax is {\tt void install(GEN (*f)(), char *name, int valence)} where {\tt f} is the (address of) the function, {\tt name} its new name, and {\tt valence} the number of its arguments, an integer between 0 and 3. For \ref{output}, there exist four different functions. First, you can use the function {\bf \ref{sor}} with the following syntax: {\tt void sor(GEN x,char format,long dec,long field);} Here format is either {\tt 'e'}, {\tt 'f'} or {\tt 'g'} corresponding to the three output formats of GP, {\tt dec} is the number of printed significant digits for real numbers, and should be put equal to -1 if all of them are wanted, and {\tt field} corresponds to the field width of GP used for printing integers. A default use of this function is to use the macro {\tt \ref{outbeaut}(GEN x)} which is equivalent to {\tt sor(x,'g',-1,0)} followed by a newline and a buffer flush. The second format corresponds to the ``raw'' format of GP (see section 2.2.5) and is obtained by using the function {\bf \ref{brute}} with the following syntax: {\tt void brute(GEN x,char format,long dec);} A default use of this function is to use the macro {\tt \ref{output}(GEN x)} which is equivalent to {\tt brute(x,'g',-1)} followed by a newline and a buffer flush. Note that during debugging, the function {\tt brute} cannot be nicely used since the buffer is not flushed, and the macro {\tt output} cannot be executed, being a macro. For this reason only, Pari provides you with the fuunction {\tt void imprimer(GEN x)} which has the same effect as {\tt output}, but can now be called under a debugger. For example, under GDB, if x is a GEN, to see the contents of x simply type {\tt p imprimer(x)}. The third format corresponds to the {\tt texprint} function of GP, and gives a \TeX{} output of the result. It is obtained by using the function {\bf texe} with the following syntax: {\tt void \ref{texe}(GEN x,char format,long dec);} Finally, you can use the \ref{hexadecimal tree output} corresponding to the GP command \bs x using the function {\bf \ref{voir}} with the following syntax: {\tt void voir(GEN x,-1);} Again the last parameter must be given and put equal to -1. In principle this last type of output is used only for debugging purposes. {\bf Remark.} All Pari output is done on the stream {\tt \ref{outfile}}, which by default is initialized to {\tt stdout}. If you want that your output be directed to another file, you should use the function {\tt \ref{switchout}(name)} where {\tt name} is a character string giving the name of the file you are going to use. The printing is {\it appended\/} at the end of the file. When you want to close the file, you simply call {\tt switchout(NULL)}. Similarly, errors are sent to the stream {\tt \ref{errfile}} ({\tt stderr} by default), and input is done on the stream {\tt infile} using the function {\tt \ref{switchin}} analogous to {\tt stdout}. The is currently no way to send a Pari output to a string. \section{Creation, destruction, and implementation of the PARI objects.} By far the most important functions which are specific to the library mode are the functions which allow the programmer to create and delete PARI objects, and the assignment statements. \subsec{Creation of PARI objects.}\sref{creation} The basic function which creates a PARI object is the function {\bf \ref{cgetg}} whose syntax is: {\tt cgetg(long length,long type);} Here length is a long specifying the number of longwords to be allocated to the object, and type is a long which is the type number of the object (see chapter 1 or below for the list of these type numbers). Note that the length is the first argument, and type the second. The precise effect of this function is as follows: it first creates on the stack a chunk of memory of size ``length'' longwords, and puts the address of the chunk as the returned value. If not enough memory is available, a message to the effect that the PARI stack overflows will be printed. Otherwise, it sets the type of the chunk to ``type'', sets the length, and sets the reference count to 1. In effect, it fills correctly and completely the first codeword (z[0] or *z) of the PARI object. Many PARI objects having also a second codeword (types 1,2,7,10,11), {\it it is the programmer's responsibility to fill this second codeword\/}, either explicitly (see below), or implicitly using an assignment statement. Note that the argument ``length'' is forced for a number of types: 3 for types 3,4,5,6,9,13,14 and 16, 4 for type 8 (Quad), and 5 for type 7 ($p$-adic). However for efficiency's sake no checking is made in the function {\tt cgetg} so disastrous results can occur if you give an incorrect length. Notes: 1) the creation of leaves, i.e. integers or reals, being so common, {\bf \ref{cgeti}} and {\bf \ref{cgetr}} should be used instead of {\tt cgetg( ,1)} and {\tt cgetg( ,2)}. 2) the macros {\tt \ref{lgetg}}, {\tt \ref{lgeti}} and {\tt \ref{lgetr}} are predefined as {\tt (long)cgetg}, {\tt (long)cgeti} and {\tt (long)cgetr} respectively. Examples: {\tt z = cgeti(6);} (or {\tt cgetg(6, 1);}) creates an integer type which can hold numbers of absolute value less than $2^{128}$ (2 codewords + 4 mantissa longwords). {\tt z = cgetg(3, 6); z[1] = lgetr(5); z[2] = lgetr(5);} creates a complex type whose real and imaginary part can hold real numbers of precision 96 bits. {\tt z = cgetg(4, 19); for(i=1; i<4; i++) z[i] = lgetg(5, 18);} creates a matrix type for $4\times 3$ matrices. One must also create space for the matrix elements themselves using a double loop. These last two examples illustrates the fact that since PARI types are recursive, all the branches of the tree must be created. The function {\tt cgetg} creates only the ``root'', and other calls to {\tt cgetg} must be made to get the whole tree. For matrices, a common mistake is to think that {\tt z = cgetg(4, 19)} (for example) will create the root of the matrix: one needs also to create the column vectors of the matrix. This is because a matrix is really nothing else than a line vector of column vectors (hence a priori not a basic type), but it has been given a special type number so that operations with matrices become possible. \subsec{Implementation of the PARI types}. Although it is a little tedious, we must go through each type and explain their implementation. Let z be a PARI object. Common to all the types is the first codeword (z[0]), which we don't have to worry about since this is taken care of by {\tt cgetg}. However we need its precise description: the first byte is the \ref{type number}, the second byte is the \ref{reference count} (used only for garbage collecting purposes), and the last two bytes (the low order word if you prefer) is the \ref{length} of the root in longwords. For instance in the example {\tt z = cgeti(6)} above, z[0] will be set to {\tt 0x01010006}. These bytes can be handled through the following functions: {\tt long \ref{typ}(GEN z);} returns the type number of z. {\tt void \ref{settyp}(GEN z, long n);} sets equal to {\tt n} the type number of z (you should not have to use this function if you use cgetg). {\tt long \ref{pere}(GEN z);} returns the reference count of z. {\tt void \ref{setpere}(GEN z, long n);} sets equal to {\tt n} the reference count of z. {\tt void \ref{incpere}(GEN z);} increments the reference count of z (with saturation at 255). {\tt long \ref{lg}(GEN z);} returns the length (in longwords) of the root of z. {\tt long \ref{setlg}(GEN z, long l);} sets equal to {\tt l} the length of z (you should not have to use this function if you use {\tt cgetg}; however, see an example of the use of this function in the matexp function given in section 4.3). Let us now look precisely at the types: {\bf Type 1 (integers)}: \sref{Integer} this type has a second codeword z[1] which contains the following information. The first byte is the sign of z, i.e. 1 if z $>$ 0, 0 if z $=$ 0, $-1$ if z $<$ 0. The second byte is unused. The low order word is the {\bf effective length} of z, i.e. the total number of significant longwords. This means the following: apart from the integer 0 (whose second codeword is equal to 2), every integer is ``normalized'', meaning that the first mantissa longword (i.e. z[2]) is non zero. However, the integer may have been created with a longer length. Hence the ``length'' which is in z[0] can be larger than the ``effective length'' which is in z[1]. In fact, it would be a disaster to try to access z[i] for i larger than or equal to the effective length. This information can be handled using the following functions: {\tt long \ref{signe}(GEN z);} returns the sign of z. {\tt void \ref{setsigne}(GEN z, long s);} sets equal to {\tt s} the sign of z. {\tt long \ref{lgef}(GEN z);} returns the \ref{effective length} of z. {\tt void \ref{setlgef}(GEN z, long l);} sets equal to {\tt l} the effective length of z. The integer 0 can be characterized either by its sign equal to 0, or by its effective length equal to 2. Apart from z $=$ 0, z[2] exists and is non zero, and the absolute value of z is (z[2],z[3],\dots,z[lgef(z)-1]) in base $2^{32}$, where as usual in this notation z[2] is the high order longword. {\bf Type 2 (real numbers)}: \sref{Real number} this type has a second codeword z[1] whose first byte is also the sign (obtained or set using the same functions as for the integers), but whose 3 low order bytes contains a biased binary exponent (i.e. the exponent plus $2^{23}$). This exponent can be handled using the following functions: {\tt long \ref{expo}(GEN z);} returns the true (unbiased) exponent of z. This is defined even when z is equal to zero, see section 1.2.6. Note also the function {\tt long \ref{gexpo}(GEN z)} which tries to return an exponent for $z$ even if $z$ is not a real number. {\tt void \ref{setexpo}(GEN z, long e);} sets equal to {\tt e} the exponent of z, of course after adding the bias. The real zero is characterized by having its sign equal to 0. However, usually the first mantissa word z[2] is defined and equal to 0. This fact must {\it never\/} be used to characterize 0. If z is not equal to 0, the first mantissa word z[2] is normalized, i.e. its first bit is 1. The mantissa is (z[2],z[3],\dots,z[lg(z]-1]) in base $2^{32}$ where z[2] is the most significant longword, and the mantissa is understood to be between 1 and 2. Thus the real number $-3.5$ is represented (if the length is 4) as {\tt 0x02010004, 0xff800001, 0xe0000000, 0}. {\bf Type 3 (integermods)}: \sref{Integermod} z[1] points to the modulus, and z[2] on the number representing the class z. Both must be of type integer. In principle z[1] $>$ 0 and 0 $<=$ z[2] $<$ z[1], but this rule does not have to be strictly obeyed by the user. Any integermod obtained after a PARI operation will in principle satisfy these conditions. {\bf Types 4 and 5 (rational numbers)}: \sref{Rational number} z[1] points to the numerator, and z[2] on the denominator. Both must be of type integer. In principle z[2] $>$ 0, but this rule does not have to be strictly obeyed either. {\bf Type 6 (complex numbers)}: \sref{Complex number}z[1] points on the real part, and z[2] on the imaginary part. A priori z[1] and z[2] can be of any type, but only certain types are useful and make sense. {\bf Type 7 ($p$-adic numbers)}: \sref{p-adic number} this type has a second codeword z[1] which contains the following information: the first {\it word\/} contains the exponent of $p$ modulo which the $p$-adic unit corresponding to z is defined (if z is not 0), i.e. one less than the number of significant $p$-adic digits. the second word contains the biased exponent of z (the bias being here equal to $2^{15}$). This information can be handled using the following functions: {\tt long \ref{precp}(GEN z);} returns the $p$-adic precision of z, i.e. one less than the number of significant $p$-adic digits. {\tt void \ref{setprecp}(GEN z, long l);} sets equal to {\tt l} the $p$-adic precision of z. {\tt long \ref{valp}(GEN z);} returns the $p$-adic valuation of z (i.e. the unbiased exponent). This is defined even if z is equal to 0, see section 1.2.6. {\tt void \ref{setvalp}(GEN z, long e);} sets equal to {\tt e} the $p$-adic valuation of z. In addition to this codeword, z[2] points to the prime $p$, z[3] points to $p^{\text{precp(z)}}$, and z[4] points to an integer representing the $p$-adic unit associated to z modulo z[3] (or points to zero if z is zero). To summarize, we have the equality: $$z=p^{\text{valp(z)}}*(z[4]+O(z[3]))=p^{\text{valp(z)}}*(z[4]+O(p^{\text{precp(z)}}))$$ {\bf Type 8 (quadratic numbers)}: \sref{Quadratic number} z[1] points on the polynomial defining the quadratic field, z[2] on the ``real part'' and z[3] on the ``imaginary part'', where this is to be taken as the coefficients of z on the ``canonical'' basis (1,w) (see section 1.2.3). Complex numbers are particular cases of quadratics but deserve a separate type. {\bf Type 9 (polymods)}: \sref{Polymod} exactly as for integermods, z[1] points on the modulus, and z[2] on a polynomial representing the class of z. Both must be of type polynomial. However, one must obey the rules explained in chapter 2 concerning the hierarchical structure of the variables of a polymod. {\bf Type 10 (polynomials)}: \sref{Polynomial} this type has a second codeword which is analogous to the one for integers. The first byte is the ``sign'', i.e. 0 if the polynomial is equal to 0, and 1 if not (see however the important remark below). The second byte is the variable number (e.g. 0 for $X$, 1 for $Y$, etc\dots). This number can be handled with the following functions: {\tt long \ref{varn}(GEN z);} returns the variable number of the object z. Note also the function {\tt long \ref{gvar}(GEN z)} which tries to return a variable number for $z$, even if $z$ is not a polynomial or power series. The variable number of a scalar type is set by definition equal to 32767. {\tt void \ref{setvarn}(GEN z,long v);} sets equal to {\tt v} the variable number of z. The least significant word of the second codeword is the effective length of the polynomial. Note that the degree of the polynomial is equal to the effective length minus three. The functions used to handle this codeword are the same as for integers. The components z[2], z[3],\dots z[lgef(z)-1] point to the coefficients of the polynomial {\it in ascending order\/}, with z[2] being the constant term and so on. {\sl Important remark\/}. A zero polynomial can be characterized by the fact that its sign is 0. However, its effective length may be equal to 2, or greater than 2. If it is greater than 2, this means that all the coefficients of the polynomial are equal to zero (as they should for a zero polynomial), but not all of these zeros are exact zeros, and more precisely the leading term z[lgef(z)-1] is not an exact zero. {\bf Type 11 (power series)}: \sref{Power series} This type also has a second codeword. Its first byte is the ``sign'', i.e. 0 if the power series is 0, and 1 if not. The second byte is the variable number, as for polynomials. The last two bytes code for a biased exponent with a bias of $2^{15}$. This information can be handled with the following functions: {\tt \ref{signe}, \ref{setsigne}, \ref{varn}, \ref{setvarn}} as for polynomials, and {\tt \ref{valp}, \ref{setvalp}} for the exponent. Beware not to use {\tt expo} and {\tt setexpo} for power series. If the power series is non zero, z[2], z[3],\dots z[lg(z)-1] point to the coefficients of z in ascending order, z[2] being the first non-zero coefficient. Note that the exponent of a power series can be negative, i.e. we can deal with Laurent series with a finite number of negative terms. {\bf Types 13 and 14 (rational functions)}: \sref{Rational function} z[1] points on the numerator, and z[2] on the denominator. Both must be of type polynomial. {\bf Type 16 (binary quadratic forms)}: \sref{Binary quadratic form} z[1], z[2], z[3] point on the three coefficients of the form. All three should be of type integer. {\bf Types 17 and 18 (vectors)}: \sref{Row vector}\sref{Column vector} z[1], z[2],\dots z[lg(z)-1] point on the components of the vector. {\bf Type 19 (matrices)}: \sref{Matrix} z[1], z[2],\dots z[lg(z)-1] point on the column vectors of z, i.e. they must be of type 18 and of the same length. \subsec{Assignment and copying of PARI objects}. The general \ref{assignment} function is the function {\bf \ref{gaffect}} with the following syntax: {\tt void gaffect(GEN x, GEN y);} Its effect is to assign the PARI object x into the (preexisting) object y. Many conditions must be met for the assignment to be possible. For instance it is allowed to assign an integer into a real number, but the converse is forbidden. For that, you must use the truncation or rounding function of your choice (see section 3.2). It can also happen that y is not large enough or does not have the proper tree structure to receive the object x. As an extreme example, assume y is the zero integer with length equal to 2. Then all assignments of a non zero integer into y will result in an error message since y is not large enough to fit a non zero integer. In general common sense will tell you what is possible, keeping in mind the PARI philosophy which says that if it makes sense it is legal (the assignment of an imprecise object into a precise one does {\it not\/} make sense. However, a change in precision of imprecise objects is allowed). It is also necessary to assign ordinary 32-bit long numbers of C into a PARI object. This is done using the function {\bf \ref{gaffsg}} with the following syntax: {\tt void gaffsg(long s, GEN y);} It is also very useful to copy a PARI object, not just by moving around a pointer, but by physically creating a copy of the whole tree structure. The function which does this is called {\bf \ref{gcopy}}, with the predefined macro {\tt \ref{lcopy}} as a synonym for {\tt (long)gcopy}. Its syntax is: {\tt GEN gcopy(GEN x);} and the effect is to create a new copy of x on the PARI stack. However, it may also be useful to create a {\sl permanent} copy of a PARI object. This will not be created on the stack but on the heap. The function which does this is called {\bf \ref{gclone}}, with the predefined macro {\tt \ref{lclone}} as a synonym for {\tt (long)gclone}. Its syntax is: {\tt GEN gclone(GEN x);} A final class of functions should be mentioned in this subsection: functions which convert from C objects to PARI objects (with creation of the objects on the stack). They are as follows: {\tt GEN \ref{stoi}(long s); GEN \ref{dbltor}(double s);} Their meaning is clear from their names. The converse functions {\tt long \ref{gtolong}(GEN z); double \ref{gtodouble}(GEN z);} also exist, but are seldom used. \subsec{Destruction of PARI objects and garbage collection}. If you want to destroy (i.e. give back the memory occupied by) the latest PARI object on the stack (e.g. the latest one obtained by a cgetg or a function), this is very simple: just use the function {\bf \ref{cgiv}} with the syntax:\sref{destruction} {\tt void cgiv(z);} where z is the object (which can be a tree), which you want to give back. Unfortunately life is not so simple, and in general you are doing a computation and want to give back accumulated garbage. A simple example is the following. Let x and y be two preexisting PARI objects and suppose that we want to compute x$*$x$+$y$*$y. This can trivially be done using the following program (we skip the necessary declarations): {\tt p1=gmul(x,x);p2=gmul(y,y);z=gadd(p1,p2);} The GEN z will indeed point on the PARI object equal to x$*$x$+$y$*$y. However, consider the stack. It contains as unnecessary garbage p1 and p2. More precisely it contains (in that order) z, p1, p2. We need a way to get rid of this garbage (in this case it causes no harm except that it occupies memory space, but in other cases it could disconnect PARI objects and this is forbidden). It would not have been possible to get rid of p1, p2 before since they are used in the final operation. It is not possible to use the function {\tt cgiv} since p1 and p2 are not at the bottom of the stack. Hence PARI provides us with a much more powerful tool whose correct handling is not always easy: the function {\bf \ref{gerepile}} (with the macro {\tt \ref{lpile}} as a synonym to {\tt (long)gerepile}). This function has the following syntax: {\tt GEN gerepile(long ltop, long lbot, GEN z);} The understanding of the behavior of this function is certainly the most difficult part of this chapter, hence we will try to go slowly in the explanation. First, the PARI stack has a current stack pointer which is a global variable (hence must not be declared in a program) called {\bf \ref{avma}} (which stands for {\bf av}ailable {\bf m}emory {\bf a}ddress, which as its name indicates points always just after the first free address on the stack (remember that the stack grows from high to low addresses). For certain reasons this variable is of type long and not of type GEN as would be natural. Now let us see the effects of the function {\tt gerepile}. For this function to work one must have lbot $<$ ltop. As a first approximation, we describe the effects of this function as follows. Give back the memory locations between lbot and ltop, and move the object z upwards so that no space is lost. Specifically, we rewrite our previous example as follows: {\tt \obeylines\parskip=0pt plus 1pt \hbox{} ltop=avma; /* keep in mind the current address of the top of the stack */ p1=gmul(x,x); p2=gmul(y,y); lbot=avma; /* keep the address of the bottom of the garbage pile */ z=gadd(p1,p2); z=gerepile(ltop,lbot,z); /* do the garbage collecting */} The last two instructions could also have been written more simply: {\tt z=gerepile(ltop,lbot,gadd(p1,p2));} As you can see, in simple conditions the use of {\tt gerepile} is not really difficult. However it is absolutely necessary to understand what has happened. When entering an instruction such as {\tt gerepile(ltop,lbot,z)}, we must have avma $<=$ lbot $<$ ltop. As we will see, the variable z is in fact used very little. The function then considers all the PARI objects between avma and lbot (i.e. the ones that we want to keep) and looks at every component of such an object which is not a codeword. This component is a pointer on an object whose address is either between avma and lbot, in which case it will be suitably updated, larger than or equal to ltop, in which case it will not change, or between lbot and ltop in which case {\tt gerepile} will scream an error message at you. If all goes well, the pointers are suitably updated, the chunk of memory from avma to lbot$-$1 is physically copied to addresses avma$+$ltop$-$lbot to ltop$-$1, avma is updated (to the value avma$+$ltop$-$lbot), and finally the displacement ltop$-$lbot is added to z and given as a result of the function. In addition, if z is not 0 (in the sense of being address 0), the PARI object pointed to by ltop is checked to see whether certain of its pointers need also to be updated. This is a horrible hack which should not be used. {\bf Important remark}: as we will see presently it is often necessary to do several {\tt gerepile} in a computation. However, the least the better. The only condition for {\tt gerepile} to work is that the garbage be connected. If the computation can be arranged so that there is a minimal number of connected pieces of garbage then it should be done. For example suppose we want to write a function of two GEN variables x and y which creates the vector [x$*$x$+$y$,$y$*$y$+$x]. Without garbage collecting, one would write: {\tt \obeylines\parskip=0pt plus 1pt \hbox{} p1=gmul(x,x);p2=gadd(p1,y);p3=gmul(y,y);p4=gadd(p3,x); z=cgetg(3,17);z[1]=(long)p2;z[2]=(long)p4;} This leaves a dirty stack with (in that order) z,p4,p3,p2,p1 and p1 and p3 being garbage. But if we compute p3 {\it before \/} p2 then the garbage becomes connected and we get the following program with garbage collecting: {\tt \obeylines\parskip=0pt plus 1pt \hbox{} ltop=avma;p1=gmul(x,x);p3=gmul(y,y);lbot=avma; z=cgetg(3,17);z[1]=ladd(p1,y);z[2]=ladd(p3,x); z=gerepile(ltop,lbot,z);} We take our next example directly from the implementation of {\tt gmul} in the file {\tt gen2.c}, case 6 times case 6. In other words we want to write a program to compute the product of two complex numbers x and y, using a method which takes only 3 multiplications instead of 4. Let z=x$*$y, and set x=xr$+$i$*$xi and similarly for y and z. The well known trick is to compute p1=xr$*$yr, p2=xi$*$yi, p3=(xr$+$xi)$*$(yr$+$yi), and then we have zr=p1$-$p2, zi=p3$-$p1$-$p2. The program is essentially as follows: {\tt \obeylines\parskip=0pt plus 1pt \hbox{} ltop=avma;p1=gmul(x[1],y[1]);p2=gmul(x[2],y[2]); p3=gmul(gadd(x[1],x[2]),gadd(y[1],y[2]));p4=gadd(p1,p2);lbot=avma; z=cgetg(3,6);z[1]=lsub(p1,p2);z[2]=lsub(p3,p4); z=gerepile(ltop,lbot,z);} Let us now look at a less trivial example where more than one {\tt gerepile} is needed in practice (in theory one can always use only one, see below, but generally at the expense of efficiency). Suppose for instance that we want to write a function which multiplies a line vector by a matrix (such a function is of course already part of {\tt gmul}, but let's ignore this for a moment). Then the most natural way is to do a {\tt cgetg} of the result immediately, and then for each coefficient of the result vector to do a {\tt gerepile} to get rid of the garbage which has accumulated for that particular coefficient. We leave the details to the reader, who can look at the answer in the file gen2.c, in the function gmul, case 17 times case 19. As was said above, it would theoretically be possible to have a single connected piece of garbage, but it would be a much less natural and unnecessarily complicated program. Finally, let us take an example which is probably the least trivial way of using {\tt gerepile}, but is unfortunately sometimes necessary. Although it is not an infrequent occurrence, we will not give a specific example but a general one. Suppose that we want to do a computation (usually inside a larger function) giving more than one PARI object as a result, say two for instance. Then even if we set up the work properly, before cleaning up we will have a stack which has the desired results z1, z2 (say), and then connected garbage from lbot to ltop. If we write {\tt z1=gerepile(ltop,lbot,z1);} then the stack will be cleaned, the pointers OK, but we will have lost the address of z2. Hence the best way in this case is to write the following, assuming that {\tt dec} has been declared as a long and not as a GEN: {\tt dec=lpile(ltop,lbot,0)/4;z1+=dec;z2+=dec;} This works because {\tt gerepile} thinks that 0 is a valid address, hence puts in {\tt dec} the displacement in bytes. Since z1 and z2 are pointers, we need to divide {\tt dec} by 4 to get a longword displacement. If you followed us this far, congratulations, and rejoice, the rest is much easier. \subsec{Some tricks and hints}. Even for the authors, the use of {\tt gerepile} was not evident at first. Hence we give some indications on how to avoid most problems connected with {\tt gerepile}, at the expense of speed. First, although it looks complicated, {\tt gerepile} has turned out to be a very flexible and fast garbage collector, which can compare very favorably with much more sophisticated methods used in other systems. A few tests that we have made indicate that the price paid for using {\tt gerepile} is never more than 10 percent, and usually around 5 or 6 percent of the total time, which is quite acceptable. Second, in many cases, in particular when the tree structure and the size of the PARI objects which will appear in a computation are under control, one can avoid {\tt gerepile} altogether by creating sufficiently large objects at the beginning (using {\tt cgetg}), and then using assignment statements and operations ending with z (such as {\tt gaddz}). Coming back to our first example, note that if we know that x and y are of type real and of length less than or equal to 5, we can program without using {\tt gerepile} at all: {\tt z=cgetr(5);ltop=avma;p1=gmul(x,x);p2=gmul(y,y);gaddz(p1,p2,z);avma=ltop;} Here the cleaning up is done using the simple assignment avma=ltop which takes essentially no time at all. Third, the philosophy of {\tt gerepile} is the following: keep the value of the stack pointer avma at the beginning, and just {\bf before} the last operation. Afterwards, it would be too late since the garbage address would be lost. Finally, if everything seems hopeless, at the expense of speed you can do the following: after keeping in ltop the value of avma, perform your computation as you wish, in any order, leaving a messy stack. Let z be your result (let's assume you have only one. If you have more than one, you can always create a vector which is a single PARI object whose components are your results; after all, we are already cheating, we can cheat some more). Then write the following: {\tt lbot=avma;z=gerepile(ltop,lbot,gcopy(z));} The trick is to force a copy of z at the bottom of the stack, hence all the rest including the initial z becomes connected garbage. Note that this may not work in some cases where you have created new universal objects between lbot and ltop, for instance moduli for mod $p$ or $p$-adic computations. In that case you must also copy these objects on top of the stack and modify the components of your results which point to these objects. This in fact involves rewriting a simplified version of a part of {\tt gerepile}, hence is not recommended. If you are in this situation and still want to use this trick, we strongly advise you to create all your new universal objects before starting your computation, hence before the instruction {\tt ltop=avma;}. Then you should have no problems. Another solution is to use clones, i.e. the function {\tt gclone} (see 4.2.3), which creates a permanent object on the heap, and not on the stack. \section{A complete program.} Now that the preliminaries are out of the way, the best way to learn how to use the library mode is to work through a detailed nontrivial example of a main program. We will write a program which computes the exponential of a square matrix $x$. The complete listing is given in Appendix C, but each part of the program will be written and commented here. We will use an algorithm which is not optimal but is not far from the one used for the PARI function {\tt gexp} (in fact in the function {\tt mpexp1}). This consists in calculating the sum of the series: $$e^{x/(2^n)}=\sum_{k=0}^\infty \dfrac{(x/(2^n))^k}{k!}$$ for a suitable positive integer $n$, and then computing $e^x$ by repeated squarings. First, we will need to compute the $L^2$-norm of the matrix $x$, i.e. the quantity: $$z=\|x\|_2=\sqrt{\sum_{i,j}x_{i,j}^2}.$$ We will then choose the integer $n$ such that the $L^2$-norm of $x/(2^n)$ is less than or equal to 1, i.e. $$n=\left\lceil\dfrac{\ln(z)}{\ln(2)}\right\rceil$$ if $z>=1$. Then the convergence of the series will be at least as good as the one for $e^1$, and will be easy to estimate. In fact a larger value of $n$ would be preferable, but this is slightly machine dependent and more complicated, and will be left to the reader. We will now start writing our program. So as to be able to use it in other contexts, we will structure it in the following way: a main program which will do the input and output, and a function which we shall call {\bf matexp} which does the real work. The main program is easy to write. It can be something like this: {\tt \obeylines\parskip=0pt plus 1pt \hbox{} \#include GEN matexp(); \hbox{} main() \obr \quad GEN x,y; \quad long prec,d; \quad char s[512]; \hbox{} \quad init(1000000,2); /* take a million bytes of memory for the stack */ \quad printf("precision of the computation in decimal digits? "); \quad scanf("\%d",\&d);if(d<0) prec=3;else prec=(long)(d*K1+3); \quad printf("input your matrix in GP format:\bs n"); \quad s[0]=0;while(!s[0]) gets(s);x=lisexpr(s); \quad y=matexp(x,prec); \quad sor(y,'g',d,0); \cbr } The variable {\tt prec} represents the length in longwords of the real numbers used. K1 is the constant (defined in {\tt gencom.h}) equal to $\ln(10)/(32*\ln(2))$, which allows us to convert from a number of decimal digits to a number of longwords. The function {\tt lisexpr} was mentioned earlier and avoids any trouble in the input. In fact, as was also mentioned, {\tt lisexpr} can take any legal GP expression hence can do not only input but also computations. Note that we have used the construction {\tt s[0]=0;while(!s[0]) gets(s);} to get our string, instead of {\tt scanf} which would make trouble in this case. Finally, {\tt sor} is the general output routine. We have chosen to give d significant digits since this is what was asked. Note that there is a trick hidden here: if d was given to be negative, then the computation will be done in precision 3 (i.e. about 9.7 decimal digits) and in the function {\tt sor}, giving a negative third argument outputs all the significant digits, hence nothing is wrong. Now let us attack the main course, the function matexp: {\tt \obeylines\parskip=0pt plus 1pt \hbox{} GEN matexp(x,prec) \qquad GEN x; \qquad long prec; \hbox{} \obr \quad GEN y,r,s,p1,p2; \quad long tx=typ(x),lx=lg(x),i,k,n,lbot,ltop; \hbox{} /* check that x is a square matrix */ \hbox{} \quad if(tx!=19) \obr printf("This expression is not a matrix\bs n");return(gzero);\cbr \quad if(lx!=lg((GEN)x[1])) \obr printf("Not a square matrix\bs n");return(gzero);\cbr \hbox{} /* compute the L2 norm of x */ \hbox{} \quad ltop=avma;s=gzero; \quad for(i=1;i>}. We now come to the heart of the function. We have now a GEN p1 which points to a certain matrix of which we want to take the exponential. We will want to transform this matrix into a matrix with real (or complex of real) entries before starting the computation. To do this, we simply multiply by the real number 1 in precision prec$+$1 (to be safe). Then, to compute the series we will keep three variables: a variable p2 which at stage k will contain ${\tt p1}^k/k!$, a variable y which will contain $\sum_{i=0}^k {\tt p1}^i/i!$, and the variable r which will contain ${\tt s}^k/k!$. Note that we do not use Horner's rule. This is simply because we are lazy and do not want to compute in advance the number of terms that we need. We leave this modification (and many other improvements!) to the reader. The program continues as follows: {\tt \obeylines\parskip=0pt plus 1pt \hbox{} /* initializations before the loop */ \hbox{} \quad r=cgetr(prec+1);gaffsg(1,r);p1=gmul(r,p1); \quad y=gscalmat(r,lx-1); /* this creates the scalar matrix r$*{\tt I_{\text {lx-1}}}$ */ \quad p2=p1;r=s;k=1; \quad y=gadd(y,p2); \hbox{} /* now the main loop */ \hbox{} \quad while(expo(r) >= -32*(prec-1)) \qquad \obr k++;p2=gdivgs(gmul(p2,p1),k);r=gdivgs(gmul(s,r),k);y=gadd(y,p2);\cbr \hbox{} /* now square back n times if necessary */ \hbox{} \quad if(!n)\obr lbot=avma;y=gerepile(ltop,lbot,gcopy(y));\cbr \quad else \quad \obr \qquad for(i=0;i