NAME

libPARI - Programming PARI in Library Mode

@3The User's Guide to Pari/GP gives in three chapters a general presentation of the system, of the gp calculator, and detailed explanation of high level PARI routines available through the calculator. The present manual assumes general familiarity with the contents of these chapters and the basics of ANSI C programming, and focuses on the usage of the PARI library. In this chapter, we introduce the general concepts of PARI programming and describe useful general purpose functions; the following chapters describes all public low or high-level functions, underlying or extending the GP functions seen in Chapter 3 of the User's guide.


Introduction: initializations, universal objects

To use PARI in library mode, you must write a C program and link it to the PARI library. See the installation guide or the Appendix to the User's Guide to Pari/GP on how to create and install the library and include files. A sample Makefile is presented in Appendix A, and a more elaborate one in examples/Makefile. The best way to understand how programming is done is to work through a complete example. We will write such a program in Label se:prog. Before doing this, a few explanations are in order.

First, one must explain to the outside world what kind of objects and routines we are going to use. This is doneFOOTNOTE<<< This assumes that PARI headers are installed in a directory which belongs to your compiler's search path for header files. You might need to add flags like -I/usr/local/include or modify C_INCLUDE_PATH. >>> with the directive

  #include <pari/pari.h>

In particular, this defines the fundamental type for all PARI objects: the type GEN, which is simply a pointer to long.

Before any PARI routine is called, one must initialize the system, and in particular the PARI stack which is both a scratchboard and a repository for computed objects. This is done with a call to the function

void pari_init(size_t size, ulong maxprime)

@3The first argument is the number of bytes given to PARI to work with, and the second is the upper limit on a precomputed prime number table; size should not reasonably be taken below 500000 but you may set maxprime = 0, although the system still needs to precompute all primes up to about 2^{16}. For lower-level variants allowing finer control, e.g. preventing PARI from installing its own error or signal handlers, see Label se:pari_init_tech.

@3We have now at our disposal:

@3* a PARI stack containing nothing. This is a big connected chunk of size bytes of memory, where all computations take place. In large computations, intermediate results quickly clutter up memory so some kind of garbage collecting is needed. Most systems do garbage collecting when the memory is getting scarce, and this slows down the performance. PARI takes a different approach, admittedly more demanding on the programmer: you must do your own cleaning up when the intermediate results are not needed anymore. We will see later how (and when) this is done.

@3* the following universal objects (by definition, objects which do not belong to the stack): the integers 0, 1, -1, 2 and -2 (respectively called gen_0, gen_1, gen_m1, gen_2 and gen_m2), the fraction (1)/(2) (ghalf). All of these are of type GEN.

@3* a heap which is just a linked list of permanent universal objects. For now, it contains exactly the ones listed above. You will probably very rarely use the heap yourself; and if so, only as a collection of copies of objects taken from the stack (called clones in the sequel). Thus you need not bother with its internal structure, which may change as PARI evolves. Some complex PARI functions create clones for special garbage collecting purposes, usually destroying them when returning.

@3* a table of primes (in fact of differences between consecutive primes), called diffptr, of type byteptr (pointer to unsigned char). Its use is described in Label se:primetable later. Using it directly is deprecated, high-level iterators provide a cleaner and more flexible interface, see Label se:primeiter (such iterators use the private prime table, but extend it dynamically).

@3* access to all the built-in functions of the PARI library. These are declared to the outside world when you include pari.h, but need the above things to function properly. So if you forget the call to pari_init, you will get a fatal error when running your program.


Important technical notes

Backward compatibility

The PARI function names evolved over time, and deprecated functions are eventually deleted. The file pariold.h contains macros implementing a weak form of backward compatibility. In particular, whenever the name of a documented function changes, a #define is added to this file so that the old name expands to the new one (provided the prototype didn't change also).

This file is included by pari.h, but a large section is commented out by default. Define PARI_OLD_NAMES before including pari.h to pollute your namespace with lots of obsolete names like unFOOTNOTE<<< For (long)gen_1. Since 2004 and version 2.2.9, typecasts are completely unnecessary in PARI programs. >>>: that might enable you to compile old programs without having to modify them. The preferred way to do that is to add -DPARI_OLD_NAMES to your compiler CFLAGS, so that you don't need to modify the program files themselves.

Of course, it's better to fix the program if you can!

Types

Although PARI objects all have the C type GEN, we will freely use the word type to refer to PARI dynamic subtypes: t_INT, t_REAL, etc. The declaration

    GEN x;

declares a C variable of type GEN, but its ``value'' will be said to have type t_INT, t_REAL, etc. The meaning should always be clear from the context.

Type recursivity

Conceptually, most PARI types are recursive. But the GEN type is a pointer to long, not to GEN. So special macros must be used to access GEN's components. The simplest one is gel(V, i), where el stands for element, to access component number i of the GEN V. This is a valid lvalue (may be put on the left side of an assignment), and the following two constructions are exceedingly frequent

    gel(V, i) = x;
    x = gel(V, i);

where x and V are GENs. This macro accesses and modifies directly the components of V and do not create a copy of the coefficient, contrary to all the library functions.

More generally, to retrieve the values of elements of lists of...of lists of vectors we have the gmael macros (for multidimensional array element). The syntax is gmaeln(V,a_1,...,a_n), where V is a GEN, the a_i are indexes, and n is an integer between 1 and 5. This stands for x[a_1][a_2]...[a_n], and returns a GEN. The macros gel (resp. gmael) are synonyms for gmael1 (resp. gmael2).

Finally, the macro gcoeff(M, i, j) has exactly the meaning of M[i,j] in GP when M is a matrix. Note that due to the implementation of t_MATs as horizontal lists of vertical vectors, gcoeff(x,y) is actually equivalent to gmael(y,x). One should use gcoeff in matrix context, and gmael otherwise.

Variations on basic functions

In the library syntax descriptions in Chapter 3, we have only given the basic names of the functions. For example gadd(x,y) assumes that x and y are GENs, and creates the result x+y on the PARI stack. For most of the basic operators and functions, many other variants are available. We give some examples for gadd, but the same is true for all the basic operators, as well as for some simple common functions (a complete list is given in Chapter 6):

GEN gaddgs(GEN x, long y)

GEN gaddsg(long x, GEN y)

@3In the following one, z is a preexisting GEN and the result of the corresponding operation is put into z. The size of the PARI stack does not change:

void gaddz(GEN x, GEN y, GEN z)

@3(This last form is inefficient in general and deprecated outside of PARI kernel programming.) Low level kernel functions implement these operators for specialized arguments and are also available: Level 0 deals with operations at the word level (longs and ulongs), Level 1 with t_INT and t_REAL and Level 2 with the rest (modular arithmetic, polynomial arithmetic and linear algebra). Here are some examples of Level 1 functions:

GEN addii(GEN x, GEN y): here x and y are GENs of type t_INT (this is not checked).

GEN addrr(GEN x, GEN y): here x and y are GENs of type t_REAL (this is not checked).

There also exist functions addir, addri, mpadd (whose two arguments can be of type t_INT or t_REAL), addis (to add a t_INT and a long) and so on.

The Level 1 names are self-explanatory once you know that i stands for a t_INT, r for a t_REAL, mp for i or r, s for a signed C long integer, u for an unsigned C long integer; finally the suffix z means that the result is not created on the PARI stack but assigned to a preexisting GEN object passed as an extra argument. Chapter 6 gives a description of these low-level functions.

Level 2 names are more complicated, see Label se:level2names for all the gory details, and we content ourselves with a simple example used to implement t_INTMOD arithmetic:

GEN Fp_add(GEN x, GEN y, GEN m): returns the sum of x and y modulo m. Here x, y, m are t_INTs (this is not checked). The operation is more efficient if the inputs x, y are reduced modulo m, but this is not a necessary condition.

@3Important Note. These specialized functions are of course more efficient than the generic ones, but note the hidden danger here: the types of the objects involved (which is not checked) must be severely controlled, e.g. using addii on a t_FRAC argument will cause disasters. Type mismatches may corrupt the PARI stack, though in most cases they will just immediately overflow the stack. Because of this, the PARI philosophy of giving a result which is as exact as possible, enforced for generic functions like gadd or gmul, is dropped in kernel routines of Level 1, where it is replaced by the much simpler rule: the result is a t_INT if and only if all arguments are integer types (t_INT but also C long and ulong) and a t_REAL otherwise. For instance, multiplying a t_REAL by a t_INT always yields a t_REAL if you use mulir, where gmul returns the t_INT gen_0 if the integer is 0.

Portability: 32-bit / 64-bit architectures

PARI supports both 32-bit and 64-bit based machines, but not simultaneously! The library is compiled assuming a given architecture, and some of the header files you include (through pari.h) will have been modified to match the library.

Portable macros are defined to bypass most machine dependencies. If you want your programs to run identically on 32-bit and 64-bit machines, you have to use these, and not the corresponding numeric values, whenever the precise size of your long integers might matter. Here are the most important ones:

  64-bit    32-bit  X<BITS_IN_LONG>C<BITS_IN_LONG>    64        32 
X<LONG_IS_64BIT>C<LONG_IS_64BIT>   defined   undefined 
X<DEFAULTPREC>C<DEFAULTPREC>     3         4   (C< ~ > 19 decimal digits, 
see formula below) 
X<MEDDEFAULTPREC>C<MEDDEFAULTPREC>  4         6   (C< ~ > 38 decimal digits) 
X<BIGDEFAULTPREC>C<BIGDEFAULTPREC>  5         8   (C< ~ > 57 decimal digits) 
@3For instance, suppose you call a transcendental function, such as

