Annotation of OpenXM_contrib/pari-2.2/doc/usersch4.tex, Revision 1.2
1.2 ! noro 1: % $Id: usersch4.tex,v 1.36 2002/09/05 14:36:59 karim Exp $
1.1 noro 2: % Copyright (c) 2000 The PARI Group
3: %
4: % This file is part of the PARI/GP documentation
5: %
6: % Permission is granted to copy, distribute and/or modify this document
7: % under the terms of the GNU Free Documentation License
8: \chapter{Programming PARI in Library Mode}
9:
10: \section{Introduction: initializations, universal objects}
11: \label{se:intro4}
12:
13: \noindent
14: To be able to use PARI in \idx{library mode}, you must write a C program and
15: link it to the PARI library. See the installation guide (in Appendix~A)
16: on how to create and install the library and include files. A sample Makefile
17: is presented in Appendix~B.
18:
19: Probably the best way to understand how programming is done is to work
20: through a complete example. We will write such a program in
21: \secref{se:prog}. Before doing this, a few explanations are in order.
22:
23: First, one must explain to the outside world what kind of objects and
24: routines we are going to use. This is done simply with the statement
25:
26: \bprog
27: #include <pari.h>
28: @eprog
29: \noindent
30: This file \tet{pari.h} imports all the necessary constants,
31: variables and functions, defines some important macros, and also defines the
32: fundamental type for all PARI objects: the type \teb{GEN}, which is
33: simply a pointer to \kbd{long}.
34:
35: \misctitle{Technical note}: we would have liked to define a type \kbd{GEN}
36: to be a pointer to itself. This unfortunately is not possible in C, except by
37: using structures, but then the names become unwieldy. The result of this is
38: that when we use a component of a PARI object, it will be a \kbd{long},
39: hence will need to be typecast to a \kbd{GEN} again if we want to avoid
40: warnings from the compiler. This will sometimes be quite tedious, but of
41: course is trivially done. See the discussion on typecasts in the next
42: section.
43:
44: After declaring the use of the file \kbd{pari.h}, the first executable
45: statement of a main program should be to initialize the PARI system, and in
46: particular the PARI stack which will be both a scratchboard and a repository
47: for computed objects. This is done with a call to the function
48:
1.2 ! noro 49: \fun{void}{pari_init}{long size, ulong maxprime)}
1.1 noro 50:
51: \noindent The first argument is the number of bytes given to PARI to work
52: with (it should not reasonably be taken below 500000), and the second is the
53: upper limit on a precomputed prime number table. If you don't want prime
54: numbers, just put $\tet{maxprime} = 2$. Be careful because lots a PARI
55: functions need this table (certainly all the ones of interest to number
56: theorists). If you wind up with the error message ``not enough precomputed
57: primes'', try to increase this value.
58:
59: \noindent We have now at our disposal:
60:
61: $\bullet$ a large PARI \tev{stack} containing nothing. It's a big
62: connected chunk of memory whose size you chose when invoking
63: \tet{pari_init}. All your computations are going to take place here.
64: When doing large computations, unwanted intermediate results clutter up
65: memory very fast so some kind of garbage collecting is needed. Most large
66: systems do garbage collecting when the memory is getting scarce, and this slows
67: down the performance. In PARI we have taken a different approach: you must do
68: your own cleaning up when the intermediate results are not needed anymore.
69: Special purpose routines have been written to do this; we will see later how
70: (and when, if at all) you should use them.
71:
72: $\bullet$ the following {\it universal objects\/} (by definition, objects
73: which do not belong on the stack): the integers 0, 1 and 2 (respectively
74: called \teb{gzero}, \teb{gun}, and \teb{gdeux}), the
75: fraction $\dfrac{1}{2}$ (\teb{ghalf}), the complex number $i$
76: (\teb{gi}). All of these are of type \kbd{GEN}.
77:
78: In addition, space is reserved for the polynomials $x_v$
79: \sidx{variable}
80: (\teb{polx}\kbd{[$v$]}), and the polynomials $1_v$ (\teb{polun}\kbd{[$v$]}).
81: Here, $x_v$ is the name of variable number $v$, where $0\le v\le
82: \tet{MAXVARN}$ (the exact value of which depends on your machine, at least
83: 16383 in any case). Both \kbd{polun} and \kbd{polx} are arrays of
84: \kbd{GEN}s, the index being the polynomial variable number.
85:
86: However, except for the ones corresponding to variables $0$ and \kbd{MAXVARN},
87: these polynomials are {\it not\/} created upon initialization. It
88: is the programmer's responsibility to fill them before use. We'll see how
89: this is done in \secref{se:vars} ({\it never\/} through direct assignment).
90:
91: $\bullet$ a {\it \idx{heap}\/} which is just a linked list of permanent
92: universal objects. For now, it contains exactly the ones listed above. You
93: will probably very rarely use the heap yourself; and if so, only as a
94: collection of individual copies of objects taken from the stack
95: (called \idx{clone}s in the sequel). Thus you need not bother with its
96: internal structure, which may change as PARI evolves. Some complex PARI
97: functions may create clones for special garbage collecting purposes, usually
98: destroying them when returning.
99:
100: $\bullet$ a table of primes (in fact of {\it differences\/} between
101: consecutive primes), called \teb{diffptr}, of type \kbd{byteptr}
102: (pointer to \kbd{unsigned char}). Its use is described in appendix~C.
103:
104: $\bullet$ access to all the built-in functions of the PARI library.
105: These are declared to the outside world when you include \kbd{pari.h}, but
106: need the above things to function properly. So if you forget the call to
107: \tet{pari_init}, you will immediately get a fatal error when running your
1.2 ! noro 108: program.
1.1 noro 109:
110: \section{Important technical notes}
111:
112: \subsec{Typecasts}.\label{se:typecast}\sidx{typecast}
113:
114: \noindent
115: We have seen that, due to the non-recursiveness of the PARI types from the
116: compiler's point of view, many typecasts will be necessary when programming
117: in PARI. To take an example, a vector $V$ of dimension 2 (two components)
118: will be represented by a chunk of memory pointed to by the \kbd{GEN}~\kbd{V}.
119: \kbd{V[0]} contains coded information, in particular about the type of the
120: object, its length, etc. \kbd{V[1]} and \kbd{V[2]} contain pointers to
121: the two components of \kbd{V}. Those coefficients \kbd{V[i]} themselves are in
122: chunks of memory whose complexity depends on their own types, and so on. This
123: is where typecasting will be necessary: a priori, \kbd{V[i]} (for
124: $\kbd{i}=1,2$) is a \kbd{long}, but we will want to use it as a \kbd{GEN}.
125: The following two constructions will be exceedingly frequent (\kbd{x} and
126: \kbd{V} are \kbd{GEN}s):
127: %
128: \bprog
129: V[i] = (long) x;
130: x = (GEN) V[i];
131: @eprog
132: \noindent Note that a typecast is not a valid lvalue (cannot be put on the
133: left side of an assignment), so \kbd{(GEN)V[i] = x} would be incorrect, though
134: some compilers may accept it.
135:
136: Due to this annoyance, the PARI functions and variables that occur most
137: frequently have analogues which are macros including the typecast. The complete
138: list can be found in the file \kbd{paricast.h} (which is included by
139: \kbd{pari.h} and can be found at the same place). For instance you can
140: abbreviate:
141: %
142: \bprog
143: (long) gzero -----> zero
144: (long) gun -----> un
145: (long) polx[v] -----> lpolx[v]
146: (long) gadd(x,y) -----> ladd(x,y)
147: @eprog\noindent%
148: \kbdsidx{zero}%
149: \kbdsidx{un}%
150: \kbdsidx{lpolx}
151: In general, replacing a leading \kbd{g} by an \kbd{l} in the name of a PARI
1.2 ! noro 152: function will typecast the result to \kbd{long}. Note that \kbd{ldiv} is an ANSI C
! 153: function which is hidden in PARI by a macro of the same name representing
! 154: \kbd{(long)gdiv}.
1.1 noro 155:
156: The macro $\teb{coeff}(x,m,n)$ exists with exactly the meaning of \kbd{x[m,n]}
157: under GP when \kbd{x} is a matrix. This is a purely syntactical trick
158: to reduce the number of typecasts and thus does not create a copy of the
159: coefficient (contrary to all the library {\it functions\/}). It can be put on
160: the left side of an assignment statement, and its value, of type \kbd{long}
161: integer, is a pointer to the desired coefficient object. The macro
162: \kbd{gcoeff} is a synonym for \kbd{(GEN) coeff}, hence cannot be put on
163: the left side of an assignment.
164:
165: To retrieve the values of elements of lists of \dots\ of
166: lists of vectors, without getting infuriated by gigantic lists of typecasts,
167: we have the \teb{mael} macros (for {\bf m}ultidimensional {\bf a}rray {\bf
168: el}ement). The syntax is $\key{mael}n(x,a_1,\dots,a_n)$, where $x$ is a
169: \kbd{GEN}, the $a_i$ are indexes, and $n$ is an integer between $2$ and $5$
170: (with a standalone \kbd{mael} as a synonym for \kbd{mael2}). This stands
171: for $x[a_1][a_2]\dots[a_n]$ (with all the necessary typecasts), and returns a
172: \kbd{long} (i.e.~they are valid lvalues). The $\kbd{gmael}n$ macros are
173: synonyms for $\kbd{(GEN) mael}n$. Note that due to the implementation of matrix
174: types in PARI (i.e.~as horizontal lists of vertical vectors), \kbd{coeff(x,y)}
175: is actually completely equivalent to \kbd{mael(y,x)}. It is suggested that
176: you use \kbd{coeff} in matrix context, and \kbd{mael} otherwise.
177:
178: \subsec{Variations on basic functions}.\label{se:low_level} In the library
179: syntax descriptions in Chapter~3, we have only given the basic names of the
180: functions. For example \kbd{gadd}$(x,y)$ assumes that $x$ and $y$ are PARI
181: objects (of type \kbd{GEN}), and {\it creates\/} the result $x+y$ on the PARI
182: stack. For most of the basic operators and functions, many other variants
183: are available. We give some examples for \kbd{gadd}, but the same is true for
184: all the basic operators, as well as for some simple common functions (a
185: more complete list is given in Chapter~5):
186:
187: \fun{GEN}{gaddgs}{GEN x, long y}
188:
189: \fun{GEN}{gaddsg}{long x, GEN y}
190:
191: \noindent In the following three, \kbd{z} is a preexisting \kbd{GEN} and the
192: result of the corresponding operation is put into~\kbd{z}. The size of the PARI
193: stack does not change:
194:
195: \fun{void}{gaddz}{GEN x, GEN y, GEN z}
196:
197: \fun{void}{gaddgsz}{GEN x, long y, GEN z}
198:
199: \fun{void}{gaddsgz}{GEN x, GEN y, GEN z}
200:
201: \noindent There are also low level functions which are special cases of the
202: above:
203:
204: \fun{GEN}{addii}{GEN x, GEN y}: here $x$ and $y$ are \kbd{GEN}s of type
205: \typ{INT} (this is not checked).
206:
207: \fun{GEN}{addrr}{GEN x, GEN y}: here $x$ and $y$ are \kbd{GEN} reals
208: (type \typ{REAL}).
209:
210: \noindent
211: There also exist functions \teb{addir}, \teb{addri}, \teb{mpadd} (whose
212: two arguments can be of type integer or real), \teb{addis} (to add a \typ{INT}
213: and a \kbd{long}) and so on.
214:
215: All these functions can of course be called by the user but we feel that
216: the few microseconds lost in calling more general functions (in this case
217: \kbd{gadd}) are compensated by the fact that one needs to remember a much
218: smaller number of functions, and also because there is a hidden danger here:
219: the types of the objects that you use, if they are themselves results of a
220: previous computation, are not completely predetermined. For instance the
221: multiplication of a type real \typ{REAL} by a type integer \typ{INT}
222: {\it usually\/} gives a result of type real, except when the integer is~0, in
223: which case according to the PARI philosophy the result is the exact integer~0.
224: Hence if afterwards you call a function which specifically needs a real
225: type argument, you are going to be in trouble.
226:
227: If you really want to use these functions, their names are self-explanatory
228: once you know that {\bf i} stands for a PARI integer, {\bf r} for a PARI
229: real, {\bf mp} for i or r, {\bf s} for an ordinary signed long, whereas {\bf
230: z} (as a suffix) means that the result is not created on the PARI
231: stack but assigned to a preexisting GEN object passed as an extra argument.
232:
233: For completeness, Chapter 5 gives a description of all these
234: low-level functions.
235:
236: Please note that in the present version \vers{} the names of the functions
237: are not always consistent. This will be changed. Hence anyone programming in
238: PARI must be aware that the names of almost all functions that he uses might
239: be subject to change. If the need arises (i.e.~if there really are people out
240: there who delve into the innards of PARI), updated versions with no name
241: changes will be released.
242:
243: \subsec{Portability: 32-bit / 64-bit architectures}.
244:
245: \noindent
246: PARI supports both 32-bit and 64-bit based machines, but not simultaneously!
247: The library will have been compiled assuming a given architecture (a~priori
248: following a guess by the \kbd{Configure} program, see Appendix~A), and some
249: of the header files you include (through \kbd{pari.h}) will have been modified
250: to match the library.
251:
252: Portable macros are defined to bypass most machine dependencies. If you
253: want your programs to run identically on 32-bit and 64-bit machines, you will
254: have to use these, and not the corresponding numeric values, whenever the
255: precise size of your \kbd{long} integers might matter. Here are the most
256: important ones:
257:
258: \settabs\+
259: -----------------------------&---------------&------------&\cr
260: \+
261: & 64-bit & 32-bit
262: \cr\+
263: \tet{BITS_IN_LONG} & 64 & 32
264: \cr\+
265: \tet{LONG_IS_64BIT} & defined & undefined
266: \cr\+
267: \tet{DEFAULTPREC} & 3 & 4 & ($\approx$ 19 decimal digits, %
268: see formula below)
269: \cr\+
270: \tet{MEDDEFAULTPREC}& 4 & 6 & ($\approx$ 38 decimal digits)
271: \cr\+
272: \tet{BIGDEFAULTPREC}& 5 & 8 & ($\approx$ 57 decimal digits)
273: \cr
274: \noindent
275: For instance, suppose you call a transcendental function, such as
276:
277: \kbd{GEN \key{gexp}(GEN x, long prec)}.
278:
279: \noindent The last argument \kbd{prec} is only used if \kbd{x} is an exact
280: object (otherwise the relative precision is determined by the precision
281: of~\kbd{x}). But since \kbd{prec} sets the size of the inexact result counted
282: in (\kbd{long}) {\it words\/} (including codewords), the same value of
283: \kbd{prec} will yield different results on 32-bit and 64-bit machines. Real
284: numbers have two codewords (see~\secref{se:impl}), so the formula for
285: computing the bit accuracy is
286: $$ \tet{bit_accuracy}(\kbd{prec}) = (\kbd{prec} - 2) * \tet{BITS_IN_LONG}$$
287: (this is actually the definition of a macro). The corresponding accuracy
288: expressed in decimal digits would be
289: %
290: $$ \kbd{bit\_accuracy(prec)} * \log(2) / \log(10).$$
291: %
292: For example if the value of \kbd{prec} is 5, the corresponding accuracy for
293: 32-bit machines is $(5-2)*\log(2^{32})/\log(10)\approx 28$ decimal digits,
294: while for 64-bit machines it is $(5-2)*\log(2^{64})/\log(10)\approx 57$
295: decimal digits.
296:
297: Thus, you must take care to change the \kbd{prec} parameter you are supplying
298: according to the bit size, either using the default precisions given by the
299: various \kbd{DEFAULTPREC}s, or by using conditional constructs of the form:
300: %
301: \bprog
302: #ifndef LONG_IS_64BIT
303: prec = 4;
304: #else
305: prec = 6;
306: #endif
307: @eprog
308: \noindent which is in this case equivalent to the statement
309: \kbd{prec = MEDDEFAULTPREC;}.
310:
311: Note that for parity reasons, half the accuracies available on 32-bit
312: architectures (the odd ones) have no precise equivalents on 64-bit machines.
313:
314: \section{Creation of PARI objects, assignments, conversions}
315:
316: \subsec{Creation of PARI objects}.\sidx{creation}
317: The basic function which creates a PARI object is the function
318: \teb{cgetg} whose prototype is:
319:
320: \kbd{GEN \key{cgetg}(long length, long type)}.
321:
322: \noindent
323: Here \kbd{length} specifies the number of longwords to be allocated to the
324: object, and type is the type number of the object, preferably in symbolic
325: form (see \secref{se:impl} for the list of these). The precise effect of
326: this function is as follows: it first creates on the PARI {\it stack\/} a
327: chunk of memory of size \kbd{length} longwords, and saves the address of the
328: chunk which it will in the end return. If the stack has been used up, a
329: message to the effect that ``the PARI stack overflows'' will be printed,
330: and an error raised. Otherwise, it sets the type and length of the PARI object.
331: In effect, it fills correctly and completely its first codeword (\kbd{z[0]} or
332: \kbd{*z}). Many PARI objects also have a second codeword (types \typ{INT},
333: \typ{REAL}, \typ{PADIC}, \typ{POL}, and \typ{SER}). In case you want to
334: produce one of those from scratch (this should be exceedingly rare), {\it it
335: is your responsibility to fill this second codeword}, either explicitly (using
336: the macros described in \secref{se:impl}), or implicitly using an assignment
337: statement (using \kbd{gaffect}).
338:
339: Note that the argument \kbd{length} is predetermined for a number of types:
340: 3 for types \typ{INTMOD}, \typ{FRAC}, \typ{FRACN}, \typ{COMPLEX},
341: \typ{POLMOD}, \typ{RFRAC} and \typ{RFRACN}, 4 for type \typ{QUAD}
342: and \typ{QFI}, and 5 for type \typ{PADIC} and \typ{QFR}. However for the sake
343: of efficiency, no checking is done in the function \kbd{cgetg}, so
344: disasters will occur if you give an incorrect length.
345:
346: \misctitle{Notes}:
1.2 ! noro 347: 1) The main use of this function is to prepare for later assignments
1.1 noro 348: (see \secref{se:assign}). Most of the time you will use \kbd{GEN} objects
349: as they are created and returned by PARI functions. In this case you don't need
350: to use \kbd{cgetg} to create space to hold them.
351:
352: \noindent 2) For the creation of leaves, i.e.~integers or reals, which is
353: very common,
354:
355: \fun{GEN}{cgeti}{long length}
356:
357: \fun{GEN}{cgetr}{long length}
358:
359: \noindent should be used instead of \kbd{cgetg(length, t\_INT)} and
360: \kbd{cgetg(length, t\_REAL)} respectively.
361:
362: \noindent 3) The macros \teb{lgetg}, \teb{lgeti}, \teb{lgetr} are
363: predefined as \kbd{(long)cgetg}, \kbd{(long)cgeti}, \kbd{(long)cgetr},
364: respectively.
365:
366: \misctitle{Examples}: 1) \kbd{z = cgeti(DEFAULTPREC)} and
367: \kbd{cgetg(DEFAULTPREC, t\_INT)} create an integer object whose ``precision''
368: is \kbd{bit\_accuracy(DEFAULTPREC)} = 64. This means \kbd{z} can hold rational
369: integers of absolute value less than $2^{64}$. Note that in both cases, the
370: second codeword will {\it not\/} be filled. Of course we could use numerical
371: values, e.g.~\kbd{cgeti(4)}, but this would have different meanings on
372: different machines as \kbd{bit\_accuracy(4)} equals 64 on 32-bit machines,
373: but 128 on 64-bit machines.
374:
375: \noindent 2) The following creates a type ``complex'' object whose real and
376: imaginary parts can hold real numbers of precision
377: $\kbd{bit\_accuracy(MEDDEFAULTPREC)} = 96\hbox{ bits:}$
378: %
379: \bprog
380: z = cgetg(3, t_COMPLEX);
381: z[1] = lgetr(MEDDEFAULTPREC);
382: z[2] = lgetr(MEDDEFAULTPREC);
383: @eprog
384:
385: \noindent 3) To create a matrix object for $4\times 3$ matrices:
386: %
387: \bprog
388: z = cgetg(4, t_MAT);
389: for(i=1; i<4; i++) z[i] = lgetg(5, t_COL);
390: @eprog
391: %
392: \noindent If one wishes to create space for the matrix elements themselves,
393: one has to follow this with a double loop to fill each column vector.
394:
395: These last two examples illustrate the fact that since PARI types are
396: recursive, all the branches of the tree must be created. The function
397: \teb{cgetg} creates only the ``root'', and other calls to \kbd{cgetg} must be
398: made to produce the whole tree. For matrices, a common mistake is to think
399: that \kbd{z = cgetg(4, t\_MAT)} (for example) will create the root of the
400: matrix: one needs also to create the column vectors of the matrix (obviously,
401: since we specified only one dimension in the first \kbd{cgetg}!). This is
402: because a matrix is really just a row vector of column vectors (hence a
403: priori not a basic type), but it has been given a special type number so that
404: operations with matrices become possible.
405:
406: \subsec{Assignments}.
407: Firstly, if \kbd{x} and \kbd{y} are both declared as \kbd{GEN} (i.e.~pointers
408: to something), the ordinary C assignment \kbd{y = x} makes perfect sense: we
409: are just moving a pointer around. However, physically modifying either
410: \kbd{x} or \kbd{y} (for instance, \kbd{x[1] = 0}) will also change the other
411: one, which is usually not desirable. \label{se:assign}
412:
413: \misctitle{Very important note}: Using the functions described in this
414: paragraph is very inefficient and often awkward: one of the \tet{gerepile}
415: functions (see~\secref{se:garbage}) should be preferred. See the paragraph
416: end for some exceptions to this rule.
417:
418: \noindent
419: The general PARI \idx{assignment} function is the function \teb{gaffect} with
420: the following syntax:
421:
422: \fun{void}{gaffect}{GEN x, GEN y}
423:
424: \noindent
425: Its effect is to assign the PARI object \kbd{x} into the {\it preexisting\/}
426: object \kbd{y}. This copies the whole structure of \kbd{x} into \kbd{y} so
427: many conditions must be met for the assignment to be possible. For instance
428: it is allowed to assign an integer into a real number, but the converse is
429: forbidden. For that, you must use the truncation or rounding function of
430: your choice (see section 3.2). It can also happen that \kbd{y} is not large
431: enough or does not have the proper tree structure to receive the object
432: \kbd{x}. As an extreme example, assume \kbd{y} is the zero integer with length
433: equal to 2. Then all assignments of a non-zero integer into \kbd{y} will
434: result in an error message since \kbd{y} is not large enough to accommodate
435: a non-zero integer. In general common sense will tell you what is possible,
436: keeping in mind the PARI philosophy which says that if it makes sense it is
437: legal. For instance, the assignment of an imprecise object into a precise one
438: does {\it not\/} make sense. However, a change in precision of imprecise
439: objects is allowed.
440:
441: All functions ending in ``\kbd{z}'' such as \teb{gaddz}
442: (see~\secref{se:low_level}) implicitly use this function. In fact what they
443: exactly do is record {\teb{avma}} (see~\secref{se:garbage}), perform the
444: required operation, \teb{gaffect} the result to the last operand, then
445: restore the initial \kbd{avma}.
446:
447: You can assign ordinary C long integers into a PARI object (not necessarily
448: of type \typ{INT}). Use the function \teb{gaffsg} with the following
449: syntax:
450:
451: \fun{void}{gaffsg}{long s, GEN y}
452:
1.2 ! noro 453: \misctitle{Note}: due to the requirements mentioned above, it's usually
1.1 noro 454: a bad idea to use \tet{gaffect} statements. Two exceptions:
455:
456: $\bullet$ for simple objects (e.g.~leaves) whose size is controlled, they can
457: be easier to use than gerepile, and about as efficient.
458:
459: $\bullet$ to coerce an inexact object to a given precision. For instance
460: \bprog
461: gaffect(x, (tmp=cgetr(3))); x = tmp;
462: @eprog
463: \noindent at the beginning of a routine where precision can be kept to a
464: minimum (otherwise the precision of \kbd{x} will be used in all subsequent
465: computations, which will be a disaster if \kbd{x} is known to thousands of
466: digits).
467:
468: \subsec{Copy}. It is also very useful to \idx{copy} a PARI object, not
469: just by moving around a pointer as in the \kbd{y = x} example, but by
470: creating a copy of the whole tree structure, without pre-allocating a
471: possibly complicated \kbd{y} to use with \kbd{gaffect}. The function which
472: does this is called \teb{gcopy}, with the predefined macro
473: \teb{lcopy} as a synonym for \kbd{(long)gcopy}. Its syntax is:
474:
475: \fun{GEN}{gcopy}{GEN x}
476:
477: \noindent and the effect is to create a new copy of x on the PARI stack.
478: Beware that universal objects which occur in specific components of certain
479: types (mainly moduli for types \typ{INTMOD} and \typ{PADIC}) are not
480: copied, as they are assumed to be permanent. In this case, \kbd{gcopy} only
481: copies the pointer. Use \fun{GEN}{forcecopy}{GEN x} if you want a complete
482: copy.
483:
484: Please be sure at this point that you really understand the difference between
485: \kbd{y = x}, \kbd{y = gcopy(x)}, and \kbd{gaffect(x,y)}: this will save you
486: from many ``obvious'' mistakes later on.
487:
488: \subsec{Clones}.\sidx{clone}
489: Sometimes, it may be more efficient to create a {\it permanent\/} copy of a
490: PARI object. This will not be created on the stack but on the heap. The
491: function which does this is called \teb{gclone}, with the predefined macro
492: \teb{lclone} as a synonym for \kbd{(long)gclone}. Its syntax is:
493:
494: \fun{GEN}{gclone}{GEN x}
495:
496: A clone can be removed from the heap (thus destroyed) using
497:
498: \fun{void}{gunclone}{GEN x}
499:
500: \noindent No PARI object should keep references to a clone which has been
501: destroyed. If you want to copy a clone back to the stack then delete it, use
502: \tet{forcecopy} and not \tet{gcopy}, otherwise some components might not be
503: copied (moduli of \typ{INTMOD}s and \typ{POLMOD}s for instance).
504:
505: \subsec{Conversions}.\sidx{conversions}
506: The following functions convert C objects to PARI objects (creating them on
507: the stack as usual):
508:
509: \fun{GEN}{stoi}{long s}: C \kbd{long} integer (``small'') to PARI integer
510: (\typ{INT})
511:
512: \fun{GEN}{dbltor}{double s}: C \kbd{double} to PARI real (\typ{REAL}). The
513: accuracy of the result is 19 decimal digits, i.e.~a type \typ{REAL} of
514: length \kbd{DEFAULTPREC}, although on 32-bit machines only 16 of them will
515: be significant.
516:
517: \noindent We also have the converse functions:
518:
519: \fun{long}{itos}{GEN x}: \kbd{x} must be of type \typ{INT},
520:
521: \fun{double}{rtodbl}{GEN x}: \kbd{x} must be of type \typ{REAL},
522:
523: \noindent as well as the more general ones:
524:
525: \fun{long}{gtolong}{GEN x},
526:
527: \fun{double}{gtodouble}{GEN x}.
528:
529: \section{Garbage collection}\label{se:garbage}\sidx{garbage collecting}
530:
531: \subsec{Why and how}.
532:
533: \noindent
534: As we have seen, the \kbd{pari\_init} routine allocates a big range of
535: addresses (the {\it stack\/}) that are going to be used throughout. Recall that
536: all PARI objects are pointers. But for universal objects, they will all point
537: at some part of the stack.
538:
539: The stack starts at the address \kbd{bot} and ends just before \kbd{top}. This
540: means that the quantity
541: %
542: $$ (\kbd{top} - \kbd{bot})\,/\,\kbd{sizeof(long)} $$
543: %
544: is equal to the \kbd{size} argument of \kbd{pari\_init}. The PARI
545: stack also has a ``current stack pointer'' called \teb{avma}, which
546: stands for {\bf av}ailable {\bf m}emory {\bf a}ddress. These three variables
547: are global (declared for you by \kbd{pari.h}). For historical reasons they
548: are of type \tet{ulong} (unsigned long) and not of type \kbd{GEN} as would
549: seem more natural.
550:
551: The stack is oriented upside-down: the more recent an object, the closer to
552: \kbd{bot}. Accordingly, initially \kbd{avma} = \kbd{top}, and \kbd{avma} gets
553: {\it decremented\/} as new objects are created. As its name indicates,
554: \kbd{avma} always points just {\it after\/} the first free address on the
555: stack, and \kbd{(GEN)avma} is always (a pointer to) the latest created object.
556: When \kbd{avma} reaches \kbd{bot}, the stack overflows, aborting all
557: computations, and an error message is issued. To avoid this {\it you\/} will
558: need to clean up the stack from time to time, when some bunch of intermediate
559: objects will not be needed anymore. This is called ``{\it garbage
560: collecting}.''
561:
562: We are now going to describe briefly how this is done. We will see many
563: concrete examples in the next subsection.
564:
565: \noindent$\bullet$
566: First, PARI routines will do their own garbage collecting, which
567: means that whenever a documented function from the library returns, only its
568: result(s) will have been added to the stack (non-documented ones may not do
569: this, for greater speed). In particular, a PARI function that does not return
570: a \kbd{GEN} does not clutter the stack. Thus, if your computation is small
571: enough (i.e.~you call few PARI routines, or most of them return \kbd{long}
572: integers), then you don't need to do any garbage collecting. This will probably
573: be the case in many of your subroutines. Of course the objects that were on
574: the stack {\it before\/} the function call are left alone. Except for the ones
575: listed below, PARI functions only collect their own garbage.
576:
577: \noindent$\bullet$
578: It may happen that all objects that were created after a certain point can
579: be deleted~--- for instance, if the final result you need is not a
580: \kbd{GEN}, or if some search proved futile. Then, it is enough to record
581: the value of \kbd{avma} just {\it before\/} the first garbage is created,
582: and restore it upon exit:
583:
584: \bprog
585: ulong ltop = avma; /*@Ccom record initial avma */
586:
587: garbage ...
588: avma = ltop; /*@Ccom restore it */
589: @eprog
590: \noindent All objects created in the \kbd{garbage} zone will eventually
591: be overwritten: they should not be accessed anymore once \kbd{avma} has been
592: restored.
593:
594: \noindent$\bullet$
595: If you want to destroy (i.e.~give back the memory occupied by) the
596: latest PARI object on the stack (e.g.~the latest one obtained from a function
597: call), and the above method is not available (because the initial value of
598: \kbd{avma} has been lost), just use the function\sidx{destruction}%
599: \vadjust{\penalty500}%discourage page break
600:
601: \fun{void}{cgiv}{GEN z}
602:
603: \noindent where \kbd{z} is the object you want to give back.
604:
605: \noindent$\bullet$
606: Unfortunately life is not so simple, and sometimes you will want
607: to give back accumulated garbage {\it during\/} a computation without losing
608: recent data. For this you need the \kbd{gerepile} function (or one of its
609: variants described hereafter):
610:
611: \fun{GEN}{gerepile}{ulong ltop, ulong lbot, GEN q}
612:
613: \noindent
614: This function cleans up the stack between \kbd{ltop} and \kbd{lbot}, where
615: $\kbd{lbot} < \kbd{ltop}$, and returns the updated object \kbd{q}. This means:
616:
617: 1) we translate (copy) all the objects in the interval
618: $[\kbd{avma}, \kbd{lbot}[$, so that its right extremity abuts the address
619: \kbd{ltop}. Graphically
620:
621: \vbox{\bprog
622: bot avma lbot ltop top
623: End of stack |-------------[++++++[-/-/-/-/-/-/-|++++++++| Start
624: free memory garbage
625: @eprog
626: \noindent becomes:
627: \bprog
628: bot avma ltop top
629: End of stack |---------------------------[++++++[++++++++| Start
630: free memory
631: @eprog
632: }
633: \noindent where \kbd{++} denote significant objects, \kbd{--} the unused part
634: of the stack, and \kbd{-/-} the garbage we remove.
635:
636: 2) The function then inspects all the PARI objects between \kbd{avma} and
637: \kbd{lbot} (i.e.~the ones that we want to keep and that have been translated)
638: and looks at every component of such an object which is not a codeword. Each
639: such component is a pointer to an object whose address is either
640:
641: --- between \kbd{avma} and \kbd{lbot}, in which case it will be suitably
642: updated,
643:
644: --- larger than or equal to \kbd{ltop}, in which case it will not change, or
645:
646: --- between \kbd{lbot} and \kbd{ltop} in which case \kbd{gerepile} will
647: scream an error message at you (``significant pointers lost in gerepile'').
648:
649: 3) \key{avma} is updated (we add $\kbd{ltop} - \kbd{lbot}$ to the old value).
650:
651: 4) We return the (possibly updated) object \kbd{q}: if \kbd{q} initially
652: pointed between \kbd{avma} and \kbd{lbot}, we return the translated
653: address, as in~2). If not, the original address is still valid (and we return
654: it!).
655:
656: As stated above, no component of the remaining objects (in particular
657: \kbd{q}) should belong to the erased segment [\kbd{lbot}, \kbd{ltop}[, and
658: this is checked within \kbd{gerepile}. But beware as well that the addresses
659: of all the objects in the translated zone will have changed after a call to
660: \kbd{gerepile}: every pointer you may have kept around elsewhere, outside the
661: stack objects, which previously pointed into the zone below \kbd{ltop}
662: must be discarded. If you need to recover more than one object, use one of
663: the \kbd{gerepilemany} functions below.
664:
665: As a consequence of the preceding explanation, we must now state the most
666: important law about programming in PARI:
667:
668: {\bf If a given PARI object is to be relocated by \hbox{gerepile} then,
669: apart from universal objects, the chunks of memory used by its components
670: should be in consecutive memory locations}. All \kbd{GEN}s created by
671: documented PARI function are guaranteed to satisfy this.
672:
673: This is because the \kbd{gerepile} function knows only about {\it two
674: connected zones\/}: the garbage that will be erased (between \kbd{lbot} and
675: \kbd{ltop}) and the significant pointers that will be copied and updated.
676: If there is garbage interspersed with your objects, disasters will occur when
677: we try to update them and consider the corresponding ``pointers''. So be
678: {\it very\/} wary when you allow objects to become disconnected. Have a look
679: at the examples, it's not as complicated as it seems.
680:
681: \noindent In practice this is achieved by the following programming idiom:
682: \bprog
683: ltop=avma; garbage(); lbot=avma; q=anything();
684: return gerepile(ltop, lbot, q); /*@Ccom returns the updated q */
685: @eprog
686:
687: \noindent Beware that
688: \bprog
689: ltop=avma; garbage();
690: return gerepile(ltop, avma, anything())
691: @eprog
692:
693: \noindent might work, but should be frowned upon. We can't predict whether
694: \kbd{avma} is going to be evaluated after or before the call to
695: \kbd{anything()}: it depends on the compiler. If we are out of luck, it will
696: be {\it after\/} the call, so the result will belong to the garbage zone and
697: the \kbd{gerepile} statement becomes equivalent to \kbd{avma = ltop}. Thus we
698: would return a pointer to random garbage.
699:
700: \noindent$\bullet$ A simple variant is
701:
702: \fun{GEN}{gerepileupto}{ulong ltop, GEN q}
703:
704: \noindent which cleans the stack between \kbd{ltop} and the {\it connected\/}
705: object \kbd{q} and returns \kbd{q} updated. For this to work, \kbd{q} must
706: have been created {\it before\/} all its components, otherwise they would
707: belong to the garbage zone! Documented PARI functions guarantee this. If you
708: stumble upon one that does not, consider it a bug worth reporting.
709:
710: \noindent$\bullet$ Another variant (a special case of \tet{gerepilemany}
711: below, where $n=1$) is
712:
713: \fun{GEN}{gerepilecopy}{ulong ltop, GEN x)}
714:
715: \noindent which is functionnally equivalent to \kbd{gerepileupto(ltop,
716: gcopy(x))} but more efficient. In this case, the \kbd{GEN} parameter \kbd{x}
717: need not satisfy any property before the garbage collection (it may be
718: disconnected, components created before the root and so on). Of course, this
719: is less efficient than either \tet{gerepileupto} or \tet{gerepile}, because
720: \kbd{x} has to be copied to a clean stack zone first.
721:
722: \noindent$\bullet$
723: To cope with complicated cases where many objects have to be
724: preserved, you can use
725:
1.2 ! noro 726: \fun{void}{gerepileall}{ulong ltop, int n, ...}
1.1 noro 727:
1.2 ! noro 728: \noindent where the routine expects $n$ further arguments, which are the {\it
! 729: addresses} of the \kbd{GEN}s you want to preserve. It cleans up the most
! 730: recent part of the stack (between \kbd{ltop} and \kbd{avma}), updating all
! 731: the \kbd{GEN}s added to the argument list. A copy is done just before the
! 732: cleaning to preserve them, so they don't need to be connected before the
! 733: call. With \kbd{gerepilecopy}, this is the most robust of the \kbd{gerepile}
! 734: functions (the less prone to user error), but also the slowest.
! 735:
! 736: An alternative, unnecessarily complicated, syntax is given by
! 737:
! 738: \fun{void}{gerepilemany}{ulong ltop, GEN *gptr[], int n}
! 739:
! 740: \noindent which works exactly as above, except that the preserved \kbd{GEN}s
! 741: are the elements of the array \kbd{gptr} (of length $n$): \kbd{gptr[0]},
! 742: \kbd{gptr[1]}, \dots, \kbd{gptr[$n$-1]}.
1.1 noro 743:
744: \noindent$\bullet$ More efficient, but tricky to use is
745:
1.2 ! noro 746: \fun{void}{gerepilemanysp}{ulong ltop, ulong lbot, GEN *gptr[], int n}
1.1 noro 747:
748: \noindent which cleans the stack between \kbd{lbot} and \kbd{ltop} and
749: updates the \kbd{GEN}s pointed at by the elements of \kbd{gptr} without doing
750: any copying. This is subject to the same restrictions as \kbd{gerepile}, the
751: only difference being that more than one address gets updated.
752:
753: \subsec{Examples}.
754:
755: \noindent
756: Let \kbd{x} and \kbd{y} be two preexisting PARI objects and suppose that we
757: want to compute $\kbd{x}^2+ \kbd{y}^2$. This can trivially be done using the
758: following program (we skip the necessary declarations; everything in sight is
759: a \kbd{GEN}):
760: \bprog
761: p1 = gsqr(x);
762: p2 = gsqr(y); z = gadd(p1,p2);
763: @eprog\noindent
764: The \kbd{GEN} \kbd{z} indeed points at the desired quantity. However,
765: consider the stack: it contains as unnecessary garbage \kbd{p1} and \kbd{p2}.
766: More precisely it contains (in this order) \kbd{z}, \kbd{p2}, \kbd{p1}
767: (recall that, since the stack grows downward from the top, the most recent
768: object comes first). We need a way to get rid of this garbage (in this case
769: it causes no harm except that it occupies memory space, but in other cases
770: it could disconnect other PARI objects and this is dangerous).
771:
772: It would not have been possible to get rid of \kbd{p1}, \kbd{p2} before
773: \kbd{z} is computed, since they are used in the final operation. We cannot
774: record \kbd{avma} before \kbd{p1} is computed and restore it later, since
775: this would destroy \kbd{z} as well. It is not possible either to use the
776: function \kbd{cgiv} since \kbd{p1} and \kbd{p2} are not at the bottom of the
777: stack and we don't want to give back~\kbd{z}.
778:
779: But using \kbd{gerepile}, we can give back the memory locations corresponding
780: to \kbd{p1}, \kbd{p2}, and move the object \kbd{z} upwards so that no
781: space is lost. Specifically:
782: \bprog
783: ltop = avma; /*@Ccom remember the current address of the top of the stack */
784: p1 = gsqr(x); p2 = gsqr(y);
785: lbot = avma; /*@Ccom keep the address of the bottom of the garbage pile */
786: z = gadd(p1, p2); /*@Ccom z is now the last object on the stack */
787: z = gerepile(ltop, lbot, z); /*@Ccom garbage collecting */
788: @eprog
789: \noindent Of course, the last two instructions could also have been
790: written more simply:
791:
792: \kbd{z = gerepile(ltop, lbot, gadd(p1,p2));}
793:
794: \noindent In fact \kbd{gerepileupto} is even simpler to use, because the
795: result of \kbd{gadd} will be the last object on the stack and \kbd{gadd} is
796: guaranteed to return an object suitable for \kbd{gerepileupto}:
797: \bprog
798: ltop = avma;
799: z = gerepileupto(ltop, gadd(gsqr(x), gsqr(y)));
800: @eprog
801: \noindent
802: As you can see, in simple conditions the use of \kbd{gerepile} is not really
803: difficult. However make sure you understand exactly what has happened before
804: you go on (use the figure from the preceding section).
805:
806: \misctitle{Important remark}: as we will see presently it is often
807: necessary to do several \kbd{gerepile}s during a computation. However, the
808: fewer the better. The only condition for \kbd{gerepile} to work is that the
809: garbage be connected. If the computation can be arranged so that there is a
810: minimal number of connected pieces of garbage, then it should be done that
811: way.
812:
813: For example suppose we want to write a function of two \kbd{GEN} variables
814: \kbd{x} and \kbd{y} which creates the vector $\kbd{[x}^2+\kbd{y},
815: \kbd{y}^2+\kbd{x]}$. Without garbage collecting, one would write:
816: %
817: \bprog
818: p1 = gsqr(x); p2 = gadd(p1, y);
819: p3 = gsqr(y); p4 = gadd(p3, x);
820: z = cgetg(3,t_VEC);
821: z[1] = (long)p2;
822: z[2] = (long)p4;
823: @eprog
824: \noindent
825: This leaves a dirty stack containing (in this order) \kbd{z}, \kbd{p4},
826: \kbd{p3}, \kbd{p2}, \kbd{p1}. The garbage here consists of \kbd{p1} and
827: \kbd{p3}, which are separated by \kbd{p2}. But if we compute \kbd{p3}
828: {\it before\/} \kbd{p2} then the garbage becomes connected, and we get the
829: following program with garbage collecting:
830: %
831: \bprog
832: ltop = avma; p1 = gsqr(x); p3 = gsqr(y); lbot = avma;
833: z = cgetg(3,t_VEC);
834: z[1] = ladd(p1,y);
835: z[2] = ladd(p3,x);
836: z = gerepile(ltop,lbot,z);
837: @eprog
838:
839: \noindent Finishing by \kbd{z = gerepileupto(ltop, z)} would be ok as well.
840: But when you have the choice, it's usually clearer to brace the garbage
841: between \kbd{ltop}~/ \kbd{lbot} pairs.
842:
843: \noindent Beware that
844: \bprog
845: ltop = avma; p1 = gadd(gsqr(x), y); p3 = gadd(gsqr(y), x);
846: z = cgetg(3,t_VEC);
847: z[1] = (long)p1;
848: z[2] = (long)p3
849: z = gerepileupto(ltop,z); /*@Ccom WRONG !!! */
850: @eprog
851: \noindent would be a disaster since \kbd{p1} and \kbd{p3} would be created
852: before \kbd{z}, so the call to \kbd{gerepileupto} would overwrite them,
853: leaving \kbd{z[1]} and \kbd{z[2]} pointing at random data!
854:
855: We next want to write a program to compute the product of two complex numbers
856: $x$ and $y$, using a method which takes only 3 multiplications instead of 4.
857: Let $z = x*y$, and set $x = x_r + i*x_i$ and similarly for $y$ and $z$. The
858: well-known trick is to compute $p_1 = x_r*y_r$, $p_2=x_i*y_i$,
859: $p_3=(x_r+x_i)*(y_r+y_i)$, and then we have $z_r=p_1-p_2$,
860: $z_i=p_3-(p_1+p_2)$. The program is essentially as follows:
861: %
862: \bprog
863: ltop = avma;
864: p1 = gmul(x[1],y[1]);
865: p2 = gmul(x[2],y[2]);
866: p3 = gmul(gadd(x[1],x[2]), gadd(y[1],y[2]));
867: p4 = gadd(p1,p2); lbot = avma;
868: z = cgetg(3,t_COMPLEX);
869: z[1] = lsub(p1,p2);
870: z[2] = lsub(p3,p4);
871: z = gerepile(ltop,lbot,z);
872: @eprog
873: \noindent
874: ``Essentially'', because for instance \kbd{x[1]} is a \kbd{long} and not a
875: \kbd{GEN}, so we need to insert many annoying typecasts:
876: \kbd{p1 = gmul((GEN)x[1], (GEN)y[1])} and so on.
877:
878: Let us now look at a less trivial example where more than one \kbd{gerepile}
879: is needed in practice (at the expense of efficiency, one can always use only
880: one using \kbd{gerepilecopy}). Suppose that we want to write a function
881: which multiplies a line vector by a matrix (such a function is of course
882: already part of \kbd{gmul}, but let's ignore this for a moment). Then the
883: most natural way is to do a \kbd{cgetg} of the result immediately, and then a
884: \kbd{gerepile} for each coefficient of the result vector to get rid of the
885: garbage which has accumulated while this particular coefficient was computed.
886: We leave the details to the reader, who can look at the answer in the file
887: \kbd{basemath/gen1.c}, in the function \kbd{gmul}, case \typ{VEC} times case
888: \typ{MAT}. It would theoretically be possible to have a single connected
889: piece of garbage, but it would be a much less natural and unnecessarily
890: complicated program, which would in fact be slower.
891:
892: Let us now take an example which is probably the least trivial way of using
893: \kbd{gerepile}, but is unfortunately sometimes necessary. Although it is not
894: an infrequent occurrence, we will not give a specific example but a general
895: one: suppose that we want to do a computation (usually inside a larger
896: function) producing more than one PARI object as a result, say two for
897: instance. Then even if we set up the work properly, before cleaning up we
898: will have a stack which has the desired results \kbd{z1}, \kbd{z2} (say),
899: and then connected garbage from lbot to ltop. If we write
900:
901: \kbd{z1 = gerepile(ltop,lbot,z1);}
902:
903: \noindent
904: then the stack will be cleaned, the pointers fixed up, but we will have lost
905: the address of \kbd{z2}. This is where we need one of the \idx{gerepilemany}
906: functions: we declare
907: \bprog
908: GEN *gptr[2]; /*@Ccom Array of pointers to GENs */
909: gptr[0] = &z1; gptr[1] = &z2;
910: @eprog
911: \noindent and now the call \kbd{gerepilemany(ltop, gptr, 2)} copies \kbd{z1}
912: and \kbd{z2} to new locations, cleans the stack from \kbd{ltop} to the old
913: \kbd{avma}, and updates the pointers \kbd{z1} and \kbd{z2}. Here we don't
914: assume anything about the stack: the garbage can be disconnected and
915: \kbd{z1}, \kbd{z2} need not be at the bottom of the stack. If all of these
916: assumptions are in fact satisfied, then we can call \kbd{gerepilemanysp}
917: instead, which will usually be faster since we don't need the initial copy
1.2 ! noro 918: (on the other hand, it is less cache friendly). The whole thing could (and
! 919: should) be abbreviated as
! 920: \bprog
! 921: gerepileall(ltop, 2, &z1, &z2)
! 922: @eprog
! 923: so that the array \kbd{gptr} is not needed anymore.
1.1 noro 924:
925: Another important usage is ``random'' garbage collection during loops
926: whose size requirements we cannot (or don't bother to) control in advance:
927: \bprog
928: ulong ltop = avma, limit = (avma+bot)/2;
929: GEN x, y;
930:
931: while (...)
932: {
933: garbage(); x = anything();
934: garbage(); y = anything()
935: garbage();
936: if (avma < limit) /*@Ccom memory is running low (half spent since entry) */
1.2 ! noro 937: gerepileall(ltop, 2, &x, &y);
1.1 noro 938: }
939: @eprog
940: \noindent Here we assume that only \kbd{x} and \kbd{y} are needed from one
941: iteration to the next. As it would be too costly to call gerepile once for
942: each iteration, we only do it when it seems to have become necessary. Of
943: course, when the need arises, you can use bigger \kbd{gptr} arrays: in the
944: PARI library source code, we once needed to preserve up to 10 objects at a
945: time (in a variant of the LLL algorithm)!
946:
947: \misctitle{Technical note:} the statement \kbd{limit = (avma+bot)/2} is
948: dangerous since the addition can overflow, which would result in \kbd{limit}
949: being always smaller than \kbd{avma}. This will prevent garbage collection in
950: the loop. To avoid this problem, we provide a robust macro
951: \tet{stack_lim}\kbd{(avma,$n$)}, which denotes an address where $2^{n-1} /
952: (2^{n-1}+1)$ of the total stack space is exhausted ($1/2$ for $n=1$, $2/3$
953: for $n=2$). Hence, the above snippet should be written as
954:
955: \bprog
956: ulong ltop = avma, limit = stack_lim(avma,1);
957: @dots
958: @eprog
959:
960: \subsec{Some hints and tricks}. In this section, we give some indications
961: on how to avoid most problems connected with garbage collecting:
962:
963: First, although it looks complicated, \kbd{gerepile} has turned out to be a
964: very flexible and fast garbage collector, which compares very favorably
965: with much more sophisticated methods used in other systems. A few tests that
966: we have done indicate that the price paid for using \kbd{gerepile}, when
967: properly used, is usually around 1 or 2 percents of the total running time,
968: which is quite acceptable.
969:
970: Secondly, in many cases, in particular when the tree structure and the size of
971: the PARI objects which will appear in a computation are under control, one
972: can avoid \kbd{gerepile} altogether by creating sufficiently large objects at
973: the beginning (using \teb{cgetg}), and then using assignment statements and
974: operations ending with z (such as \teb{gaddz}). Coming back to our first
975: example, note that if we know that x and y are of type real and of length
976: less than or equal to 5, we can program without using \kbd{gerepile} at all:
977:
978: \bprog
979: z = cgetr(5); ltop = avma;
980: p1 = gsqr(x); p2 = gsqr(y); gaddz(p1,p2,z);
981: avma = ltop;
982: @eprog
983: \noindent This practice will usually be {\it slower\/} than a craftily used
984: \kbd{gerepile} though, and is certainly more cumbersome to use. As a rule,
985: assignment statements should generally be avoided.
986: \smallskip
987:
988: Thirdly, the philosophy of \kbd{gerepile} is the following: keep the value of
989: the stack pointer \kbd{avma} at the beginning, and just {\it before\/} the
990: last operation. Afterwards, it would be too late since the lower end address
991: of the garbage zone would have been lost. Of course you can always use
992: \kbd{gerepileupto}, but you will have to assume that the object was created
993: {\it before\/} its components.
994:
995: Lastly, if all seems lost, just use \tet{gerepilecopy} (or
996: \tet{gerepilemany}) to fix up the stack for you.
997:
998: \smallskip If you followed us this far, congratulations, and rejoice: the
999: rest is much easier.
1000:
1001: \section{Implementation of the PARI types}
1002: \label{se:impl}
1003:
1004: \noindent
1005: Although it is a little tedious, we now go through each type and explain its
1006: implementation. Let \kbd{z} be a \kbd{GEN}, pointing at a PARI object. In
1007: the following paragraphs, we will constantly mix two points of view: on the
1008: one hand, \kbd{z} will be treated as the C pointer it is (in the context of
1009: program fragments like \kbd{z[1]}), on the other, as PARI's handle on (the
1010: internal representation of) some mathematical entity, so we will shamelessly
1011: write $\kbd{z} \ne 0$ to indicate that the {\it value\/} thus represented
1012: is nonzero (in which case the {\it pointer\/}~\kbd{z} certainly will be
1013: non-\kbd{NULL}). We offer no apologies for this style. In fact, you had
1014: better feel comfortable juggling both views simultaneously in your mind if
1015: you want to write correct PARI programs.
1016:
1017: Common to all the types is the
1018: first codeword \kbd{z[0]}, which we don't have to worry about since this is
1019: taken care of by \kbd{cgetg}. Its precise structure will depend on the
1020: machine you are using, but it always contain the following data: the
1021: {\it internal \idx{type number}\/} associated to the symbolic type name, the
1022: {\it\idx{length}\/} of the root in longwords, and a technical bit which
1023: indicates whether the object is a clone (see below) or not. This last one is
1024: used by GP for internal garbage collecting, you won't have to worry about it.
1025:
1026: \noindent These data can be handled through the following {\it macros\/}:
1027:
1028: \fun{long}{typ}{GEN z} returns the type number of \kbd{z}.
1029:
1030: \fun{void}{settyp}{GEN z, long n} sets the type number of \kbd{z} to
1031: \kbd{n} (you should not have to use this function if you use \kbd{cgetg}).
1032:
1033: \fun{long}{lg}{GEN z} returns the length (in longwords) of the root of \kbd{z}.
1034:
1035: \fun{long}{setlg}{GEN z, long l} sets the length of \kbd{z} to \kbd{l} (you
1036: should not have to use this function if you use \kbd{cgetg}; however, see
1037: an advanced example in \secref{se:prog}).
1038:
1039: \noindent
1040: (If you know enough about PARI to need to access the ``clone'' bit, then you'll
1041: be able to find usage hints in the code (esp. \kbd{killbloc()} and
1042: \kbd{matrix\_block()}). It {\it is\/} technical after all.)
1043:
1044: These macros are written in such a way that you don't need to worry about
1045: type casts when using them: i.e.~if \kbd{z} is a \kbd{GEN}, \kbd{typ(z[2])}
1046: will be accepted by your compiler, without your having to explicitly type
1047: \kbd{typ((GEN)z[2])}. Note that for the sake of efficiency, none of the
1048: codeword-handling macros check the types of their arguments even when there
1049: are stringent restrictions on their use.
1050:
1051: The possible second codeword is used differently by the different types, and
1052: we will describe it as we now consider each of them in turn:
1053:
1054:
1055: \subsec{Type \typ{INT} (integer):}\sidx{integer}\kbdsidx{t_INT} this type has
1056: a second codeword \kbd{z[1]} which contains the following information:
1057:
1058: the sign of \kbd{z}: coded as $1$, $0$ or $-1$ if $\kbd{z} > 0$, $\kbd{z} = 0$,
1059: $\kbd{z} < 0$ respectively.
1060:
1061: the {\it effective length\/} of \kbd{z}, i.e.~the total number of significant
1062: longwords. This means the following: apart from the integer 0, every integer
1063: is ``normalized'', meaning that the first mantissa longword (i.e.~\kbd{z[2]})
1064: is non-zero. However, the integer may have been created with a longer length.
1065: Hence the ``length'' which is in \kbd{z[0]} can be larger than the
1066: ``effective length'' which is in \kbd{z[1]}. Accessing \kbd{z[i]} for \kbd{i}
1067: larger than or equal to the effective length will yield random results.
1068:
1069: \noindent This information is handled using the following macros:
1070:
1071: \fun{long}{signe}{GEN z} returns the sign of \kbd{z}.
1072:
1073: \fun{void}{setsigne}{GEN z, long s} sets the sign of \kbd{z} to \kbd{s}.
1074:
1075: \fun{long}{lgefint}{GEN z} returns the \idx{effective length} of \kbd{z}.
1076:
1077: \fun{void}{setlgefint}{GEN z, long l} sets the effective length
1078: of \kbd{z} to \kbd{l}.
1079:
1080: The integer 0 can be recognized either by its sign being~0, or by its
1081: effective length being equal to~2. When $\kbd{z} \ne 0$, the word
1082: \kbd{z[2]} exists and is non-zero, and the absolute value of \kbd{z}
1083: is (\kbd{z[2]},\kbd{z[3]},\dots,\kbd{z[lgefint(z)-1]}) in base
1084: \kbd{2\pow BITS\_IN\_LONG}, where as usual in this notation \kbd{z[2]} is
1085: the highest order longword.
1086:
1087: \noindent The following further macros are available:
1088:
1089: \fun{long}{mpodd}{GEN x} which is 1 if \kbd{x} is odd, and 0 otherwise.
1090:
1091: \fun{long}{mod2}{GEN x}, \fun{}{mod4}{x}, and so on up to \fun{}{mod64}{x},
1092: which give the residue class of \kbd{x} modulo the corresponding power of
1.2 ! noro 1093: 2, for {\it positive\/}~\kbd{x}. By definition, $\kbd{mod}n(x) :=
! 1094: \kbd{mod}n(|x|)$ for $x < 0$ (the macros disregard the sign), and the
! 1095: result is undefined if $x = 0$.
1.1 noro 1096:
1097: These macros directly access the binary data and are thus much faster than
1098: the generic modulo functions. Besides, they return long integers instead of
1099: \kbd{GEN}s, so they don't clutter up the stack.
1100:
1101: \subsec{Type \typ{REAL} (real number):}\kbdsidx{t_REAL}\sidx{real number}
1102: this type has a second codeword z[1] which also encodes its sign (obtained
1.2 ! noro 1103: or set using the same functions as for the integers), and a binary
! 1104: exponent. This exponent can be handled using the following macros:
1.1 noro 1105:
1.2 ! noro 1106: \fun{long}{expo}{GEN z} returns the exponent of \kbd{z}.
1.1 noro 1107: This is defined even when \kbd{z} is equal to zero, see
1108: \secref{se:whatzero}.
1109:
1.2 ! noro 1110: \fun{void}{setexpo}{GEN z, long e} sets the exponent of \kbd{z} to \kbd{e}.
1.1 noro 1111:
1112: \noindent Note the functions:
1113:
1114: \fun{long}{gexpo}{GEN z} which tries to return an exponent for \kbd{z},
1115: even if \kbd{z} is not a real number.
1116:
1117: \fun{long}{gsigne}{GEN z} which returns a sign for \kbd{z}, even when
1118: \kbd{z} is neither real nor integer (a rational number for instance).
1119:
1.2 ! noro 1120: The real zero is characterized by having its sign equal to 0. If \kbd{z} is
! 1121: not equal to~0, then is is represented as $2^e M$, where $e$ is the exponent,
! 1122: and $M\in [1, 2[$ is the mantissa of $z$, whose digits are stored in
! 1123: $\kbd{z[2]},\dots, \kbd{z[lg(z)-1]}$.
! 1124:
! 1125: More precisely, let $m$ be the integer (\kbd{z[2]},\dots, \kbd{z[lg(z)-1]})
! 1126: in base \kbd{2\pow BITS\_IN\_LONG}; here, \kbd{z[2]} is the most significant
! 1127: longword and is normalized, i.e.~its most significant bit is~1. Then we have
! 1128: $M := m \cdot 2^{1 - \kbd{bit\_accuracy(lg(z))}}$.
! 1129:
! 1130: Thus, the real number $3.5$ to accuracy \kbd{bit\_accuracy(lg(z))} is
! 1131: represented as \kbd{z[0]} (encoding $\kbd{type} = \typ{REAL}$, \kbd{lg(z)}),
! 1132: \kbd{z[1]} (encoding $\kbd{sign} = 1$, $\kbd{expo} = 1$), $\kbd{z[2]} =
! 1133: \kbd{0xe0000000}$, $\kbd{z[3]} =\dots = \kbd{z[lg(z)-1]} = \kbd{0x0}$.
1.1 noro 1134:
1135: \subsec{Type \typ{INTMOD} (integermod):}\kbdsidx{t_INTMOD}\sidx{integermod}
1136: \kbd{z[1]} points to the modulus, and \kbd{z[2]} at the number representing
1137: the class \kbd{z}. Both are separate \kbd{GEN} objects, and both must be of
1138: type integer, satisfying the inequality $0 \le \kbd{z[2]} < \kbd{z[1]}$.
1139:
1140: It is good practice to keep the modulus object on the heap, so that new
1141: integermods resulting from operations can point at this common object,
1142: instead of carrying along their own copies of it on the stack. The library
1143: functions implement this practice almost by default.
1144:
1145: \subsec{Type \typ{FRAC} and \typ{FRACN} (rational number):}%
1146: \kbdsidx{t_FRAC}\kbdsidx{t_FRACN}\sidx{rational number}
1147: \kbd{z[1]} points to the numerator, and \kbd{z[2]} to the denominator. Both
1148: must be of type integer. In principle $\kbd{z[2]} > 0$, but this rule does not
1149: have to be strictly obeyed. Note that a type \typ{FRACN} rational number can be
1150: converted to irreducible form using the function \fun{GEN}{gred}{GEN x}.
1151:
1152: \subsec{Type \typ{COMPLEX} (complex number):}%
1153: \kbdsidx{t_COMPLEX}\sidx{complex number}
1154: \kbd{z[1]} points to the real part, and \kbd{z[2]} to the imaginary part. A
1155: priori \kbd{z[1]} and \kbd{z[2]} can be of any type, but only certain types
1156: are useful and make sense.
1157:
1158: \subsec{Type \typ{PADIC} ($p$-adic numbers):}%
1159: \sidx{p-adic number}\kbdsidx{t_PADIC} this type has a second codeword
1160: \kbd{[1]} which contains the following information: the $p$-adic precision
1161: (the exponent of $p$ modulo which the $p$-adic unit corresponding to
1162: \kbd{z} is defined if \kbd{z} is not~0), i.e.~one less than the number of
1.2 ! noro 1163: significant $p$-adic digits, and the exponent of \kbd{z}. This information
! 1164: can be handled using the following functions:
1.1 noro 1165:
1166: \fun{long}{precp}{GEN z} returns the $p$-adic precision of \kbd{z}.
1167:
1168: \fun{void}{setprecp}{GEN z, long l} sets the $p$-adic precision of \kbd{z}
1169: to \kbd{l}.
1170:
1171: \fun{long}{valp}{GEN z} returns the $p$-adic valuation of \kbd{z} (i.e. the
1.2 ! noro 1172: exponent). This is defined even if \kbd{z} is equal to~0, see
1.1 noro 1173: \secref{se:whatzero}.
1174:
1175: \fun{void}{setvalp}{GEN z, long e} sets the $p$-adic valuation of \kbd{z}
1176: to \kbd{e}.
1177:
1178: In addition to this codeword, \kbd{z[2]} points to the prime $p$,
1179: \kbd{z[3]} points to $p^{\text{precp(z)}}$, and \kbd{z[4]} points to an
1180: integer representing the $p$-adic unit associated to \kbd{z} modulo
1181: \kbd{z[3]} (and points to zero if \kbd{z} is zero). To summarize, if $z\neq
1182: 0$, we have the equality:
1183: $$ \kbd{z} = p^{\text{valp(z)}} * (\kbd{z[4]} + O(\kbd{z[3]})) =
1184: p^{\text{valp(z)}} * (\kbd{z[4]} + O(p^{\text{precp(z)}})) $$
1185:
1186: \subsec{Type \typ{QUAD} (quadratic number):}\sidx{quadratic
1187: number}\kbdsidx{t_QUAD} \kbd{z[1]} points to the polynomial defining the
1188: quadratic field, \kbd{z[2]} to the ``real part'' and \kbd{z[3]} to the
1189: ``imaginary part'', which are to be taken as the coefficients of \kbd{z}
1190: with respect to the ``canonical'' basis $(1,w)$, see~\secref{se:compquad}.
1191: Complex numbers are a particular case of quadratics but deserve a separate
1192: type.
1193:
1194: \subsec{Type \typ{POLMOD} (polmod):}\kbdsidx{t_POLMOD}\sidx{polmod} exactly as
1195: for integermods, \kbd{z[1]} points to the modulus, and \kbd{z[2]} to a
1196: polynomial representing the class of~\kbd{z}. Both must be of type
1197: \typ{POL} in the same variable. However, \kbd{z[2]} is allowed to be a
1198: simplification of such a polynomial, e.g a scalar. This is quite tricky
1199: considering the hierarchical structure of the variables; in particular, a
1.2 ! noro 1200: polynomial in variable of \var{lesser} priority (see \secref{se:priority})
! 1201: than the modulus variable is valid, since it can be considered as the
! 1202: constant term of a polynomial of degree 0 in the correct variable. On the
! 1203: other hand a variable of \var{greater} priority would not be acceptable.
1.1 noro 1204:
1205: \subsec{Type \typ{POL} (polynomial):}\kbdsidx{t_POL}\sidx{polynomial} this
1206: type has a second codeword which is analogous to the one for integers. It
1207: contains a ``{\it sign\/}'': 0 if the polynomial is equal to~0, and 1 if
1208: not (see however the important remark below), a {\it variable number\/}
1209: (e.g.~0 for $x$, 1 for $y$, etc\dots), and an {\it effective length}.
1210:
1211: \noindent These data can be handled with the following macros:
1212:
1213: \teb{signe} and \teb{setsigne} as for reals and integers.
1214:
1215: \fun{long}{lgef}{GEN z} returns the \idx{effective length} of \kbd{z}.
1216:
1217: \fun{void}{setlgef}{GEN z, long l} sets the effective length of \kbd{z} to
1218: \kbd{l}.
1219:
1220: \fun{long}{varn}{GEN z} returns the variable number of the object \kbd{z}.
1221:
1222: \fun{void}{setvarn}{GEN z, long v} sets the variable number of \kbd{z} to
1223: \kbd{v}.
1224:
1.2 ! noro 1225: The variable numbers encode the relative priorities of variables as discussed
! 1226: in \secref{se:priority}. We will give more details in \secref{se:vars}.
1.1 noro 1227: Note also the function \fun{long}{gvar}{GEN z} which tries to return a
1228: \idx{variable number} for \kbd{z}, even if \kbd{z} is not a polynomial or
1229: power series. The variable number of a scalar type is set by definition
1230: equal to \tet{BIGINT}.
1231:
1232: The components \kbd{z[2]}, \kbd{z[3]},\dots \kbd{z[lgef(z)-1]} point to the
1233: coefficients of the polynomial {\it in ascending order}, with \kbd{z[2]}
1234: being the constant term and so on. Note that the {\it \idx{degree}\/} of the
1235: polynomial is equal to its effective length minus three. The function
1236:
1237: \fun{long}{degree}{GEN x} returns the degree of \kbd{x} with respect to its
1238: main variable even when \kbd{x} is not a polynomial (a rational function
1239: for instance). By convention, the degree of $0$ is~$-1$.
1240:
1241: \misctitle{Important remark}. A zero polynomial can be characterized by the
1242: fact that its sign is~0. However, its effective length may be equal to 2, or
1243: greater than 2. If it is greater than 2, this means that all the coefficients
1244: of the polynomial are equal to zero (as they should for a zero polynomial),
1245: but not all of these zeros are exact zeros, and more precisely the leading
1246: term \kbd{z[lgef(z)-1]} is not an exact zero.
1247:
1248: \subsec{Type \typ{SER} (power series):}\kbdsidx{t_SER}\sidx{power series} This
1249: type also has a second codeword, which encodes a ``{\it sign\/}'', i.e.~0
1250: if the power series is 0, and 1 if not, a {\it variable number\/} as for
1.2 ! noro 1251: polynomials, and an {\it exponent\/}. This information can be handled with
! 1252: the following functions: \teb{signe}, \teb{setsigne}, \teb{varn},
! 1253: \teb{setvarn} as for polynomials, and \teb{valp}, \teb{setvalp} for the
! 1254: exponent as for $p$-adic numbers. Beware: do {\it not\/} use \teb{expo} and
! 1255: \teb{setexpo} on power series.
1.1 noro 1256:
1257: If the power series is non-zero, \kbd{z[2]}, \kbd{z[3]},\dots
1258: \kbd{z[lg(z)-1]} point to the coefficients of \kbd{z} in ascending order,
1259: \kbd{z[2]} being the first non-zero coefficient. Note that the exponent of a
1260: power series can be negative, i.e.~we are then dealing with a Laurent series
1261: (with a finite number of negative terms).
1262:
1263: \subsec{Type \typ{RFRAC} and \typ{RFRACN} (rational function):}%
1264: \kbdsidx{t_RFRAC}\kbdsidx{t_RFRACN}\sidx{rational function}
1265: \kbd{z[1]} points to the numerator, and \kbd{z[2]} on the denominator. The
1266: denominator must be of type polynomial. Note that a type \typ{RFRACN}
1267: rational function can be converted to irreducible form using the function
1268: \teb{gred}.
1269:
1270: \subsec{Type \typ{QFR} (indefinite binary quadratic form):}%
1271: \kbdsidx{t_QFR}\sidx{indefinite binary quadratic form}
1272: \kbd{z[1]}, \kbd{z[2]}, \kbd{z[3]} point to the three coefficients of the
1273: form and should be of type integer. \kbd{z[4]} is Shanks's distance
1274: function, and should be of type real.
1275:
1276: \subsec{Type \typ{QFI} (definite binary quadratic form):}%
1277: \kbdsidx{t_QFI}\sidx{definite binary quadratic form}
1278: \kbd{z[1]}, \kbd{z[2]}, \kbd{z[3]} point to the three coefficients of the
1279: form. All three should be of type integer.
1280:
1281: \subsec{Type \typ{VEC} and \typ{COL} (vector):}%
1282: \kbdsidx{t_VEC}\kbdsidx{t_COL}\sidx{row vector}\sidx{column vector}
1283: \kbd{z[1]}, \kbd{z[2]},\dots \kbd{z[lg(z)-1]} point to the components of
1284: the vector.
1285:
1286: \subsec{Type \typ{MAT} (matrix):}\kbdsidx{t_MAT}\sidx{matrix}
1287: \kbd{z[1]}, \kbd{z[2]},\dots \kbd{z[lg(z)-1]} point to the column vectors
1288: of \kbd{z}, i.e.~they must be of type \typ{COL} and of the same length.
1289:
1290: \noindent The last two were introduced for specific GP use, and you'll be
1291: much better off using the standard malloc'ed C constructs when programming
1292: in library mode. We quote them just for completeness (advising you not to
1293: use them):
1294:
1295: \subsec{Type \typ{LIST} (list):}\kbdsidx{t_LIST}\sidx{list}
1296: This one has a second codeword which contains an effective length (handled
1297: through \teb{lgef}~/ \teb{setlgef}). \kbd{z[2]},\dots, \kbd{z[lgef(z)-1]}
1298: contain the components of the list.
1299:
1300: \subsec{Type \typ{STR} (character string):}%
1301: \kbdsidx{t_STR}\sidx{character string}
1302: \fun{char *}{GSTR}{z} (= \kbd{(z+1)}) points to the first character of the
1303: (\kbd{NULL}-terminated) string.
1304:
1.2 ! noro 1305: \misctitle{Implementation note:} for the types including an exponent (or a
! 1306: valuation), we actually store a biased non-negative exponent (bit-ORing the
! 1307: biased exponent to the codeword), obtained by adding a constant to the true
! 1308: exponent: either \kbd{HIGHEXPOBIT} (for \typ{REAL}) or \kbd{HIGHVALPBIT} (for
! 1309: \typ{PADIC} and \typ{SER}). Of course, this is encapsulated by the
! 1310: exponent/valuation-handling macros and need not concern the library user.
! 1311:
1.1 noro 1312: \section{PARI variables}\label{se:vars}
1.2 ! noro 1313: \subsec{Multivariate objects}\sidx{variable (priority)}
1.1 noro 1314:
1315: \noindent
1.2 ! noro 1316: We now consider variables and formal computations, and give the technical
! 1317: details corresponding to the general discussion in \secref{se:priority}. As
! 1318: we have seen in \secref{se:impl}, the codewords for types \typ{POL} and
! 1319: \typ{SER} encode a ``variable number''. This is an integer, ranging from $0$
! 1320: to \kbd{MAXVARN}. The lower it is, the higher the variable priority.
1.1 noro 1321:
1322: In fact, the way an object will be considered in formal computations depends
1323: entirely on its ``principal variable number'' which is given by the function
1324:
1325: \fun{long}{gvar}{GEN z}
1326:
1327: \noindent which returns a \idx{variable number} for \kbd{z}, even if \kbd{z}
1328: is not a polynomial or power series. The variable number of a scalar type is
1329: set by definition equal to \tet{BIGINT} which is bigger than any legal
1330: variable number. The variable number of a recursive type which is not a
1331: polynomial or power series is the minimal variable number of its components.
1.2 ! noro 1332: But for polynomials and power series only the ``outermost'' number counts (we
! 1333: directly access $\tet{varn}(x)$ in the codewords): the representation is not
! 1334: symmetrical at all.
1.1 noro 1335:
1336: Under GP, one need not worry too much since the interpreter will define
1.2 ! noro 1337: the variables as it sees them\footnote{*}{
! 1338: More precisely, the first time a given identifier is read by the GP parser
! 1339: (and is not immediately interpreted as a function) a new variable is created,
! 1340: and it is assigned a strictly lower priority than any variable in use at this
! 1341: point. On startup, before any user input has taken place, 'x' is defined in
! 1342: this way and will thus always have maximal priority (and variable number $0$).}
! 1343: %
! 1344: and do the right thing with the polynomials produced (however, have a look at
! 1345: the remark in \secref{se:rempolmod}).
! 1346:
! 1347: But in library mode, they are tricky objects if you intend to build
! 1348: polynomials yourself (and not just let PARI functions produce them, which is
! 1349: usually less efficient). For instance, it does not make sense to have a
! 1350: variable number occur in the components of a polynomial whose main variable
! 1351: has a higher number (lower priority), even though there's nothing PARI can do
! 1352: to prevent you from doing it.
1.1 noro 1353:
1354: \subsec{Creating variables}
1355: A basic difficulty is to ``create'' a variable. As we have seen in
1356: \secref{se:intro4}, a plethora of objects is associated to variable
1357: number~$v$. Here is the complete list: \teb{polun}$[v]$ and
1358: \teb{polx}$[v]$, which you can use in library mode and which represent,
1359: respectively, the monic monomials of degrees 0 and 1 in~$v$;
1360: \teb{varentries}$[v]$, and \teb{polvar}$[v]$. The latter two are only
1361: meaningful to GP, but they have to be set nevertheless. All of them must be
1362: properly defined before you can use a given integer as a variable number.
1363:
1364: Initially, this is done for $0$ (the variable \kbd{x} under GP), and
1365: \tet{MAXVARN}, which is there to address the need for a ``temporary'' new
1.2 ! noro 1366: variable in library mode and cannot be input under GP.
! 1367: No documented library function can create from scratch an object involving
! 1368: \tet{MAXVARN} (of course, if the operands originally involve \kbd{MAXVARN},
! 1369: the function will abide). We call the latter type a ``temporary variable''.
! 1370: The regular variables meant to be used in regular objects, are called ``user
1.1 noro 1371: variables\sidx{variable (user)}''.
1372:
1373: \subsubsec{User variables}\sidx{variable (user)}:
1374: When the program starts, \kbd{x} is the only user variable (number~$0$). To
1375: define new ones, use
1376:
1377: \fun{long}{fetch_user_var}{char *$s$}
1378:
1379: \noindent which inspects the user variable named $s$ (creating it if
1380: needed), and returns its variable number.
1381: \bprog
1382: long v = fetch_user_var("y");
1383: GEN gy = polx[v];
1384: @eprog
1385: This function raises an error if $s$ is already known as a function name to
1386: the interpreter.
1387:
1388: \misctitle{Caveat:} it is possible to use \teb{flissexpr}
1389: (see~\secref{se:flisexpr}) to execute a GP command and create GP variables
1390: on the fly as needed:
1391:
1392: \bprog
1393: GEN gy = flissexpr("y"); /*@Ccom supposedly returns polx[$v$], for some $v$ */
1394: long v = gvar(gy);
1395: @eprog
1396:
1397: \noindent This is dangerous, especially when programming functions that
1398: will be used under GP. The code above reads the value of \kbd{y}, as it is
1399: currently known by the GP interpreter (possibly creating it in the
1400: process). All is well and good if \kbd{y} hasn't been tampered with in
1401: previous GP commands. But if \kbd{y} has been modified (e.g \kbd {y = 1}),
1402: then the value of \kbd{gy} is not what you expected it to be and corresponds
1403: instead to the current value of the GP variable (e.g \kbd{gun}).
1404:
1405: \subsubsec{Temporary variables}\sidx{variable (temporary)}:
1406: \kbd{MAXVARN} is available, but is better left to pari internal functions
1407: (some of which don't check that \kbd{MAXVARN} is free for them to use,
1408: which can be considered a bug). You can create more temporary variables
1409: using
1410:
1411: \fun{long}{fetch_var}{}\label{se:fetch_var}
1412:
1413: \noindent
1414: This returns a variable number which is guaranteed to be unused by the
1415: library at the time you get it and as long as you do not delete it (we'll see
1416: how to do that shortly). This has {\it lower\/} number (i.e.~{\it higher\/}
1417: priority) than any temporary variable produced so far (\kbd{MAXVARN} is
1418: assumed to be the first such). This call updates all the aforementioned
1419: internal arrays. In particular, after the statement \kbd{v = fetch\_var()},
1420: you can use \kbd{polun[v]} and \kbd{polx[v]}. The variables created in this
1421: way have no identifier assigned to them though, and they will be printed as
1422: \kbd{\#<\text{number}>}, except for \kbd{MAXVARN} which will be printed
1423: as~\kbd{\#}. You can assign a name to a temporary variable, after creating
1424: it, by calling the function
1425:
1426: \fun{void}{name_var}{long n, char *s}
1427:
1428: \noindent after which the output machinery will use the name \kbd{s} to
1429: represent the variable number~\kbd{n}. The GP parser will {\it not\/}
1430: recognize it by that name, however, and calling this on a variable known
1431: to~GP will raise an error. Temporary variables are meant to be used as free
1432: variables, and you should never assign values or functions to them as you
1433: would do with variables under~GP. For that, you need a user variable.
1434:
1435: All objects created by \kbd{fetch\_var} are on the heap and not on the stack,
1436: thus they are not subject to standard garbage collecting (they won't be
1437: destroyed by a \kbd{gerepile} or \kbd{avma = ltop} statement). When you don't
1438: need a variable number anymore, you can delete it using
1439:
1440: \fun{long}{delete_var}{}
1441:
1442: \noindent which deletes the {\it latest\/} temporary variable created and
1443: returns the variable number of the previous one (or simply returns 0 if you
1444: try, in vain, to delete \kbd{MAXVARN}). Of course you should make sure that
1445: the deleted variable does not appear anywhere in the objects you use later
1446: on. Here is an example:
1447:
1448: \bprog
1449: {
1450: long first = fetch_var();
1451: long n1 = fetch_var();
1452: long n2 = fetch_var(); /*@Ccom prepare three variables for internal use */
1453: ...
1454: /*@Ccom delete all variables before leaving */
1455: do { num = delete_var(); } while (num && num <= first);
1456: }@eprog
1457:
1458: \noindent
1459: The (dangerous) statement
1460:
1461: \bprog
1462: while (delete_var()) /*@Ccom empty */;
1463: @eprog
1464:
1465: \noindent removes all temporary variables that were in use, except
1466: \kbd{MAXVARN} which cannot be deleted.
1467:
1468: \section{Input and output}
1469:
1470: \noindent
1471: Two important aspects have not yet been explained which are specific to
1472: library mode: input and output of PARI objects.
1473:
1474: \subsec{Input}.
1475:
1476: \noindent
1477: For \idx{input}, PARI provides you with two powerful high level functions
1478: which enables you to input your objects as if you were under GP. In fact,
1479: the second one {\it is\/} essentially the GP syntactical parser, hence you
1480: can use it not only for input but for (most) computations that you can do
1481: under GP. These functions are called \teb{flisexpr} and \teb{flisseq}. The
1482: first one has the following syntax:\label{se:flisexpr}
1483:
1484: \fun{GEN}{flisexpr}{char *s}
1485:
1486: \noindent
1487: Its effect is to analyze the input string s and to compute the result as in
1488: GP. However it is limited to one expression. If you want to read and
1489: evaluate a sequence of expressions, use
1490:
1491: \fun{GEN}{flisseq}{char *s}
1492:
1493: \noindent\sidx{filter}
1494: In fact these two functions start by {\it filtering\/} out all spaces and
1495: comments in the input string (that's what the initial \kbd{f} stands for).
1496: They then call the underlying basic functions, the GP parser proper:
1497: \fun{GEN}{lisexpr}{char *s} and \fun{GEN}{lisseq}{char *s}, which are
1498: slightly faster but which you probably don't need.
1499:
1500: To read a \kbd{GEN} from a file, you can use the simpler interface
1501:
1502: \fun{GEN}{lisGEN}{FILE *file}
1503:
1504: which reads a character string of arbitrary length from the stream \kbd{file}
1505: (up to the first newline character), applies \kbd{flisexpr} to it, and
1506: returns the resulting \kbd{GEN}. This way, you won't have to worry about
1507: allocating buffers to hold the string. To interactively input an expression,
1508: use \kbd{lisGEN(stdin)}. This function returns \kbd{NULL} if \kbd{EOF} is
1509: encountered before a complete expression could be read.
1510:
1511: Once in a while, it may be necessary to evaluate a GP expression sequence
1512: involving a call to a function you have defined in~C. This is easy using
1513: \teb{install} which allows you to manipulate quite an arbitrary function (GP
1514: knows about pointers!). The syntax is
1515:
1516: \fun{void}{install}{void *f, char *name, char *code}
1517:
1518: \noindent where \kbd{f} is the (address of) the function (cast to the C type
1519: \kbd{void*}), \kbd{name} is the name by which you want to access your
1520: function from within your GP expressions, and \kbd{code} is a character
1521: string describing the function call prototype (see~\secref{se:gp.interface}
1522: for the precise description of prototype strings). In case the function
1523: returns a \kbd{GEN}, it should satisfy \kbd{gerepileupto} assumptions (see
1524: \secref{se:garbage}).
1525:
1526: \subsec{Output}.
1527:
1528: \noindent
1529: For \idx{output}, there exist essentially three different functions (with
1530: variants), corresponding to the three main GP output formats (as described in
1531: \secref{se:output}), plus three extra ones, respectively devoted to
1532: \TeX\ output, string output, and (advanced) debugging.
1533:
1534: \noindent $\bullet$ ``raw'' format, obtained by using the function
1535: \teb{brute} with the following syntax:
1536:
1537: \fun{void}{brute}{GEN obj, char x, long n}
1538:
1539: \noindent
1540: This prints the PARI object \kbd{obj} in \idx{format} \kbd{x0.n}, using the
1541: notations from \secref{se:format}. Recall that here \kbd{x} is either
1542: \kbd{'e'}, \kbd{'f'} or \kbd{'g'} corresponding to the three numerical output
1543: formats, and \kbd{n} is the number of printed significant digits, and should
1544: be set to $-1$ if all of them are wanted (these arguments only affect the
1545: printing of real numbers). Usually you won't need that much flexibility, so
1546: most of the time you will get by with the function
1547:
1548: \fun{void}{outbrute}{GEN obj}, which is equivalent to \kbd{brute(x,'g',-1)},
1549:
1550: \noindent or even better, with
1551:
1552: \fun{void}{output}{GEN obj} which is equivalent to \kbd{outbrute(obj)}
1553: followed by a newline and a buffer flush. This is especially nice during
1554: debugging. For instance using \kbd{dbx} or \kbd{gdb}, if \kbd{obj} is a
1555: \kbd{GEN}, typing \kbd{print output(obj)} will enable you to see the
1556: content of \kbd{obj} (provided the optimizer has not put it into a
1557: register, but it's rarely a good idea to debug optimized code).
1558:
1559: \noindent $\bullet$ ``prettymatrix'' format: this format is identical to the
1560: preceding one except for matrices. The relevant functions are:
1561:
1562: \fun{void}{matbrute}{GEN obj, char x, long n}
1563:
1564: \fun{void}{outmat}{GEN obj}, which is followed by a newline and a buffer flush.
1565:
1566: \noindent $\bullet$ ``prettyprint'' format: the basic function has an
1567: additional parameter \kbd{m}, corresponding to the (minimum) field width
1568: used for printing integers:
1569:
1570: \fun{void}{sor}{GEN obj, char x, long n, long m}
1571:
1572: \noindent The simplified version is
1573:
1574: \fun{void}{outbeaut}{GEN obj} which is equivalent to
1575: \kbd{sor(obj,'g',-1,0)} followed by a newline and a buffer flush.
1576:
1577: \noindent $\bullet$ The first extra format corresponds to the \teb{texprint}
1578: function of GP, and gives a \TeX\ output of the result. It is obtained by
1579: using:
1580:
1581: \fun{void}{exe}{GEN obj, char x, long n}
1582:
1583: \noindent $\bullet$ The second one is the function \teb{GENtostr} which
1584: converts a PARI \kbd{GEN} to an ASCII string. The syntax is
1585:
1586: \fun{char*}{GENtostr}{GEN obj}, wich returns a \kbd{malloc}'ed character
1587: string (which you should \kbd{free} after use).
1588:
1589: \noindent $\bullet$ The third and final one outputs the \idx{hexadecimal tree}
1590: corresponding to the GP command \kbd{\b x} using the function
1591:
1592: \fun{void}{voir}{GEN obj, long nb}, which will only output the first
1593: \kbd{nb} words corresponding to leaves (very handy when you have a look at
1594: big recursive structures). If you set this parameter to $-1$ all
1595: significant words will be printed. Usually this last type of output would
1596: only be used for debugging purposes.
1597:
1598: \misctitle{Remark}. Apart from \teb{GENtostr}, all PARI output is done on
1599: the stream \teb{outfile}, which by default is initialized to \teb{stdout}. If
1600: you want that your output be directed to another file, you should use the
1601: function \fun{void}{switchout}{char *name} where \kbd{name} is a
1602: character string giving the name of the file you are going to use. The
1603: output will be {\it appended\/} at the end of the file. In order to close
1604: the file, simply call \kbd{switchout(NULL)}.
1605:
1606: Similarly, errors are sent to the stream \teb{errfile} (\teb{stderr}
1607: by default), and input is done on the stream \teb{infile}, which you can change
1608: using the function \teb{switchin} which is analogous to \teb{switchout}.
1609:
1610: \misctitle{(Advanced) Remark}. All output is done according to the values
1611: of the \teb{pariOut}~/ \teb{pariErr} global variables which are pointers to
1612: structs of pointer to functions. If you really intend to use these, this
1613: probably means you are rewriting GP. In that case, have a look at the code in
1614: \kbd{language/es.c} (\kbd{init80()} or \kbd{GENtostr()} for instance).
1615:
1616: \subsec{Errors}.\sidx{error}\kbdsidx{talker}
1617:
1618: \noindent
1619: If you want your functions to issue error messages, you can use the general
1620: error handling routine \teb{err}. The basic syntax is
1621: %
1622: \bprog
1623: err(talker, "error message");
1624: @eprog
1625:
1626: \noindent
1627: This will print the corresponding error message and exit the program (in
1628: library mode; go back to the GP prompt otherwise).\label{se:err} You can
1629: also use it in the more versatile guise
1630: \bprog
1631: err(talker, format, ...);
1632: @eprog\noindent
1633: where \kbd{format} describes the format to use to write the remaining
1634: operands, as in the \teb{printf} function (however, see the next section).
1635: The simple syntax above is just a special case with a constant format and no
1636: remaining arguments.
1637:
1638: \noindent
1639: The general syntax is
1640:
1641: \fun{void}{err}{numerr,...}
1642:
1643: \noindent where \kbd{numerr} is a codeword which indicates what to do with
1644: the remaining arguments and what message to print. The list of valid keywords
1645: is in \kbd{language/errmessages.c} together with the basic corresponding
1646: message. For instance, \kbd{err(typeer,"matexp")} will print the message:
1647:
1648: \kbd{ *** incorrect type in matexp.}
1649:
1650: \noindent
1651: Among the codewords are {\it warning\/} keywords (all those which start with
1652: the prefix \kbd{warn}). In that case, \teb{err} does {\it not\/} abort the
1653: computation, just print the requested message and go on. The basic example is
1654:
1655: \kbd{err(warner, "Strategy 1 failed. Trying strategy 2")}
1656:
1657: \noindent which is the exact equivalent of \kbd{err(talker,...)} except that
1658: you certainly don't want to stop the program at this point, just inform the
1.2 ! noro 1659: user that something important has occurred (in particular, this output would be
1.1 noro 1660: suitably highlighted under GP, whereas a simple \kbd{printf} would not).
1661:
1662: \subsec{Debugging output}.\sidx{debugging}\sidx{format}\label{se:dbg_output}
1663:
1664: \noindent
1665: The global variables \teb{DEBUGLEVEL} and \teb{DEBUGMEM} (corresponding
1666: to the default \teb{debug} and \teb{debugmem}, see \secref{se:defaults})
1667: are used throughout the PARI code to govern the amount of diagnostic and
1668: debugging output, depending on their values. You can use them to debug your
1669: own functions, especially after having made them accessible under GP through
1670: the command \teb{install} (see \secref{se:install}).
1671:
1672: For debugging output, you can use \kbd{printf} and the standard output
1673: functions (\teb{brute} or \teb{output} mainly), but also some special purpose
1674: functions which embody both concepts, the main one being
1675:
1676: \fun{void}{fprintferr}{char *pariformat, ...}
1677:
1678: \noindent
1679: Now let's define what a PARI format is. It is a character string, similar
1680: to the one \kbd{printf} uses, where \kbd{\%} characters have a special
1681: meaning. It describes the format to use when printing the remaining operands.
1682: But, in addition to the standard format types, you can use \kbd{\%Z} to
1683: denote a \kbd{GEN} object (we would have liked to pick \kbd{\%G} but it was
1684: already in use!). For instance you could write:
1685: \bprog
1686: err(talker, "x[%d] = %Z is not invertible!", i, x[i])
1687: @eprog
1688: \noindent since the \teb{err} function accepts PARI formats. Here \kbd{i}
1689: is an \kbd{int}, \kbd{x} a \kbd{GEN} which is not a leaf and this would
1690: insert in raw format the value of the \kbd{GEN} \kbd{x[i]}.
1691:
1692: \subsec{Timers and timing output}.
1693:
1694: \noindent
1695: To profile your functions, you can use the PARI timer. The functions
1696: \fun{long}{timer}{} and \fun{long}{timer2}{} return the elapsed time since
1697: the last call of the same function (in milliseconds). Two different
1698: functions (identical except for their independent time-of-last-call
1699: memories!) are provided so you can have both global timing and fine tuned
1700: profiling.
1701:
1702: You can also use \fun{void}{msgtimer}{char *format,...}, which prints
1703: prints \kbd{Time}, then the remaining arguments as specified by
1704: \kbd{format} (which is a PARI format), then the output of \kbd{timer2}.
1705:
1706: This mechanism is simple to use but not foolproof. If some other function
1.2 ! noro 1707: uses these timers, and many PARI functions do use \kbd{timer2} when
! 1708: \tet{DEBUGLEVEL} is high enough, the timings will be meaningless. To handle
! 1709: timing in a reentrant way, PARI defines a dedicated datatype,
! 1710: \tet{pari_timer}. The functions
! 1711:
! 1712: \fun{long}{TIMER}{pari\_timer *T}
! 1713:
! 1714: \fun{long}{msgTIMER}{pari\_timer *T, char *format,...}
! 1715:
! 1716: \noindent are equivalent to \kbd{timer} and \kbd{msgtimer} respectively,
! 1717: except they use a unique timer \kbd{T} containing all the information needed,
! 1718: so that no other function can mess with your timings. They are used as
! 1719: follows:
! 1720: \bprog
! 1721: pari_timer T;
! 1722:
! 1723: (void)TIMER(&T); /* initialize timer */
! 1724: ...
! 1725: printf("Total time: %ld\n", TIMER(&T));
! 1726: @eprog
! 1727: or
! 1728: \bprog
! 1729: pari_timer T;
! 1730: long i;
! 1731: GEN L;
! 1732:
! 1733: (void)TIMER(&T); /* initialize timer */
! 1734: for (i = 1; i < 10; i++) {
! 1735: ...
! 1736: msgTIMER(&T, "for i = %ld (L[i] = %Z)", i, L[i]);
! 1737: }
! 1738: @eprog
1.1 noro 1739:
1740: \section{A complete program}
1741: \label{se:prog}
1742:
1743: \noindent
1744: Now that the preliminaries are out of the way, the best way to learn how to
1745: use the library mode is to work through a detailed non-trivial example of a
1746: main program. We will write a program which computes the exponential of a
1747: square matrix~$x$. The complete listing is given in Appendix~B, but each
1748: part of the program will be produced and explained here. We will use an
1749: algorithm which is not optimal but is not far from the one used for the PARI
1750: function \teb{gexp} (in fact embodied in the function \kbd{mpexp1}). This
1751: consists in calculating the sum of the series:
1752: $$e^{x/(2^n)}=\sum_{k=0}^\infty \dfrac{(x/(2^n))^k}{k!}$$
1753: for a suitable positive integer $n$, and then computing $e^x$ by repeated
1754: squarings. First, we will need to compute the $L^2$-norm of the matrix~$x$,
1755: i.e.~the quantity:
1756: $$z=\|x\|_2=\sqrt{\sum x_{i,j}^2}.$$
1757: We will then choose the integer $n$ such that the $L^2$-norm of $x/(2^n)$ is
1758: less than or equal to~1, i.e.
1759: $$ n = \left\lceil{\ln(z)}\big/{\ln(2)}\right\rceil $$
1760: if $z\ge1$, and~$n=0$ otherwise. Then the series will converge at least as
1761: fast as the usual one for $e^1$, and the cutoff error will be easy to
1762: estimate. In fact a larger value of $n$ would be preferable, but this is
1763: slightly machine dependent and more complicated, and will be left to the
1764: reader.
1765:
1766: Let us start writing our program. So as to be able to use it in other
1767: contexts, we will structure it in the following way: a main program which
1768: will do the input and output, and a function which we shall call \kbd{matexp}
1769: which does the real work. The main program is easy to write. It can be
1770: something like this:
1771:
1772: \bprog
1773: #include <pari.h>
1774:
1775: GEN matexp(GEN x, long prec);
1776:
1777: int
1778: main()
1779: {
1780: long d, prec = 3;
1781: GEN x;
1782:
1783: /*@Ccom take a stack of $10^6$ bytes, no prime table */
1784: pari_init(1000000, 2);
1785: printf("precision of the computation in decimal digits:\n");
1786: d = itos(lisGEN(stdin));
1787: if (d > 0) prec = (long) (d*pariK1+3);
1788:
1789: printf("input your matrix in GP format:\n");
1790: x = matexp(lisGEN(stdin), prec);
1791:
1792: sor(x, 'g', d, 0);
1793: exit(0);
1794: }@eprog
1795: \noindent
1796: The variable \kbd{prec} represents the length in longwords of the real
1797: numbers used. \teb{pariK1} is a constant (defined in \kbd{paricom.h}) equal
1798: to $\ln(10) / (\ln(2) * \kbd{BITS\_IN\_LONG})$, which allows us to convert
1799: from a number of decimal digits to a number of longwords, independently of
1800: the actual bit size of your long integers. The function \teb{lisGEN} reads an
1801: expression (here from standard input) and converts it to a \kbd{GEN}, like
1802: the GP parser itself would. This means it takes care of whitespace etc.\ in
1803: the input, and can do computations (e.g.~\kbd{matid(2)} or \kbd{[1,0; 0,1]}
1804: are equally valid inputs).
1805:
1806: Finally, \kbd{sor} is the general output routine. We have chosen to give
1807: \kbd{d} significant digits since this is what was asked for. Note that there
1808: is a trick hidden here: if a negative \kbd{d} was input, then the computation
1809: will be done in precision 3 (i.e.~about 9.7 decimal digits for 32-bit
1810: machines and 19.4 for 64-bit machines) and in the function \kbd{sor}, giving
1811: a negative third argument outputs all the significant digits, which is entirely
1812: appropriate. Now let us attack the main course, the function \kbd{matexp}:
1813: %
1814: \bprog
1815: GEN
1816: matexp(GEN x, long prec)
1817: {
1818: long lx=lg(x),i,k,n,lbot, ltop = avma;
1819: GEN y,r,s,p1,p2;
1820:
1821: /*@Ccom check that x is a square matrix */
1822: if (typ(x) != t_MAT) err(talker,"this expression is not a matrix");
1823: if (lx == 1) return cgetg(1, t_MAT);
1824: if (lx != lg(x[1])) err(talker,"not a square matrix");
1825:
1826: /*@Ccom compute the $L_2$ norm of x */
1827: s = gzero;
1828: for (i=1; i<lx; i++)
1829: s = gadd(s, gnorml2((GEN)x[i]));
1830: if (typ(s) == t_REAL) setlg(s,3);
1831: s = gsqrt(s,3); /*@Ccom we do not need much precision on s */
1832:
1833: /*@Ccom if $s<1$, we are happy */
1834: k = expo(s);
1835: if (k < 0) { n = 0; p1 = x; }
1836: else { n = k+1; p1 = gmul2n(x,-n); setexpo(s,-1); }@eprog
1837: \noindent
1838: Before continuing, several remarks are in order.
1839:
1840: First, before starting this computation which will produce garbage on the
1841: stack, we have carefully saved the value of the stack pointer \kbd{avma} in
1842: \kbd{ltop}. Note that we are going to assume throughout that the garbage does
1843: not overflow the currently available stack. If it ever did, we would
1844: have several options~--- allocate a larger stack in the main program (for
1845: instance change 1000000 into 2000000), do some \kbd{gerepile}ing along the
1846: way, or (if you know what you are doing) use \tet{allocatemoremem}.
1847:
1848: Secondly, the \teb{err} function is the general error handler for the
1849: PARI library. This will abort the program after printing the required
1850: message.
1851:
1852: Thirdly, notice how we handle the special case $\kbd{lx} = 1$ (empty matrix)
1853: {\it before\/} accessing \kbd{lx(x[1])}. Doing it the other way round could
1854: produce a fatal error (a segmentation fault or a bus error, most probably).
1855: Indeed, if \kbd{x} is of length 1, then \kbd{x[1]} is not a component of
1856: \kbd{x}. It is just the contents of the memory cell which happens to follow
1857: the one pointed to by \kbd{x}, and thus has no reason to be a valid \kbd{GEN}.
1858: Now recall that none of the codeword handling macros do any kind of type
1859: checking (see \secref{se:impl}), thus \teb{lg} would consider \kbd{x[1]}
1860: as a valid address, and try to access \kbd{*((GEN)x[1])} (the first codeword)
1861: which is unlikely to be a legal memory address.
1862:
1863: In the fourth place, to compute the square of the $L^2$-norm of \kbd{x} we
1864: just add the squares of the $L^2$-norms of the column vectors which we obtain
1865: using the library function \teb{gnorml2}. Had this function not existed, the
1866: norm computation would of course have been just as easy to write, but we would
1867: have needed a double loop.
1868:
1869: We then take the square root of \kbd{s}, in precision~3 (the smallest
1870: possible). The \kbd{prec} argument of transcendental functions (here~$3$) is
1871: only taken into account when the arguments are {\it exact\/} objects, and thus
1872: no a priori precision can be determined from the objects themselves. To cater
1873: for this possibility, if \kbd{s} is of type \typ{REAL}, we use the function
1874: \teb{setlg} which effectively sets the precision of \kbd{s} to the required
1875: value. Note that here, since we are using a numeric value for a \kbd{cget}
1876: function, the program will run slightly differently on 32-bit and 64-bit
1877: machines: we want to use the smallest possible bit accuracy, and this is equal
1878: to \kbd{BITS\_IN\_LONG}.
1879:
1880: Note that the matrix \kbd{x} is allowed to have complex entries, but the
1881: function \kbd{gnorml2} guarantees that \kbd{s} is a non-negative real number
1882: (not necessarily of type \typ{REAL} of course). If we had not known this fact,
1883: we would simply have added the instruction \kbd{s = greal(s);} just after the
1884: \kbd{for} loop.
1885:
1886: Note also that the function \kbd{gnorml2} works as desired on matrices, so we
1887: really did not need this loop at all (\kbd{s = gnorml2(x)} would have been
1888: enough), but we wanted to give examples of function usage. Similarly, it is of
1889: course not necessary to take the square root for testing whether the norm
1890: exceeds~$1$.
1891:
1892: In the fifth place, note that we initialized the sum \kbd{s} to \teb{gzero},
1893: which is an exact zero. This is logical, but has some disadvantages: if all
1894: the entries of the matrix are integers (or rational numbers), the computation
1895: will take rather long, about twice as long as with real numbers of the same
1896: length. It would be better to initialize \kbd{s} to a real zero, using for
1897: instance the instructions:
1898:
1899: \kbd{s = cgetr(prec+1); gaffsg(0,s);}
1900:
1901: \noindent
1902: This raises the question: which real zero does this produce (have a look
1903: at~\secref{se:whatzero})? In fact, the following choice has been made: it will
1904: give you the zero with exponent equal to $-\kbd{BITS\_IN\_LONG}$ times the
1905: number of longwords in the mantissa, i.e.~$-\kbd{bit\_accuracy(lg(s))}$.
1906: Instead of the above idiom, you can also use the function
1907: \fun{GEN}{realzero}{long prec}, which simply returns a real zero to
1908: accuracy $-\tet{bit_accuracy}(\var{prec})$.
1909:
1910: The sixth remark here is about how to determine the approximate size
1911: of a real number. The fastest way to do this is to look at its binary
1912: exponent. Hence we need to have \kbd{s} actually represented as a real number,
1913: and not as an integer or a rational number. The result of transcendental
1914: functions is guaranteed to be of type \typ{REAL}, or complex with \typ{REAL}
1915: components, thus this is indeed the case after the call to \kbd{gsqrt} since
1916: its argument is a nonnegative (real) number.
1917:
1918: Finally, note the use of the function \teb{gmul2n}. It has the following
1919: syntax:
1920:
1921: \fun{GEN}{gmul2n}{GEN x, long n}
1922:
1923: \noindent and the effect is simply to multiply \kbd{x} by $2^{\text n}$,
1924: where \kbd{n} can be positive or negative. This is much faster than
1925: \kbd{gmul} or \kbd{gmulgs}.
1926:
1927: There is another function \teb{gshift} with exactly the same syntax.
1928: When \kbd{n} is non-negative, the effects of these two functions are the same.
1929: However, when \kbd{n} is negative, \kbd{gshift} acts like a right shift of
1.2 ! noro 1930: \kbd{-n}, hence does not normally perform an exact division on integers. The
1.1 noro 1931: function \kbd{gshift} is the PARI analogue of the C or GP operators \kbd{<<}
1932: and~\kbd{>>}.
1933:
1934: We now come to the heart of the function. We have a \kbd{GEN} \kbd{p1} which
1935: points to a certain matrix of which we want to take the exponential. We will
1936: want to transform this matrix into a matrix with real (or complex of real)
1937: entries before starting the computation. To do this, we simply multiply by
1938: the real number 1 in precision $\kbd{prec}+1$ (to be on the side of safety).
1939: To sum the series, we will use three variables: a variable \kbd{p2} which at
1940: stage $k$ will contain $\kbd{p1}^k/k!$, a variable \kbd{y} which will contain
1941: $\sum_{i=0}^k \kbd{p1}^i/i!$, and a variable \kbd{r} which will contain the
1942: size estimate $\kbd{s}^k/k!$. Note that we do not use Horner's rule. This is
1943: simply because we are lazy and do not want to compute in advance the number
1944: of terms that we need. We leave this modification (and many other
1945: improvements!) to the reader. The program continues as follows:
1946:
1947: \bprog
1948: /*@Ccom initializations before the loop */
1949: r = cgetr(prec+1); gaffsg(1,r); p1 = gmul(r,p1);
1950: y = gscalmat(r,lx-1); /*@Ccom creates scalar matrix with r on diagonal */
1951: p2 = p1; r = s; k = 1;
1952: y = gadd(y,p2);
1953:
1954: /*@Ccom now the main loop */
1955: while (expo(r) >= -BITS_IN_LONG*(prec-1))
1956: {
1957: k++; p2 = gdivgs(gmul(p2,p1),k);
1958: r = gdivgs(gmul(s,r),k); y = gadd(y,p2);
1959: }
1960:
1961: /*@Ccom now square back n times if necessary */
1962: if (!n) { lbot = avma; y = gcopy(y); }
1963: else
1964: {
1965: for (i=0; i<n; i++) { lbot = avma; y = gsqr(y); }
1966: }
1967: return gerepile(ltop,lbot,y);
1968: }@eprog
1969: \noindent
1970: A few remarks once again. First note the use of the function
1971: \teb{gscalmat} with the following syntax:
1972:
1973: \fun{GEN}{gscalmat}{GEN x, long m}
1974:
1975: \noindent
1976: The effect of this function is to create the $\kbd{m}\times\kbd{m}$ scalar
1977: matrix whose diagonal entries are~\kbd{x}. Hence the length of the matrix
1978: including the codeword will in fact be \kbd{m+1}. There is a corresponding
1979: function \teb{gscalsmat} which takes a long as a first argument.
1980:
1981: If we refer to what has been said above, the main loop should be self-evident.
1982:
1983: When we do the final squarings, according to the fundamental dogma on the
1984: use of \teb{gerepile}, we keep the value of \kbd{avma} in \kbd{lbot} just {\it
1985: before\/} the squaring, so that if it is the last one, \kbd{lbot} will indeed
1986: be the bottom address of the garbage pile, and \kbd{gerepile} will work. Note
1987: that it takes a completely negligible time to do this in each loop compared
1988: to a matrix squaring. However, when \kbd{n} is initially equal to 0, no
1989: squaring has to be done, and we have our final result ready but we lost the
1990: address of the bottom of the garbage pile. Hence we use the trick of copying
1991: \kbd{y} again to the top of the stack. This is inefficient, but does the trick.
1992: If we wanted to avoid this using only \kbd{gerepile}, the best thing to do
1993: would be to put the instruction \kbd{lbot=avma} just before both occurrences
1994: of the instruction \kbd{y=gadd(p2,y)}. Of course, we could also rewrite the
1995: last block as follows:
1996:
1997: \bprog
1998: /*@Ccom now square back n times */
1999: for (i=0; i<n; i++) y = gsqr(y);
2000: return gerepileupto(ltop,y);
2001: @eprog
2002: \noindent
2003: because it does not matter to \teb{gerepileupto} that we have lost the address
2004: just before the final result (note that the loop is not executed if \kbd{n}
2005: is~0). It is safe to use \kbd{gerepileupto} here as \kbd{y} will have been
2006: created by either \kbd{gsqr} or \kbd{gadd}, both of which are guaranteed to
2007: return suitable objects.
2008:
2009: \misctitle{Remarks}. As such, the program should work most of the time if
2010: \kbd{x} is a square matrix with real or complex entries. Indeed, since
2011: essentially the first thing that we do is to multiply by the real number~1,
2012: the program should work for integer, real, rational, complex or quadratic
2013: entries. This is in accordance with the behavior of transcendental functions.
2014:
2015: Furthermore, since this program is intended to be only an illustrative
2016: example, it has been written a little sloppily. In particular many error
2017: checks have been omitted, and the efficiency is far from optimal. An evident
2018: improvement would be the use of \kbd{gerepileupto} mentioned above. Another
2019: improvement is to multiply the matrix x by the real number 1 right at the
2020: beginning, speeding up the computation of the $L^2$-norm in many cases. These
2021: improvements are included in the version given in Appendix~B. Still another
2022: improvement would come from a better choice of~\kbd{n}. If the reader takes a
2023: look at the implementation of the function \kbd{mpexp1} in the file
2024: \kbd{basemath/trans1.c}, he can make the necessary changes himself. Finally,
2025: there exist other algorithms of a different nature to compute the exponential
2026: of a matrix.
2027:
2028: \section{Adding functions to PARI}
2029: \subsec{Nota Bene}.
2030: %
2031: As mentioned in the \kbd{COPYING} file, modified versions of the PARI package
2032: can be distributed under the conditions of the GNU General Public License. If
2033: you do modify PARI, however, it is certainly for a good reason, hence we
2034: would like to know about it, so that everyone can benefit from it. There is
2035: then a good chance that the modifications that you have made will be
2036: incorporated into the next release.
2037:
2038: (Recall the e-mail address: \kbd{pari@math.u-bordeaux.fr}, or use the mailing
2039: lists).
2040:
2041: Roughly four types of modifications can be made. The first type includes all
2042: improvements to the documentation, in a broad sense. This includes correcting
1.2 ! noro 2043: typos or inaccurracies of course, but also items which are not really covered
1.1 noro 2044: in this document, e.g.~if you happen to write a tutorial, or pieces of code
2045: exemplifying some fine points that you think were unduly omitted.
2046:
2047: The second type is to expand or modify the configuration routines and skeleton
2048: files (the \kbd{Configure} script and anything in the \kbd{config/}
2049: subdirectory) so that compilation is possible (or easier, or more efficient)
2050: on an operating system previously not catered for. This includes discovering
2051: and removing idiosyncrasies in the code that would hinder its portability.
2052:
2053: The third type is to modify existing (mathematical) code, either to correct
2054: bugs, to add new functionalities to existing functions, or to improve their
2055: efficiency.
2056:
2057: Finally the last type is to add new functions to PARI. We explain here how
2058: to do this, so that in particular the new function can be called from GP.
2059:
2060: \subsec{The calling interface from GP, parser codes}.
2061: \label{se:gp.interface}
2062: A \idx{parser code} is a character string describing all the GP parser
2063: needs to know about the function prototype. It contains a sequence of the
2064: following atoms:
2065:
2066: \settabs\+\indent&\kbd{Dxxx}\quad&\cr
2067: \noindent $\bullet$ Syntax requirements, used by functions like
2068: \kbd{for}, \kbd{sum}, etc.:
2069: %
2070: \+& \kbd{=} & separator \kbd{=} required at this point (between two
2071: arguments)\cr
2072:
2073: \noindent$\bullet$ Mandatory arguments, appearing in the same order as the
2074: input arguments they describe:
2075: %
2076: \+& \kbd{G} & \kbd{GEN}\cr
2077: \+& \kbd{\&}& \kbd{*GEN}\cr
2078: \+& \kbd{L} & long {\rm (we implicitly identify \kbd{int} with \kbd{long})}\cr
2079: \+& \kbd{S} & symbol (i.e.~GP identifier name). Function expects a
2080: \kbd{*entree}\cr
2081: \+& \kbd{V} & variable (as \kbd{S}, but rejects symbols associated to
2082: functions)\cr
2083: \+& \kbd{n} & variable, expects a \idx{variable number} (a \kbd{long}, not an
2084: \kbd{*entree})\cr
2085: \+& \kbd{I} & string containing a sequence of GP statements (a \var{seq}), %
2086: to be processed by \kbd{lisseq}\cr
2087: \+&&(useful for control statements)\cr
2088: \+& \kbd{E} & string containing a {\it single} GP statement (an %
2089: \var{expr}), to be processed by \kbd{lisexpr}\cr
2090: \+& \kbd{r} & raw input (treated as a string without quotes). Quoted %
2091: args are copied as strings\cr
2092: \+&&\quad Stops at first unquoted \kbd{')'} or \kbd{','}. Special chars can
2093: be quoted using '\kbd{\bs}'\cr
2094: \+&&\quad Example: \kbd{aa"b\bs n)"c} yields the string \kbd{"aab\bs{n})c"}\cr
2095: \+& \kbd{s} & expanded string. Example: \kbd{Pi"x"2} yields \kbd{"3.142x2"}\cr
2096: \+&& Unquoted components can be of any PARI type (converted following current
2097: output\cr
2098: \+&& format)\cr
2099:
2100: \noindent$\bullet$ Optional arguments:
2101: %
2102: \+& \kbd{s*} & any number of strings, possibly 0 (see \kbd{s})\cr
2103: \+& \kbd{D\var{xxx}} & argument has a default value\cr
2104:
1.2 ! noro 2105: The \kbd{s*} code is technical and you probably do not need it, but we give
! 2106: its description for completeness. It reads all remaining arguments in
! 2107: \tev{string context} (see \secref{se:strings}), and sends a
! 2108: (\kbd{NULL}-terminated) list of \kbd{GEN*} pointing to these. The automatic
! 2109: concatenation rules in string context are implemented so that adjacent strings
! 2110: are read as different arguments, as if they had been comma-separated. For
! 2111: instance, if the remaining argument sequence is: \kbd{"xx" 1, "yy"}, the
! 2112: \kbd{s*} atom will send a \kbd{GEN *g = \obr \&a, \&b, \&c, NULL\cbr}, where
! 2113: $a$, $b$, $c$ are \kbd{GEN}s of type \typ{STR} (content \kbd{xx}), \typ{INT}
! 2114: and \typ{STR} (content \kbd{yy}).
! 2115:
1.1 noro 2116: The format to indicate a default value (atom starts with a \kbd{D}) is
2117: ``\kbd{D\var{value},\var{type},}'', where \var{type} is the code for any
1.2 ! noro 2118: mandatory atom (previous group), \var{value} is any valid GP expression
1.1 noro 2119: which is converted according to \var{type}, and the ending comma is
2120: mandatory. For instance \kbd{D0,L,} stands for ``this optional argument will
2121: be converted to a \kbd{long}, and is \kbd{0} by default''. So if the
2122: user-given argument reads \kbd{1 + 3} at this point, \kbd{(long)4} is sent to
2123: the function (via \tet{itos}$()$); and \kbd{(long)0} if the argument is
1.2 ! noro 2124: omitted. The following special syntaxes are available:
1.1 noro 2125:
2126: \settabs\+\indent\indent&\kbd{Dxxx}\quad& optional \kbd{*GEN},&\cr
2127: \+&\kbd{DG}& optional \kbd{GEN}, & send \kbd{NULL} if argument omitted.\cr
2128:
2129: \+&\kbd{D\&}& optional \kbd{*GEN}, send \kbd{NULL} if argument omitted.\cr
2130:
2131: \+&\kbd{DV}& optional \kbd{*entree}, send \kbd{NULL} if argument omitted.\cr
2132:
2133: \+&\kbd{DI}& optional \kbd{*char}, send \kbd{NULL} if argument omitted.\cr
2134:
2135: \+&\kbd{Dn}& optional variable number, $-1$ if omitted.\cr
2136:
2137: \noindent$\bullet$ Automatic arguments:
2138: %
2139: \+& \kbd{f} & Fake \kbd{*long}. C function requires a pointer but we
2140: don't use the resulting \kbd{long}\cr
2141: \+& \kbd{p} & real precision (default \kbd{realprecision})\cr
2142: \+& \kbd{P} & series precision (default \kbd{seriesprecision},
2143: global variable \kbd{precdl} for the library)\cr
2144:
2145: \noindent $\bullet$ Return type: \kbd{GEN} by default, otherwise the
2146: following can appear at the start of the code string:
2147: %
2148: \+& \kbd{l} & return \kbd{long}\cr
2149: \+& \kbd{v} & return \kbd{void}\cr
2150:
2151: No more than 8 arguments can be given (syntax requirements and return types
2152: are not considered as arguments). This is currently hardcoded but can
2153: trivially be changed by modifying the definition of \kbd{argvec} in
2154: \kbd{anal.c:identifier()}. This limitation should disappear in future
2155: versions.
2156:
2157: When the function is called under GP, the prototype is scanned and each time
2158: an atom corresponding to a mandatory argument is met, a user-given argument
2159: is read (GP outputs an error message it the argument was missing). Each time
2160: an optional atom is met, a default value is inserted if the user omits the
2161: argument. The ``automatic'' atoms fill in the argument list transparently,
2162: supplying the current value of the corresponding variable (or a dummy
2163: pointer).
2164:
2165: For instance, here is how you would code the following prototypes (which
2166: don't involve default values):
2167: \bprog
2168: GEN name(GEN x, GEN y, long prec) ----> "GGp"
2169: void name(GEN x, GEN y, long prec) ----> "vGGp"
2170: void name(GEN x, long y, long prec) ----> "vGLp"
2171: long name(GEN x) ----> "lG"
2172: @eprog
2173:
2174: If you want more examples, GP gives you easy access to the parser codes
2175: associated to all GP functions: just type \kbd{\b{h} \var{function}}. You
2176: can then compare with the C prototypes as they stand in the code.
2177:
2178: \misctitle{Remark}: If you need to implement complicated control statements
2179: (probably for some improved summation functions), you'll need to know about
2180: the \teb{entree} type, which is not documented. Check the comment before
2181: the function list at the end of \kbd{language/init.c} and the source code
2182: in \kbd{language/sumiter.c}. You should be able to make something of it.
2183: \smallskip
2184:
2185: \subsec{Coding guidelines}.
2186: \noindent
2187: Code your function in a file of its own, using as a guide other functions
2188: in the PARI sources. One important thing to remember is to clean the stack
2189: before exiting your main function (usually using \kbd{gerepile}), since
2190: otherwise successive calls to the function will clutter the stack with
2191: unnecessary garbage, and stack overflow will occur sooner. Also, if it
2192: returns a \kbd{GEN} and you want it to be accessible to GP, you have to
2193: make sure this \kbd{GEN} is suitable for \kbd{gerepileupto} (see
2194: \secref{se:garbage}).
2195:
2196: If error messages are to be generated in your function, use the general
2197: error handling routine \kbd{err} (see \secref{se:err}). Recall that, apart
2198: from the \kbd{warn} variants, this function does not return but ends with
2199: a \kbd{longjmp} statement. As well, instead of explicit \kbd{printf}~/
2200: \kbd{fprintf} statements, use the following encapsulated variants:
2201:
2202: \fun{void}{pariputs}{char *s}: write \kbd{s} to the GP output stream.
2203:
2204: \fun{void}{fprintferr}{char *s}: write \kbd{s} to the GP error stream (this
2205: function is in fact much more versatile, see \secref{se:dbg_output}).
2206:
2207: Declare all public functions in an appropriate header file, if you
2208: expect somebody will access them from C. For example, if dynamic
2209: loading is not available, you may need to modify PARI to access these
2210: functions, so put them in \kbd{paridecl.h}. The other functions should
2211: be declared \kbd{static} in your file.
2212:
2213: Your function is now ready to be used in library mode after compilation and
2214: creation of the library. If possible, compile it as a shared library (see
2215: the \kbd{Makefile} coming with the \kbd{matexp} example in the
2216: distribution). It is however still inaccessible from GP.\smallskip
2217:
2218: \subsec{Integration with GP as a shared module}
2219:
2220: To tell GP about your function, you must do the following. First, find a
2221: name for it. It does not have to match the one used in library mode, but
2222: consistency is nice. It has to be a valid GP identifier, i.e.~use only
2223: alphabetic characters, digits and the underscore character (\kbd{\_}), the
2224: first character being alphabetic.
2225:
2226: Then you have to figure out the correct \idx{parser code} corresponding to
2227: the function prototype. This has been explained above
2228: (\secref{se:gp.interface}).
2229:
2230: Now, assuming your Operating System is supported by \tet{install}, simply
2231: write a GP script like the following:
2232:
2233: \bprog
2234: install(libname, code, gpname, library)
2235: addhelp(gpname, "some help text")
2236: @eprog
2237: \noindent(see \secref{se:addhelp} and~\ref{se:install}). The \idx{addhelp}
2238: part is not mandatory, but very useful if you want others to use your
2239: module. \kbd{libname} is how the function is named in the library,
2240: usually the same name as one visible from C.
2241:
2242: Read that file from your GP session (from your \idx{preferences file} for
2243: instance, see \secref{se:gprc}), and that's it, you can use the new
2244: function \var{gpname} under GP (and we would very much like to hear about
2245: it!).
2246:
2247: \subsec{Integration the hard way}
2248:
2249: If \tet{install} is not available for your Operating System, it's more
2250: complicated: you have to hardcode your function in the GP binary (or
2251: install \idx{Linux}). Here's what needs to be done:
2252:
2253: In the definition of \kbd{functions\_basic} (file \kbd{language/init.c}),
2254: add your entry in exact alphabetical order by its GP name (note that digits
2255: come before letters), in a line of the form:
2256: \bprog
2257: { "gpname", V, (void*)libname, secno, "code" }
2258: @eprog
2259:
2260: \noindent where
2261:
2262: \kbd{libname} is the name of your function in library mode,
2263:
2264: \kbd{gpname} the name that you have chosen to call it under GP,
2265:
2266: \kbd{secno} is the section number of Chapter~3 in which this function would
2267: belong (type \kbd{?} in GP to see the list),
2268:
2269: \kbd{V} is a number between 0 and 99. Right now, for PARI there are only two
2270: significant values: zero means that it's possible to call the function
2271: without argument, and non-zero means it needs at least one argument.
2272: A binding of PARI to an external language (such as \kbd{Math::Pari}
2273: Perl module) may actually distinguish between different non-zero values.
2274: Better use 99 if you want a non-zero value which will not confuse anybody.
2275:
2276: \kbd{code} is the parser code.
2277:
2278: Once this has been done, in the file \kbd{language/helpmessages.c} add in
2279: exact alphabetical order a short message describing the effect of your
2280: function:
2281: \kbd{"name(x,y,...)=short descriptive message",}
2282:
2283: The message must be a single line, of arbitrary length. Do not use
2284: \kbd{\bs{n}}; the necessary newlines will be inserted by GP's online help
2285: functions. Optional arguments should be shown between braces (see the other
2286: messages for comparison).\smallskip
2287:
2288: Now, you can recompile GP.
2289:
2290: \subsec{Example}.
2291: %
2292: A complete description could look like this:
2293:
2294: \bprog
2295: {
2296: install(bnfinit0, GD0,L,DGp, ClassGroupInit, "libpari.so");
2297: addhelp(ClassGroupInit, "ClassGroupInit(P,{flag=0},{data=[]}):
2298: compute the necessary data for ...");
2299: }
2300: @eprog
2301: which means we have a function \kbd{ClassGroupInit} under GP, which calls
2302: the library function \kbd{bnfinit0} . The function has one mandatory
2303: argument, and possibly two more (two \kbd{'D'} in the code), plus the
2304: current real precision. More precisely, the first argument is a \kbd{GEN},
2305: the second one is converted to a \kbd{long} using \kbd{itos} (\kbd{0} is
2306: passed if it is omitted), and the third one is also a \kbd{GEN}, but we
2307: pass \kbd{NULL} if no argument was supplied by the user. This matches the C
2308: prototype (from \kbd{paridecl.h}):
2309: %
2310: \bprog
2311: GEN bnfinit0(GEN P, long flag, GEN data, long prec)
2312: @eprog
2313:
2314: This function is in fact coded in \kbd{basemath/buch2.c}, and will in this
2315: case be completely identical to the GP function \kbd{bnfinit} but GP does
2316: not need to know about this, only that it can be found somewhere in the
2317: shared library \kbd{libpari.so}.
2318:
2319: \misctitle{Important note:} You see in this example that it is the
2320: function's responsibility to correctly interpret its operands: \kbd{data =
2321: NULL} is interpreted {\it by the function\/} as an empty vector. Note that
2322: since \kbd{NULL} is never a valid \kbd{GEN} pointer, this trick always
2323: enables you to distinguish between a default value and actual input: the
2324: user could explicitly supply an empty vector!
2325:
2326: \misctitle{Note:} If \kbd{install} were not available we would have to
2327: modify \kbd{language/helpmessages.c}, and \kbd{language/init.c} and
2328: recompile GP. The entry in \kbd{functions\_basic} corresponding to the
2329: function above is actually
2330: \bprog
2331: { "bnfinit", 91, (void*)bnfinit0, 6, "GD0,L,DGp" }
2332: @eprog
2333: \vfill\eject
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to usersch4.tex CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari-2.2 / doc