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.

Creation of PARI objects.creation The basic function which creates a PARI object is the function [*] whose syntax is:

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 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 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, [*] and [*] should be used instead of cgetg( ,1) and cgetg( ,2).

2) the macros [*], [*] and [*] are predefined as (long)cgetg, (long)cgeti and (long)cgetr respectively.

Examples:

z = cgeti(6); (or cgetg(6, 1);) creates an integer type which can hold numbers of absolute value less than 2128 (2 codewords + 4 mantissa longwords).

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.

z = cgetg(4, 19); for(i=1; i<4; i++) z[i] = lgetg(5, 18); creates a matrix type for 4×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 cgetg creates only the ``root'', and other calls to cgetg must be made to get the whole tree. For matrices, a common mistake is to think that 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.

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 cgetg. However we need its precise description: the first byte is the [*], the second byte is the [*] (used only for garbage collecting purposes), and the last two bytes (the low order word if you prefer) is the [*] of the root in longwords. For instance in the example z = cgeti(6) above, z[0] will be set to 0x01010006.

These bytes can be handled through the following functions:

long [*](GEN z); returns the type number of z.

void [*](GEN z, long n); sets equal to n the type number of z (you should not have to use this function if you use cgetg).

long [*](GEN z); returns the reference count of z.

void [*](GEN z, long n); sets equal to n the reference count of z.

void [*](GEN z); increments the reference count of z (with saturation at 255).

long [*](GEN z); returns the length (in longwords) of the root of z.

long [*](GEN z, long l); sets equal to l the length of z (you should not have to use this function if you use 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:

Type 1 (integers): 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 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:

long [*](GEN z); returns the sign of z.

void [*](GEN z, long s); sets equal to s the sign of z.

long [*](GEN z); returns the [*] of z.

void [*](GEN z, long l); sets equal to 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],...,z[lgef(z)-1]) in base 232, where as usual in this notation z[2] is the high order longword.

Type 2 (real numbers): 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 223). This exponent can be handled using the following functions:

long [*](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 long [*](GEN z) which tries to return an exponent for z even if z is not a real number.

void [*](GEN z, long e); sets equal to 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 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],...,z[lg(z]-1]) in base 232 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 0x02010004, 0xff800001, 0xe0000000, 0.

Type 3 (integermods): 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.

Types 4 and 5 (rational numbers): 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.

Type 6 (complex numbers): Complex numberz[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.

Type 7 (p-adic numbers): p-adic number this type has a second codeword z[1] which contains the following information: the first 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 215). This information can be handled using the following functions:

long [*](GEN z); returns the p-adic precision of z, i.e. one less than the number of significant p-adic digits.

void [*](GEN z, long l); sets equal to l the p-adic precision of z.

long [*](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.

void [*](GEN z, long e); sets equal to e the p-adic valuation of z.

In addition to this codeword, z[2] points to the prime p, z[3] points to pprecp(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 = pvalp(z)*(z[4] + O(z[3])) = pvalp(z)*(z[4] + O(pprecp(z)))

Type 8 (quadratic numbers): 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.

Type 9 (polymods): 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.

Type 10 (polynomials): 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...). This number can be handled with the following functions:

long [*](GEN z); returns the variable number of the object z.

Note also the function long [*](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.

void [*](GEN z,long v); sets equal to 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],...z[lgef(z)-1] point to the coefficients of the polynomial in ascending order, with z[2] being the constant term and so on.

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.

Type 11 (power series): 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 215. This information can be handled with the following functions:

[*], [*], [*], [*] as for polynomials, and [*], [*] for the exponent. Beware not to use expo and setexpo for power series.

If the power series is non zero, z[2], z[3],...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.

Types 13 and 14 (rational functions): Rational function z[1] points on the numerator, and z[2] on the denominator. Both must be of type polynomial.

Type 16 (binary quadratic forms): Binary quadratic form z[1], z[2], z[3] point on the three coefficients of the form. All three should be of type integer.

Types 17 and 18 (vectors): Row vectorColumn vector z[1], z[2],...z[lg(z)-1] point on the components of the vector.

Type 19 (matrices): Matrix z[1], z[2],...z[lg(z)-1] point on the column vectors of z, i.e. they must be of type 18 and of the same length.

Assignment and copying of PARI objects.

The general [*] function is the function [*] with the following syntax:

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 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 [*] with the following syntax:

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 [*], with the predefined macro [*] as a synonym for (long)gcopy. Its syntax is:

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 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 [*], with the predefined macro [*] as a synonym for (long)gclone. Its syntax is:

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:

GEN [*](long s);

GEN [*](double s);

Their meaning is clear from their names. The converse functions

long [*](GEN z);

double [*](GEN z);

also exist, but are seldom used.

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 [*] with the syntax:destruction

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):

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 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 [*] (with the macro [*] as a synonym to (long)gerepile). This function has the following syntax:

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 [*] (which stands for available memory address, 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 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:

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:

z=gerepile(ltop,lbot,gadd(p1,p2));

As you can see, in simple conditions the use of gerepile is not really difficult. However it is absolutely necessary to understand what has happened. When entering an instruction such as 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 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.

Important remark: as we will see presently it is often necessary to do several gerepile in a computation. However, the least the better. The only condition for 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:

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 before p2 then the garbage becomes connected and we get the following program with garbage collecting:

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 gmul in the file 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:

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 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 gmul, but let's ignore this for a moment). Then the most natural way is to do a cgetg of the result immediately, and then for each coefficient of the result vector to do a 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 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

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 dec has been declared as a long and not as a GEN:

dec=lpile(ltop,lbot,0)/4;z1+=dec;z2+=dec;

This works because gerepile thinks that 0 is a valid address, hence puts in dec the displacement in bytes. Since z1 and z2 are pointers, we need to divide dec by 4 to get a longword displacement.

If you followed us this far, congratulations, and rejoice, the rest is much easier.

Some tricks and hints. Even for the authors, the use of gerepile was not evident at first. Hence we give some indications on how to avoid most problems connected with gerepile, at the expense of speed.

First, although it looks complicated, 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 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 gerepile altogether by creating sufficiently large objects at the beginning (using cgetg), and then using assignment statements and operations ending with z (such as 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 gerepile at all:

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 gerepile is the following: keep the value of the stack pointer avma at the beginning, and just 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:

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 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 ltop=avma;. Then you should have no problems.

Another solution is to use clones, i.e. the function gclone (see 4.2.3), which creates a permanent object on the heap, and not on the stack.