GEN gexp(GEN x, long prec).

@3The last argument prec is an integer >= 3, corresponding to the default floating point precision required. It is only used if x is an exact object, otherwise the relative precision is determined by the precision of x. Since the parameter prec sets the size of the inexact result counted in (long) words (including codewords), the same value of prec will yield different results on 32-bit and 64-bit machines. Real numbers have two codewords (see Label se:impl), so the formula for computing the bit accuracy is

   bit_accuracy(prec) = (prec - 2) * BITS_IN_LONG

(this is actually the definition of an inline function). The corresponding accuracy expressed in decimal digits would be

   bit_accuracy(prec) * log (2) / log (10).

For example if the value of prec is 5, the corresponding accuracy for 32-bit machines is (5-2)* log (2^{32})/ log (10) ~ 28 decimal digits, while for 64-bit machines it is (5-2)* log (2^{64})/ log (10) ~ 57 decimal digits.

Thus, you must take care to change the prec parameter you are supplying according to the bit size, either using the default precisions given by the various DEFAULTPRECs, or by using conditional constructs of the form:

  #ifndef LONG_IS_64BIT
    prec = 4;
  #else
    prec = 6;
  #endif

@3which is in this case equivalent to the statement prec = MEDDEFAULTPREC;.

Note that for parity reasons, half the accuracies available on 32-bit architectures (the odd ones) have no precise equivalents on 64-bit machines.

Using malloc / free

You should make use of the PARI stack as much as possible, and avoid allocating objects using the customary functions. If you do, you should use, or at least have a very close look at, the following wrappers:

void* pari_malloc(size_t size) calls malloc to allocate size bytes and returns a pointer to the allocated memory. If the request fails, an error is raised. The SIGINT signal is blocked until malloc returns, to avoid leaving the system stack in an inconsistent state.

void* pari_realloc(void* ptr, size_t size) as pari_malloc but calls realloc instead of malloc.

void* pari_calloc(size_t size) as pari_malloc, setting the memory to zero.

void pari_free(void* ptr) calls free to liberate the memory space pointed to by ptr, which must have been allocated by malloc (pari_malloc) or realloc (pari_realloc). The SIGINT signal is blocked until free returns.

If you use the standard libc functions instead of our wrappers, then your functions will be subtly incompatible with the gp calculator: when the user tries to interrupt a computation, the calculator may crash (if a system call is interrupted at the wrong time).


Garbage collection

Why and how

As we have seen, pari_init allocates a big range of addresses, the stack, that are going to be used throughout. Recall that all PARI objects are pointers. Except for a few universal objects, they all point at some part of the stack.

The stack starts at the address bot and ends just before top. This means that the quantity

   (top - bot)/sizeof(long)

is (roughly) equal to the size argument of pari_init. The PARI stack also has a ``current stack pointer'' called avma, which stands for available memory address. These three variables are global (declared by pari.h). They are of type pari_sp, which means pari stack pointer.

The stack is oriented upside-down: the more recent an object, the closer to bot. Accordingly, initially avma = top, and avma gets decremented as new objects are created. As its name indicates, avma always points just after the first free address on the stack, and (GEN)avma is always (a pointer to) the latest created object. When avma reaches bot, the stack overflows, aborting all computations, and an error message is issued. To avoid this you need to clean up the stack from time to time, when intermediate objects are not needed anymore. This is called ``garbage collecting.''

We are now going to describe briefly how this is done. We will see many concrete examples in the next subsection.

@3* First, PARI routines do their own garbage collecting, which means that whenever a documented function from the library returns, only its result(s) have been added to the stack, possibly up to a very small overhead (non-documented ones may not do this). In particular, a PARI function that does not return a GEN does not clutter the stack. Thus, if your computation is small enough (e.g. you call few PARI routines, or most of them return long integers), then you do not need to do any garbage collecting. This is probably the case in many of your subroutines. Of course the objects that were on the stack before the function call are left alone. Except for the ones listed below, PARI functions only collect their own garbage.

@3* It may happen that all objects that were created after a certain point can be deleted --- for instance, if the final result you need is not a GEN, or if some search proved futile. Then, it is enough to record the value of avma just before the first garbage is created, and restore it upon exit:

  pari_sp av = avma; /* record initial avma */
  garbage ...
  avma = av; /* restore it */

@3All objects created in the garbage zone will eventually be overwritten: they should no longer be accessed after avma has been restored.

@3* 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 from a function call), you can use the function

void cgiv(GEN z)

@3where z is the object you want to give back. This is equivalent to the above where the initial av is computed from z.

@3* Unfortunately life is not so simple, and sometimes you will want to give back accumulated garbage during a computation without losing recent data. We shall start with the lowest level function to get a feel for the underlying mechanisms, we shall describe simpler variants later:

GEN gerepile(pari_sp ltop, pari_sp lbot, GEN q). This function cleans up the stack between ltop and lbot, where lbot < ltop, and returns the updated object q. This means:

