Testing latest pari + WASM + node.js... and it works?! Wow.
License: GPL3
ubuntu2004
% Copyright (c) 2000 The PARI Group1%2% This file is part of the PARI/GP documentation3%4% Permission is granted to copy, distribute and/or modify this document5% under the terms of the GNU General Public License6\chapter{Programming PARI in Library Mode}78\noindent The \emph{User's Guide to Pari/GP} gives in three chapters a9general presentation of the system, of the \kbd{gp} calculator, and detailed10explanation of high level PARI routines available through the calculator. The11present manual assumes general familiarity with the contents of these12chapters and the basics of ANSI C programming, and focuses on the usage of13the PARI library. In this chapter, we introduce the general concepts of PARI14programming and describe useful general purpose functions; the following15chapters describes all public low or high-level functions, underlying or16extending the GP functions seen in Chapter 3 of the User's guide.1718\section{Introduction: initializations, universal objects}19\label{se:intro4}2021\noindent22To use PARI in \idx{library mode}, you must write a C program and link it to23the PARI library. See the installation guide or the Appendix to the24\emph{User's Guide to Pari/GP} on how to create and install the library and25include files. A sample Makefile is presented in Appendix~A, and a more26elaborate one in \kbd{examples/Makefile}. The best way to understand how27programming is done is to work through a complete example. We will write such28a program in~\secref{se:prog}. Before doing this, a few explanations are in29order.3031First, one must explain to the outside world what kind of objects and32routines we are going to use. This is done\footnote{*}{This assumes that PARI33headers are installed in a directory which belongs to your compiler's search34path for header files. You might need to add flags like35\kbd{-I/usr/local/include} or modify \tet{C_INCLUDE_PATH}.}36with the directive3738\bprog39#include <pari/pari.h>40@eprog41\noindent42In particular, this defines the fundamental type for all PARI objects: the43type \teb{GEN}, which is simply a pointer to \kbd{long}.4445Before any PARI routine is called, one must initialize the system, and in46particular the PARI stack which is both a scratchboard and a repository for47computed objects. This is done with a call to the function4849\fun{void}{pari_init}{size_t size, ulong maxprime}5051\noindent The first argument is the number of bytes given to PARI to work52with, and the second is the upper limit on a precomputed prime number table;53\kbd{size} should not reasonably be taken below $500000$ but you may set54$\tet{maxprime} = 0$, although the system still needs to precompute all55primes up to about $2^{16}$. For lower-level variants allowing finer56control, e.g.~preventing PARI from installing its own error or signal57handlers, see~\secref{se:pari_init_tech}.5859\noindent We have now at our disposal:6061\item a PARI \tev{stack} containing nothing. This is a big62connected chunk of \kbd{size} bytes of memory, where all computations63take place. In large computations, intermediate results quickly64clutter up memory so some kind of garbage collecting is needed. Most65systems do garbage collecting when the memory is getting scarce, and this66slows down the performance. PARI takes a different approach, admittedly more67demanding on the programmer: you must do your own cleaning up when the68intermediate results are not needed anymore. We will see later how (and when)69this is done.7071\item the following \emph{universal objects} (by definition, objects72which do not belong to the stack): the integers $0$, $1$, $-1$, $2$ and73$-2$ (respectively called \tet{gen_0}, \tet{gen_1}, \tet{gen_m1},74\tet{gen_2} and \tet{gen_m2}), the fraction $\dfrac{1}{2}$ (\tet{ghalf}).75All of these are of type \kbd{GEN}.7677\item a \tev{heap} which is just a linked list of permanent78universal objects. For now, it contains exactly the ones listed above. You79will probably very rarely use the heap yourself; and if so, only as a80collection of copies of objects taken from the stack (called \idx{clone}s in81the sequel). Thus you need not bother with its internal structure, which may82change as PARI evolves. Some complex PARI functions create clones for special83garbage collecting purposes, usually destroying them when returning.8485\item a table of primes (in fact of \emph{differences} between86consecutive primes), called \teb{diffptr}, of type \kbd{byteptr}87(pointer to \kbd{unsigned char}). Its use is described in88\secref{se:primetable} later. Using it directly is deprecated,89high-level iterators provide a cleaner and more flexible interface, see90\secref{se:primeiter} (such iterators use the private prime table, but extend91it dynamically).9293\item access to all the built-in functions of the PARI library.94These are declared to the outside world when you include \kbd{pari.h}, but95need the above things to function properly. So if you forget the call to96\tet{pari_init}, you will get a fatal error when running your program.9798\section{Important technical notes}99100\subsec{Backward compatibility} The PARI function names evolved over time,101and deprecated functions are eventually deleted. The file \kbd{pariold.h}102contains macros implementing a weak form of backward compatibility.103In particular, whenever the name of a documented function changes, a104\kbd{\#define} is added to this file so that the old name expands to the new105one (provided the prototype didn't change also).106107This file is included by \kbd{pari.h}, but a large section is commented out108by default. Define \tet{PARI_OLD_NAMES} before including \kbd{pari.h} to109pollute your namespace with lots of obsolete names like110\kbd{un}\footnote{*}{For \kbd{(long)gen\_1}. Since 2004 and version 2.2.9,111typecasts are completely unnecessary in PARI programs.}: that might enable112you to compile old programs without having to modify them. The preferred way113to do that is to add \kbd{-DPARI\_OLD\_NAMES} to your compiler \kbd{CFLAGS},114so that you don't need to modify the program files themselves.115116Of course, it's better to fix the program if you can!117118\subsec{Types}119120\noindent121Although PARI objects all have the C type \kbd{GEN}, we will freely use122the word \teb{type} to refer to PARI dynamic subtypes: \typ{INT}, \typ{REAL},123etc. The declaration124\bprog125GEN x;126@eprog\noindent127declares a C variable of type \kbd{GEN}, but its ``value'' will be said to128have type \typ{INT}, \typ{REAL}, etc. The meaning should always be clear from129the context.130131\subsec{Type recursivity}132133\noindent134Conceptually, most PARI types are recursive. But the \kbd{GEN} type is a135pointer to \kbd{long}, not to \kbd{GEN}. So special macros must be used to136access \kbd{GEN}'s components. The simplest one is \tet{gel}$(V, i)$, where137\key{el} stands for \key{el}ement, to access component number $i$ of the138\kbd{GEN} $V$. This is a valid \kbd{lvalue} (may be put on the left side of139an assignment), and the following two constructions are exceedingly frequent140%141\bprog142gel(V, i) = x;143x = gel(V, i);144@eprog\noindent145where \kbd{x} and \kbd{V} are \kbd{GEN}s. This macro accesses and modifies146directly the components of $V$ and do not create a copy of the coefficient,147contrary to all the library \emph{functions}.148149More generally, to retrieve the values of elements of lists of \dots\ of150lists of vectors we have the \tet{gmael} macros (for {\bf m}ultidimensional151{\bf a}rray {\bf el}ement). The syntax is $\key{gmael}n(V,a_1,\dots,a_n)$,152where $V$ is a \kbd{GEN}, the $a_i$ are indexes, and $n$ is an integer153between $1$ and $5$. This stands for $x[a_1][a_2]\dots[a_n]$, and returns a154\kbd{GEN}. The macros \tet{gel} (resp.~\tet{gmael}) are synonyms for155\tet{gmael1} (resp.~\kbd{gmael2}).156157Finally, the macro $\tet{gcoeff}(M, i, j)$ has exactly the meaning of158\kbd{M[i,j]} in GP when \kbd{M} is a matrix. Note that due to the159implementation of \typ{MAT}s as horizontal lists of vertical vectors,160\kbd{gcoeff(x,y)} is actually equivalent to \kbd{gmael(y,x)}. One should use161\kbd{gcoeff} in matrix context, and \kbd{gmael} otherwise.162163\subsec{Variations on basic functions}\label{se:low_level} In the library164syntax descriptions in Chapter~3, we have only given the basic names of the165functions. For example \kbd{gadd}$(x,y)$ assumes that $x$ and $y$ are166\kbd{GEN}s, and \emph{creates} the result $x+y$ on the PARI stack. For most167of the basic operators and functions, many other variants are available. We168give some examples for \kbd{gadd}, but the same is true for all the basic169operators, as well as for some simple common functions (a complete list170is given in Chapter~6):171172\fun{GEN}{gaddgs}{GEN x, long y}173174\fun{GEN}{gaddsg}{long x, GEN y}175176\noindent In the following one, \kbd{z} is a preexisting \kbd{GEN} and the177result of the corresponding operation is put into~\kbd{z}. The size of the PARI178stack does not change:179180\fun{void}{gaddz}{GEN x, GEN y, GEN z}181182\noindent (This last form is inefficient in general and deprecated outside of183PARI kernel programming.) Low level kernel functions implement these184operators for specialized arguments and are also available: Level 0 deals185with operations at the word level (\kbd{long}s and \kbd{ulong}s), Level 1186with \typ{INT} and \typ{REAL} and Level 2 with the rest (modular arithmetic,187polynomial arithmetic and linear algebra). Here are some examples of Level188$1$ functions:189190\fun{GEN}{addii}{GEN x, GEN y}: here $x$ and $y$ are \kbd{GEN}s of type191\typ{INT} (this is not checked).192193\fun{GEN}{addrr}{GEN x, GEN y}: here $x$ and $y$ are \kbd{GEN}s of194type \typ{REAL} (this is not checked).195196\noindent197There also exist functions \teb{addir}, \teb{addri}, \teb{mpadd} (whose198two arguments can be of type \typ{INT} or \typ{REAL}), \teb{addis} (to add a199\typ{INT} and a \kbd{long}) and so on.200201The Level $1$ names are self-explanatory once you know that {\bf i} stands for a202\typ{INT}, {\bf r} for a \typ{REAL}, {\bf mp} for i or r, {\bf s} for a signed C203long integer, {\bf u} for an unsigned C long integer; finally the suffix {\bf z}204means that the result is not created on the PARI stack but assigned to a205preexisting GEN object passed as an extra argument. Chapter 6 gives a206description of these low-level functions.207208Level $2$ names are more complicated, see \secref{se:level2names} for all the209gory details, and we content ourselves with a simple example used to implement210\typ{INTMOD} arithmetic:211212\fun{GEN}{Fp_add}{GEN x, GEN y, GEN m}: returns the sum of $x$ and $y$ modulo213$m$. Here $x, y, m$ are \typ{INT}s (this is not checked). The operation is214more efficient if the inputs $x$, $y$ are reduced modulo $m$, but this is not215a necessary condition.216217\misctitle{Important Note} These specialized functions are of course more218efficient than the generic ones, but note the hidden danger here: the types219of the objects involved (which is not checked) must be severely controlled,220e.g.~using \kbd{addii} on a \typ{FRAC} argument will cause disasters. Type221mismatches may corrupt the PARI stack, though in most cases they will just222immediately overflow the stack. Because of this, the PARI philosophy of223giving a result which is as exact as possible, enforced for generic functions224like \kbd{gadd} or \kbd{gmul}, is dropped in kernel routines of Level $1$,225where it is replaced by the much simpler rule: the result is a \typ{INT} if226and only if all arguments are integer types (\typ{INT} but also C \kbd{long}227and \kbd{ulong}) and a \typ{REAL} otherwise. For instance, multiplying a228\typ{REAL} by a \typ{INT} always yields a \typ{REAL} if you use \kbd{mulir},229where \kbd{gmul} returns the \typ{INT} \kbd{gen\_0} if the integer is $0$.230231\subsec{Portability: 32-bit / 64-bit architectures}232233\noindent234PARI supports both 32-bit and 64-bit based machines, but not simultaneously!235The library is compiled assuming a given architecture, and some236of the header files you include (through \kbd{pari.h}) will have been237modified to match the library.238239Portable macros are defined to bypass most machine dependencies. If you want240your programs to run identically on 32-bit and 64-bit machines, you have to241use these, and not the corresponding numeric values, whenever the precise242size of your \kbd{long} integers might matter. Here are the most important243ones:244245\settabs\+ -----------------------------&---------------&------------&\cr \+246& 64-bit & 32-bit \cr\+ \tet{BITS_IN_LONG} & 64 & 32 \cr\+247\tet{LONG_IS_64BIT} & defined & undefined \cr\+248\tet{DEFAULTPREC} & 3 & 4 & ($\approx$ 19 decimal digits, %249see formula below) \cr\+250\tet{MEDDEFAULTPREC}& 4 & 6 & ($\approx$ 38 decimal digits) \cr\+251\tet{BIGDEFAULTPREC}& 5 & 8 & ($\approx$ 57 decimal digits) \cr252\noindent For instance, suppose you call a transcendental function, such as253254\kbd{GEN \key{gexp}(GEN x, long prec)}.255256\noindent The last argument \kbd{prec} is an integer $\geq 3$, corresponding257to the default floating point precision required. It is \emph{only} used if258\kbd{x} is an exact object, otherwise the relative precision is determined by259the precision of~\kbd{x}. Since the parameter \kbd{prec} sets the size of the260inexact result counted in (\kbd{long}) \emph{words} (including codewords),261the same value of \kbd{prec} will yield different results on 32-bit and26264-bit machines. Real numbers have two codewords (see~\secref{se:impl}), so263the formula for computing the bit accuracy is264$$ \tet{bit_accuracy}(\kbd{prec}) = (\kbd{prec} - 2) * \tet{BITS_IN_LONG}$$265(this is actually the definition of an inline function). The corresponding266accuracy expressed in decimal digits would be267%268$$ \kbd{bit\_accuracy(prec)} * \log(2) / \log(10).$$269%270For example if the value of \kbd{prec} is 5, the corresponding accuracy for27132-bit machines is $(5-2)*\log(2^{32})/\log(10)\approx 28$ decimal digits,272while for 64-bit machines it is $(5-2)*\log(2^{64})/\log(10)\approx 57$273decimal digits.274275Thus, you must take care to change the \kbd{prec} parameter you are supplying276according to the bit size, either using the default precisions given by the277various \kbd{DEFAULTPREC}s, or by using conditional constructs of the form:278%279\bprog280#ifndef LONG_IS_64BIT281prec = 4;282#else283prec = 6;284#endif285@eprog286\noindent which is in this case equivalent to the statement287\kbd{prec = MEDDEFAULTPREC;}.288289Note that for parity reasons, half the accuracies available on 32-bit290architectures (the odd ones) have no precise equivalents on 64-bit machines.291292\subsec{Using \kbd{malloc} / \kbd{free}}293You should make use of the PARI stack as much as possible, and avoid294allocating objects using the customary functions. If you do, you should295use, or at least have a very close look at, the following wrappers:296297\fun{void*}{pari_malloc}{size_t size} calls \kbd{malloc} to allocate298\kbd{size} bytes and returns a pointer to the allocated memory. If the299request fails, an error is raised. The \kbd{SIGINT} signal is blocked until300\kbd{malloc} returns, to avoid leaving the system stack in an inconsistent301state.302303\fun{void*}{pari_realloc}{void* ptr, size_t size} as \tet{pari_malloc} but304calls \kbd{realloc} instead of \kbd{malloc}.305306\fun{void}{pari_realloc_ip}{void** ptr, size_t size}307equivalent to \kbd{*ptr= realloc(*ptr, size)}, while blocking \kbd{SIGINT}.308309\fun{void*}{pari_calloc}{size_t size} as \tet{pari_malloc}, setting the310memory to zero.311312\fun{void}{pari_free}{void* ptr} calls \kbd{free} to liberate the memory313space pointed to by \kbd{ptr}, which must have been allocated by \kbd{malloc}314(\tet{pari_malloc}) or \kbd{realloc} (\tet{pari_realloc}). The \kbd{SIGINT}315signal is blocked until \kbd{free} returns.316317If you use the standard \kbd{libc} functions instead of our wrappers, then318your functions will be subtly incompatible with the \kbd{gp} calculator: when319the user tries to interrupt a computation, the calculator may crash320(if a system call is interrupted at the wrong time).321322\section{Garbage collection}\label{se:garbage}\sidx{garbage collecting}323324\subsec{Why and how}325326\noindent327As we have seen, \kbd{pari\_init} allocates a big range of328addresses, the \tev{stack}, that are going to be used throughout. Recall329that all PARI objects are pointers. Except for a few universal objects,330they all point at some part of the stack.331332The stack starts at the address \kbd{bot} and ends just before \kbd{top}. This333means that the quantity334%335$$ (\kbd{top} - \kbd{bot})\,/\,\kbd{sizeof(long)} $$336%337is (roughly) equal to the \kbd{size} argument of \kbd{pari\_init}. The PARI338stack also has a ``current stack pointer'' called \teb{avma}, which stands339for {\bf av}ailable {\bf m}emory {\bf a}ddress. These three variables are340global (declared by \kbd{pari.h}). They are of type \tet{pari_sp}, which341means \emph{pari stack pointer}.342343The stack is oriented upside-down: the more recent an object, the closer to344\kbd{bot}. Accordingly, initially \kbd{avma} = \kbd{top}, and \kbd{avma} gets345\emph{decremented} as new objects are created. As its name indicates,346\kbd{avma} always points just \emph{after} the first free address on the347stack, and \kbd{(GEN)avma} is always (a pointer to) the latest created object.348When \kbd{avma} reaches \kbd{bot}, the stack overflows, aborting all349computations, and an error message is issued. To avoid this \emph{you}350need to clean up the stack from time to time, when intermediate objects are351not needed anymore. This is called ``\emph{garbage collecting}.''352353We are now going to describe briefly how this is done. We will see many354concrete examples in the next subsection.355356\noindent\item357First, PARI routines do their own garbage collecting, which means that358whenever a documented function from the library returns, only its result(s)359have been added to the stack, possibly up to a very small overhead360(undocumented ones may not do this). In361particular, a PARI function that does not return a \kbd{GEN} does not clutter362the stack. Thus, if your computation is small enough (e.g.~you call few PARI363routines, or most of them return \kbd{long} integers), then you do not need364to do any garbage collecting. This is probably the case in many of your365subroutines. Of course the objects that were on the stack \emph{before} the366function call are left alone. Except for the ones listed below, PARI367functions only collect their own garbage.368369\noindent\item370It may happen that all objects that were created after a certain point can371be deleted~--- for instance, if the final result you need is not a372\kbd{GEN}, or if some search proved futile. Then, it is enough to record373the value of \kbd{avma} just \emph{before} the first garbage is created,374and restore it upon exit:375376\bprog377pari_sp av = avma; /*@Ccom record initial avma */378379garbage ...380avma = av; /*@Ccom restore it */381@eprog382\noindent All objects created in the \kbd{garbage} zone will eventually383be overwritten: they should no longer be accessed after \kbd{avma} has been384restored.385386\noindent\item387If you want to destroy (i.e.~give back the memory occupied by) the388\emph{latest} PARI object on the stack (e.g.~the latest one obtained from a389function call), you can use the function\sidx{destruction}%390\vadjust{\penalty500}%discourage page break391392\fun{void}{cgiv}{GEN z}393394\noindent where \kbd{z} is the object you want to give back. This is395equivalent to the above where the initial \kbd{av} is computed from \kbd{z}.396397\noindent\item398Unfortunately life is not so simple, and sometimes you will want399to give back accumulated garbage \emph{during} a computation without losing400recent data. We shall start with the lowest level function to get a feel for401the underlying mechanisms, we shall describe simpler variants later:402403\fun{GEN}{gerepile}{pari_sp ltop, pari_sp lbot, GEN q}. This function cleans404up the stack between \kbd{ltop} and \kbd{lbot}, where $\kbd{lbot} <405\kbd{ltop}$, and returns the updated object \kbd{q}. This means:4064071) we translate (copy) all the objects in the interval408$[\kbd{avma}, \kbd{lbot}[$, so that its right extremity abuts the address409\kbd{ltop}. Graphically410411\vbox{\bprog412bot avma lbot ltop top413End of stack |-------------[++++++[-/-/-/-/-/-/-|++++++++| Start414free memory garbage415@eprog416\noindent becomes:417\bprog418bot avma ltop top419End of stack |---------------------------[++++++[++++++++| Start420free memory421@eprog422}423\noindent where \kbd{++} denote significant objects, \kbd{--} the unused part424of the stack, and \kbd{-/-} the garbage we remove.4254262) The function then inspects all the PARI objects between \kbd{avma} and427\kbd{lbot} (i.e.~the ones that we want to keep and that have been translated)428and looks at every component of such an object which is not a codeword. Each429such component is a pointer to an object whose address is either430431--- between \kbd{avma} and \kbd{lbot}, in which case it is suitably updated,432433--- larger than or equal to \kbd{ltop}, in which case it does not change, or434435--- between \kbd{lbot} and \kbd{ltop} in which case \kbd{gerepile}436raises an error (``significant pointers lost in gerepile'').4374383) \key{avma} is updated (we add $\kbd{ltop} - \kbd{lbot}$ to the old value).4394404) We return the (possibly updated) object \kbd{q}: if \kbd{q} initially441pointed between \kbd{avma} and \kbd{lbot}, we return the updated address, as442in~2). If not, the original address is still valid, and is returned!443444As stated above, no component of the remaining objects (in particular445\kbd{q}) should belong to the erased segment [\kbd{lbot}, \kbd{ltop}[, and446this is checked within \kbd{gerepile}. But beware as well that the addresses447of the objects in the translated zone change after a call to \kbd{gerepile},448so you must not access any pointer which previously pointed into the zone449below \kbd{ltop}. If you need to recover more than one object, use the450\kbd{gerepileall} function below.451452\misctitle{Remark}453As a consequence of the preceding explanation, if a PARI object is to be454relocated by \hbox{gerepile} then, apart from universal objects, the chunks455of memory used by its components should be in consecutive memory locations.456All \kbd{GEN}s created by documented PARI functions are guaranteed to satisfy457this. This is because the \kbd{gerepile} function knows only about \emph{two458connected zones}: the garbage that is erased (between \kbd{lbot} and459\kbd{ltop}) and the significant pointers that are copied and updated. If460there is garbage interspersed with your objects, disaster occurs when we try461to update them and consider the corresponding ``pointers''. In most cases of462course the said garbage is in fact a bunch of other \kbd{GEN}s, in which case463we simply waste time copying and updating them for nothing. But be wary when464you allow objects to become disconnected.465466\noindent In practice this is achieved by the following programming idiom:467\bprog468ltop = avma; garbage(); lbot = avma; q = anything();469return gerepile(ltop, lbot, q); /*@Ccom returns the updated q */470@eprog\noindent or directly471\bprog472ltop = avma; garbage(); lbot = avma;473return gerepile(ltop, lbot, anything());474@eprog\noindent475Beware that476\bprog477ltop = avma; garbage();478return gerepile(ltop, avma, anything())479@eprog480481\noindent might work, but should be frowned upon. We cannot predict whether482\kbd{avma} is evaluated after or before the call to \kbd{anything()}: it483depends on the compiler. If we are out of luck, it is \emph{after} the484call, so the result belongs to the garbage zone and the \kbd{gerepile}485statement becomes equivalent to \kbd{avma = ltop}. Thus we return a486pointer to random garbage.487488\subsec{Variants}489490\fun{GEN}{gerepileupto}{pari_sp ltop, GEN q}. Cleans the stack between491\kbd{ltop} and the \emph{connected} object \kbd{q} and returns \kbd{q}492updated. For this to work, \kbd{q} must have been created \emph{before} all493its components, otherwise they would belong to the garbage zone! Unless494mentioned otherwise, documented PARI functions guarantee this.495496\fun{GEN}{gerepilecopy}{pari_sp ltop, GEN x}. Functionally equivalent to,497but more efficient than498\bprog499gerepileupto(ltop, gcopy(x))500@eprog\noindent In this case, the \kbd{GEN} parameter \kbd{x} need not501satisfy any property before the garbage collection: it may be disconnected,502components created before the root, and so on. Of course, this is about503twice slower than either \tet{gerepileupto} or \tet{gerepile}, because504\kbd{x} has to be copied to a clean stack zone first. This function is a505special case of \tet{gerepileall} below, where $n=1$.506507\fun{void}{gerepileall}{pari_sp ltop, int n, ...}.508To cope with complicated cases where many objects have to be preserved. The509routine expects $n$ further arguments, which are the \emph{addresses} of510the \kbd{GEN}s you want to preserve:511\bprog512pari_sp ltop = avma;513...; y = ...; ... x = ...; ...;514gerepileall(ltop, 2, &x, &y);515@eprog\noindent516It cleans up the most recent part of the517stack (between \kbd{ltop} and \kbd{avma}), updating all the \kbd{GEN}s added518to the argument list. A copy is done just before the cleaning to preserve519them, so they do not need to be connected before the call. With520\kbd{gerepilecopy}, this is the most robust of the \kbd{gerepile} functions521(the less prone to user error), hence the slowest.522523\fun{void}{gerepileallsp}{pari_sp ltop, pari_sp lbot, int n, ...}.524More efficient, but trickier than \kbd{gerepileall}. Cleans the stack between525\kbd{lbot} and \kbd{ltop} and updates the \kbd{GEN}s pointed at by the526elements of \kbd{gptr} without any further copying. This is subject to the527same restrictions as \kbd{gerepile}, the only difference being that more than528one address gets updated.529530\subsec{Examples}531532\subsubsec{gerepile}533534Let \kbd{x} and \kbd{y} be two preexisting PARI objects and suppose that we535want to compute $\kbd{x}^2 + \kbd{y}^2$. This is done using the following536program:537\bprog538GEN x2 = gsqr(x);539GEN y2 = gsqr(y), z = gadd(x2,y2);540@eprog\noindent541The \kbd{GEN} \kbd{z} indeed points at the desired quantity. However,542consider the stack: it contains as unnecessary garbage \kbd{x2} and \kbd{y2}.543More precisely it contains (in this order) \kbd{z}, \kbd{y2}, \kbd{x2}.544(Recall that, since the stack grows downward from the top, the most recent545object comes first.)546547It is not possible to get rid of \kbd{x2}, \kbd{y2} before \kbd{z} is548computed, since they are used in the final operation. We cannot record549\kbd{avma} before \kbd{x2} is computed and restore it later, since this would550destroy \kbd{z} as well. It is not possible either to use the function551\kbd{cgiv} since \kbd{x2} and \kbd{y2} are not at the bottom of the stack and552we do not want to give back~\kbd{z}.553554But using \kbd{gerepile}, we can give back the memory locations corresponding555to \kbd{x2}, \kbd{y2}, and move the object \kbd{z} upwards so that no556space is lost. Specifically:557\bprog558pari_sp ltop = avma; /*@Ccom remember the current top of the stack */559GEN x2 = gsqr(x);560GEN y2 = gsqr(y);561pari_sp lbot = avma; /*@Ccom the bottom of the garbage pile */562GEN z = gadd(x2, y2); /*@Ccom z is now the last object on the stack */563z = gerepile(ltop, lbot, z);564@eprog565\noindent Of course, the last two instructions could also have been566written more simply:567\bprog568z = gerepile(ltop, lbot, gadd(x2,y2));569@eprog\noindent In fact \kbd{gerepileupto} is even simpler to use, because570the result of \kbd{gadd} is the last object on the stack and \kbd{gadd}571is guaranteed to return an object suitable for \kbd{gerepileupto}:572\bprog573ltop = avma;574z = gerepileupto(ltop, gadd(gsqr(x), gsqr(y)));575@eprog\noindent576Make sure you understand exactly what has happened before you go on!577578\misctitle{Remark on assignments and gerepile} When the tree structure and579the size of the PARI objects which will appear in a computation are under580control, one may allocate sufficiently large objects at the beginning,581use assignment statements, then simply restore \kbd{avma}. Coming back to the582above example, note that \emph{if} we know that x and y are of type real583fitting into \kbd{DEFAULTPREC} words, we can program without using584\kbd{gerepile} at all:585\bprog586z = cgetr(DEFAULTPREC); ltop = avma;587gaffect(gadd(gsqr(x), gsqr(y)), z);588avma = ltop;589@eprog\noindent This is often \emph{slower} than a craftily used590\kbd{gerepile} though, and certainly more cumbersome to use. As a rule,591assignment statements should generally be avoided.592593\misctitle{Variations on a theme} it is often necessary to do several594\kbd{gerepile}s during a computation. However, the fewer the better. The only595condition for \kbd{gerepile} to work is that the garbage be connected. If the596computation can be arranged so that there is a minimal number of connected597pieces of garbage, then it should be done that way.598599For example suppose we want to write a function of two \kbd{GEN} variables600\kbd{x} and \kbd{y} which creates the vector $\kbd{[x}^2+\kbd{y},601\kbd{y}^2+\kbd{x]}$. Without garbage collecting, one would write:602%603\bprog604p1 = gsqr(x); p2 = gadd(p1, y);605p3 = gsqr(y); p4 = gadd(p3, x);606z = mkvec2(p2, p4); /* not suitable for gerepileupto! */607@eprog\noindent608This leaves a dirty stack containing (in this order) \kbd{z}, \kbd{p4},609\kbd{p3}, \kbd{p2}, \kbd{p1}. The garbage here consists of \kbd{p1} and610\kbd{p3}, which are separated by \kbd{p2}. But if we compute \kbd{p3}611\emph{before} \kbd{p2} then the garbage becomes connected, and we get the612following program with garbage collecting:613%614\bprog615ltop = avma; p1 = gsqr(x); p3 = gsqr(y);616lbot = avma; z = cgetg(3, t_VEC);617gel(z, 1) = gadd(p1,y);618gel(z, 2) = gadd(p3,x); z = gerepile(ltop,lbot,z);619@eprog\noindent Finishing by \kbd{z = gerepileupto(ltop, z)} would be ok as620well. Beware that621\bprog622ltop = avma; p1 = gadd(gsqr(x), y); p3 = gadd(gsqr(y), x);623z = cgetg(3, t_VEC);624gel(z, 1) = p1;625gel(z, 2) = p3; z = gerepileupto(ltop,z); /*@Ccom WRONG */626@eprog\noindent627is a disaster since \kbd{p1} and \kbd{p3} are created before628\kbd{z}, so the call to \kbd{gerepileupto} overwrites them, leaving629\kbd{gel(z, 1)} and \kbd{gel(z, 2)} pointing at random data! The following630does work:631\bprog632ltop = avma; p1 = gsqr(x); p3 = gsqr(y);633lbot = avma; z = mkvec2(gadd(p1,y), gadd(p3,x));634z = gerepile(ltop,lbot,z);635@eprog\noindent but is very subtly wrong in the sense that636\kbd{z = gerepileupto(ltop, z)} would \emph{not} work. The reason being637that \kbd{mkvec2} creates the root \kbd{z} of the vector \emph{after}638its arguments have been evaluated, creating the components of \kbd{z}639too early; \kbd{gerepile} does not care, but the created \kbd{z} is a time640bomb which will explode on any later \kbd{gerepileupto}.641On the other hand642\bprog643ltop = avma; z = cgetg(3, t_VEC);644gel(z, 1) = gadd(gsqr(x), y);645gel(z, 2) = gadd(gsqr(y), x); z = gerepileupto(ltop,z); /*@Ccom INEFFICIENT */646@eprog\noindent647leaves the results of \kbd{gsqr(x)} and \kbd{gsqr(y)} on the stack (and648lets \kbd{gerepileupto} update them for naught). Finally, the most elegant649and efficient version (with respect to time and memory use) is as follows650\bprog651z = cgetg(3, t_VEC);652ltop = avma; gel(z, 1) = gerepileupto(ltop, gadd(gsqr(x), y));653ltop = avma; gel(z, 2) = gerepileupto(ltop, gadd(gsqr(y), x));654@eprog\noindent655which avoids updating the container \kbd{z} and cleans up its components656individually, as soon as they are computed.657658\misctitle{One last example} Let us compute the product of two complex659numbers $x$ and $y$, using the $3M$ method which requires 3 multiplications660instead of the obvious 4. Let $z = x*y$, and set $x = x_r + i*x_i$ and661similarly for $y$ and $z$. We compute $p_1 = x_r*y_r$, $p_2=x_i*y_i$,662$p_3=(x_r+x_i)*(y_r+y_i)$, and then we have $z_r=p_1-p_2$,663$z_i=p_3-(p_1+p_2)$. The program is as follows:664%665\bprog666ltop = avma;667p1 = gmul(gel(x,1), gel(y,1));668p2 = gmul(gel(x,2), gel(y,2));669p3 = gmul(gadd(gel(x,1), gel(x,2)), gadd(gel(y,1), gel(y,2)));670p4 = gadd(p1,p2);671lbot = avma; z = cgetg(3, t_COMPLEX);672gel(z, 1) = gsub(p1,p2);673gel(z, 2) = gsub(p3,p4); z = gerepile(ltop,lbot,z);674@eprog675676\misctitle{Exercise} Write a function which multiplies a matrix by a column677vector. Hint: start with a \kbd{cgetg} of the result, and use \kbd{gerepile}678whenever a coefficient of the result vector is computed. You can look at the679answer in \kbd{src/basemath/RgV.c:RgM\_RgC\_mul()}.680681\subsubsec{gerepileall}682683Let us now see why we may need the \kbd{gerepileall} variants. Although it684is not an infrequent occurrence, we do not give a specific example but a685general one: suppose that we want to do a computation (usually inside a686larger function) producing more than one PARI object as a result, say two for687instance. Then even if we set up the work properly, before cleaning up we688have a stack which has the desired results \kbd{z1}, \kbd{z2} (say), and689then connected garbage from lbot to ltop. If we write690\bprog691z1 = gerepile(ltop, lbot, z1);692@eprog\noindent693then the stack is cleaned, the pointers fixed up, but we have lost the694address of \kbd{z2}. This is where we need the \idx{gerepileall}695function:696\bprog697gerepileall(ltop, 2, &z1, &z2)698@eprog699\noindent copies \kbd{z1} and \kbd{z2} to new locations, cleans the stack700from \kbd{ltop} to the old \kbd{avma}, and updates the pointers \kbd{z1} and701\kbd{z2}. Here we do not assume anything about the stack: the garbage can be702disconnected and \kbd{z1}, \kbd{z2} need not be at the bottom of the stack.703If all of these assumptions are in fact satisfied, then we can call704\kbd{gerepilemanysp} instead, which is usually faster since we do not need705the initial copy (on the other hand, it is less cache friendly).706707A most important usage is ``random'' garbage collection during loops708whose size requirements we cannot (or do not bother to) control in advance:709\bprog710pari_sp av = avma;711GEN x, y;712while (...)713{714garbage(); x = anything();715garbage(); y = anything(); garbage();716if (gc_needed(av,1)) /*@Ccom memory is running low (half spent since entry) */717gerepileall(av, 2, &x, &y);718}719@eprog720\noindent Here we assume that only \kbd{x} and \kbd{y} are needed from one721iteration to the next. As it would be costly to call gerepile once for each722iteration, we only do it when it seems to have become necessary.723724More precisely, the macro \tet{stack_lim}\kbd{(av,$n$)} denotes an address725where $2^{n-1} / (2^{n-1}+1)$ of the remaining stack space since reference726point \kbd{av} is exhausted ($1/2$ for $n=1$, $2/3$ for $n=2$). The test727\tet{gc_needed}\kbd{(av,$n$)} becomes true whenever \kbd{avma} drops below728that address.729730\subsec{Comments}731732First, \kbd{gerepile} has turned out to be a flexible and fast garbage733collector for number-theoretic computations, which compares favorably with734more sophisticated methods used in other systems. Our benchmarks indicate735that the price paid for using \kbd{gerepile} and \kbd{gerepile}-related736copies, when properly used, is usually less than 1\% of the total737running time, which is quite acceptable!738739Second, it is of course harder on the programmer, and quite error-prone740if you do not stick to a consistent PARI programming style. If all seems741lost, just use \tet{gerepilecopy} (or \tet{gerepileall}) to fix up the stack742for you. You can always optimize later when you have sorted out exactly which743routines are crucial and what objects need to be preserved and their usual744sizes.745746\smallskip If you followed us this far, congratulations, and rejoice: the747rest is much easier.748749\section{Creation of PARI objects, assignments, conversions}750751\subsec{Creation of PARI objects}\sidx{creation}752The basic function which creates a PARI object is753754\fun{GEN}{cgetg}{long l, long t}755$l$ specifies the number of longwords to be allocated to the756object, and $t$ is the type of the object, in symbolic757form (see \secref{se:impl} for the list of these). The precise effect of758this function is as follows: it first creates on the PARI \emph{stack} a759chunk of memory of size \kbd{length} longwords, and saves the address of the760chunk which it will in the end return. If the stack has been used up, a761message to the effect that ``the PARI stack overflows'' is printed,762and an error raised. Otherwise, it sets the type and length of the PARI object.763In effect, it fills its first codeword (\kbd{z[0]}). Many PARI764objects also have a second codeword (types \typ{INT}, \typ{REAL},765\typ{PADIC}, \typ{POL}, and \typ{SER}). In case you want to produce one of766those from scratch, which should be exceedingly rare, \emph{it is your767responsibility to fill this second codeword}, either explicitly (using the768macros described in \secref{se:impl}), or implicitly using an assignment769statement (using \kbd{gaffect}).770771Note that the length argument $l$ is predetermined for a number of types:7723 for types \typ{INTMOD}, \typ{FRAC}, \typ{COMPLEX}, \typ{POLMOD},773\typ{RFRAC}, 4 for type \typ{QUAD}, and 5 for type \typ{PADIC} and \typ{QFB}.774However for the sake of efficiency, \kbd{cgetg} does not check this: disasters775will occur if you give an incorrect length for those types.776777\misctitle{Notes} 1) The main use of this function is create efficiently778a constant object, or to prepare for later assignments (see779\secref{se:assign}). Most of the time you will use \kbd{GEN} objects as they780are created and returned by PARI functions. In this case you do not need to781use \kbd{cgetg} to create space to hold them.782783\noindent 2) For the creation of leaves, i.e.~\typ{INT} or \typ{REAL},784785\fun{GEN}{cgeti}{long length}786787\fun{GEN}{cgetr}{long length}788789\noindent should be used instead of \kbd{cgetg(length, t\_INT)} and790\kbd{cgetg(length, t\_REAL)} respectively. Finally791792\fun{GEN}{cgetc}{long prec}793794\noindent creates a \typ{COMPLEX} whose real and imaginary part are795\typ{REAL}s allocated by \kbd{cgetr(prec)}.796797\misctitle{Examples} 1) Both \kbd{z = cgeti(DEFAULTPREC)} and798\kbd{cgetg(DEFAULTPREC, t\_INT)} create a \typ{INT} whose ``precision'' is799\kbd{bit\_accuracy(DEFAULTPREC)} = 64. This means \kbd{z} can hold rational800integers of absolute value less than $2^{64}$. Note that in both cases, the801second codeword is \emph{not} filled. Of course we could use numerical802values, e.g.~\kbd{cgeti(4)}, but this would have different meanings on803different machines as \kbd{bit\_accuracy(4)} equals 64 on 32-bit machines,804but 128 on 64-bit machines.805806\noindent 2) The following creates a \emph{complex number} whose real and807imaginary parts can hold real numbers of precision808$\kbd{bit\_accuracy(MEDDEFAULTPREC)} = 96\hbox{ bits:}$809%810\bprog811z = cgetg(3, t_COMPLEX);812gel(z, 1) = cgetr(MEDDEFAULTPREC);813gel(z, 2) = cgetr(MEDDEFAULTPREC);814@eprog\noindent815or simply \kbd{z = cgetc(MEDDEFAULTPREC)}.816817\noindent 3) To create a matrix object for $4\times 3$ matrices:818%819\bprog820z = cgetg(4, t_MAT);821for(i=1; i<4; i++) gel(z, i) = cgetg(5, t_COL);822@eprog\noindent823or simply \kbd{z = zeromatcopy(4, 3)}, which further initializes all entries824to \kbd{gen\_0}.825826These last two examples illustrate the fact that since PARI types are827recursive, all the branches of the tree must be created. The function828\teb{cgetg} creates only the ``root'', and other calls to \kbd{cgetg} must be829made to produce the whole tree. For matrices, a common mistake is to think830that \kbd{z = cgetg(4, t\_MAT)} (for example) creates the root of the831matrix: one needs also to create the column vectors of the matrix (obviously,832since we specified only one dimension in the first \kbd{cgetg}!). This is833because a matrix is really just a row vector of column vectors (hence a834priori not a basic type), but it has been given a special type number so that835operations with matrices become possible.836837Finally, to facilitate input of constant objects when speed is not paramount,838there are four \tet{varargs} functions:839840\fun{GEN}{mkintn}{long n, ...}841returns the nonnegative \typ{INT} whose development in base $2^{32}$842is given by the following $n$ 32bit-words (\kbd{unsigned int}).843\bprog844mkintn(3, a2, a1, a0);845@eprog846\noindent returns $a_2 2^{64} + a_1 2^{32} + a_0$.847848\fun{GEN}{mkpoln}{long n, ...}849Returns the \typ{POL} whose $n$ coefficients (\kbd{GEN}) follow, in order of850decreasing degree.851\bprog852mkpoln(3, gen_1, gen_2, gen_0);853@eprog854\noindent returns the polynomial $X^2 + 2X$ (in variable $0$, use855\tet{setvarn} if you want other variable numbers). Beware that $n$ is the856number of coefficients, hence \emph{one more} than the degree.857858\fun{GEN}{mkvecn}{long n, ...}859returns the \typ{VEC} whose $n$ coefficients (\kbd{GEN}) follow.860861\fun{GEN}{mkcoln}{long n, ...}862returns the \typ{COL} whose $n$ coefficients (\kbd{GEN}) follow.863864\misctitle{Warning} Contrary to the policy of general PARI functions, the865latter three functions do \emph{not} copy their arguments, nor do they produce866an object a priori suitable for \tet{gerepileupto}. For instance867\bprog868/*@Ccom gerepile-safe: components are universal objects */869z = mkvecn(3, gen_1, gen_0, gen_2);870871/*@Ccom not OK for gerepileupto: stoi(3) creates component before root */872z = mkvecn(3, stoi(3), gen_0, gen_2);873874/*@Ccom NO! First vector component \kbd{x} is destroyed */875x = gclone(gen_1);876z = mkvecn(3, x, gen_0, gen_2);877gunclone(x);878@eprog879880\noindent The following function is also available as a special case of881\tet{mkintn}:882883\fun{GEN}{uu32toi}{ulong a, ulong b}884885Returns the \kbd{GEN} equal to $2^{32} a + b$, \emph{assuming} that886$a,b < 2^{32}$. This does not depend on \kbd{sizeof(long)}: the behavior is887as above on both $32$ and $64$-bit machines.888889\subsec{Sizes}890891\fun{long}{gsizeword}{GEN x} returns the total number of \B-bit words occupied892by the tree representing~\kbd{x}.893894\fun{long}{gsizebyte}{GEN x} returns the total number of bytes occupied895by the tree representing~\kbd{x}, i.e.~\kbd{gsizeword(x)} multiplied by896\kbd{sizeof(long)}. This is normally useless since PARI functions use897a number of \emph{words} as input for lengths and precisions.898899\subsec{Assignments}900Firstly, if \kbd{x} and \kbd{y} are both declared as \kbd{GEN} (i.e.~pointers901to something), the ordinary C assignment \kbd{y = x} makes perfect sense: we902are just moving a pointer around. However, physically modifying either903\kbd{x} or \kbd{y} (for instance, \kbd{x[1] = 0}) also changes the other904one, which is usually not desirable. \label{se:assign}905906\misctitle{Very important note} Using the functions described in this907paragraph is inefficient and often awkward: one of the \tet{gerepile}908functions (see~\secref{se:garbage}) should be preferred. See the paragraph909end for one exception to this rule.910911\noindent912The general PARI \idx{assignment} function is the function \teb{gaffect} with913the following syntax:914915\fun{void}{gaffect}{GEN x, GEN y}916917\noindent918Its effect is to assign the PARI object \kbd{x} into the \emph{preexisting}919object \kbd{y}. Both \kbd{x} and \kbd{y} must be \emph{scalar} types. For920convenience, vector or matrices of scalar types are also allowed.921922This copies the whole structure of \kbd{x} into \kbd{y} so many conditions923must be met for the assignment to be possible. For instance it is allowed to924assign a \typ{INT} into a \typ{REAL}, but the converse is forbidden. For925that, you must use the truncation or rounding function of your choice,926e.g.\kbd{mpfloor}.927928It can also happen that \kbd{y} is not large enough or does not have the proper929tree structure to receive the object \kbd{x}. For instance, let \kbd{y} the zero930integer with length equal to 2; then \kbd{y} is too small to accommodate any931nonzero \typ{INT}. In general common sense tells you what is possible,932keeping in mind the PARI philosophy which says that if it makes sense it is933valid. For instance, the assignment of an imprecise object into a precise one934does \emph{not} make sense. However, a change in precision of imprecise935objects is allowed, even if it \emph{increases} its accuracy: we complement936the ``mantissa'' with infinitely many $0$ digits in this case. (Mantissa937between quotes, because this is not restricted to \typ{REAL}s, it also938applies for $p$-adics for instance.)939940All functions ending in ``\kbd{z}'' such as \teb{gaddz}941(see~\secref{se:low_level}) implicitly use this function. In fact what they942exactly do is record {\teb{avma}} (see~\secref{se:garbage}), perform the943required operation, \teb{gaffect} the result to the last operand, then944restore the initial \kbd{avma}.945946You can assign ordinary C long integers into a PARI object (not necessarily947of type \typ{INT}) using948949\fun{void}{gaffsg}{long s, GEN y}950951\misctitle{Note} Due to the requirements mentioned above, it is usually952a bad idea to use \tet{gaffect} statements. There is one exception: for simple953objects (e.g.~leaves) whose size is controlled, they can be easier to use than954\kbd{gerepile}, and about as efficient.955956\misctitle{Coercion} It is often useful to coerce an inexact object to a957given precision. For instance at the beginning of a routine where precision958can be kept to a minimum; otherwise the precision of the input is used in all959subsequent computations, which is inefficient if the latter is known to960thousands of digits. One may use the \kbd{gaffect} function for this, but it961is easier and more efficient to call962963\fun{GEN}{gtofp}{GEN x, long prec} converts the complex number~\kbd{x}964(\typ{INT}, \typ{REAL}, \typ{FRAC}, \typ{QUAD} or \typ{COMPLEX}) to either965a \typ{REAL} or \typ{COMPLEX} whose components are \typ{REAL} of length966\kbd{prec}.967968\subsec{Copy} It is also very useful to \idx{copy} a PARI object, not969just by moving around a pointer as in the \kbd{y = x} example, but by970creating a copy of the whole tree structure, without pre-allocating a971possibly complicated \kbd{y} to use with \kbd{gaffect}. The function which972does this is called \teb{gcopy}. Its syntax is:973974\fun{GEN}{gcopy}{GEN x}975976\noindent and the effect is to create a new copy of x on the PARI stack.977978Sometimes, on the contrary, a quick copy of the skeleton of \kbd{x} is979enough, leaving pointers to the original data in \kbd{x} for the sake of980speed instead of making a full recursive copy. Use981\fun{GEN}{shallowcopy}{GEN x} for this. Note that the result is not suitable982for \tet{gerepileupto} !983984Make sure at this point that you understand the difference between \kbd{y =985x}, \kbd{y = gcopy(x)}, \kbd{y = shallowcopy(x)} and \kbd{gaffect(x,y)}.986987\subsec{Clones}\sidx{clone}\label{se:clone}988Sometimes, it is more efficient to create a \emph{persistent} copy of a PARI989object. This is not created on the stack but on the heap, hence unaffected by990\tet{gerepile} and friends. The function which does this is called991\teb{gclone}. Its syntax is:992993\fun{GEN}{gclone}{GEN x}994995A clone can be removed from the heap (thus destroyed) using996997\fun{void}{gunclone}{GEN x}998999\noindent No PARI object should keep references to a clone which has been1000destroyed!10011002\subsec{Conversions}\sidx{conversions}1003The following functions convert C objects to PARI objects (creating them on1004the stack as usual):10051006\fun{GEN}{stoi}{long s}: C \kbd{long} integer (``small'') to \typ{INT}.10071008\fun{GEN}{dbltor}{double s}: C \kbd{double} to \typ{REAL}. The accuracy of1009the result is 19 decimal digits, i.e.~a type \typ{REAL} of length1010\kbd{DEFAULTPREC}, although on 32-bit machines only 16 of them are1011significant.10121013\noindent We also have the converse functions:10141015\fun{long}{itos}{GEN x}: \kbd{x} must be of type \typ{INT},10161017\fun{double}{rtodbl}{GEN x}: \kbd{x} must be of type \typ{REAL},10181019\noindent as well as the more general ones:10201021\fun{long}{gtolong}{GEN x},10221023\fun{double}{gtodouble}{GEN x}.10241025\section{Implementation of the PARI types}1026\label{se:impl}10271028\noindent1029We now go through each type and explain its implementation. Let \kbd{z} be a1030\kbd{GEN}, pointing at a PARI object. In the following paragraphs, we will1031constantly mix two points of view: on the one hand, \kbd{z} is treated as the1032C pointer it is, on the other, as PARI's handle on some mathematical entity,1033so we will shamelessly write $\kbd{z} \ne 0$ to indicate that the1034\emph{value} thus represented is nonzero (in which case the1035\emph{pointer}~\kbd{z} is certainly not \kbd{NULL}). We offer no apologies1036for this style. In fact, you had better feel comfortable juggling both views1037simultaneously in your mind if you want to write correct PARI programs.10381039Common to all the types is the first codeword \kbd{z[0]}, which we do not1040have to worry about since this is taken care of by \kbd{cgetg}. Its precise1041structure depends on the machine you are using, but it always contains the1042following data: the \emph{internal type number}\sidx{type number} attached1043to the symbolic type name, the \emph{length} of the root in longwords, and a1044technical bit which indicates whether the object is a clone or not (see1045\secref{se:clone}). This last one is used by \kbd{gp} for internal garbage1046collecting, you will not have to worry about it.10471048Some types have a second codeword, different for each type, which we will1049soon describe as we will shortly consider each of them in turn.10501051\noindent The first codeword is handled through the following \emph{macros}:10521053\fun{long}{typ}{GEN z} returns the type number of \kbd{z}.10541055\fun{void}{settyp}{GEN z, long n} sets the type number of \kbd{z} to1056\kbd{n} (you should not have to use this function if you use \kbd{cgetg}).10571058\fun{long}{lg}{GEN z} returns the length (in longwords) of the root of \kbd{z}.10591060\fun{long}{setlg}{GEN z, long l} sets the length of \kbd{z} to \kbd{l}; you1061should not have to use this function if you use \kbd{cgetg}.10621063\fun{void}{lg_increase}{GEN z} increase the length of \kbd{z} by 1; you1064should not have to use this function if you use \kbd{cgetg}.10651066\fun{long}{isclone}{GEN z} is \kbd{z} a clone?10671068\fun{void}{setisclone}{GEN z} sets the \emph{clone} bit.10691070\fun{void}{unsetisclone}{GEN z} clears the \emph{clone} bit.10711072\misctitle{Important remark} For the sake of efficiency, none of the1073codeword-handling macros check the types of their arguments even when there1074are stringent restrictions on their use. It is trivial to create invalid1075objects, or corrupt one of the ``universal constants'' (e.g. setting the sign1076of \kbd{gen\_0} to $1$), and they usually provide negligible savings.1077Use higher level functions whenever possible.10781079\misctitle{Remark} The clone bit is there so that \kbd{gunclone} can check1080it is deleting an object which was allocated by \kbd{gclone}. Miscellaneous1081vector entries are often cloned by \kbd{gp} so that a GP statement like1082\kbd{v[1] = x} does not involve copying the whole of \kbd{v}: the component1083\kbd{v[1]} is deleted if its clone bit is set, and is replaced by a clone of1084\kbd{x}. Don't set/unset yourself the clone bit unless you know what you are1085doing: in particular \emph{never} set the clone bit of a vector component1086when the said vector is scheduled to be uncloned. Hackish code may abuse the1087clone bit to tag objects for reasons unrelated to the above instead of using1088proper data structures. Don't do that.10891090\subsec{Type \typ{INT} (integer)}1091\sidx{integer}\kbdsidx{t_INT}this type has a second codeword \kbd{z[1]} which1092contains the following information:10931094the sign of \kbd{z}: coded as $1$, $0$ or $-1$ if $\kbd{z} > 0$, $\kbd{z} = 0$,1095$\kbd{z} < 0$ respectively.10961097the \emph{effective length} of \kbd{z}, i.e.~the total number of significant1098longwords. This means the following: apart from the integer 0, every integer1099is ``normalized'', meaning that the most significant mantissa longword is1100nonzero. However, the integer may have been created with a longer length.1101Hence the ``length'' which is in \kbd{z[0]} can be larger than the1102``effective length'' which is in \kbd{z[1]}.11031104\noindent This information is handled using the following macros:11051106\fun{long}{signe}{GEN z} returns the sign of \kbd{z}.11071108\fun{void}{setsigne}{GEN z, long s} sets the sign of \kbd{z} to \kbd{s}.11091110\fun{long}{lgefint}{GEN z} returns the \idx{effective length} of \kbd{z}.11111112\fun{void}{setlgefint}{GEN z, long l} sets the effective length1113of \kbd{z} to \kbd{l}.11141115The integer 0 can be recognized either by its sign being~0, or by its1116effective length being equal to~2. Now assume that $\kbd{z} \ne 0$, and let1117$$ |z| = \sum_{i = 0}^n z_i B^i,1118\quad\text{where}~z_n\ne 0~\text{and}~B = 2^{\kbd{BITS\_IN\_LONG}}.1119$$1120With these notations, $n$ is \kbd{lgefint(z) - 3}, and the mantissa of1121$\kbd{z}$ may be manipulated via the following interface:11221123\fun{GEN}{int_MSW}{GEN z} returns a pointer to the most significant word of1124\kbd{z}, $z_n$.11251126\fun{GEN}{int_LSW}{GEN z} returns a pointer to the least significant word of1127\kbd{z}, $z_0$.11281129\fun{GEN}{int_W}{GEN z, long i} returns the $i$-th significant word of1130\kbd{z}, $z_i$. Accessing the $i$-th significant word for $i > n$1131yields unpredictable results.11321133\fun{GEN}{int_W_lg}{GEN z, long i, long lz} returns the $i$-th significant1134word of \kbd{z}, $z_i$, assuming \kbd{lgefint(z)} is \kbd{lz} ($= n + 3$).1135Accessing the $i$-th significant word for $i > n$ yields unpredictable1136results.11371138\fun{GEN}{int_precW}{GEN z} returns the previous (less significant) word of1139\kbd{z}, $z_{i-1}$ assuming \kbd{z} points to $z_i$.11401141\fun{GEN}{int_nextW}{GEN z} returns the next (more significant) word of \kbd{z},1142$z_{i+1}$ assuming \kbd{z} points to $z_i$.11431144Unnormalized integers, such that $z_n$ is possibly $0$, are explicitly1145forbidden. To enforce this, one may write an arbitrary mantissa then call11461147\fun{void}{int_normalize}{GEN z, long known0}11481149\noindent normalizes in place a nonnegative integer (such that $z_n$ is1150possibly $0$), assuming at least the first \kbd{known0} words are zero.11511152\noindent For instance a binary \kbd{and} could be implemented in the1153following way:1154\bprog1155GEN AND(GEN x, GEN y) {1156long i, lx, ly, lout;1157long *xp, *yp, *outp; /* mantissa pointers */1158GEN out;11591160if (!signe(x) || !signe(y)) return gen_0;1161lx = lgefint(x); xp = int_LSW(x);1162ly = lgefint(y); yp = int_LSW(y); lout = min(lx,ly); /* > 2 */11631164out = cgeti(lout); out[1] = evalsigne(1) | evallgefint(lout);1165outp = int_LSW(out);1166for (i=2; i < lout; i++)1167{1168*outp = (*xp) & (*yp);1169outp = int_nextW(outp);1170xp = int_nextW(xp);1171yp = int_nextW(yp);1172}1173if ( !*int_MSW(out) ) out = int_normalize(out, 1);1174return out;1175}1176@eprog11771178\noindent This low-level interface is mandatory in order to write portable1179code since PARI can be compiled using various multiprecision kernels, for1180instance the native one or GNU MP, with incompatible internal structures1181(for one thing, the mantissa is oriented in different directions).11821183\subsec{Type \typ{REAL} (real number)}1184\kbdsidx{t_REAL}\sidx{real number}this type has a second codeword z[1] which1185also encodes its sign, obtained or set using the same functions as for a1186\typ{INT}, and a binary exponent. This exponent is handled using the1187following macros:11881189\fun{long}{expo}{GEN z} returns the exponent of \kbd{z}.1190This is defined even when \kbd{z} is equal to zero.11911192\fun{void}{setexpo}{GEN z, long e} sets the exponent of \kbd{z} to \kbd{e}.11931194\noindent Note the functions:11951196\fun{long}{gexpo}{GEN z} which tries to return an exponent for \kbd{z},1197even if \kbd{z} is not a real number.11981199\fun{long}{gsigne}{GEN z} which returns a sign for \kbd{z}, even when1200\kbd{z} is a real number of type \typ{INT}, \typ{FRAC} or \typ{REAL},1201an infinity (\typ{INFINITY}) or a \typ{QUAD} of positive discriminant.12021203The real zero is characterized by having its sign equal to 0. If \kbd{z} is1204not equal to~0, then it is represented as $2^e M$, where $e$ is the exponent,1205and $M\in [1, 2[$ is the mantissa of $z$, whose digits are stored in1206$\kbd{z[2]},\dots, \kbd{z[lg(z)-1]}$. For historical reasons, the \kbd{prec}1207parameter attached to floating point functions is measured in \B-bit words1208and is equal to the length of \kbd{x}: yes, this includes the two code words1209and depends on \kbd{sizeof(long)}. For clarity we advise to use1210\kbd{bit\_accuracy}, which computes the true length of the mantissa in bits,1211and convert between bits and \kbd{prec} using the \kbd{prec2nbits} and1212\kbd{nbits2prec} macros. But keep in mind that the accuracy of \typ{REAL}1213actually increases by increments of \B bits.12141215More precisely, let $m$ be the integer (\kbd{z[2]},\dots, \kbd{z[lg(z)-1]})1216in base \kbd{2\pow BITS\_IN\_LONG}; here, \kbd{z[2]} is the most significant1217longword and is normalized, i.e.~its most significant bit is~1. Then we have1218$M := m / 2^{\kbd{bit\_accuracy(lg(z))} - 1 - \kbd{expo}(z)}$.12191220\fun{GEN}{mantissa_real}{GEN z, long *e} returns the mantissa $m$ of $z$, and1221sets \kbd{*e} to the exponent $\kbd{bit\_accuracy(lg(z))}-1-\kbd{expo}(z)$,1222so that $z = m / 2^e$.12231224Thus, the real number $3.5$ to accuracy \kbd{bit\_accuracy(lg(z))} is1225represented as \kbd{z[0]} (encoding $\kbd{type} = \typ{REAL}$, \kbd{lg(z)}),1226\kbd{z[1]} (encoding $\kbd{sign} = 1$, $\kbd{expo} = 1$), $\kbd{z[2]} =1227\kbd{0xe0000000}$, $\kbd{z[3]} =\dots = \kbd{z[lg(z)-1]} = \kbd{0x0}$.12281229\subsec{Type \typ{INTMOD}}\kbdsidx{t_INTMOD}1230\kbd{z[1]} points to the modulus, and \kbd{z[2]} at the number representing1231the class \kbd{z}. Both are separate \kbd{GEN} objects, and both must be1232\typ{INT}s, satisfying the inequality $0 \le \kbd{z[2]} < \kbd{z[1]}$.12331234\subsec{Type \typ{FRAC} (rational number)}%1235\kbdsidx{t_FRAC}\sidx{rational number}1236\kbd{z[1]} points to the numerator $n$, and \kbd{z[2]} to the denominator1237$d$. Both must be of type \typ{INT} such that $n\neq 0$, $d > 0$ and1238$(n,d) = 1$.12391240\subsec{Type \typ{FFELT} (finite field element)}%1241\kbdsidx{t_FFELT}\sidx{finite field element} (Experimental)12421243Components of this type should normally not be accessed directly. Instead,1244finite field elements should be created using \kbd{ffgen}.12451246\noindent The second codeword \kbd{z[1]} determines the storage format of the1247element, among12481249\item \tet{t_FF_FpXQ}: \kbd{A=z[2]} and \kbd{T=z[3]} are \kbd{FpX},1250\kbd{p=z[4]} is a \typ{INT}, where $p$ is a prime number, $T$ is irreducible1251modulo $p$, and $\deg A < \deg T$.1252This represents the element $A\pmod{T}$ in $\F_p[X]/T$.12531254\item \tet{t_FF_Flxq}: \kbd{A=z[2]} and \kbd{T=z[3]} are \kbd{Flx},1255\kbd{l=z[4]} is a \typ{INT}, where $l$ is a prime number, $T$ is irreducible1256modulo $l$, and $\deg A < \deg T$ This represents the element $A\pmod{T}$ in1257$\F_l[X]/T$.12581259\item \tet{t_FF_F2xq}: \kbd{A=z[2]} and \kbd{T=z[3]} are \kbd{F2x},1260\kbd{l=z[4]} is the \typ{INT} $2$, $T$ is irreducible modulo $2$, and1261$\deg A < \deg T$. This represents the element $A\pmod{T}$ in $\F_2[X]/T$.12621263\subsec{Type \typ{COMPLEX} (complex number)}%1264\kbdsidx{t_COMPLEX}\sidx{complex number}1265\kbd{z[1]} points to the real part, and \kbd{z[2]} to the imaginary part.1266The components \kbd{z[1]} and \kbd{z[2]} must be of type1267\typ{INT}, \typ{REAL} or \typ{FRAC}. For historical reasons \typ{INTMOD}1268and \typ{PADIC} are also allowed (the latter for $p = 2$ or1269congruent to 3 mod 4 only), but one should rather use the more general1270\typ{POLMOD} construction.12711272\subsec{Type \typ{PADIC} ($p$-adic numbers)}%1273\sidx{p-adic number}\kbdsidx{t_PADIC} this type has a second codeword1274\kbd{z[1]} which contains the following information: the $p$-adic precision1275(the exponent of $p$ modulo which the $p$-adic unit corresponding to1276\kbd{z} is defined if \kbd{z} is not~0), i.e.~one less than the number of1277significant $p$-adic digits, and the exponent of \kbd{z}. This information1278can be handled using the following functions:12791280\fun{long}{precp}{GEN z} returns the $p$-adic precision of \kbd{z}. This is1281$0$ if $\kbd{z} = 0$.12821283\fun{void}{setprecp}{GEN z, long l} sets the $p$-adic precision of \kbd{z}1284to \kbd{l}.12851286\fun{long}{valp}{GEN z} returns the $p$-adic valuation of \kbd{z}1287(i.e. the exponent). This is defined even if \kbd{z} is equal to~0.12881289\fun{void}{setvalp}{GEN z, long e} sets the $p$-adic valuation of \kbd{z}1290to \kbd{e}.12911292In addition to this codeword, \kbd{z[2]} points to the prime $p$,1293\kbd{z[3]} points to $p^{\text{precp(z)}}$, and \kbd{z[4]} points to1294a\typ{INT} representing the $p$-adic unit attached to \kbd{z} modulo1295\kbd{z[3]} (and to zero if \kbd{z} is zero). To summarize, if $z\neq12960$, we have the equality:1297$$ \kbd{z} = p^{\text{valp(z)}} * (\kbd{z[4]} + O(\kbd{z[3]})),\quad1298\text{where}\quad \kbd{z[3]} = O(p^{\text{precp(z)}}). $$12991300\subsec{Type \typ{QUAD} (quadratic number)}1301\sidx{quadratic number}\kbdsidx{t_QUAD}\kbd{z[1]} points to the canonical1302polynomial $P$ defining the quadratic field (as output by \tet{quadpoly}),1303\kbd{z[2]} to the ``real part'' and \kbd{z[3]} to the ``imaginary part''. The1304latter are of type \typ{INT}, \typ{FRAC}, \typ{INTMOD}, or \typ{PADIC} and1305are to be taken as the coefficients of \kbd{z} with respect to the canonical1306basis $(1,X)$ of $\Q[X]/(P(X))$. Exact complex numbers may be implemented as1307quadratics, but \typ{COMPLEX} is in general more versatile (\typ{REAL}1308components are allowed) and more efficient.13091310Operations involving a \typ{QUAD} and \typ{COMPLEX} are implemented by1311converting the \typ{QUAD} to a \typ{REAL} (or \typ{COMPLEX} with \typ{REAL}1312components) to the accuracy of the \typ{COMPLEX}. As a consequence,1313operations between \typ{QUAD} and \emph{exact} \typ{COMPLEX}s are not allowed.13141315\subsec{Type \typ{POLMOD} (polmod)}\kbdsidx{t_POLMOD}\sidx{polmod}1316as for \typ{INTMOD}s, \kbd{z[1]} points to the modulus, and \kbd{z[2]}1317to a polynomial representing the class of~\kbd{z}. Both must be of type1318\typ{POL} in the same variable, satisfying the inequality $\deg \kbd{z[2]}1319< \deg \kbd{z[1]}$. However, \kbd{z[2]} is allowed to be a simplification1320of such a polynomial, e.g.~a scalar. This is tricky considering the1321hierarchical structure of the variables; in particular, a polynomial in1322variable of \emph{lesser} priority (see \secref{se:vars}) than the1323modulus variable is valid, since it is considered as the constant term of1324a polynomial of degree 0 in the correct variable. On the other hand a1325variable of \emph{greater} priority is not acceptable.13261327\subsec{Type \typ{POL} (polynomial)}\kbdsidx{t_POL}\sidx{polynomial} this1328type has a second codeword. It contains a ``\emph{sign}'': 0 if the1329polynomial is equal to~0, and 1 if not (see however the important remark1330below) and a \emph{variable number} (e.g.~0 for $x$, 1 for $y$, etc\dots).13311332\noindent These data can be handled with the following macros: \teb{signe}1333and \teb{setsigne} as for \typ{INT} and \typ{REAL},13341335\fun{long}{varn}{GEN z} returns the variable number of the object \kbd{z},13361337\fun{void}{setvarn}{GEN z, long v} sets the variable number of \kbd{z} to1338\kbd{v}.13391340The variable numbers encode the relative priorities of variables, we will1341give more details in \secref{se:vars}. Note also the function1342\fun{long}{gvar}{GEN z} which tries to return a \idx{variable number} for1343\kbd{z}, even if \kbd{z} is not a polynomial or power series. The variable1344number of a scalar type is set by definition equal to \tet{NO_VARIABLE},1345which has lower priority than any other variable number.13461347The components \kbd{z[2]}, \kbd{z[3]},\dots \kbd{z[lg(z)-1]} point to the1348coefficients of the polynomial \emph{in ascending order}, with \kbd{z[2]}1349being the constant term and so on.13501351For a \typ{POL} of nonzero sign, \tet{degpol}, \tet{leading_coeff},1352\tet{constant_coeff}, return its degree, and a pointer to the leading,1353resp. constant, coefficient with respect to the main variable. Note that no1354copy is made on the PARI stack so the returned value is not safe for a basic1355\kbd{gerepile} call. Applied to any other type than \typ{POL}, the result is1356unspecified. Those three functions are still defined when the sign is $0$,1357see \secref{se:accessors} and \secref{se:polynomials}.13581359\fun{long}{degree}{GEN x} returns the degree of \kbd{x} with respect to its1360main variable even when \kbd{x} is not a polynomial (a rational function for1361instance). By convention, the degree of a zero polynomial is~$-1$.13621363\misctitle{Important remark}1364The leading coefficient of a \typ{POL} may be equal to zero:13651366\item it is not allowed to be an exact rational $0$, such as \tet{gen_0};13671368\item an exact nonrational $0$, like \kbd{Mod(0,2)}, is possible for1369constant polynomials, i.e. of length $3$ and no other coefficient: this1370carries information about the base ring for the polynomial;13711372\item an inexact $0$, like \kbd{0.E-38} or \kbd{O(3\pow 5)}, is always1373possible. Inexact zeroes do not correspond to an actual $0$, but to a1374very small coefficient according to some metric; we keep them to give1375information on how much cancellation occurred in previous computations.13761377A polynomial disobeying any of these rules is an invalid \emph{unnormalized}1378object. We advise \emph{not} to use low-level constructions to build a1379\typ{POL} coefficient by coefficient, such as1380\bprog1381GEN T = cgetg(4, t_POL);1382T[1] = evalvarn(0);1383gel(T, 2) = x;1384gel(T, 3) = y;1385@eprog\noindent But if you do and it is not clear whether the result will be1386normalized, call13871388\fun{GEN}{normalizepol}{GEN x} applied to an unnormalized \typ{POL}~\kbd{x}1389(with all coefficients correctly set except that \kbd{leading\_term(x)} might1390be zero), normalizes \kbd{x} correctly in place and returns~\kbd{x}. This1391functions sets \kbd{signe} (to $0$ or $1$) properly.13921393\misctitle{Caveat} A consequence of the remark above is that zero polynomials1394are characterized by the fact that their sign is~0. It is in general1395incorrect to check whether \kbd{lg(x)} is $2$ or \kbd{degpol(x)} $< 0$,1396although both tests are valid when the coefficient types are under control:1397for instance, when they are guaranteed to be \typ{INT}s or \typ{FRAC}s.1398The same remark applies to \typ{SER}s.13991400\subsec{Type \typ{SER} (power series)}1401\kbdsidx{t_SER}\sidx{power series}This type also has a second codeword, which1402encodes a ``\emph{sign}'', i.e.~0 if the power series is 0, and 1 if not, a1403\emph{variable number} as for polynomials, and an \emph{exponent}. This1404information can be handled with the following functions: \teb{signe},1405\teb{setsigne}, \teb{varn}, \teb{setvarn} as for polynomials, and \teb{valp},1406\teb{setvalp} for the exponent as for $p$-adic numbers. Beware: do \emph{not}1407use \teb{expo} and \teb{setexpo} on power series.14081409The coefficients \kbd{z[2]}, \kbd{z[3]},\dots \kbd{z[lg(z)-1]} point to1410the coefficients of \kbd{z} in ascending order. As for polynomials1411(see remark there), the sign of a \typ{SER} is $0$ if and only all1412its coefficients are equal to $0$. (The leading coefficient cannot be an1413integer $0$.) A series whose coefficients are integers equal to zero1414is represented as $O(x^n)$ (\kbd{zeroser}$(\var{vx},n)$).1415A series whose coefficients are exact zeroes, but not all of1416them integers (e.g. an \typ{INTMOD} such as \kbd{Mod(0,2)}) is1417represented as $z*x^{n-1} +O(x^n)$, where $z$ is the $0$ of the1418base ring, as per \tet{Rg_get_0}.14191420Note that the exponent of a power series can be negative, i.e.~we are1421then dealing with a Laurent series (with a finite number of negative1422terms).14231424\subsec{Type \typ{RFRAC} (rational function)}%1425\kbdsidx{t_RFRAC}\sidx{rational function} \kbd{z[1]} points to the1426numerator $n$, and \kbd{z[2]} on the denominator $d$. The denominator must be1427of type \typ{POL}, with variable of higher priority than the numerator. The1428numerator $n$ is not an exact $0$ and $(n,d) = 1$ (see \tet{gred_rfac2}).14291430\subsec{Type \typ{QFB} (binary quadratic form)}%1431\kbdsidx{t_QFB}\sidx{binary quadratic form} \kbd{z[1]}, \kbd{z[2]},1432\kbd{z[3]} point to the three coefficients of the form, and \kbd{z[4]} point to1433the form discriminant. All four are of type \typ{INT}.14341435\subsec{Type \typ{VEC} and \typ{COL} (vector)}%1436\kbdsidx{t_VEC}\kbdsidx{t_COL}\sidx{row vector}\sidx{column vector}1437\kbd{z[1]}, \kbd{z[2]},\dots \kbd{z[lg(z)-1]} point to the components of the1438vector.14391440\subsec{Type \typ{MAT} (matrix)}\kbdsidx{t_MAT}\sidx{matrix} \kbd{z[1]},1441\kbd{z[2]},\dots \kbd{z[lg(z)-1]} point to the column vectors of \kbd{z},1442i.e.~they must be of type \typ{COL} and of the same length.14431444\subsec{Type \typ{VECSMALL} (vector of small integers)}\kbdsidx{t_VECSMALL}1445\kbd{z[1]}, \kbd{z[2]},\dots \kbd{z[lg(z)-1]} are ordinary signed long1446integers. This type is used instead of a \typ{VEC} of \typ{INT}s for1447efficiency reasons, for instance to implement efficiently permutations,1448polynomial arithmetic and linear algebra over small finite fields, etc.14491450\subsec{Type \typ{STR} (character string)}%1451\kbdsidx{t_STR}\sidx{character string}14521453\fun{char *}{GSTR}{z} (= \kbd{(z+1)}) points to the first character of the1454(\kbd{NULL}-terminated) string.14551456\subsec{Type \typ{ERROR} (error context)}\kbdsidx{t_ERROR}1457This type holds error messages, as well as details about the error, as1458returned by the exception handling system. The second codeword \kbd{z[1]}1459contains the error type (an \kbd{int}, as passed to \tet{pari_err}). The1460subsequent words \kbd{z[2]},\dots \kbd{z[lg(z)-1]} are \kbd{GEN}s containing1461additional data, depending on the error type.\sidx{error context}14621463\subsec{Type \typ{CLOSURE} (closure)}\kbdsidx{t_CLOSURE}\sidx{closure}1464This type holds GP functions and closures, in compiled form. The internal detail1465of this type is subject to change each time the GP language evolves. Hence we1466do not describe it here and refer to the Developer's Guide. However1467functions to create or to evaluate \typ{CLOSURE}s are documented in1468\secref{se:closure}.14691470\fun{long}{closure_arity}{GEN C} returns the arity of the \typ{CLOSURE}.14711472\fun{long}{closure_is_variadic}{GEN C} returns $1$ if the closure \kbd{C} is1473variadic, $0$ else.14741475\subsec{Type \typ{INFINITY} (infinity)}\kbdsidx{t_INFINITY}\sidx{infinity}14761477This type has a single \typ{INT} component, which is either $1$ or $-1$,1478corresponding to $+\infty$ and $-\infty$ respectively.14791480\fun{GEN}{mkmoo}{} returns $-\infty$14811482\fun{GEN}{mkoo}{} returns $\infty$14831484\fun{long}{inf_get_sign}{GEN x} returns $1$ if $x$ is $+\infty$, and $-1$1485if $x$ is $-\infty$.14861487\subsec{Type \typ{LIST} (list)}\kbdsidx{t_LIST}\sidx{list}1488this type was introduced for specific \kbd{gp} use and is rather inefficient1489compared to a straightforward linked list implementation (it requires more1490memory, as well as many unnecessary copies). Hence we do not describe it1491here and refer to the Developer's Guide.14921493\misctitle{Implementation note} For the types including an exponent (or a1494valuation), we actually store a biased nonnegative exponent (bit-ORing the1495biased exponent to the codeword), obtained by adding a constant to the true1496exponent: either \kbd{HIGHEXPOBIT} (for \typ{REAL}) or \kbd{HIGHVALPBIT} (for1497\typ{PADIC} and \typ{SER}). Of course, this is encapsulated by the1498exponent/valuation-handling macros and needs not concern the library user.14991500\section{PARI variables}\label{se:vars}1501\subsec{Multivariate objects}1502\sidx{variable (priority)}15031504\noindent We now consider variables and formal computations. As we have seen1505in \secref{se:impl}, the codewords for types \typ{POL} and \typ{SER} encode a1506``variable number''. This is an integer, ranging from $0$ to \kbd{MAXVARN}.1507Relative priorities may be ascertained using15081509\fun{int}{varncmp}{long v, long w}15101511\noindent which is $>0$, $=0$, $<0$ whenever $v$ has lower, resp.~same,1512resp.~higher priority than $w$.15131514The way an object is considered in formal computations depends entirely on1515its ``principal variable number'' which is given by the function15161517\fun{long}{gvar}{GEN z}15181519\noindent which returns a \idx{variable number} for \kbd{z}, even if \kbd{z}1520is not a polynomial or power series. The variable number of a scalar type is1521set by definition equal to \tet{NO_VARIABLE} which has lower priority than any1522valid variable number. The variable number of a recursive type which is not a1523polynomial or power series is the variable number with highest priority among1524its components. But for polynomials and power series only the ``outermost''1525number counts (we directly access $\tet{varn}(x)$ in the codewords): the1526representation is not symmetrical at all.15271528Under \kbd{gp}, one needs not worry too much since the interpreter defines1529the variables as it sees them\footnote{*}{The first time a given identifier1530is read by the GP parser a new variable is created, and it is assigned a1531strictly lower priority than any variable in use at this point. On startup,1532before any user input has taken place, 'x' is defined in this way and has1533initially maximal priority (and variable number $0$).}1534%1535and do the right thing with the polynomials produced.15361537But in library mode, they are tricky objects if you intend to build1538polynomials yourself (and not just let PARI functions produce them, which is1539less efficient). For instance, it does not make sense to have a variable1540number occur in the components of a polynomial whose main variable has a1541lower priority, even though PARI cannot prevent you from doing it.15421543\subsec{Creating variables} A basic difficulty is to ``create'' a variable.1544Some initializations are needed before you can use a given integer $v$ as a1545variable number.15461547Initially, this is done for $0$ and $1$ (the variables \kbd{x} and1548\kbd{y} under \kbd{gp}), and $2,\dots,9$ (printed as \kbd{t2}, \dots1549\kbd{t9}), with decreasing priority.15501551\subsubsec{User variables}\sidx{variable (user)} When the program starts,1552\kbd{x} (number~$0$) and \kbd{y} (number~$1$) are the only available1553variables, numbers $2$ to $9$ (decreasing priority) are reserved for building1554polynomials with predictable priorities.15551556To define further ones, you may use15571558\fun{GEN}{varhigher}{const char *s}15591560\fun{GEN}{varlower}{const char *s}15611562to recover a monomial of degree $1$ in a new variable, which is guaranteed1563to have higer (resp.~lower) priority than all existing ones at the1564time of the function call. The variable is printed as $s$, but is not part1565of GP's interpreter: it is not a symbol bound to a value.15661567On the other hand15681569\fun{long}{fetch_user_var}{char *s}: inspects the user variable whose name is1570the string pointed to by \kbd{s}, creating it if needed, and returns its1571variable number.1572\bprog1573long v = fetch_user_var("y");1574GEN gy = pol_x(v);1575@eprog\noindent1576The function raises an exception if the name is already in use for an1577\tet{install}ed or built-in function, or an alias. This function1578is mostly useless since it returns a variable with unpredictable priority.1579Don't use it to create new variables.15801581\misctitle{Caveat} You can use \tet{gp_read_str}1582(see~\secref{se:gp_read_str}) to execute a GP command and create GP1583variables on the fly as needed:1584\bprog1585GEN gy = gp_read_str("'y"); /*@Ccom returns \kbd{pol\_x}($v$), for some $v$ */1586long v = varn(gy);1587@eprog\noindent1588But please note the quote \kbd{'y} in the above. Using \kbd{gp\_read\_str("y")}1589might work, but is dangerous, especially when programming functions to1590be used under \kbd{gp}. The latter reads the value of \kbd{y}, as1591\emph{currently} known by the \kbd{gp} interpreter, possibly creating it1592in the process. But if \kbd{y} has been modified by previous \kbd{gp}1593commands (e.g.~\kbd {y = 1}), then the value of \kbd{gy} is not what you1594expected it to be and corresponds instead to the current value of the1595\kbd{gp} variable (e.g.~\kbd{gen\_1}).15961597\fun{GEN}{fetch_var_value}{long v} returns a shallow copy of the current1598value of the variable numbered $v$. Returns \kbd{NULL} if that variable1599number is unknown to the interpreter, e.g. it is a user variable. Note1600that this may not be the same as \kbd{pol\_x(v)} if assignments have been1601performed in the interpreter.16021603\subsubsec{Temporary variables}\sidx{variable (temporary)}1604You can create temporary variables using16051606\fun{long}{fetch_var}{}1607returns a new variable with \emph{lower} priority than any variable currently1608in use.16091610\fun{long}{fetch_var_higher}{}1611returns a new variable with \emph{higher} priority than any variable1612currently in use.16131614\noindent1615After the statement \kbd{v = fetch\_var()}, you can use1616\kbd{pol\_1(v)} and \kbd{pol\_x(v)}. The variables created in this way have no1617identifier assigned to them though, and are printed as1618\kbd{t\text{number}}. You can assign a name to a temporary variable, after1619creating it, by calling the function16201621\fun{void}{name_var}{long n, char *s}16221623\noindent after which the output machinery will use the name \kbd{s} to1624represent the variable number~\kbd{n}. The GP parser will \emph{not}1625recognize it by that name, however, and calling this on a variable known1626to~\kbd{gp} raises an error. Temporary variables are meant to be used as free1627variables to build polynomials and power series, and you should never assign1628values or functions to them as you would do with variables under~\kbd{gp}.1629For that, you need a user variable.16301631All objects created by \kbd{fetch\_var} are on the heap and not on the stack,1632thus they are not subject to standard garbage collecting (they are not1633destroyed by a \kbd{gerepile} or \kbd{avma = ltop} statement). When you do1634not need a variable number anymore, you can delete it using16351636\fun{long}{delete_var}{}16371638\noindent which deletes the \emph{latest} temporary variable created and1639returns the variable number of the previous one (or simply returns 0 if none1640remain). Of course you should make sure that1641the deleted variable does not appear anywhere in the objects you use later1642on. Here is an example:16431644\bprog1645long first = fetch_var();1646long n1 = fetch_var();1647long n2 = fetch_var(); /*@Ccom prepare three variables for internal use */1648...1649/*@Ccom delete all variables before leaving */1650do { num = delete_var(); } while (num && num <= first);1651@eprog\noindent1652The (dangerous) statement16531654\bprog1655while (delete_var()) /*@Ccom empty */;1656@eprog\noindent1657removes all temporary variables in use.16581659\subsec{Comparing variables}16601661Let us go back to \kbd{varncmp}. There is an interesting corner case, when1662one of the compared variables (from \kbd{gvar}, say) is \kbd{NO\_VARIABLE}.1663In this case, \kbd{varncmp} declares it has lower priority than any other1664variable; of course, comparing \kbd{NO\_VARIABLE} with itself yields1665$0$ (same priority);16661667In addition to \kbd{varncmp} we have16681669\fun{long}{varnmax}{long v, long w} given two variable numbers (possibly1670\kbd{NO\_VARIABLE}), returns the variable with the highest priority.1671This function always returns a valid variable number unless it is comparing1672\kbd{NO\_VARIABLE} to itself.16731674\fun{long}{varnmin}{long x, long y} given two variable numbers (possibly1675\kbd{NO\_VARIABLE}), returns the variable with the lowest priority. Note1676that when comparing a true variable with \kbd{NO\_VARIABLE}, this function1677returns \kbd{NO\_VARIABLE}, which is not a valid variable number.16781679\section{Input and output}16801681\noindent1682Two important aspects have not yet been explained which are specific to1683library mode: input and output of PARI objects.16841685\subsec{Input}16861687\noindent1688For \idx{input}, PARI provides several powerful high level functions1689which enable you to input your objects as if you were under \kbd{gp}. In fact,1690it \emph{is} essentially the GP syntactical parser.16911692There are two similar functions available to parse a string:16931694\fun{GEN}{gp_read_str}{const char *s}\label{se:gp_read_str}16951696\fun{GEN}{gp_read_str_multiline}{const char *s, char *last}16971698\noindent1699Both functions read the whole string \kbd{s}. The function1700\kbd{gp\_read\_str} ignores newlines: it assumes that the input is one1701expression and returns the result of this expression.17021703The function \kbd{gp\_read\_str\_multiline} processes the text in the1704same way as the GP command \tet{read}: newlines are significant and can1705be used to separate expressions.1706The return value is that of the last nonempty expression evaluated.17071708In \kbd{gp\_read\_str\_multiline}, if \kbd{last} is not \kbd{NULL},1709then \kbd{*last} receives the last character from the \emph{filtered}1710input: this can be used to check if the last character was a semi-colon1711(to hide the output in interactive usage). If (and only if) the1712input contains no statements, then \kbd{*last} is set to \kbd{0}.17131714For both functions, \kbd{gp}'s metacommands \emph{are} recognized.17151716Two variants allow to specify a default precision while evaluating the1717string:17181719\fun{GEN}{gp_read_str_prec}{const char *s, long prec}1720As \kbd{gp\_read\_str}, but set the precision to \kbd{prec} words while evaluating $s$.17211722\fun{GEN}{gp_read_str_bitprec}{const char *s, long bitprec}1723As \kbd{gp\_read\_str}, but set the precision to \kbd{bitprec} bits while evaluating $s$.17241725\misctitle{Note} The obsolete form17261727\fun{GEN}{readseq}{char *t}17281729still exists for backward compatibility (assumes filtered input, without1730spaces or comments). Don't use it.17311732To read a \kbd{GEN} from a file, you can use the simpler interface17331734\fun{GEN}{gp_read_stream}{FILE *file}17351736\noindent which reads a character string of arbitrary length from the stream1737\kbd{file} (up to the first complete expression sequence), applies1738\kbd{gp\_read\_str} to it, and returns the resulting \kbd{GEN}. This way, you1739do not have to worry about allocating buffers to hold the string. To1740interactively input an expression, use \kbd{gp\_read\_stream(stdin)}.1741Return \kbd{NULL} when there are no more expressions to read (we reached1742EOF).17431744Finally, you can read in a whole file, as in GP's \tet{read} statement17451746\fun{GEN}{gp_read_file}{char *name}17471748\noindent As usual, the return value is that of the last nonempty expression1749evaluated. There is one technical exception: if \kbd{name} is a \emph{binary}1750file (from \tet{writebin}) containing more than one object, a \typ{VEC}1751containing them all is returned. This is because binary objects bypass the1752parser, hence reading them has no useful side effect.17531754\subsec{Output to screen or file, output to string}\sidx{output}17551756General output functions return nothing but print a character string as a1757side effect. Low level routines are available to write on PARI output stream1758\tet{pari_outfile} (\tet{stdout} by default):17591760\fun{void}{pari_putc}{char c}: write character \kbd{c} to the output stream.17611762\fun{void}{pari_puts}{char *s}: write \kbd{s} to the output stream.17631764\fun{void}{pari_flush}{}: flush output stream; most streams are buffered by1765default, this command makes sure that all characters output so are actually1766written.17671768\fun{void}{pari_printf}{const char *fmt, ...}: the most versatile such1769function. \kbd{fmt} is a character string similar to the one1770\tet{printf} uses. In there, \kbd{\%} characters have a special meaning, and1771describe how to print the remaining operands. In addition to the standard1772format types (see the GP function \tet{printf}), you can use the \emph{length1773modifier}~\kbd{P} (for PARI of course!) to specify that an argument is a1774\kbd{GEN}. For instance, the following are valid conversions for a \kbd{GEN}1775argument1776\bprog1777%Ps @com convert to \kbd{char*} (will print an arbitrary \kbd{GEN})1778%P.10s @com convert to \kbd{char*}, truncated to 10 chars1779%P.2f @com convert to floating point format with 2 decimals1780%P4d @com convert to integer, field width at least 417811782pari_printf("x[%d] = %Ps is not invertible!\n", i, gel(x,i));1783@eprog\noindent1784Here \kbd{i} is an \kbd{int}, \kbd{x} a \kbd{GEN} which is not a leaf1785(presumably a vector, or a polynomial) and this would insert the value of its1786$i$-th \kbd{GEN} component: \kbd{gel(x,i)}.17871788\noindent Simple but useful variants to \kbd{pari\_printf} are17891790\fun{void}{output}{GEN x} prints \kbd{x} in raw format, followed by a1791newline and a buffer flush. This is more or less equivalent to1792\bprog1793pari_printf("%Ps\n", x);1794pari_flush();1795@eprog17961797\fun{void}{outmat}{GEN x} as above except if $x$ is a \typ{MAT}, in which1798case a multi-line display is used to display the matrix. This is prettier for1799small dimensions, but quickly becomes unreadable and cannot be pasted and1800reused for input. If all entries of $x$ are small integers, you may use the1801recursive features of \kbd{\%Pd} and obtain the same (or better) effect with1802\bprog1803pari_printf("%Pd\n", x);1804pari_flush();1805@eprog\noindent1806A variant like \kbd{"\%5Pd"} would improve alignment by imposing18075 chars for each coefficient. Similarly if all entries are to be converted to1808floats, a format like \kbd{"\%5.1Pf"} could be useful.180918101811These functions write on (PARI's idea of) standard output, and must be used1812if you want your functions to interact nicely with \kbd{gp}. In most1813programs, this is not a concern and it is more flexible to write to an1814explicit \kbd{FILE*}, or to recover a character string:18151816\fun{void}{pari_fprintf}{FILE *file, const char *fmt, ...} writes the1817remaining arguments to stream \kbd{file} according to the format1818specification \kbd{fmt}.18191820\fun{char*}{pari_sprintf}{const char *fmt, ...} produces a string from the1821remaining arguments, according to the PARI format \kbd{fmt} (see \tet{printf}).1822This is the \kbd{libpari} equivalent of \tet{strprintf}, and returns a1823\kbd{malloc}'ed string, which must be freed by the caller. Note that contrary1824to the analogous \tet{sprintf} in the \kbd{libc} you do not provide a buffer1825(leading to all kinds of buffer overflow concerns); the function provided is1826actually closer to the GNU extension \kbd{asprintf}, although the latter has1827a different interface.18281829Simple variants of \tet{pari_sprintf} convert a \kbd{GEN} to a1830\kbd{malloc}'ed ASCII string, which you must still \kbd{free} after use:18311832\fun{char*}{GENtostr}{GEN x}, using the current default output format1833(\kbd{prettymat} by default).18341835\fun{char*}{GENtoTeXstr}{GEN x}, suitable for inclusion in a \TeX\ file.183618371838Note that we have \tet{va_list} analogs of the functions of \kbd{printf} type1839seen so far:18401841\fun{void}{pari_vprintf}{const char *fmt, va_list ap}18421843\fun{void}{pari_vfprintf}{FILE *file, const char *fmt, va_list ap}18441845\fun{char*}{pari_vsprintf}{const char *fmt, va_list ap}18461847\subsec{Errors}\sidx{error}\kbdsidx{e_MISC}18481849\noindent1850If you want your functions to issue error messages, you can use the general1851error handling routine \tet{pari_err}. The basic syntax is1852%1853\bprog1854pari_err(e_MISC, "error message");1855@eprog\noindent1856This prints the corresponding error message and exit the program (in1857library mode; go back to the \kbd{gp} prompt otherwise).\label{se:err} You can1858also use it in the more versatile guise1859\bprog1860pari_err(e_MISC, format, ...);1861@eprog\noindent1862where \kbd{format} describes the format to use to write the remaining1863operands, as in the \tet{pari_printf} function. For instance:1864\bprog1865pari_err(e_MISC, "x[%d] = %Ps is not invertible!", i, gel(x,i));1866@eprog\noindent1867The simple syntax seen above is just a special case with a constant format1868and no remaining arguments. The general syntax is18691870\fun{void}{pari_err}{numerr,...}18711872\noindent where \kbd{numerr} is a codeword which specifies the error class1873and what to do with the remaining arguments and what message to print.1874For instance, if $x$ is a \kbd{GEN} with internal type \typ{STR}, say,1875\kbd{pari\_err(e\_TYPE,"extgcd", $x$)} prints the message:1876\bprog1877*** incorrect type in extgcd (t_STR),1878@eprog\noindent See \secref{se:errors} for details. In the libpari code1879itself, the general-purpose \kbd{e\_MISC} is used sparingly: it is so1880flexible that the corresponding error contexts (\typ{ERROR}) become hard to1881use reliably. Other more rigid error types are generally more useful: for1882instance the error context attached to the \kbd{e\_TYPE} exception above is1883precisely documented and contains \kbd{"extgcd"} and $x$ (not only its type)1884as readily available components.18851886\subsec{Warnings}18871888\noindent To issue a warning, use18891890\fun{void}{pari_warn}{warnerr,...}1891In that case, of course, we do \emph{not} abort the computation, just print1892the requested message and go on. The basic example is1893%1894\bprog1895pari_warn(warner, "Strategy 1 failed. Trying strategy 2")1896@eprog\noindent1897which is the exact equivalent of \kbd{pari\_err(e\_MISC,...)} except that1898you certainly do not want to stop the program at this point, just inform the1899user that something important has occurred; in particular, this output would be1900suitably highlighted under \kbd{gp}, whereas a simple \kbd{printf} would not.19011902The valid \emph{warning} keywords are \tet{warner} (general), \tet{warnprec}1903(increasing precision), \tet{warnmem} (garbage collecting) and \tet{warnfile}1904(error in file operation), used as follows:1905\bprog1906pari_warn(warnprec, "bnfinit", newprec);1907pari_warn(warnmem, "bnfinit");1908pari_warn(warnfile, "close", "afile"); /* error when closing "afile" */1909@eprog19101911\subsec{Debugging output}\sidx{debugging}\sidx{format}\label{se:dbg_output}19121913For debugging output, you can use the standard output1914functions, \tet{output} and \tet{pari_printf} mainly. Corresponding to the1915\kbd{gp} metacommand \kbd{\b x}, you can also output the \idx{hexadecimal1916tree} attached to an object:19171918\fun{void}{dbgGEN}{GEN x, long nb = -1}, displays the recursive structure of1919\kbd{x}. If $\kbd{nb} = -1$, the full structure is printed, otherwise1920the leaves (nonrecursive components) are truncated to \kbd{nb} words.19211922\noindent The function \tet{output} is vital under debuggers, since none of1923them knows how to print PARI objects by default. Seasoned PARI developers1924add the following \kbd{gdb} macro to their \kbd{.gdbinit}:1925\bprog1926define oo1927call output((GEN)$arg0)1928end1929define xx1930call dbgGEN($arg0,-1)1931end1932@eprog\noindent1933Typing \kbd{i x} at a breakpoint in \kbd{gdb} then prints the value of the1934\kbd{GEN} \kbd{x} (provided the optimizer has not put it into a register, but1935it is rarely a good idea to debug optimized code).19361937\noindent1938The global variables \teb{DEBUGLEVEL} and \teb{DEBUGMEM} (corresponding1939to the default \teb{debug} and \teb{debugmem})1940are used throughout the PARI code to govern the amount of diagnostic and1941debugging output, depending on their values. You can use them to debug your1942own functions, especially if you \tet{install} the latter under \kbd{gp}.19431944\fun{void}{dbg_pari_heap}{void} print debugging statements about the PARI1945stack, heap, and number of variables used. Corresponds to \kbd{\bs s}1946under gp.19471948\subsec{Timers and timing output}19491950\noindent To handle timings in a reentrant way, PARI defines a dedicated data1951type, \tet{pari_timer}, together with the following methods:19521953\fun{void}{timer_start}{pari_timer *T} start (or reset) a timer.19541955\fun{long}{timer_delay}{pari_timer *T} returns the number of milliseconds1956elapsed since the timer was last reset. Resets the timer as a side effect.1957Assume $T$ was started by \kbd{timer\_start}.19581959\fun{long}{timer_get}{pari_timer *T} returns the number of milliseconds1960elapsed since the timer was last reset. Does \emph{not} reset the timer.1961Assume $T$ was started by \kbd{timer\_start}.19621963\fun{void}{walltimer_start}{pari_timer *T} start a timer, as if it1964had been started at the Unix epoch (see \kbd{getwalltime}).19651966\fun{long}{walltimer_delay}{pari_timer *T} returns the number of milliseconds1967elapsed since the timer was last checked. Assume $T$ was started by1968\kbd{walltimer\_start}.19691970\fun{long}{walltimer_get}{pari_timer *T} returns the number of milliseconds1971elapsed since the timer was last reset. Does \emph{not} reset the timer.1972Assume $T$ was started by \kbd{walltimer\_start}.19731974\fun{long}{timer_printf}{pari_timer *T, char *format,...} This diagnostics1975function is equivalent to the following code1976\bprog1977err_printf("Time ")1978... prints remaining arguments according to format ...1979err_printf(": %ld", timer_delay(T));1980@eprog\noindent Resets the timer as a side effect.19811982\noindent They are used as follows:1983\bprog1984pari_timer T;1985timer_start(&T); /* initialize timer */1986...1987printf("Total time: %ldms\n", timer_delay(&T));1988@eprog\noindent1989or1990\bprog1991pari_timer T;1992timer_start(&T);1993for (i = 1; i < 10; i++) {1994...1995timer_printf(&T, "for i = %ld (L[i] = %Ps)", i, gel(L,i));1996}1997@eprog19981999The following functions provided the same functionality, in a2000nonreentrant way, and are now deprecated.20012002\fun{long}{timer}{void}20032004\fun{long}{timer2}{void}20052006\fun{void}{msgtimer}{const char *format, ...}20072008The following function implements \kbd{gp}'s timer and should not be used in2009libpari programs:2010\fun{long}{gettime}{void} equivalent to \tet{timer_delay}$(T)$ attached to a2011private timer $T$.20122013\section{Iterators, Numerical integration, Sums, Products}2014\subsec{Iterators}2015Since it is easier to program directly simple loops in library mode, some GP2016iterators are mainly useful for GP programming. Here are the others:20172018\item \tet{fordiv} is a trivial iteration over a list produced by2019\tet{divisors}.20202021\item \tet{forell}, \tet{forqfvec} and \tet{forsubgroup} are currently not2022implemented as an iterator but as a procedure with callbacks.20232024\fun{void}{forell}{void *E, long fun(void*, GEN), GEN a, GEN b, long flag}2025goes through the same curves as \tet{forell(ell,a,b,,flag)}, calling2026\tet{fun(E, ell)} for each curve \kbd{ell}, stopping if \kbd{fun} returns a2027nonzero value.20282029\fun{void}{forqfvec}{void *E, long (*fun)(void *, GEN, GEN, double), GEN q, GEN b}:2030Evaluate \kbd{fun(E,U,v,m)} on all $v$ such that $q(U\*v)<b$, where $U$ is a2031\typ{MAT}, $v$ is a \typ{VECSMALL} and $m=q(v)$ is a C double. The function2032\kbd{fun} must return $0$, unless \kbd{forqfvec} should stop, in which case,2033it should return $1$.20342035\fun{void}{forqfvec1}{void *E, long (*fun)(void *, GEN), GEN q, GEN b}:2036Evaluate \kbd{fun(E,v)} on all $v$ such that $q(v)<b$, where $v$ is a2037\typ{COL}. The function \kbd{fun} must return $0$, unless \kbd{forqfvec}2038should stop, in which case, it should return $1$.20392040\fun{void}{forsubgroup}{void *E, long fun(void*, GEN), GEN G, GEN B}2041goes through the same subgroups as \tet{forsubgroup(H = G, B,)}, calling2042\tet{fun(E, H)} for each subgroup $H$, stopping if \kbd{fun} returns a2043nonzero value.20442045\item \tet{forprime} and \tet{forprimestep}, iterators over primes and2046primes in arithmetic progressions, for which we refer you to the2047next subsection.20482049\item \tet{forcomposite}, we provide an iterator over composite integers:20502051\fun{int}{forcomposite_init}{forcomposite_t *T, GEN a, GEN b} initialize an2052iterator $T$ over composite integers in $[a,b]$; over composites $\geq a$ if2053$b = \kbd{NULL}$. We must have $a\geq 0$. Return $0$ if the range is known to2054be empty from the start (as if $b < a$ or $b < 0$), and return $1$ otherwise.20552056\fun{GEN}{forcomposite_next}{forcomposite_t *T} returns the next composite in2057the range, assuming that $T$ was initialized by \tet{forcomposite_init}.20582059\item \tet{forvec}, for which we provide a convenient iterator. To2060initialize the analog of \kbd{forvec(X = v, ..., flag)}, call20612062\fun{int}{forvec_init}{forvec_t *T, GEN v, long flag}2063initialize an iterator $T$ over the vectors generated by2064\kbd{forvec(X = $v$,..., flag)}. This returns $0$ if this vector list is2065empty, and $1$ otherwise.20662067\fun{GEN}{forvec_next}{forvec_t *T} returns the next element in the2068\kbd{forvec} sequence, or \kbd{NULL} if we are done. The return value must be2069used immediately or copied since the next call to the iterator destroys it:2070the relevant vector is updated in place. The iterator works hard to not2071use up PARI stack, and is more efficient when all lower bounds in the2072initialization vector $v$ are integers. In that case, the cost is linear in2073the number of tuples enumerated, and you can expect to run over more than2074$10^9$ tuples per minute. If speed is critical and all integers involved2075would fit in $C$ \kbd{long}s, write a simple direct backtracking algorithm2076yourself.20772078\item \tet{forpart} is a variant of \kbd{forvec} which iterates over2079partitions. See the documentation of the \kbd{forpart} GP function for2080details. This function is available as a loop with callbacks:20812082\fun{void}{forpart}{void *data, long (*call)(void*,GEN), long k, GEN a, GEN n}20832084\noindent It is also available as an iterator:20852086\fun{void}{forpart_init}{forpart_t *T, long k, GEN a, GEN n} initializes an2087iterator over the partitions of $k$, with length restricted by $n$,2088and components restricted by $a$, either of which can be set to \kbd{NULL}2089to run without restriction.20902091\fun{GEN}{forpart_next}{forpart_t *T} returns the next partition, or2092\kbd{NULL} when all partitions have been exhausted.20932094\fun{GEN}{forpart_prev}{forpart_t *T} returns the previous partition, or2095\kbd{NULL} when all partitions have been exhausted.20962097In both cases, the partition must be used or copied before the next call2098since it is returned from a state array which will be modified in place.2099You may \emph{not} mix calls to \tet{forpart_next} and \tet{forpart_prev}:2100the first one called determines the ordering used to iterate over the2101partitions; you can not go back since the \tet{forpart_t} structure is used2102in incompatible ways.21032104\item \tet{forperm} to loop over permutations of $k$. See the documentation2105of the \kbd{forperm} GP function for details. This function is available as2106an iterator:21072108\fun{void}{forperm_init}{forperm_t *T, GEN k} initializes an iterator2109over the permutations of $k$ (\typ{INT}, \typ{VEC} or \typ{VECSMALL}).21102111\fun{GEN}{forperm_next}{forperm_t *T} returns the next permutation2112as a \typ{VECSMALL} or \kbd{NULL} whell all permutations have been2113exhausted. The permutation must be used or copied before the next call2114since it is returned from a state array which will be modified in place.21152116\item \tet{forsubset} to loop over subsets. See the documentation2117of the \kbd{forsubset} GP function for details. This function2118is available as two iterators:21192120\fun{void}{forallsubset_init}{forsubset_t *T, long n}21212122\fun{void}{forksubset_init}{forsubset_t *T, long n, long k}21232124\noindent It is also available in generic form:21252126\fun{void}{forsubset_init}{forsubset_t *T, GEN nk} where \kbd{nk} is either2127a \typ{INT} $n$ or a \typ{VEC} with two integral components $[n,k]$.21282129In all three cases, \fun{GEN}{forsubset_next}{forsubset_t *T} returns the2130next subset as a \typ{VECSMALL} or \kbd{NULL} when all subsets have been2131exhausted.21322133\subsec{Iterating over primes}\label{se:primeiter}21342135The library provides a high-level iterator, which stores its (private) data2136in a \kbd{struct} \tet{forprime_t} and runs over arbitrary ranges of primes,2137without ever overflowing.21382139The iterator has two flavors, one providing the successive primes as2140\kbd{ulong}s, the other as \kbd{GEN}. They are initialized as follows,2141where we expect to run over primes $\geq a$ and $\leq b$:21422143\fun{int}{u_forprime_init}{forprime_t *T, ulong a, ulong b} for the2144\kbd{ulong} variant, where $b = \kbd{ULONG\_MAX}$ means we will run through2145all primes representable in a \kbd{ulong} type.21462147\fun{int}{forprime_init}{forprime_t *T, GEN a, GEN b} for the \kbd{GEN}2148variant, where $b = \kbd{NULL}$ means $+\infty$.21492150\fun{int}{forprimestep_init}{forprime_t *T, GEN a, GEN b, GEN q} initialize2151an iterator $T$ over primes in an arithmetic progression, $p \geq a$ and $p2152\leq b$ (where $b = \kbd{NULL}$ means $+\infty$). The argument $q$ is2153either a \typ{INT} ($p \equiv a \pmod{q}$) or a \typ{INTMOD} \kbd{Mod(c,N)}2154and we restrict to that congruence class.21552156All variants return $1$ on success, and $0$ if the iterator would run over an2157empty interval (if $a > b$, for instance). They allocate the \tet{forprime_t}2158data structure on the PARI stack.21592160\noindent The successive primes are then obtained using21612162\fun{GEN}{forprime_next}{forprime_t *T}, returns \kbd{NULL} if no more primes2163are available in the interval and the next suitable prime as a \typ{INT}2164otherwise.21652166\fun{ulong}{u_forprime_next}{forprime_t *T}, returns $0$ if no more primes2167are available in the interval and fitting in an \kbd{ulong} and the next2168suitable prime otherwise.21692170These two functions leave alone the PARI stack, and write their state2171information in the preallocated \tet{forprime_t} struct. The typical usage is2172thus:2173\bprog2174forprime_t T;2175GEN p;2176pari_sp av = avma, av2;21772178forprime_init(&T, gen_2, stoi(1000));2179av2 = avma;2180while ( (p = forprime_next(&T)) )2181{2182...2183if ( prime_is_OK(p) ) break;2184avma = av2; /* delete garbage accumulated in this iteration */2185}2186avma = av; /* delete all */2187@eprog\noindent Of course, the final \kbd{avma = av} could be replaced2188by a \kbd{gerepile} call. Beware that swapping the2189\kbd{av2 = avma} and \tet{forprime_init} call would be incorrect: the2190first \kbd{avma = av2} would delete the \tet{forprime_t} structure!21912192\subsec{Parallel iterators}21932194Theses iterators loops over the values of a \typ{CLOSURE}2195taken at some data, where the evaluations are done in parallel.21962197\item \tet{parfor}. To initialize the analog of2198\kbd{parfor(i = a, b, ...)}, call21992200\fun{void}{parfor_init}{parfor_t *T, GEN a, GEN b, GEN code}2201initialize an iterator over the evaluation of \kbd{code} on the integers2202between $a$ and $b$.22032204\fun{GEN}{parfor_next}{parfor_t *T} returns a \typ{VEC} \kbd{[i,code(i)]} where2205$i$ is one of the integers and \kbd{code(i)} is the evaluation, \kbd{NULL} when2206all data have been exhausted. Once it happens, \kbd{parfor\_next} must not be2207called anymore with the same initialization.22082209\fun{void}{parfor_stop}{parfor_t *T} needs to be called when leaving the2210iterator before \kbd{parfor\_next} returned \kbd{NULL}.22112212The following returns an integer $1\leq i\leq N$ such that \kbd{fun(i)} is not2213zero, or \kbd{NULL}.2214\bprog2215GEN2216parfirst(GEN fun, GEN N)2217{2218parfor_t T;2219GEN e;2220parfor_init(&T, gen_1, N, fun);2221while ((e = parfor_next(&T)))2222{2223GEN i = gel(e,1), funi = gel(e,2);2224if (!gequal0(funi))2225{ /* found: stop the iterator and return the index */2226parfor_stop(&T);2227return i;2228}2229}2230return NULL; /* not found */2231}2232@eprog22332234\item \tet{parforeach}. To initialize the analog of2235\kbd{parforeach(V, X, ...)}, call22362237\fun{void}{parforeach_init}{parforeach_t *T, GEN V, GEN code}2238initialize an iterator over the evaluation of \kbd{code} on the components2239of $V$.22402241\fun{GEN}{parforeach_next}{parforeach_t *T} returns a \typ{VEC}2242\kbd{[V[i],code(V[i])]} where $V[i]$ is one of the components of $V$ and \kbd{code(V[i])} is2243the evaluation, \kbd{NULL} when all data have been exhausted. Once it happens,2244\kbd{parforprime\_next}2245must not be called anymore with the same initialization.22462247\fun{void}{parforeach_stop}{parforeach_t *T} needs to be called when leaving2248the iterator before \kbd{parforeach\_next} returned \kbd{NULL}.22492250\item \tet{parforprime}. To initialize the analog of2251\kbd{parforprime(p = a, b, ...)}, call22522253\fun{void}{parforprime_init}{parforprime_t *T, GEN a, GEN b, GEN code}2254initialize an iterator over the evaluation of \kbd{code} on the prime numbers2255between $a$ and $b$.22562257\item \tet{parforprimestep}. To initialize the analog of2258\kbd{parforprimestep(p = a, b, q, ...)}, call22592260\fun{void}{parforprimestep_init}{parforprime_t *T, GEN a, GEN b, GEN q, GEN code}2261initialize an iterator over the evaluation of \kbd{code} on the prime numbers2262between $a$ and $b$ in the congruence class defined by $q$.22632264\fun{GEN}{parforprime_next}{parforprime_t *T} returns a \typ{VEC}2265\kbd{[p,code(p)]} where $p$ is one of the prime numbers and \kbd{code(p)} is2266the evaluation, \kbd{NULL} when all data have been exhausted. Once it happens,2267\kbd{parforprime\_next}2268must not be called anymore with the same initialization.22692270\fun{void}{parforprime_stop}{parforprime_t *T} needs to be called when leaving2271the iterator before \kbd{parforprime\_next} returned \kbd{NULL}.22722273\item \tet{parforvec}. To initialize the analog of2274\kbd{parforvec(X = V, ..., flag)}, call22752276\fun{void}{parforvec_init}{parforvec_t *T, GEN V, GEN code, long flag}2277initialize an iterator over the evaluation of \kbd{code} on the vectors2278specified by \kbd{V} and \kbd{flag}, see \kbd{forvec} for detail.22792280\fun{GEN}{parforvec_next}{parforvec_t *T} returns a \typ{VEC} \kbd{[v,code(v)]}2281where $v$ is one of the vectors and \kbd{code(v)} is the evaluation, \kbd{NULL}2282when all data have been exhausted. Once it happens, \kbd{parforvec\_next}2283must not be called anymore with the same initialization.22842285\fun{void}{parforvec_stop}{parforvec_t *T} needs to be called when leaving the2286iterator before \kbd{parforvec\_next} returned \kbd{NULL}.22872288\subsec{Numerical analysis}22892290Numerical routines code a function (to be integrated, summed, zeroed, etc.)2291with two parameters named2292\bprog2293void *E;2294GEN (*eval)(void*, GEN)2295@eprog\noindent2296The second is meant to contain all auxiliary data needed by your function.2297The first is such that \kbd{eval(x, E)} returns your function evaluated at2298\kbd{x}. For instance, one may code the family of functions2299$f_t: x \to (x+t)^2$ via2300\bprog2301GEN fun(void *t, GEN x) { return gsqr(gadd(x, (GEN)t)); }2302@eprog\noindent2303One can then integrate $f_1$ between $a$ and $b$ with the call2304\bprog2305intnum((void*)stoi(1), &fun, a, b, NULL, prec);2306@eprog\noindent2307Since you can set \kbd{E} to a pointer to any \kbd{struct} (typecast to2308\kbd{void*}) the above mechanism handles arbitrary functions. For simple2309functions without extra parameters, you may set \kbd{E = NULL} and ignore2310that argument in your function definition.23112312\section{Catching exceptions}23132314\subsec{Basic use}23152316PARI provides a mechanism to trap exceptions generated via \kbd{pari\_err}2317using the \tet{pari_CATCH} construction. The basic usage is as follows2318\bprog2319pari_CATCH(err_code) {2320@com recovery branch2321}2322pari_TRY {2323@com main branch2324}2325pari_ENDCATCH2326@eprog\noindent This fragment executes the main branch, then the recovery2327branch \emph{if} exception \kbd{err\_code} is thrown, e.g. \kbd{e\_TYPE}.2328See \secref{se:errors} for the description of all error classes.2329The special error code \tet{CATCH_ALL} is available to catch all errors.23302331One can replace the \tet{pari_TRY} keyword by \tet{pari_RETRY}, in which case2332once the recovery branch is run, we run the main branch again, still catching2333the same exceptions.23342335\misctitle{Restrictions}23362337\item Such constructs can be nested without adverse effect, the innermost2338handler catching the exception.23392340\item It is \emph{valid} to leave either branch using \tet{pari_err}.23412342\item It is \emph{invalid} to use C flow control instructions (\kbd{break},2343\kbd{continue}, \kbd{return}) to directly leave either branch without seeing2344the \tet{pari_ENDCATCH} keyword. This would leave an invalid structure in the2345exception handler stack, and the next exception would crash.23462347\item In order to leave using \kbd{break}, \kbd{continue} or \kbd{return},2348one must precede the keyword by a call to23492350\fun{void}{pari_CATCH_reset}{} disable the current handler, allowing to leave2351without adverse effect.23522353\subsec{Advanced use}23542355In the recovery branch, the exception context can be examined via the2356following helper routines:23572358\fun{GEN}{pari_err_last}{} returns the exception context, as a \typ{ERROR}.2359The exception $E$ returned by \tet{pari_err_last} can be rethrown, using2360\bprog2361pari_err(0, E);2362@eprog23632364\fun{long}{err_get_num}{GEN E} returns the error symbolic name. E.g2365\kbd{e\_TYPE}.23662367\fun{GEN}{err_get_compo}{GEN E, long i} error $i$-th component, as documented2368in \secref{se:errors}.23692370\noindent For instance2371\bprog2372pari_CATCH(CATCH_ALL) { /* catch everything */2373GEN x, E = pari_err_last();2374long code = err_get_num(E);2375if (code != e_INV) pari_err(0, E); /* unexpected error, rethrow */2376x = err_get_compo(E, 2);2377/* e_INV has two components, 1: function name 2: noninvertible x */2378if (typ(x) != t_INTMOD) pari_err(0, E); /* unexpected type, rethrow */2379pari_CATCH_reset();2380return x; /* leave ! */2381@com @dots2382} pari_TRY {2383@com main branch2384}2385pari_ENDCATCH2386@eprog23872388\section{A complete program}2389\label{se:prog}23902391\noindent2392Now that the preliminaries are out of the way, the best way to learn how to2393use the library mode is to study a detailed example. We want to write a2394program which computes the gcd of two integers, together with the Bezout2395coefficients. We shall use the standard quadratic algorithm which is not2396optimal but is not too far from the one used in the PARI function2397\teb{bezout}.23982399Let $x,y$ two integers and initially2400$ \pmatrix{s_x & s_y \cr t_x & t_y } =2401\pmatrix{1 & 0 \cr 0 & 1}$, so that2402$$ \pmatrix{s_x & s_y \cr2403t_x & t_y }2404\pmatrix{x \cr y } =2405\pmatrix{x \cr y }.2406$$2407To apply the ordinary Euclidean algorithm to the right hand side,2408multiply the system from the left by2409$ \pmatrix{0 & 1 \cr 1 & -q }$,2410with $q = \kbd{floor}(x / y)$. Iterate until $y = 0$ in the right hand side,2411then the first line of the system reads2412$$ s_x x + s_y y = \gcd(x,y).$$2413In practice, there is no need to update $s_y$ and $t_y$ since2414$\gcd(x,y)$ and $s_x$ are enough to recover $s_y$. The following program2415is now straightforward. A couple of new functions appear in there, whose2416description can be found in the technical reference manual in Chapter 5,2417but whose meaning should be clear from their name and the context.24182419This program can be found in \kbd{examples/extgcd.c} together with a2420proper \kbd{Makefile}. You may ignore the first comment2421\bprog2422/*2423GP;install("extgcd", "GG&&", "gcdex", "./libextgcd.so");2424*/2425@eprog\noindent which instruments the program so that \kbd{gp2c-run extgcd.c}2426can import the \kbd{extgcd()} routine into an instance of the \kbd{gp}2427interpreter (under the name \kbd{gcdex}). See the \kbd{gp2c} manual for2428details.2429\newpage24302431\bprogfile{../examples/extgcd.c}24322433\noindent For simplicity, the inner loop does not include any garbage2434collection, hence memory use is quadratic in the size of the inputs instead2435of linear. Here is a better version of that loop:2436\bprog2437pari_sp av = avma;2438...2439while (!gequal0(b))2440{2441GEN r, q = dvmdii(a, b, &r), v = vx;24422443vx = subii(ux, mulii(q, vx));2444ux = v; a = b; b = r;2445if (gc_needed(av,1))2446gerepileall(av, 4, &a, &b, &ux, &vx);2447}2448@eprog2449\newpage245024512452