1) we translate (copy) all the objects in the interval [avma, lbot[, so that its right extremity abuts the address ltop. Graphically

               bot           avma   lbot          ltop     top
  End of stack  |-------------[++++++[-/-/-/-/-/-/-|++++++++|  Start
                  free memory            garbage

@3becomes:

               bot                         avma   ltop     top
  End of stack  |---------------------------[++++++[++++++++|  Start
                         free memory

@3where ++ denote significant objects, -- the unused part of the stack, and -/- the garbage we remove.

2) The function then inspects all the PARI objects between avma and lbot (i.e. the ones that we want to keep and that have been translated) and looks at every component of such an object which is not a codeword. Each such component is a pointer to an object whose address is either

--- between avma and lbot, in which case it is suitably updated,

--- larger than or equal to ltop, in which case it does not change, or

--- between lbot and ltop in which case gerepile raises an error (``significant pointers lost in gerepile'').

3) avma is updated (we add ltop - lbot to the old value).

4) We return the (possibly updated) object q: if q initially pointed between avma and lbot, we return the updated address, as in 2). If not, the original address is still valid, and is returned!

As stated above, no component of the remaining objects (in particular q) should belong to the erased segment [lbot, ltop[, and this is checked within gerepile. But beware as well that the addresses of the objects in the translated zone change after a call to gerepile, so you must not access any pointer which previously pointed into the zone below ltop. If you need to recover more than one object, use the gerepileall function below.

@3Remark. As a consequence of the preceding explanation, if a PARI object is to be relocated by gerepile then, apart from universal objects, the chunks of memory used by its components should be in consecutive memory locations. All GENs created by documented PARI functions are guaranteed to satisfy this. This is because the gerepile function knows only about two connected zones: the garbage that is erased (between lbot and ltop) and the significant pointers that are copied and updated. If there is garbage interspersed with your objects, disaster occurs when we try to update them and consider the corresponding ``pointers''. In most cases of course the said garbage is in fact a bunch of other GENs, in which case we simply waste time copying and updating them for nothing. But be wary when you allow objects to become disconnected.

@3In practice this is achieved by the following programming idiom:

    ltop = avma; garbage(); lbot = avma; q = anything();
    return gerepile(ltop, lbot, q); /* returns the updated q */

@3or directly

    ltop = avma; garbage(); lbot = avma;
    return gerepile(ltop, lbot, anything());

Beware that

    ltop = avma; garbage();
    return gerepile(ltop, avma, anything())

@3might work, but should be frowned upon. We cannot predict whether avma is evaluated after or before the call to anything(): it depends on the compiler. If we are out of luck, it is after the call, so the result belongs to the garbage zone and the gerepile statement becomes equivalent to avma = ltop. Thus we return a pointer to random garbage.

Variants

GEN gerepileupto(pari_sp ltop, GEN q). Cleans the stack between ltop and the connected object q and returns q updated. For this to work, q must have been created before all its components, otherwise they would belong to the garbage zone! Unless mentioned otherwise, documented PARI functions guarantee this.

GEN gerepilecopy(pari_sp ltop, GEN x). Functionally equivalent to, but more efficient than

    gerepileupto(ltop, gcopy(x))

@3In this case, the GEN parameter x need not satisfy any property before the garbage collection: it may be disconnected, components created before the root, and so on. Of course, this is about twice slower than either gerepileupto or gerepile, because x has to be copied to a clean stack zone first. This function is a special case of gerepileall below, where n = 1.

void gerepileall(pari_sp ltop, int n, ...). To cope with complicated cases where many objects have to be preserved. The routine expects n further arguments, which are the addresses of the GENs you want to preserve:

    pari_sp ltop = avma;
    ...; y = ...; ... x = ...; ...;
    gerepileall(ltop, 2, &x, &y);

It cleans up the most recent part of the stack (between ltop and avma), updating all the GENs added to the argument list. A copy is done just before the cleaning to preserve them, so they do not need to be connected before the call. With gerepilecopy, this is the most robust of the gerepile functions (the less prone to user error), hence the slowest.

void gerepileallsp(pari_sp ltop, pari_sp lbot, int n, ...). More efficient, but trickier than gerepileall. Cleans the stack between lbot and ltop and updates the GENs pointed at by the elements of gptr without any further copying. This is subject to the same restrictions as gerepile, the only difference being that more than one address gets updated.

Examples

gerepile

Let x and y be two preexisting PARI objects and suppose that we want to compute x^2 + y^2. This is done using the following program:

    GEN x2 = gsqr(x);
    GEN y2 = gsqr(y), z = gadd(x2,y2);

The GEN z indeed points at the desired quantity. However, consider the stack: it contains as unnecessary garbage x2 and y2. More precisely it contains (in this order) z, y2, x2. (Recall that, since the stack grows downward from the top, the most recent object comes first.)

It is not possible to get rid of x2, y2 before z is computed, since they are used in the final operation. We cannot record avma before x2 is computed and restore it later, since this would destroy z as well. It is not possible either to use the function cgiv since x2 and y2 are not at the bottom of the stack and we do not want to give back z.

But using gerepile, we can give back the memory locations corresponding to x2, y2, and move the object z upwards so that no space is lost. Specifically:

    pari_sp ltop = avma;  /* remember the current top of the stack */
    GEN x2 = gsqr(x);
    GEN y2 = gsqr(y);
    pari_sp lbot = avma;  /* the bottom of the garbage pile */
    GEN z = gadd(x2, y2); /* z is now the last object on the stack */
    z = gerepile(ltop, lbot, z);

@3Of course, the last two instructions could also have been written more simply:

    z = gerepile(ltop, lbot, gadd(x2,y2));

@3In fact gerepileupto is even simpler to use, because the result of gadd is the last object on the stack and gadd is guaranteed to return an object suitable for gerepileupto:

    ltop = avma;
    z = gerepileupto(ltop, gadd(gsqr(x), gsqr(y)));

Make sure you understand exactly what has happened before you go on!

@3Remark on assignments and gerepile. When the tree structure and the size of the PARI objects which will appear in a computation are under control, one may allocate sufficiently large objects at the beginning, use assignment statements, then simply restore avma. Coming back to the above example, note that if we know that x and y are of type real fitting into DEFAULTPREC words, we can program without using gerepile at all:

    z = cgetr(DEFAULTPREC); ltop = avma;
    gaffect(gadd(gsqr(x), gsqr(y)), z);
    avma = ltop;

@3This is often slower than a craftily used gerepile though, and certainly more cumbersome to use. As a rule, assignment statements should generally be avoided.

@3Variations on a theme. it is often necessary to do several gerepiles during a computation. However, the fewer 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 that way.

For example suppose we want to write a function of two GEN variables x and y which creates the vector [x^2+y, y^2+x]. Without garbage collecting, one would write:

    p1 = gsqr(x); p2 = gadd(p1, y);
    p3 = gsqr(y); p4 = gadd(p3, x);
    z = mkvec2(p2, p4);  /* not suitable for gerepileupto! */

This leaves a dirty stack containing (in this order) z, p4, p3, p2, p1. The garbage here consists of p1 and p3, which are separated by p2. But if we compute p3 before p2 then the garbage becomes connected, and we get the following program with garbage collecting:

    ltop = avma; p1 = gsqr(x); p3 = gsqr(y);
    lbot = avma; z = cgetg(3, t_VEC);
    gel(z, 1) = gadd(p1,y);
    gel(z, 2) = gadd(p3,x); z = gerepile(ltop,lbot,z);

@3Finishing by z = gerepileupto(ltop, z) would be ok as well. Beware that

    ltop = avma; p1 = gadd(gsqr(x), y); p3 = gadd(gsqr(y), x);
    z = cgetg(3, t_VEC);
    gel(z, 1) = p1;
    gel(z, 2) = p3; z = gerepileupto(ltop,z); /* WRONG */

is a disaster since p1 and p3 are created before z, so the call to gerepileupto overwrites them, leaving gel(z, 1) and gel(z, 2) pointing at random data! The following does work:

    ltop = avma; p1 = gsqr(x); p3 = gsqr(y);
    lbot = avma; z = mkvec2(gadd(p1,y), gadd(p3,x));
    z = gerepile(ltop,lbot,z);

@3but is very subtly wrong in the sense that z = gerepileupto(ltop, z) would not work. The reason being that mkvec2 creates the root z of the vector after its arguments have been evaluated, creating the components of z too early; gerepile does not care, but the created z is a time bomb which will explode on any later gerepileupto. On the other hand

    ltop = avma; z = cgetg(3, t_VEC);
    gel(z, 1) = gadd(gsqr(x), y);
    gel(z, 2) = gadd(gsqr(y), x); z = gerepileupto(ltop,z); /* INEFFICIENT */

leaves the results of gsqr(x) and gsqr(y) on the stack (and lets gerepileupto update them for naught). Finally, the most elegant and efficient version (with respect to time and memory use) is as follows

    z = cgetg(3, t_VEC);
    ltop = avma; gel(z, 1) = gerepileupto(ltop, gadd(gsqr(x), y));
    ltop = avma; gel(z, 2) = gerepileupto(ltop, gadd(gsqr(y), x));

which avoids updating the container z and cleans up its components individually, as soon as they are computed.

@3One last example. Let us compute the product of two complex numbers x and y, using the 3M method which requires 3 multiplications instead of the obvious 4. Let z = x*y, and set x = x_r + i*x_i and similarly for y and z. We compute p_1 = x_r*y_r, p_2 = x_i*y_i, p_3 = (x_r+x_i)*(y_r+y_i), and then we have z_r = p_1-p_2, z_i = p_3-(p_1+p_2). The program is as follows:

  ltop = avma;
  p1 = gmul(gel(x,1), gel(y,1));
  p2 = gmul(gel(x,2), gel(y,2));
  p3 = gmul(gadd(gel(x,1), gel(x,2)), gadd(gel(y,1), gel(y,2)));
  p4 = gadd(p1,p2);
  lbot = avma; z = cgetg(3, t_COMPLEX);
  gel(z, 1) = gsub(p1,p2);
  gel(z, 2) = gsub(p3,p4); z = gerepile(ltop,lbot,z);

@3Exercise. Write a function which multiplies a matrix by a column vector. Hint: start with a cgetg of the result, and use gerepile whenever a coefficient of the result vector is computed. You can look at the answer in src/basemath/RgV.c:RgM_RgC_mul().

gerepileall

Let us now see why we may need the gerepileall variants. Although it is not an infrequent occurrence, we do not give a specific example but a general one: suppose that we want to do a computation (usually inside a larger function) producing 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 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 is cleaned, the pointers fixed up, but we have lost the address of z2. This is where we need the gerepileall function:

    gerepileall(ltop, 2, &z1, &z2)

@3copies z1 and z2 to new locations, cleans the stack from ltop to the old avma, and updates the pointers z1 and z2. Here we do not assume anything about the stack: the garbage can be disconnected and z1, z2 need not be at the bottom of the stack. If all of these assumptions are in fact satisfied, then we can call gerepilemanysp instead, which is usually faster since we do not need the initial copy (on the other hand, it is less cache friendly).

A most important usage is ``random'' garbage collection during loops whose size requirements we cannot (or do not bother to) control in advance:

    pari_sp ltop = avma, limit = stack_lim(avma, 1);
    GEN x, y;
    while (...)
    {
      garbage(); x = anything();
      garbage(); y = anything(); garbage();
      if (avma < limit) /* memory is running low (half spent since entry) */
        gerepileall(ltop, 2, &x, &y);
    }

@3Here we assume that only x and y are needed from one iteration to the next. As it would be costly to call gerepile once for each iteration, we only do it when it seems to have become necessary. The macro stack_lim(avma,n) denotes an address where 2^{n-1} / (2^{n-1}+1) of the remaining stack space is exhausted (1/2 for n = 1, 2/3 for n = 2).

Comments

First, gerepile has turned out to be a flexible and fast garbage collector for number-theoretic computations, which compares favorably with more sophisticated methods used in other systems. Our benchmarks indicate that the price paid for using gerepile and gerepile-related copies, when properly used, is usually less than 1% of the total running time, which is quite acceptable!

Second, it is of course harder on the programmer, and quite error-prone if you do not stick to a consistent PARI programming style. If all seems lost, just use gerepilecopy (or gerepileall) to fix up the stack for you. You can always optimize later when you have sorted out exactly which routines are crucial and what objects need to be preserved and their usual sizes.

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


Creation of PARI objects, assignments, conversions

Creation of PARI objects

The basic function which creates a PARI object is

GEN cgetg(long l, long t) l specifies the number of longwords to be allocated to the object, and t is the type of the object, in symbolic form (see Label se:impl for the list of these). The precise effect of this function is as follows: it first creates on the PARI stack a chunk of memory of size length longwords, and saves the address of the chunk which it will in the end return. If the stack has been used up, a message to the effect that ``the PARI stack overflows'' is printed, and an error raised. Otherwise, it sets the type and length of the PARI object. In effect, it fills its first codeword (z[0]). Many PARI objects also have a second codeword (types t_INT, t_REAL, t_PADIC, t_POL, and t_SER). In case you want to produce one of those from scratch, which should be exceedingly rare, it is your responsibility to fill this second codeword, either explicitly (using the macros described in Label se:impl), or implicitly using an assignment statement (using gaffect).

Note that the length argument l is predetermined for a number of types: 3 for types t_INTMOD, t_FRAC, t_COMPLEX, t_POLMOD, t_RFRAC, 4 for type t_QUAD and t_QFI, and 5 for type t_PADIC and t_QFR. However for the sake of efficiency, cgetg does not check this: disasters will occur if you give an incorrect length for those types.

@3Notes. 1) The main use of this function is create efficiently a constant object, or to prepare for later assignments (see Label se:assign). Most of the time you will use GEN objects as they are created and returned by PARI functions. In this case you do not need to use cgetg to create space to hold them.

@32) For the creation of leaves, i.e. t_INT or t_REAL,

GEN cgeti(long length)

GEN cgetr(long length)

@3should be used instead of cgetg(length, t_INT) and cgetg(length, t_REAL) respectively. Finally

GEN cgetc(long prec)

@3creates a t_COMPLEX whose real and imaginary part are t_REALs allocated by cgetr(prec).

@3Examples. 1) Both z = cgeti(DEFAULTPREC) and cgetg(DEFAULTPREC, t_INT) create a t_INT whose ``precision'' is bit_accuracy(DEFAULTPREC) = 64. This means z can hold rational integers of absolute value less than 2^{64}. Note that in both cases, the second codeword is not filled. Of course we could use numerical values, e.g. cgeti(4), but this would have different meanings on different machines as bit_accuracy(4) equals 64 on 32-bit machines, but 128 on 64-bit machines.

@32) The following creates a complex number whose real and imaginary parts can hold real numbers of precision bit_accuracy(MEDDEFAULTPREC) = 96 bits:

    z = cgetg(3, t_COMPLEX);
    gel(z, 1) = cgetr(MEDDEFAULTPREC);
    gel(z, 2) = cgetr(MEDDEFAULTPREC);

or simply z = cgetc(MEDDEFAULTPREC).

@33) To create a matrix object for 4 x 3 matrices:

    z = cgetg(4, t_MAT);
    for(i=1; i<4; i++) gel(z, i) = cgetg(5, t_COL);

or simply z = zeromatcopy(4, 3), which further initializes all entries to gen_0.

These last two examples illustrate 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 produce the whole tree. For matrices, a common mistake is to think that z = cgetg(4, t_MAT) (for example) creates the root of the matrix: one needs also to create the column vectors of the matrix (obviously, since we specified only one dimension in the first cgetg!). This is because a matrix is really just a row 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.

Finally, to facilitate input of constant objects when speed is not paramount, there are four varargs functions:

GEN mkintn(long n, ...) returns the non-negative t_INT whose development in base 2^{32} is given by the following n words (unsigned long). It is assumed that all such arguments are less than 2^{32} (the actual sizeof(long) is irrelevant, the behavior is also as above on 64-bit machines).

    mkintn(3, a2, a1, a0);

@3returns a_2 2^{64} + a_1 2^{32} + a_0.

GEN mkpoln(long n, ...) Returns the t_POL whose n coefficients (GEN) follow, in order of decreasing degree.

    mkpoln(3, gen_1, gen_2, gen_0);

@3returns the polynomial X^2 + 2X (in variable 0, use setvarn if you want other variable numbers). Beware that n is the number of coefficients, hence one more than the degree.

GEN mkvecn(long n, ...) returns the t_VEC whose n coefficients (GEN) follow.

GEN mkcoln(long n, ...) returns the t_COL whose n coefficients (GEN) follow.

@3Warning. Contrary to the policy of general PARI functions, the latter three functions do not copy their arguments, nor do they produce an object a priori suitable for gerepileupto. For instance

    /* gerepile-safe: components are universal objects */
    z = mkvecn(3, gen_1, gen_0, gen_2);
    /* not OK for gerepileupto: stoi(3) creates component before root */
    z = mkvecn(3, stoi(3), gen_0, gen_2);
    /* NO! First vector component C<x> is destroyed */
    x = gclone(gen_1);
    z = mkvecn(3, x, gen_0, gen_2);
    gunclone(x);

@3The following function is also available as a special case of mkintn:

GEN uu32toi(ulong a, ulong b)

Returns the GEN equal to 2^{32} a + b, assuming that a,b < 2^{32}. This does not depend on sizeof(long): the behavior is as above on both 32 and 64-bit machines.

Sizes

long gsizeword(GEN x) returns the total number of BIL-bit words occupied by the tree representing x.

long gsizebyte(GEN x) returns the total number of bytes occupied by the tree representing x, i.e. gsizeword(x) multiplied by sizeof(long). This is normally useless since PARI functions use a number of words as input for lengths and precisions.

Assignments

Firstly, if x and y are both declared as GEN (i.e. pointers to something), the ordinary C assignment y = x makes perfect sense: we are just moving a pointer around. However, physically modifying either x or y (for instance, x[1] = 0) also changes the other one, which is usually not desirable.

@3Very important note. Using the functions described in this paragraph is inefficient and often awkward: one of the gerepile functions (see Label se:garbage) should be preferred. See the paragraph end for one exception to this rule.

The general PARI assignment function is the function gaffect with the following syntax:

void gaffect(GEN x, GEN y)

Its effect is to assign the PARI object x into the preexisting object y. Both x and y must be scalar types. For convenience, vector or matrices of scalar types are also allowed.

This copies the whole structure of x into y so many conditions must be met for the assignment to be possible. For instance it is allowed to assign a t_INT into a t_REAL, but the converse is forbidden. For that, you must use the truncation or rounding function of your choice, e.g.mpfloor.

It can also happen that y is not large enough or does not have the proper tree structure to receive the object x. For instance, let y the zero integer with length equal to 2; then y is too small to accommodate any non-zero t_INT. In general common sense tells you what is possible, keeping in mind the PARI philosophy which says that if it makes sense it is valid. For instance, the assignment of an imprecise object into a precise one does not make sense. However, a change in precision of imprecise objects is allowed, even if it increases its accuracy: we complement the ``mantissa'' with infinitely many 0 digits in this case. (Mantissa between quotes, because this is not restricted to t_REALs, it also applies for p-adics for instance.)

All functions ending in ``z'' such as gaddz (see Label se:low_level) implicitly use this function. In fact what they exactly do is record {avma} (see Label se:garbage), perform the required operation, gaffect the result to the last operand, then restore the initial avma.

You can assign ordinary C long integers into a PARI object (not necessarily of type t_INT) using

void gaffsg(long s, GEN y)

@3Note. Due to the requirements mentioned above, it is usually a bad idea to use gaffect statements. There is one exception: for simple objects (e.g. leaves) whose size is controlled, they can be easier to use than gerepile, and about as efficient.

@3Coercion. It is often useful to coerce an inexact object to a given precision. For instance at the beginning of a routine where precision can be kept to a minimum; otherwise the precision of the input is used in all subsequent computations, which is inefficient if the latter is known to thousands of digits. One may use the gaffect function for this, but it is easier and more efficient to call

GEN gtofp(GEN x, long prec) converts the complex number x (t_INT, t_REAL, t_FRAC, t_QUAD or t_COMPLEX) to either a t_REAL or t_COMPLEX whose components are t_REAL of length prec.

Copy

It is also very useful to copy a PARI object, not just by moving around a pointer as in the y = x example, but by creating a copy of the whole tree structure, without pre-allocating a possibly complicated y to use with gaffect. The function which does this is called gcopy. Its syntax is:

GEN gcopy(GEN x)

@3and the effect is to create a new copy of x on the PARI stack.

Sometimes, on the contrary, a quick copy of the skeleton of x is enough, leaving pointers to the original data in x for the sake of speed instead of making a full recursive copy. Use GEN shallowcopy(GEN x) for this. Note that the result is not suitable for gerepileupto !

Make sure at this point that you understand the difference between y = x, y = gcopy(x), y = shallowcopy(x) and gaffect(x,y).

Clones

Sometimes, it is more efficient to create a persistent copy of a PARI object. This is not created on the stack but on the heap, hence unaffected by gerepile and friends. The function which does this is called gclone. Its syntax is:

GEN gclone(GEN x)

A clone can be removed from the heap (thus destroyed) using

void gunclone(GEN x)

@3No PARI object should keep references to a clone which has been destroyed!

Conversions

The following functions convert C objects to PARI objects (creating them on the stack as usual):

GEN stoi(long s): C long integer (``small'') to t_INT.

GEN dbltor(double s): C double to t_REAL. The accuracy of the result is 19 decimal digits, i.e. a type t_REAL of length DEFAULTPREC, although on 32-bit machines only 16 of them are significant.

@3We also have the converse functions:

long itos(GEN x): x must be of type t_INT,

double rtodbl(GEN x): x must be of type t_REAL,

@3as well as the more general ones:

long gtolong(GEN x),

double gtodouble(GEN x).


Implementation of the PARI types

We now go through each type and explain its implementation. Let z be a GEN, pointing at a PARI object. In the following paragraphs, we will constantly mix two points of view: on the one hand, z is treated as the C pointer it is, on the other, as PARI's handle on some mathematical entity, so we will shamelessly write z != 0 to indicate that the value thus represented is nonzero (in which case the pointer z is certainly non-NULL). We offer no apologies for this style. In fact, you had better feel comfortable juggling both views simultaneously in your mind if you want to write correct PARI programs.

Common to all the types is the first codeword z[0], which we do not have to worry about since this is taken care of by cgetg. Its precise structure depends on the machine you are using, but it always contains the following data: the internal type number associated to the symbolic type name, the length of the root in longwords, and a technical bit which indicates whether the object is a clone or not (see Label se:clone). This last one is used by gp for internal garbage collecting, you will not have to worry about it.

Some types have a second codeword, different for each type, which we will soon describe as we will shortly consider each of them in turn.

@3The first codeword is handled through the following macros:

long typ(GEN z) returns the type number of z.

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

long lg(GEN z) returns the length (in longwords) of the root of z.

long setlg(GEN z, long l) sets the length of z to l (you should not have to use this function if you use cgetg; however, see an advanced example in Label se:prog).

long isclone(GEN z) is z a clone?

void setisclone(GEN z) sets the clone bit.

void unsetisclone(GEN z) clears the clone bit.

@3Important remark. For the sake of efficiency, none of the codeword-handling macros check the types of their arguments even when there are stringent restrictions on their use. It is trivial to create invalid objects, or corrupt one of the ``universal constants'' (e.g. setting the sign of gen_0 to 1), and they usually provide negligible savings. Use higher level functions whenever possible.

@3Remark. The clone bit is there so that gunclone can check it is deleting an object which was allocated by gclone. Miscellaneous vector entries are often cloned by gp so that a GP statement like v[1] = x does not involve copying the whole of v: the component v[1] is deleted if its clone bit is set, and is replaced by a clone of x. Don't set/unset yourself the clone bit unless you know what you are doing: in particular never set the clone bit of a vector component when the said vector is scheduled to be uncloned. Hackish code may abuse the clone bit to tag objects for reasons unrelated to the above instead of using proper data structures. Don't do that.

Type t_INT (integer)

this type has a second codeword z[1] which contains the following information:

the sign of z: coded as 1, 0 or -1 if z > 0, z = 0, z < 0 respectively.

the effective length of z, i.e. the total number of significant longwords. This means the following: apart from the integer 0, every integer is ``normalized'', meaning that the most significant mantissa longword 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].

@3This information is handled using the following macros:

long signe(GEN z) returns the sign of z.

void setsigne(GEN z, long s) sets the sign of z to s.

long lgefint(GEN z) returns the effective length of z.

void setlgefint(GEN z, long l) sets the effective length of z to l.

The integer 0 can be recognized either by its sign being 0, or by its effective length being equal to 2. Now assume that z != 0, and let

   |z |= sum_{i = 0}^n z_i B^i, {where} z_n != 0 {and} B = 2^{BITS_IN_LONG}.

With these notations, n is lgefint(z) - 3, and the mantissa of z may be manipulated via the following interface:

GEN int_MSW(GEN z) returns a pointer to the most significant word of z, z_n.

GEN int_LSW(GEN z) returns a pointer to the least significant word of z, z_0.

GEN int_W(GEN z, long i) returns the i-th significant word of z, z_i. Accessing the i-th significant word for i > n yields unpredictable results.

GEN int_W_lg(GEN z, long i, long lz) returns the i-th significant word of z, z_i, assuming lgefint(z) is lz ( = n + 3). Accessing the i-th significant word for i > n yields unpredictable results.

GEN int_precW(GEN z) returns the previous (less significant) word of z, z_{i-1} assuming z points to z_i.

GEN int_nextW(GEN z) returns the next (more significant) word of z, z_{i+1} assuming z points to z_i.

Unnormalized integers, such that z_n is possibly 0, are explicitly forbidden. To enforce this, one may write an arbitrary mantissa then call

void int_normalize(GEN z, long known0)

@3normalizes in place a non-negative integer (such that z_n is possibly 0), assuming at least the first known0 words are zero.

@3For instance a binary and could be implemented in the following way:

  GEN AND(GEN x, GEN y) {
    long i, lx, ly, lout;
    long *xp, *yp, *outp; /* mantissa pointers */
    GEN out;
    if (!signe(x) || !signe(y)) return gen_0;
    lx = lgefint(x); xp = int_LSW(x);
    ly = lgefint(y); yp = int_LSW(y); lout = min(lx,ly); /* > 2 */
    out = cgeti(lout); out[1] = evalsigne(1) | evallgefint(lout);
    outp = int_LSW(out);
    for (i=2; i < lout; i++)
    {
      *outp = (*xp) & (*yp);
      outp  = int_nextW(outp);
      xp    = int_nextW(xp);
      yp    = int_nextW(yp);
    }
    if ( !*int_MSW(out) ) out = int_normalize(out, 1);
    return out;
  }

@3This low-level interface is mandatory in order to write portable code since PARI can be compiled using various multiprecision kernels, for instance the native one or GNU MP, with incompatible internal structures (for one thing, the mantissa is oriented in different directions).

@3The following further functions are available:

int mpodd(GEN x) which is 1 if x is odd, and 0 otherwise.

long mod2(GEN x)

long mod4(GEN x)

long mod8(GEN x)

long mod16(GEN x)

long mod32(GEN x)

long mod64(GEN x) give the residue class of x modulo the corresponding power of 2, for positive x. By definition, modn(x) := modn(|x|) for x < 0 (the functions disregard the sign), and the result is undefined if x = 0. As well,

ulong mod2BIL(GEN x) returns the least significant word of |x|, still assuming that x != 0.

These functions directly access the binary data and are thus much faster than the generic modulo functions. Besides, they return long integers instead of GENs, so they do not clutter up the stack.

Type t_REAL (real number)

this type has a second codeword z[1] which also encodes its sign, obtained or set using the same functions as for a t_INT, and a binary exponent. This exponent is handled using the following macros:

long expo(GEN z) returns the exponent of z. This is defined even when z is equal to zero, see Label se:whatzero.

void setexpo(GEN z, long e) sets the exponent of z to e.

@3Note the functions:

long gexpo(GEN z) which tries to return an exponent for z, even if z is not a real number.

long gsigne(GEN z) which returns a sign for z, even when z is neither real nor integer (a rational number for instance).

The real zero is characterized by having its sign equal to 0. If z is not equal to 0, then is is represented as 2^e M, where e is the exponent, and M belongs to [1, 2[ is the mantissa of z, whose digits are stored in z[2],..., z[lg(z)-1].

More precisely, let m be the integer (z[2],..., z[lg(z)-1]) in base 2^BITS_IN_LONG; here, z[2] is the most significant longword and is normalized, i.e. its most significant bit is 1. Then we have M := m / 2^{bit_accuracy(lg(z)) - 1 - expo(z)}.

GEN mantissa_real(GEN z, long *e) returns the mantissa m of z, and sets *e to the exponent bit_accuracy(lg(z))-1-expo(z), so that z = m / 2^e.

Thus, the real number 3.5 to accuracy bit_accuracy(lg(z)) is represented as z[0] (encoding type = t_REAL, lg(z)), z[1] (encoding sign = 1, expo = 1), z[2] = 0xe0000000, z[3] = .. .= z[lg(z)-1] = 0x0.

Type t_INTMOD

z[1] points to the modulus, and z[2] at the number representing the class z. Both are separate GEN objects, and both must be t_INTs, satisfying the inequality 0 <= z[2] < z[1].

Type t_FRAC (rational number)

z[1] points to the numerator n, and z[2] to the denominator d. Both must be of type t_INT such that n != 0, d > 0 and (n,d) = 1.

Type t_FFELT (finite field element)

(Experimental)

Components of this type should normally not be accessed directly. Instead, finite field elements should be created using ffgen.

@3The second codeword z[1] determines the storage format of the element, among

@3* t_FF_FpXQ: A = z[2] and T = z[3] are FpX, p = z[4] is a t_INT, where p is a prime number, T is irreducible modulo p, and deg A < deg T. This represents the element A (mod T) in F_p[X]/T.

@3* t_FF_Flxq: A = z[2] and T = z[3] are Flx, l = z[4] is a t_INT, where l is a prime number, T is irreducible modulo l, and deg A < deg T This represents the element A (mod T) in F_l[X]/T.

@3* t_FF_F2xq: A = z[2] and T = z[3] are F2x, l = z[4] is the t_INT 2, T is irreducible modulo 2, and deg A < deg T. This represents the element A (mod T) in F_2[X]/T.

Type t_COMPLEX (complex number)

z[1] points to the real part, and z[2] to the imaginary part. The components z[1] and z[2] must be of type t_INT, t_REAL or t_FRAC. For historical reasons t_INTMOD and t_PADIC are also allowed (the latter for p = 2 or congruent to 3 mod 4 only), but one should rather use the more general t_POLMOD construction.

Type t_PADIC (p-adic numbers)

this type has a second codeword z[1] which contains the following information: the p-adic precision (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, and the exponent of z. This information can be handled using the following functions:

long precp(GEN z) returns the p-adic precision of z. This is 0 if z = 0.

void setprecp(GEN z, long l) sets the p-adic precision of z to l.

long valp(GEN z) returns the p-adic valuation of z (i.e. the exponent). This is defined even if z is equal to 0, see Label se:whatzero.

void setvalp(GEN z, long e) sets the p-adic valuation of z to e.

In addition to this codeword, z[2] points to the prime p, z[3] points to p^{{precp(z)}}, and z[4] points to at_INT representing the p-adic unit associated to z modulo z[3] (and to zero if z is zero). To summarize, if z != 0, we have the equality:

   z = p^{{valp(z)}} * (z[4] + O(z[3])), {where} z[3] = O(p^{{precp(z)}}).

Type t_QUAD (quadratic number)

z[1] points to the canonical polynomial P defining the quadratic field (as output by quadpoly), z[2] to the ``real part'' and z[3] to the ``imaginary part''. The latter are of type t_INT, t_FRAC, t_INTMOD, or t_PADIC and are to be taken as the coefficients of z with respect to the canonical basis (1,X) or Q[X]/(P(X)), see Label se:compquad. Exact complex numbers may be implemented as quadratics, but t_COMPLEX is in general more versatile (t_REAL components are allowed) and more efficient.

Operations involving a t_QUAD and t_COMPLEX are implemented by converting the t_QUAD to a t_REAL (or t_COMPLEX with t_REAL components) to the accuracy of the t_COMPLEX. As a consequence, operations between t_QUAD and exact t_COMPLEXs are not allowed.

Type t_POLMOD (polmod)

as for t_INTMODs, z[1] points to the modulus, and z[2] to a polynomial representing the class of z. Both must be of type t_POL in the same variable, satisfying the inequality deg z[2] < deg z[1]. However, z[2] is allowed to be a simplification of such a polynomial, e.g. a scalar. This is tricky considering the hierarchical structure of the variables; in particular, a polynomial in variable of lesser priority (see Label se:priority) than the modulus variable is valid, since it is considered as the constant term of a polynomial of degree 0 in the correct variable. On the other hand a variable of greater priority is not acceptable; see Label se:priority for the problems which may arise.

Type t_POL (polynomial)

this type has a second codeword. It contains a ``sign'': 0 if the polynomial is equal to 0, and 1 if not (see however the important remark below) and a variable number (e.g. 0 for x, 1 for y, etc...).

@3These data can be handled with the following macros: signe and setsigne as for t_INT and t_REAL,

long varn(GEN z) returns the variable number of the object z,

void setvarn(GEN z, long v) sets the variable number of z to v.

The variable numbers encode the relative priorities of variables as discussed in Label se:priority. We will give more details in Label se:vars. Note also the function long 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 NO_VARIABLE, which has lower priority than any other variable number.

The components z[2], z[3],...z[lg(z)-1] point to the coefficients of the polynomial in ascending order, with z[2] being the constant term and so on.

For a t_POL of non-zero sign, degpol, leading_term, constant_term, return its degree, and a pointer to the leading, resp. constant, coefficient with respect to the main variable. Note that no copy is made on the PARI stack so the returned value is not safe for a basic gerepile call. Applied to any other type than t_POL, the result is unspecified. Those three functions are still defined when the sign is 0, see Label se:accessors and Label se:polynomials.

long degree(GEN x) returns the degree of x with respect to its main variable even when x is not a polynomial (a rational function for instance). By convention, the degree of a zero polynomial is -1.

@3Important remark. The leading coefficient of a t_POL may be equal to zero:

@3* it is not allowed to be an exact rational 0, such as gen_0;

@3* an exact non-rational 0, like Mod(0,2), is possible for constant polynomials, i.e. of length 3 and no other coefficient: this carries information about the base ring for the polynomial;

@3* an inexact 0, like 0.E-38 or O(3^5), is always possible. Inexact zeroes do not correspond to an actual 0, but to a very small coefficient according to some metric; we keep them to give information on how much cancellation occurred in previous computations.

A polynomial disobeying any of these rules is an invalid unnormalized object. We advise not to use low-level constructions to build a t_POL coefficient by coefficient, such as

    GEN T = cgetg(4, t_POL);
    T[1] = evalvarn(0);
    gel(T, 2) = x;
    gel(T, 3) = y;

@3But if you do and it is not clear whether the result will be normalized, call

GEN normalizepol(GEN x) applied to an unnormalized t_POL x (with all coefficients correctly set except that leading_term(x) might be zero), normalizes x correctly in place and returns x. This functions sets signe (to 0 or 1) properly.

@3Caveat. A consequence of the remark above is that zero polynomials are characterized by the fact that their sign is 0. It is in general incorrect to check whether lg(x) is 2 or degpol(x) < 0, although both tests are valid when the coefficient types are under control: for instance, when they are guaranteed to be t_INTs or t_FRACs. The same remark applies to t_SERs.

Type t_SER (power series)

This type also has a second codeword, which encodes a ``sign'', i.e. 0 if the power series is 0, and 1 if not, a variable number as for polynomials, and an exponent. This information can be handled with the following functions: signe, setsigne, varn, setvarn as for polynomials, and valp, setvalp for the exponent as for p-adic numbers. Beware: do not use expo and setexpo on power series.

The coefficients z[2], z[3],...z[lg(z)-1] point to the coefficients of z in ascending order. As for polynomials (see remark there), the sign of a t_SER is 0 if and only all its coefficients are equal to 0. (The leading coefficient cannot be an integer 0.)

Note that the exponent of a power series can be negative, i.e. we are then dealing with a Laurent series (with a finite number of negative terms).

Type t_RFRAC (rational function)

z[1] points to the numerator n, and z[2] on the denominator d. The denominator must be of type t_POL, with variable of higher priority than the numerator. The numerator n is not an exact 0 and (n,d) = 1 (see gred_rfac2).

Type t_QFR (indefinite binary quadratic form)

z[1], z[2], z[3] point to the three coefficients of the form and are of type t_INT. z[4] is Shanks's distance function, and must be of type t_REAL.

Type t_QFI (definite binary quadratic form)

z[1], z[2], z[3] point to the three coefficients of the form. All three are of type t_INT.

Type t_VEC and t_COL (vector)

z[1], z[2],...z[lg(z)-1] point to the components of the vector.

Type t_MAT (matrix)

z[1], z[2],...z[lg(z)-1] point to the column vectors of z, i.e. they must be of type t_COL and of the same length.

Type t_VECSMALL (vector of small integers)

z[1], z[2],...z[lg(z)-1] are ordinary signed long integers. This type is used instead of a t_VEC of t_INTs for efficiency reasons, for instance to implement efficiently permutations, polynomial arithmetic and linear algebra over small finite fields, etc.

Type t_STR (character string)

char * GSTR(z) ( = (z+1)) points to the first character of the (NULL-terminated) string.

Type t_ERROR (error context)error

context> This type holds error messages, as well as details about the error, as returned by the exception handling system. The second codeword z[1] contains the error type (an int, as passed to pari_err). The subsequent words z[2],...z[lg(z)-1] are GENs containing additional data, depending on the error type.

Type t_CLOSURE (closure)

This type holds GP functions and closures, in compiled form. The internal detail of this type is subject to change each time the GP language evolves. Hence we do not describe it here and refer to the Developer's Guide. However functions to create or to evaluate t_CLOSUREs are documented in Label se:closure.

long closure_arity(GEN C) returns the arity of the t_CLOSURE.

Type t_LIST (list)

this type was introduced for specific gp use and is rather inefficient compared to a straightforward linked list implementation (it requires more memory, as well as many unnecessary copies). Hence we do not describe it here and refer to the Developer's Guide.

@3Implementation note. For the types including an exponent (or a valuation), we actually store a biased non-negative exponent (bit-ORing the biased exponent to the codeword), obtained by adding a constant to the true exponent: either HIGHEXPOBIT (for t_REAL) or HIGHVALPBIT (for t_PADIC and t_SER). Of course, this is encapsulated by the exponent/valuation-handling macros and needs not concern the library user.


PARI variables

=head2 Multivariate objects

@3We now consider variables and formal computations, and give the technical details corresponding to the general discussion in Label se:priority. As we have seen in Label se:impl, the codewords for types t_POL and t_SER encode a ``variable number''. This is an integer, ranging from 0 to MAXVARN. Relative priorities may be ascertained using

int varncmp(long v, long w)

@3which is > 0, = 0, < 0 whenever v has lower, resp. same, resp. higher priority than w.

The way an object is considered in formal computations depends entirely on its ``principal variable number'' which is given by the function

long gvar(GEN z)

@3which returns 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 NO_VARIABLE which has lower priority than any valid variable number. The variable number of a recursive type which is not a polynomial or power series is the variable number with highest priority among its components. But for polynomials and power series only the ``outermost'' number counts (we directly access varn(x) in the codewords): the representation is not symmetrical at all.

Under gp, one needs not worry too much since the interpreter defines the variables as it sees themFOOTNOTE<<< The first time a given identifier is read by the GP parser a new variable is created, and it is assigned a strictly lower priority than any variable in use at this point. On startup, before any user input has taken place, 'x' is defined in this way and has initially maximal priority (and variable number 0). >>>

and do the right thing with the polynomials produced (however, have a look at the remark in Label se:rempolmod).

But in library mode, they are tricky objects if you intend to build polynomials yourself (and not just let PARI functions produce them, which is less efficient). For instance, it does not make sense to have a variable number occur in the components of a polynomial whose main variable has a lower priority, even though PARI cannot prevent you from doing it; see Label se:priority for a discussion of possible problems in a similar situation.

Creating variables

A basic difficulty is to ``create'' a variable. Some initializations are needed before you can use a given integer v as a variable number.

Initially, this is done for 0 (the variable x under gp), and MAXVARN, which is there to address the need for a ``temporary'' new variable in library mode and cannot be input under gp. No documented library function can create from scratch an object involving MAXVARN (of course, if the operands originally involve MAXVARN, the function abides). We call the latter type a ``temporary variable''. The regular variables meant to be used in regular objects, are called ``user variables''.

User variables When the program starts,

x is the only user variable (number 0). To define new ones, use

long fetch_user_var(char *s): inspects the user variable whose name is the string pointed to by s, creating it if needed, and returns its variable number.

  long v = fetch_user_var("y");
  GEN gy = pol_x(v);

The function raises an exception if the name is already in use for an installed or built-in function, or an alias.

@3Caveat. You can use gp_read_str (see Label se:gp_read_str) to execute a GP command and create GP variables on the fly as needed:

  GEN gy = gp_read_str("'y"); /* returns C<pol_x>(v), for some v */
  long v = varn(gy);

But please note the quote 'y in the above. Using gp_read_str("y") might work, but is dangerous, especially when programming functions to be used under gp. The latter reads the value of y, as currently known by the gp interpreter, possibly creating it in the process. But if y has been modified by previous gp commands (e.g. y = 1), then the value of gy is not what you expected it to be and corresponds instead to the current value of the gp variable (e.g. gen_1).

GEN fetch_var_value(long v) returns a shallow copy of the current value of the variable numbered v. Returns NULL if that variable number is unknown to the interpreter, e.g. it is a user variable. Note that this may not be the same as pol_x(v) if assignments have been performed in the interpreter.

Temporary variables

MAXVARN is available, but is better left to PARI internal functions (some of which do not check that MAXVARN is free for them to use, which can be considered a bug). You can create more temporary variables using

long fetch_var()

This returns a variable number which is guaranteed to be unused by the library at the time you get it and as long as you do not delete it (we will see how to do that shortly). This has higher priority than any temporary variable produced so far (MAXVARN is assumed to be the first such). After the statement v = fetch_var(), you can use pol_1(v) and pol_x(v). The variables created in this way have no identifier assigned to them though, and are printed as # < {number} > , except for MAXVARN which is printed as #. You can assign a name to a temporary variable, after creating it, by calling the function

void name_var(long n, char *s)

@3after which the output machinery will use the name s to represent the variable number n. The GP parser will not recognize it by that name, however, and calling this on a variable known to gp raises an error. Temporary variables are meant to be used as free variables, and you should never assign values or functions to them as you would do with variables under gp. For that, you need a user variable.

All objects created by fetch_var are on the heap and not on the stack, thus they are not subject to standard garbage collecting (they are not destroyed by a gerepile or avma = ltop statement). When you do not need a variable number anymore, you can delete it using

long delete_var()

@3which deletes the latest temporary variable created and returns the variable number of the previous one (or simply returns 0 if you try, in vain, to delete MAXVARN). Of course you should make sure that the deleted variable does not appear anywhere in the objects you use later on. Here is an example:

    long first = fetch_var();
    long n1 = fetch_var();
    long n2 = fetch_var(); /* prepare three variables for internal use */
    ...
    /* delete all variables before leaving */
    do { num = delete_var(); } while (num && num <= first);

The (dangerous) statement

    while (delete_var()) /* empty */;

removes all temporary variables in use, except MAXVARN which cannot be deleted.


Input and output

Two important aspects have not yet been explained which are specific to library mode: input and output of PARI objects.

Input

For input, PARI provides a powerful high level function which enables you to input your objects as if you were under gp. In fact, it is essentially the GP syntactical parser, hence you can use it not only for input but for (most) computations that you can do under gp. It has the following syntax:

GEN gp_read_str(const char *s)

Note that gp's metacommands are not recognized.

@3Note. The obsolete form

GEN readseq(char *t)

still exists for backward compatibility (assumes filtered input, without spaces or comments). Don't use it.

To read a GEN from a file, you can use the simpler interface

GEN gp_read_stream(FILE *file)

@3which reads a character string of arbitrary length from the stream file (up to the first complete expression sequence), applies gp_read_str to it, and returns the resulting GEN. This way, you do not have to worry about allocating buffers to hold the string. To interactively input an expression, use gp_read_stream(stdin).

Finally, you can read in a whole file, as in GP's read statement

GEN gp_read_file(char *name)

@3As usual, the return value is that of the last non-empty expression evaluated. There is one technical exception: if name is a binary file (from writebin) containing more than one object, a t_VEC containing them all is returned. This is because binary objects bypass the parser, hence reading them has no useful side effect.

Output to screen or file, output to string

General output functions return nothing but print a character string as a side effect. Low level routines are available to write on PARI output stream pari_outfile (stdout by default):

void pari_putc(char c): write character c to the output stream.

void pari_puts(char *s): write s to the output stream.

void pari_flush(): flush output stream; most streams are buffered by default, this command makes sure that all characters output so are actually written.

void pari_printf(const char *fmt, ...): the most versatile such function. fmt is a character string similar to the one printf uses. In there, % characters have a special meaning, and describe how to print the remaining operands. In addition to the standard format types (see the GP function printf), you can use the length modifier P (for PARI of course!) to specify that an argument is a GEN. For instance, the following are valid conversions for a GEN argument

      %Ps      convert to C<char*> (will print an arbitrary C<GEN>)
      %P.10s   convert to C<char*>, truncated to 10 chars
      %P.2f    convert to floating point format with 2 decimals
      %P4d     convert to integer, field width at least 4
      pari_printf("x[%d] = %Ps is not invertible!\n", i, gel(x,i));

@3 Here i is an int, x a GEN which is not a leaf (presumably a vector, or a polynomial) and this would insert the value of its i-th GEN component: gel(x,i).

@3Simple but useful variants to pari_printf are

void output(GEN x) prints x in raw format, followed by a newline and a buffer flush. This is more or less equivalent to

      pari_printf("%Ps\n", x);
      pari_flush();

void outmat(GEN x) as above except if x is a t_MAT, in which case a multi-line display is used to display the matrix. This is prettier for small dimensions, but quickly becomes unreadable and cannot be pasted and reused for input. If all entries of x are small integers, you may use the recursive features of %Pd and obtain the same (or better) effect with

      pari_printf("%Pd\n", x);
      pari_flush();

A variant like "%5Pd" would improve alignment by imposing 5 chars for each coefficient. Similarly if all entries are to be converted to floats, a format like "%5.1Pf" could be useful.

These functions write on (PARI's idea of) standard output, and must be used if you want your functions to interact nicely with gp. In most programs, this is not a concern and it is more flexible to write to an explicit FILE*, or to recover a character string:

void pari_fprintf(FILE *file, const char *fmt, ...) writes the remaining arguments to stream file according to the format specification fmt.

char* pari_sprintf(const char *fmt, ...) produces a string from the remaining arguments, according to the PARI format fmt (see printf). This is the libpari equivalent of Strprintf, and returns a malloc'ed string, which must be freed by the caller. Note that contrary to the analogous sprintf in the libc you do not provide a buffer (leading to all kinds of buffer overflow concerns); the function provided is actually closer to the GNU extension asprintf, although the latter has a different interface.

Simple variants of pari_sprintf convert a GEN to a malloc'ed ASCII string, which you must still free after use:

char* GENtostr(GEN x), using the current default output format (prettymat by default).

char* GENtoTeXstr(GEN x), suitable for inclusion in a TeX file.

Note that we have va_list analogs of the functions of printf type seen so far:

void pari_vprintf(const char *fmt, va_list ap)

void pari_vfprintf(FILE *file, const char *fmt, va_list ap)

char* pari_vsprintf(const char *fmt, va_list ap)

Errors

If you want your functions to issue error messages, you can use the general error handling routine pari_err. The basic syntax is

    pari_err(e_MISC, "error message");

This prints the corresponding error message and exit the program (in library mode; go back to the gp prompt otherwise). You can also use it in the more versatile guise

    pari_err(e_MISC, format, ...);

where format describes the format to use to write the remaining operands, as in the pari_printf function. For instance:

    pari_err(e_MISC, "x[%d] = %Ps is not invertible!", i, gel(x,i));

The simple syntax seen above is just a special case with a constant format and no remaining arguments. The general syntax is

void pari_err(numerr,...)

@3where numerr is a codeword which specifies the error class and what to do with the remaining arguments and what message to print. For instance, if x is a GEN with internal type t_STR, say, pari_err(e_TYPE,"extgcd", x) prints the message:

      ***   incorrect type in extgcd (t_STR),

@3See Label se:errors for details. In the libpari code itself, the general-purpose e_MISC is used sparingly: it is so flexible that the corresponding error contexts (t_ERROR) become hard to use reliably. Other more rigid error types are generally more useful: for instance the error context associated to the e_TYPE exception above is precisely documented and contains "extgcd" and x (not only its type) as readily available components.

Warnings

@3To issue a warning, use

void pari_warn(warnerr,...) In that case, of course, we do not abort the computation, just print the requested message and go on. The basic example is

      pari_warn(warner, "Strategy 1 failed. Trying strategy 2")

which is the exact equivalent of pari_err(e_MISC,...) except that you certainly do not want to stop the program at this point, just inform the user that something important has occurred; in particular, this output would be suitably highlighted under gp, whereas a simple printf would not.

The valid warning keywords are warner (general), warnprec (increasing precision), warnmem (garbage collecting) and warnfile (error in file operation), used as follows:

      pari_warn(warnprec, "bnfinit", newprec);
      pari_warn(warnmem,  "bnfinit");
      pari_warn(warnfile, "close", "afile");  /* error when closing "afile" */

Debugging output

For debugging output, you can use the standard output functions, output and pari_printf mainly. Corresponding to the gp metacommand \b x, you can also output the hexadecimal tree associated to an object:

void dbgGEN(GEN x, long nb = -1), displays the recursive structure of x. If nb = -1, the full structure is printed, otherwise the leaves (non-recursive components) are truncated to nb words.

@3The function output is vital under debuggers, since none of them knows how to print PARI objects by default. Seasoned PARI developers add the following gdb macro to their .gdbinit:

    define i
      call output((GEN)$arg0)
    end

Typing i x at a breakpoint in gdb then prints the value of the GEN x (provided the optimizer has not put it into a register, but it is rarely a good idea to debug optimized code).

The global variables DEBUGLEVEL and DEBUGMEM (corresponding to the default debug and debugmem, see Label se:defaults) are used throughout the PARI code to govern the amount of diagnostic and debugging output, depending on their values. You can use them to debug your own functions, especially if you install the latter under gp (see Label se:install).

void dbg_pari_heap(void) print debugging statements about the PARI stack, heap, and number of variables used. Corresponds to \s under gp.

Timers and timing output

@3To handle timings in a reentrant way, PARI defines a dedicated data type, pari_timer, together with the following methods:

void timer_start(pari_timer *T) start (or reset) a timer.

long timer_delay(pari_timer *T) returns the number of milliseconds elapsed since the timer was last reset. Resets the timer as a side effect.

long timer_get(pari_timer *T) returns the number of milliseconds elapsed since the timer was last reset. Does not reset the timer.

long timer_printf(pari_timer *T, char *format,...) This diagnostics function is equivalent to the following code

    err_printf("Time ")
    ... prints remaining arguments according to format ...
    err_printf(": %ld", timer_delay(T));

@3Resets the timer as a side effect.

@3They are used as follows:

    pari_timer T;
    timer_start(&T); /* initialize timer */
    ...
    printf("Total time: %ldms\n", timer_delay(&T));

or

    pari_timer T;
    timer_start(&T);
    for (i = 1; i < 10; i++) {
      ...
      timer_printf(&T, "for i = %ld (L[i] = %Ps)", i, gel(L,i));
    }

The following functions provided the same functionality, in a non-reentrant way, and are now deprecated.

long timer(void)

long timer2(void)

void msgtimer(const char *format, ...)

The following function implements gp's timer and should not be used in libpari programs: long gettime(void) equivalent to timer_delay(T) associated to a private timer T.


Iterators, Numerical integration, Sums, Products

Iterators

Since it is easier to program directly simple loops in library mode, some GP iterators are mainly useful for GP programming. Here are the others:

@3* fordiv is a trivial iteration over a list produced by divisors.

@3* forell and forsubgroup are currently not implemented as an iterator but as a procedure with callbacks.

void forell(void *E, long fun(void*, GEN), GEN a, GEN b) goes through the same curves as forell(ell,a,b,), calling fun(E, ell) for each curve ell, stopping if fun returns a non-zero value.

void forsubgroup(void *E, long fun(void*, GEN), GEN G, GEN B) goes through the same subgroups as forsubgroup(H = G, B,), calling fun(E, H) for each subgroup H, stopping if fun returns a non-zero value.

@3* forprime, for which we refer you to the next subsection.

@3* forcomposite, we provide an iterator over composite integers:

int forcomposite(forcomposite_t *T, GEN a, GEN b) initialize an iterator T over composite integers in [a,b]; over composites >= a if b = NULL. Return 0 if the range is known to be empty from the start (as if b < a or b < 0), and return 1 otherwise.

GEN forcomposite_next(forcomposite_t *T) returns the next composite in the range, assuming that T was initialized by forcomposite_init.

@3* forvec, for which we provide a convenient iterator. To initialize the analog of forvec(X = v, ..., flag), call

int forvec_init(forvec_t *T, GEN v, long flag) initialize an iterator T over the vectors generated by forvec(X = v,..., flag). This returns 0 if this vector list is empty, and 1 otherwise.

GEN forvec_next(forvec_t *T) returns the next element in the forvec sequence, or NULL if we are done. The return value must be used immediately or copied since the next call to the iterator destroys it: the relevant vector is updated in place. The iterator works hard to not use up PARI stack, and is more efficient when all lower bounds in the initialization vector v are integers. In that case, the cost is linear in the number of tuples enumerated, and you can expect to run over more than 10^9 tuples per minute. If speed is critical and all integers involved would fit in C longs, write a simple direct backtracking algorithm yourself.

@3* forpart is a variant of forvec which iterates over partitions. See the documentation of the forpart GP function for details. This function is available as a loop with callbacks:

void forpart(void *data, long (*call)(void*,GEN), long k, GEN a, GEN n)

@3It is also available as an iterator:

void forpart_init(forpart_t *T, long k, GEN a, GEN n) initializes an iterator over the partitions of k, with length restricted by n, and components restricted by a, either of which can be set to NULL to run without restriction.

GEN forpart_next(forpart_t *T) returns the next partition, or NULL when all partitions have been exhausted.

GEN forpart_prev(forpart_t *T)returns the previous partition, or NULL when all partitions have been exhausted.

You may not mix calls to forpart_next and forpart_prev: the first one called determines the ordering used to iterate over the partitions; you can not go back since the forpart_t structure is used in incompatible ways.

Iterating over primes

The library provides a high-level iterator, which stores its (private) data in a struct forprime_t and runs over arbitrary ranges of primes, without ever overflowing.

The iterator has two flavors, one providing the successive primes as ulongs, the other as GEN. They are initialized as follows, where we expect to run over primes >= a and <= b:

int forprime_init(forprime_t *T, GEN a, GEN b) for the GEN variant, where b = NULL means + oo .

int u_forprime_init(forprime_t *T, ulong a, ulong b) for the ulong variant, where b = ULONG_MAX means we will run through all primes representable in a ulong type.

Both variant return 1 on success, and 0 if the iterator would run over an empty interval (if a > b, for instance). They allocate the forprime_t data structure on the PARI stack.

@3The successive primes are then obtained using

GEN forprime_next(forprime_t *T), returns NULL if no more primes are available in the interval.

ulong u_forprime_next(forprime_t *T), returns 0 if no more primes are available in the interval.

These two functions leave alone the PARI stack, and write their state information in the preallocated forprime_t struct. The typical usage is thus:

    forprime_t T;
    GEN p;
    pari_sp av = avma, av2;
    forprime_init(&T, gen_2, stoi(1000));
    av2 = avma;
    while ( (p = forprime_next(&T)) )
    {
      ...
      if ( prime_is_OK(p) ) break;
      avma = av2; /* delete garbage accumulated in this iteration */
    }
    avma = av; /* delete all */

@3Of course, the final avma = av could be replaced by a gerepile call. Beware that swapping the av2 = avma and forprime_init call would be incorrect: the first avma = av2 would delete the forprime_t structure!

Numerical analysis

Numerical routines code a function (to be integrated, summed, zeroed, etc.) with two parameters named

    void *E;
    GEN (*eval)(void*, GEN)

The second is meant to contain all auxiliary data needed by your function. The first is such that eval(x, E) returns your function evaluated at x. For instance, one may code the family of functions f_t: x \to (x+t)^2 via

  GEN fun(void *t, GEN x) { return gsqr(gadd(x, (GEN)t)); }

One can then integrate f_1 between a and b with the call

  intnum((void*)stoi(1), &fun, a, b, NULL, prec);

Since you can set E to a pointer to any struct (typecast to void*) the above mechanism handles arbitrary functions. For simple functions without extra parameters, you may set E = NULL and ignore that argument in your function definition.


Catching exceptions

Basic use

PARI provides a mechanism to trap exceptions generated via pari_err using the pari_CATCH construction. The basic usage is as follows

   pari_CATCH(err_code) {
      recovery branch
   }
   pari_TRY {
      main branch
   }
   pari_ENDCATCH

@3This fragment executes the main branch, then the recovery branch if exception err_code is thrown, e.g. e_TYPE. See Label se:errors for the description of all error classes. The special error code CATCH_ALL is available to catch all errors.

One can replace the pari_TRY keyword by pari_RETRY, in which case once the recovery branch is run, we run the main branch again, still catching the same exceptions.

@3Restrictions.

@3* Such constructs can be nested without adverse effect, the innermost handler catching the exception.

@3* It is valid to leave either branch using pari_err.

@3* It is invalid to use C flow control instructions (break, continue, return) to directly leave either branch without seeing the pari_ENDCATCH keyword. This would leave an invalid structure in the exception handler stack, and the next exception would crash.

@3* In order to leave using break, continue or return, one must precede the keyword by a call to

void pari_CATCH_reset() disable the current handler, allowing to leave without adverse effect.

Advanced use

In the recovery branch, the exception context can be examined via the following helper routines:

GEN pari_err_last() returns the exception context, as a t_ERROR. The exception E returned by pari_err_last can be rethrown, using

    pari_err(0, E);

long err_get_num(GEN E) returns the error symbolic name. E.g e_TYPE.

GEN err_get_compo(GEN E, long i) error i-th component, as documented in Label se:errors.

@3For instance

   pari_CATCH(CATCH_ALL) { /* catch everything */
      GEN x, E = pari_err_last();
      long code = err_get_num(E);
      if (code != e_INV) pari_err(0, E); /* unexpected error, rethrow */
      x = err_get_compo(E, 2);
      /* e_INV has two components, 1: function name 2: non-invertible x */
      if (typ(x) != t_INTMOD) pari_err(0, E); /* unexpected type, rethrow */
      pari_CATCH_reset();
      return x; /* leave ! */
       ...
   } pari_TRY {
      main branch
   }
   pari_ENDCATCH


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 study a detailed example. We want to write a program which computes the gcd of two integers, together with the Bezout coefficients. We shall use the standard quadratic algorithm which is not optimal but is not too far from the one used in the PARI function bezout.

Let x,y two integers and initially \pmatrix{s_x s_y t_x t_y } = \pmatrix{1 0 0 1}, so that

   \pmatrix{s_x s_y t_x t_y } \pmatrix{x y } = \pmatrix{x y }.

To apply the ordinary Euclidean algorithm to the right hand side, multiply the system from the left by \pmatrix{0 1 1 -q }, with q = floor(x / y). Iterate until y = 0 in the right hand side, then the first line of the system reads

   s_x x + s_y y = gcd(x,y).

In practice, there is no need to update s_y and t_y since gcd(x,y) and s_x are enough to recover s_y. The following program is now straightforward. A couple of new functions appear in there, whose description can be found in the technical reference manual in Chapter 5, but whose meaning should be clear from their name and the context.

This program can be found in examples/extgcd.c together with a proper Makefile. You may ignore the first comment

  /*
  GP;install("extgcd", "GG&&", "gcdex", "./libextgcd.so");
  */

@3which instruments the program so that gp2c-run extgcd.c can import the extgcd() routine into an instance of the gp interpreter (under the name gcdex). See the gp2c manual for details. \newpage

file{../examples/extgcd.c}

@3For simplicity, the inner loop does not include any garbage collection, hence memory use is quadratic in the size of the inputs instead of linear. Here is a better version of that loop:

    pari_sp av = avma, lim = stack_lim(av,1);
    ...
    while (!gequal0(b))
    {
      GEN r, q = dvmdii(a, b, &r), v = vx;
      vx = subii(ux, mulii(q, vx));
      ux = v; a = b; b = r;
      if (low_stack(lim, stack_lim(av,1)))
        gerepileall(av, 4, &a, &b, &ux, &vx);
    }

\newpage