Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download

Testing latest pari + WASM + node.js... and it works?! Wow.

28495 views
License: GPL3
ubuntu2004
1
% Copyright (c) 2000 The PARI Group
2
%
3
% This file is part of the PARI/GP documentation
4
%
5
% Permission is granted to copy, distribute and/or modify this document
6
% under the terms of the GNU General Public License
7
\chapter{Programming PARI in Library Mode}
8
9
\noindent The \emph{User's Guide to Pari/GP} gives in three chapters a
10
general presentation of the system, of the \kbd{gp} calculator, and detailed
11
explanation of high level PARI routines available through the calculator. The
12
present manual assumes general familiarity with the contents of these
13
chapters and the basics of ANSI C programming, and focuses on the usage of
14
the PARI library. In this chapter, we introduce the general concepts of PARI
15
programming and describe useful general purpose functions; the following
16
chapters describes all public low or high-level functions, underlying or
17
extending the GP functions seen in Chapter 3 of the User's guide.
18
19
\section{Introduction: initializations, universal objects}
20
\label{se:intro4}
21
22
\noindent
23
To use PARI in \idx{library mode}, you must write a C program and link it to
24
the PARI library. See the installation guide or the Appendix to the
25
\emph{User's Guide to Pari/GP} on how to create and install the library and
26
include files. A sample Makefile is presented in Appendix~A, and a more
27
elaborate one in \kbd{examples/Makefile}. The best way to understand how
28
programming is done is to work through a complete example. We will write such
29
a program in~\secref{se:prog}. Before doing this, a few explanations are in
30
order.
31
32
First, one must explain to the outside world what kind of objects and
33
routines we are going to use. This is done\footnote{*}{This assumes that PARI
34
headers are installed in a directory which belongs to your compiler's search
35
path for header files. You might need to add flags like
36
\kbd{-I/usr/local/include} or modify \tet{C_INCLUDE_PATH}.}
37
with the directive
38
39
\bprog
40
#include <pari/pari.h>
41
@eprog
42
\noindent
43
In particular, this defines the fundamental type for all PARI objects: the
44
type \teb{GEN}, which is simply a pointer to \kbd{long}.
45
46
Before any PARI routine is called, one must initialize the system, and in
47
particular the PARI stack which is both a scratchboard and a repository for
48
computed objects. This is done with a call to the function
49
50
\fun{void}{pari_init}{size_t size, ulong maxprime}
51
52
\noindent The first argument is the number of bytes given to PARI to work
53
with, and the second is the upper limit on a precomputed prime number table;
54
\kbd{size} should not reasonably be taken below $500000$ but you may set
55
$\tet{maxprime} = 0$, although the system still needs to precompute all
56
primes up to about $2^{16}$. For lower-level variants allowing finer
57
control, e.g.~preventing PARI from installing its own error or signal
58
handlers, see~\secref{se:pari_init_tech}.
59
60
\noindent We have now at our disposal:
61
62
\item a PARI \tev{stack} containing nothing. This is a big
63
connected chunk of \kbd{size} bytes of memory, where all computations
64
take place. In large computations, intermediate results quickly
65
clutter up memory so some kind of garbage collecting is needed. Most
66
systems do garbage collecting when the memory is getting scarce, and this
67
slows down the performance. PARI takes a different approach, admittedly more
68
demanding on the programmer: you must do your own cleaning up when the
69
intermediate results are not needed anymore. We will see later how (and when)
70
this is done.
71
72
\item the following \emph{universal objects} (by definition, objects
73
which do not belong to the stack): the integers $0$, $1$, $-1$, $2$ and
74
$-2$ (respectively called \tet{gen_0}, \tet{gen_1}, \tet{gen_m1},
75
\tet{gen_2} and \tet{gen_m2}), the fraction $\dfrac{1}{2}$ (\tet{ghalf}).
76
All of these are of type \kbd{GEN}.
77
78
\item a \tev{heap} which is just a linked list of permanent
79
universal objects. For now, it contains exactly the ones listed above. You
80
will probably very rarely use the heap yourself; and if so, only as a
81
collection of copies of objects taken from the stack (called \idx{clone}s in
82
the sequel). Thus you need not bother with its internal structure, which may
83
change as PARI evolves. Some complex PARI functions create clones for special
84
garbage collecting purposes, usually destroying them when returning.
85
86
\item a table of primes (in fact of \emph{differences} between
87
consecutive primes), called \teb{diffptr}, of type \kbd{byteptr}
88
(pointer to \kbd{unsigned char}). Its use is described in
89
\secref{se:primetable} later. Using it directly is deprecated,
90
high-level iterators provide a cleaner and more flexible interface, see
91
\secref{se:primeiter} (such iterators use the private prime table, but extend
92
it dynamically).
93
94
\item access to all the built-in functions of the PARI library.
95
These are declared to the outside world when you include \kbd{pari.h}, but
96
need the above things to function properly. So if you forget the call to
97
\tet{pari_init}, you will get a fatal error when running your program.
98
99
\section{Important technical notes}
100
101
\subsec{Backward compatibility} The PARI function names evolved over time,
102
and deprecated functions are eventually deleted. The file \kbd{pariold.h}
103
contains macros implementing a weak form of backward compatibility.
104
In particular, whenever the name of a documented function changes, a
105
\kbd{\#define} is added to this file so that the old name expands to the new
106
one (provided the prototype didn't change also).
107
108
This file is included by \kbd{pari.h}, but a large section is commented out
109
by default. Define \tet{PARI_OLD_NAMES} before including \kbd{pari.h} to
110
pollute your namespace with lots of obsolete names like
111
\kbd{un}\footnote{*}{For \kbd{(long)gen\_1}. Since 2004 and version 2.2.9,
112
typecasts are completely unnecessary in PARI programs.}: that might enable
113
you to compile old programs without having to modify them. The preferred way
114
to do that is to add \kbd{-DPARI\_OLD\_NAMES} to your compiler \kbd{CFLAGS},
115
so that you don't need to modify the program files themselves.
116
117
Of course, it's better to fix the program if you can!
118
119
\subsec{Types}
120
121
\noindent
122
Although PARI objects all have the C type \kbd{GEN}, we will freely use
123
the word \teb{type} to refer to PARI dynamic subtypes: \typ{INT}, \typ{REAL},
124
etc. The declaration
125
\bprog
126
GEN x;
127
@eprog\noindent
128
declares a C variable of type \kbd{GEN}, but its ``value'' will be said to
129
have type \typ{INT}, \typ{REAL}, etc. The meaning should always be clear from
130
the context.
131
132
\subsec{Type recursivity}
133
134
\noindent
135
Conceptually, most PARI types are recursive. But the \kbd{GEN} type is a
136
pointer to \kbd{long}, not to \kbd{GEN}. So special macros must be used to
137
access \kbd{GEN}'s components. The simplest one is \tet{gel}$(V, i)$, where
138
\key{el} stands for \key{el}ement, to access component number $i$ of the
139
\kbd{GEN} $V$. This is a valid \kbd{lvalue} (may be put on the left side of
140
an assignment), and the following two constructions are exceedingly frequent
141
%
142
\bprog
143
gel(V, i) = x;
144
x = gel(V, i);
145
@eprog\noindent
146
where \kbd{x} and \kbd{V} are \kbd{GEN}s. This macro accesses and modifies
147
directly the components of $V$ and do not create a copy of the coefficient,
148
contrary to all the library \emph{functions}.
149
150
More generally, to retrieve the values of elements of lists of \dots\ of
151
lists of vectors we have the \tet{gmael} macros (for {\bf m}ultidimensional
152
{\bf a}rray {\bf el}ement). The syntax is $\key{gmael}n(V,a_1,\dots,a_n)$,
153
where $V$ is a \kbd{GEN}, the $a_i$ are indexes, and $n$ is an integer
154
between $1$ and $5$. This stands for $x[a_1][a_2]\dots[a_n]$, and returns a
155
\kbd{GEN}. The macros \tet{gel} (resp.~\tet{gmael}) are synonyms for
156
\tet{gmael1} (resp.~\kbd{gmael2}).
157
158
Finally, the macro $\tet{gcoeff}(M, i, j)$ has exactly the meaning of
159
\kbd{M[i,j]} in GP when \kbd{M} is a matrix. Note that due to the
160
implementation of \typ{MAT}s as horizontal lists of vertical vectors,
161
\kbd{gcoeff(x,y)} is actually equivalent to \kbd{gmael(y,x)}. One should use
162
\kbd{gcoeff} in matrix context, and \kbd{gmael} otherwise.
163
164
\subsec{Variations on basic functions}\label{se:low_level} In the library
165
syntax descriptions in Chapter~3, we have only given the basic names of the
166
functions. For example \kbd{gadd}$(x,y)$ assumes that $x$ and $y$ are
167
\kbd{GEN}s, and \emph{creates} the result $x+y$ on the PARI stack. For most
168
of the basic operators and functions, many other variants are available. We
169
give some examples for \kbd{gadd}, but the same is true for all the basic
170
operators, as well as for some simple common functions (a complete list
171
is given in Chapter~6):
172
173
\fun{GEN}{gaddgs}{GEN x, long y}
174
175
\fun{GEN}{gaddsg}{long x, GEN y}
176
177
\noindent In the following one, \kbd{z} is a preexisting \kbd{GEN} and the
178
result of the corresponding operation is put into~\kbd{z}. The size of the PARI
179
stack does not change:
180
181
\fun{void}{gaddz}{GEN x, GEN y, GEN z}
182
183
\noindent (This last form is inefficient in general and deprecated outside of
184
PARI kernel programming.) Low level kernel functions implement these
185
operators for specialized arguments and are also available: Level 0 deals
186
with operations at the word level (\kbd{long}s and \kbd{ulong}s), Level 1
187
with \typ{INT} and \typ{REAL} and Level 2 with the rest (modular arithmetic,
188
polynomial arithmetic and linear algebra). Here are some examples of Level
189
$1$ functions:
190
191
\fun{GEN}{addii}{GEN x, GEN y}: here $x$ and $y$ are \kbd{GEN}s of type
192
\typ{INT} (this is not checked).
193
194
\fun{GEN}{addrr}{GEN x, GEN y}: here $x$ and $y$ are \kbd{GEN}s of
195
type \typ{REAL} (this is not checked).
196
197
\noindent
198
There also exist functions \teb{addir}, \teb{addri}, \teb{mpadd} (whose
199
two arguments can be of type \typ{INT} or \typ{REAL}), \teb{addis} (to add a
200
\typ{INT} and a \kbd{long}) and so on.
201
202
The Level $1$ names are self-explanatory once you know that {\bf i} stands for a
203
\typ{INT}, {\bf r} for a \typ{REAL}, {\bf mp} for i or r, {\bf s} for a signed C
204
long integer, {\bf u} for an unsigned C long integer; finally the suffix {\bf z}
205
means that the result is not created on the PARI stack but assigned to a
206
preexisting GEN object passed as an extra argument. Chapter 6 gives a
207
description of these low-level functions.
208
209
Level $2$ names are more complicated, see \secref{se:level2names} for all the
210
gory details, and we content ourselves with a simple example used to implement
211
\typ{INTMOD} arithmetic:
212
213
\fun{GEN}{Fp_add}{GEN x, GEN y, GEN m}: returns the sum of $x$ and $y$ modulo
214
$m$. Here $x, y, m$ are \typ{INT}s (this is not checked). The operation is
215
more efficient if the inputs $x$, $y$ are reduced modulo $m$, but this is not
216
a necessary condition.
217
218
\misctitle{Important Note} These specialized functions are of course more
219
efficient than the generic ones, but note the hidden danger here: the types
220
of the objects involved (which is not checked) must be severely controlled,
221
e.g.~using \kbd{addii} on a \typ{FRAC} argument will cause disasters. Type
222
mismatches may corrupt the PARI stack, though in most cases they will just
223
immediately overflow the stack. Because of this, the PARI philosophy of
224
giving a result which is as exact as possible, enforced for generic functions
225
like \kbd{gadd} or \kbd{gmul}, is dropped in kernel routines of Level $1$,
226
where it is replaced by the much simpler rule: the result is a \typ{INT} if
227
and only if all arguments are integer types (\typ{INT} but also C \kbd{long}
228
and \kbd{ulong}) and a \typ{REAL} otherwise. For instance, multiplying a
229
\typ{REAL} by a \typ{INT} always yields a \typ{REAL} if you use \kbd{mulir},
230
where \kbd{gmul} returns the \typ{INT} \kbd{gen\_0} if the integer is $0$.
231
232
\subsec{Portability: 32-bit / 64-bit architectures}
233
234
\noindent
235
PARI supports both 32-bit and 64-bit based machines, but not simultaneously!
236
The library is compiled assuming a given architecture, and some
237
of the header files you include (through \kbd{pari.h}) will have been
238
modified to match the library.
239
240
Portable macros are defined to bypass most machine dependencies. If you want
241
your programs to run identically on 32-bit and 64-bit machines, you have to
242
use these, and not the corresponding numeric values, whenever the precise
243
size of your \kbd{long} integers might matter. Here are the most important
244
ones:
245
246
\settabs\+ -----------------------------&---------------&------------&\cr \+
247
& 64-bit & 32-bit \cr\+ \tet{BITS_IN_LONG} & 64 & 32 \cr\+
248
\tet{LONG_IS_64BIT} & defined & undefined \cr\+
249
\tet{DEFAULTPREC} & 3 & 4 & ($\approx$ 19 decimal digits, %
250
see formula below) \cr\+
251
\tet{MEDDEFAULTPREC}& 4 & 6 & ($\approx$ 38 decimal digits) \cr\+
252
\tet{BIGDEFAULTPREC}& 5 & 8 & ($\approx$ 57 decimal digits) \cr
253
\noindent For instance, suppose you call a transcendental function, such as
254
255
\kbd{GEN \key{gexp}(GEN x, long prec)}.
256
257
\noindent The last argument \kbd{prec} is an integer $\geq 3$, corresponding
258
to the default floating point precision required. It is \emph{only} used if
259
\kbd{x} is an exact object, otherwise the relative precision is determined by
260
the precision of~\kbd{x}. Since the parameter \kbd{prec} sets the size of the
261
inexact result counted in (\kbd{long}) \emph{words} (including codewords),
262
the same value of \kbd{prec} will yield different results on 32-bit and
263
64-bit machines. Real numbers have two codewords (see~\secref{se:impl}), so
264
the formula for computing the bit accuracy is
265
$$ \tet{bit_accuracy}(\kbd{prec}) = (\kbd{prec} - 2) * \tet{BITS_IN_LONG}$$
266
(this is actually the definition of an inline function). The corresponding
267
accuracy expressed in decimal digits would be
268
%
269
$$ \kbd{bit\_accuracy(prec)} * \log(2) / \log(10).$$
270
%
271
For example if the value of \kbd{prec} is 5, the corresponding accuracy for
272
32-bit machines is $(5-2)*\log(2^{32})/\log(10)\approx 28$ decimal digits,
273
while for 64-bit machines it is $(5-2)*\log(2^{64})/\log(10)\approx 57$
274
decimal digits.
275
276
Thus, you must take care to change the \kbd{prec} parameter you are supplying
277
according to the bit size, either using the default precisions given by the
278
various \kbd{DEFAULTPREC}s, or by using conditional constructs of the form:
279
%
280
\bprog
281
#ifndef LONG_IS_64BIT
282
prec = 4;
283
#else
284
prec = 6;
285
#endif
286
@eprog
287
\noindent which is in this case equivalent to the statement
288
\kbd{prec = MEDDEFAULTPREC;}.
289
290
Note that for parity reasons, half the accuracies available on 32-bit
291
architectures (the odd ones) have no precise equivalents on 64-bit machines.
292
293
\subsec{Using \kbd{malloc} / \kbd{free}}
294
You should make use of the PARI stack as much as possible, and avoid
295
allocating objects using the customary functions. If you do, you should
296
use, or at least have a very close look at, the following wrappers:
297
298
\fun{void*}{pari_malloc}{size_t size} calls \kbd{malloc} to allocate
299
\kbd{size} bytes and returns a pointer to the allocated memory. If the
300
request fails, an error is raised. The \kbd{SIGINT} signal is blocked until
301
\kbd{malloc} returns, to avoid leaving the system stack in an inconsistent
302
state.
303
304
\fun{void*}{pari_realloc}{void* ptr, size_t size} as \tet{pari_malloc} but
305
calls \kbd{realloc} instead of \kbd{malloc}.
306
307
\fun{void}{pari_realloc_ip}{void** ptr, size_t size}
308
equivalent to \kbd{*ptr= realloc(*ptr, size)}, while blocking \kbd{SIGINT}.
309
310
\fun{void*}{pari_calloc}{size_t size} as \tet{pari_malloc}, setting the
311
memory to zero.
312
313
\fun{void}{pari_free}{void* ptr} calls \kbd{free} to liberate the memory
314
space pointed to by \kbd{ptr}, which must have been allocated by \kbd{malloc}
315
(\tet{pari_malloc}) or \kbd{realloc} (\tet{pari_realloc}). The \kbd{SIGINT}
316
signal is blocked until \kbd{free} returns.
317
318
If you use the standard \kbd{libc} functions instead of our wrappers, then
319
your functions will be subtly incompatible with the \kbd{gp} calculator: when
320
the user tries to interrupt a computation, the calculator may crash
321
(if a system call is interrupted at the wrong time).
322
323
\section{Garbage collection}\label{se:garbage}\sidx{garbage collecting}
324
325
\subsec{Why and how}
326
327
\noindent
328
As we have seen, \kbd{pari\_init} allocates a big range of
329
addresses, the \tev{stack}, that are going to be used throughout. Recall
330
that all PARI objects are pointers. Except for a few universal objects,
331
they all point at some part of the stack.
332
333
The stack starts at the address \kbd{bot} and ends just before \kbd{top}. This
334
means that the quantity
335
%
336
$$ (\kbd{top} - \kbd{bot})\,/\,\kbd{sizeof(long)} $$
337
%
338
is (roughly) equal to the \kbd{size} argument of \kbd{pari\_init}. The PARI
339
stack also has a ``current stack pointer'' called \teb{avma}, which stands
340
for {\bf av}ailable {\bf m}emory {\bf a}ddress. These three variables are
341
global (declared by \kbd{pari.h}). They are of type \tet{pari_sp}, which
342
means \emph{pari stack pointer}.
343
344
The stack is oriented upside-down: the more recent an object, the closer to
345
\kbd{bot}. Accordingly, initially \kbd{avma} = \kbd{top}, and \kbd{avma} gets
346
\emph{decremented} as new objects are created. As its name indicates,
347
\kbd{avma} always points just \emph{after} the first free address on the
348
stack, and \kbd{(GEN)avma} is always (a pointer to) the latest created object.
349
When \kbd{avma} reaches \kbd{bot}, the stack overflows, aborting all
350
computations, and an error message is issued. To avoid this \emph{you}
351
need to clean up the stack from time to time, when intermediate objects are
352
not needed anymore. This is called ``\emph{garbage collecting}.''
353
354
We are now going to describe briefly how this is done. We will see many
355
concrete examples in the next subsection.
356
357
\noindent\item
358
First, PARI routines do their own garbage collecting, which means that
359
whenever a documented function from the library returns, only its result(s)
360
have been added to the stack, possibly up to a very small overhead
361
(undocumented ones may not do this). In
362
particular, a PARI function that does not return a \kbd{GEN} does not clutter
363
the stack. Thus, if your computation is small enough (e.g.~you call few PARI
364
routines, or most of them return \kbd{long} integers), then you do not need
365
to do any garbage collecting. This is probably the case in many of your
366
subroutines. Of course the objects that were on the stack \emph{before} the
367
function call are left alone. Except for the ones listed below, PARI
368
functions only collect their own garbage.
369
370
\noindent\item
371
It may happen that all objects that were created after a certain point can
372
be deleted~--- for instance, if the final result you need is not a
373
\kbd{GEN}, or if some search proved futile. Then, it is enough to record
374
the value of \kbd{avma} just \emph{before} the first garbage is created,
375
and restore it upon exit:
376
377
\bprog
378
pari_sp av = avma; /*@Ccom record initial avma */
379
380
garbage ...
381
avma = av; /*@Ccom restore it */
382
@eprog
383
\noindent All objects created in the \kbd{garbage} zone will eventually
384
be overwritten: they should no longer be accessed after \kbd{avma} has been
385
restored.
386
387
\noindent\item
388
If you want to destroy (i.e.~give back the memory occupied by) the
389
\emph{latest} PARI object on the stack (e.g.~the latest one obtained from a
390
function call), you can use the function\sidx{destruction}%
391
\vadjust{\penalty500}%discourage page break
392
393
\fun{void}{cgiv}{GEN z}
394
395
\noindent where \kbd{z} is the object you want to give back. This is
396
equivalent to the above where the initial \kbd{av} is computed from \kbd{z}.
397
398
\noindent\item
399
Unfortunately life is not so simple, and sometimes you will want
400
to give back accumulated garbage \emph{during} a computation without losing
401
recent data. We shall start with the lowest level function to get a feel for
402
the underlying mechanisms, we shall describe simpler variants later:
403
404
\fun{GEN}{gerepile}{pari_sp ltop, pari_sp lbot, GEN q}. This function cleans
405
up the stack between \kbd{ltop} and \kbd{lbot}, where $\kbd{lbot} <
406
\kbd{ltop}$, and returns the updated object \kbd{q}. This means:
407
408
1) we translate (copy) all the objects in the interval
409
$[\kbd{avma}, \kbd{lbot}[$, so that its right extremity abuts the address
410
\kbd{ltop}. Graphically
411
412
\vbox{\bprog
413
bot avma lbot ltop top
414
End of stack |-------------[++++++[-/-/-/-/-/-/-|++++++++| Start
415
free memory garbage
416
@eprog
417
\noindent becomes:
418
\bprog
419
bot avma ltop top
420
End of stack |---------------------------[++++++[++++++++| Start
421
free memory
422
@eprog
423
}
424
\noindent where \kbd{++} denote significant objects, \kbd{--} the unused part
425
of the stack, and \kbd{-/-} the garbage we remove.
426
427
2) The function then inspects all the PARI objects between \kbd{avma} and
428
\kbd{lbot} (i.e.~the ones that we want to keep and that have been translated)
429
and looks at every component of such an object which is not a codeword. Each
430
such component is a pointer to an object whose address is either
431
432
--- between \kbd{avma} and \kbd{lbot}, in which case it is suitably updated,
433
434
--- larger than or equal to \kbd{ltop}, in which case it does not change, or
435
436
--- between \kbd{lbot} and \kbd{ltop} in which case \kbd{gerepile}
437
raises an error (``significant pointers lost in gerepile'').
438
439
3) \key{avma} is updated (we add $\kbd{ltop} - \kbd{lbot}$ to the old value).
440
441
4) We return the (possibly updated) object \kbd{q}: if \kbd{q} initially
442
pointed between \kbd{avma} and \kbd{lbot}, we return the updated address, as
443
in~2). If not, the original address is still valid, and is returned!
444
445
As stated above, no component of the remaining objects (in particular
446
\kbd{q}) should belong to the erased segment [\kbd{lbot}, \kbd{ltop}[, and
447
this is checked within \kbd{gerepile}. But beware as well that the addresses
448
of the objects in the translated zone change after a call to \kbd{gerepile},
449
so you must not access any pointer which previously pointed into the zone
450
below \kbd{ltop}. If you need to recover more than one object, use the
451
\kbd{gerepileall} function below.
452
453
\misctitle{Remark}
454
As a consequence of the preceding explanation, if a PARI object is to be
455
relocated by \hbox{gerepile} then, apart from universal objects, the chunks
456
of memory used by its components should be in consecutive memory locations.
457
All \kbd{GEN}s created by documented PARI functions are guaranteed to satisfy
458
this. This is because the \kbd{gerepile} function knows only about \emph{two
459
connected zones}: the garbage that is erased (between \kbd{lbot} and
460
\kbd{ltop}) and the significant pointers that are copied and updated. If
461
there is garbage interspersed with your objects, disaster occurs when we try
462
to update them and consider the corresponding ``pointers''. In most cases of
463
course the said garbage is in fact a bunch of other \kbd{GEN}s, in which case
464
we simply waste time copying and updating them for nothing. But be wary when
465
you allow objects to become disconnected.
466
467
\noindent In practice this is achieved by the following programming idiom:
468
\bprog
469
ltop = avma; garbage(); lbot = avma; q = anything();
470
return gerepile(ltop, lbot, q); /*@Ccom returns the updated q */
471
@eprog\noindent or directly
472
\bprog
473
ltop = avma; garbage(); lbot = avma;
474
return gerepile(ltop, lbot, anything());
475
@eprog\noindent
476
Beware that
477
\bprog
478
ltop = avma; garbage();
479
return gerepile(ltop, avma, anything())
480
@eprog
481
482
\noindent might work, but should be frowned upon. We cannot predict whether
483
\kbd{avma} is evaluated after or before the call to \kbd{anything()}: it
484
depends on the compiler. If we are out of luck, it is \emph{after} the
485
call, so the result belongs to the garbage zone and the \kbd{gerepile}
486
statement becomes equivalent to \kbd{avma = ltop}. Thus we return a
487
pointer to random garbage.
488
489
\subsec{Variants}
490
491
\fun{GEN}{gerepileupto}{pari_sp ltop, GEN q}. Cleans the stack between
492
\kbd{ltop} and the \emph{connected} object \kbd{q} and returns \kbd{q}
493
updated. For this to work, \kbd{q} must have been created \emph{before} all
494
its components, otherwise they would belong to the garbage zone! Unless
495
mentioned otherwise, documented PARI functions guarantee this.
496
497
\fun{GEN}{gerepilecopy}{pari_sp ltop, GEN x}. Functionally equivalent to,
498
but more efficient than
499
\bprog
500
gerepileupto(ltop, gcopy(x))
501
@eprog\noindent In this case, the \kbd{GEN} parameter \kbd{x} need not
502
satisfy any property before the garbage collection: it may be disconnected,
503
components created before the root, and so on. Of course, this is about
504
twice slower than either \tet{gerepileupto} or \tet{gerepile}, because
505
\kbd{x} has to be copied to a clean stack zone first. This function is a
506
special case of \tet{gerepileall} below, where $n=1$.
507
508
\fun{void}{gerepileall}{pari_sp ltop, int n, ...}.
509
To cope with complicated cases where many objects have to be preserved. The
510
routine expects $n$ further arguments, which are the \emph{addresses} of
511
the \kbd{GEN}s you want to preserve:
512
\bprog
513
pari_sp ltop = avma;
514
...; y = ...; ... x = ...; ...;
515
gerepileall(ltop, 2, &x, &y);
516
@eprog\noindent
517
It cleans up the most recent part of the
518
stack (between \kbd{ltop} and \kbd{avma}), updating all the \kbd{GEN}s added
519
to the argument list. A copy is done just before the cleaning to preserve
520
them, so they do not need to be connected before the call. With
521
\kbd{gerepilecopy}, this is the most robust of the \kbd{gerepile} functions
522
(the less prone to user error), hence the slowest.
523
524
\fun{void}{gerepileallsp}{pari_sp ltop, pari_sp lbot, int n, ...}.
525
More efficient, but trickier than \kbd{gerepileall}. Cleans the stack between
526
\kbd{lbot} and \kbd{ltop} and updates the \kbd{GEN}s pointed at by the
527
elements of \kbd{gptr} without any further copying. This is subject to the
528
same restrictions as \kbd{gerepile}, the only difference being that more than
529
one address gets updated.
530
531
\subsec{Examples}
532
533
\subsubsec{gerepile}
534
535
Let \kbd{x} and \kbd{y} be two preexisting PARI objects and suppose that we
536
want to compute $\kbd{x}^2 + \kbd{y}^2$. This is done using the following
537
program:
538
\bprog
539
GEN x2 = gsqr(x);
540
GEN y2 = gsqr(y), z = gadd(x2,y2);
541
@eprog\noindent
542
The \kbd{GEN} \kbd{z} indeed points at the desired quantity. However,
543
consider the stack: it contains as unnecessary garbage \kbd{x2} and \kbd{y2}.
544
More precisely it contains (in this order) \kbd{z}, \kbd{y2}, \kbd{x2}.
545
(Recall that, since the stack grows downward from the top, the most recent
546
object comes first.)
547
548
It is not possible to get rid of \kbd{x2}, \kbd{y2} before \kbd{z} is
549
computed, since they are used in the final operation. We cannot record
550
\kbd{avma} before \kbd{x2} is computed and restore it later, since this would
551
destroy \kbd{z} as well. It is not possible either to use the function
552
\kbd{cgiv} since \kbd{x2} and \kbd{y2} are not at the bottom of the stack and
553
we do not want to give back~\kbd{z}.
554
555
But using \kbd{gerepile}, we can give back the memory locations corresponding
556
to \kbd{x2}, \kbd{y2}, and move the object \kbd{z} upwards so that no
557
space is lost. Specifically:
558
\bprog
559
pari_sp ltop = avma; /*@Ccom remember the current top of the stack */
560
GEN x2 = gsqr(x);
561
GEN y2 = gsqr(y);
562
pari_sp lbot = avma; /*@Ccom the bottom of the garbage pile */
563
GEN z = gadd(x2, y2); /*@Ccom z is now the last object on the stack */
564
z = gerepile(ltop, lbot, z);
565
@eprog
566
\noindent Of course, the last two instructions could also have been
567
written more simply:
568
\bprog
569
z = gerepile(ltop, lbot, gadd(x2,y2));
570
@eprog\noindent In fact \kbd{gerepileupto} is even simpler to use, because
571
the result of \kbd{gadd} is the last object on the stack and \kbd{gadd}
572
is guaranteed to return an object suitable for \kbd{gerepileupto}:
573
\bprog
574
ltop = avma;
575
z = gerepileupto(ltop, gadd(gsqr(x), gsqr(y)));
576
@eprog\noindent
577
Make sure you understand exactly what has happened before you go on!
578
579
\misctitle{Remark on assignments and gerepile} When the tree structure and
580
the size of the PARI objects which will appear in a computation are under
581
control, one may allocate sufficiently large objects at the beginning,
582
use assignment statements, then simply restore \kbd{avma}. Coming back to the
583
above example, note that \emph{if} we know that x and y are of type real
584
fitting into \kbd{DEFAULTPREC} words, we can program without using
585
\kbd{gerepile} at all:
586
\bprog
587
z = cgetr(DEFAULTPREC); ltop = avma;
588
gaffect(gadd(gsqr(x), gsqr(y)), z);
589
avma = ltop;
590
@eprog\noindent This is often \emph{slower} than a craftily used
591
\kbd{gerepile} though, and certainly more cumbersome to use. As a rule,
592
assignment statements should generally be avoided.
593
594
\misctitle{Variations on a theme} it is often necessary to do several
595
\kbd{gerepile}s during a computation. However, the fewer the better. The only
596
condition for \kbd{gerepile} to work is that the garbage be connected. If the
597
computation can be arranged so that there is a minimal number of connected
598
pieces of garbage, then it should be done that way.
599
600
For example suppose we want to write a function of two \kbd{GEN} variables
601
\kbd{x} and \kbd{y} which creates the vector $\kbd{[x}^2+\kbd{y},
602
\kbd{y}^2+\kbd{x]}$. Without garbage collecting, one would write:
603
%
604
\bprog
605
p1 = gsqr(x); p2 = gadd(p1, y);
606
p3 = gsqr(y); p4 = gadd(p3, x);
607
z = mkvec2(p2, p4); /* not suitable for gerepileupto! */
608
@eprog\noindent
609
This leaves a dirty stack containing (in this order) \kbd{z}, \kbd{p4},
610
\kbd{p3}, \kbd{p2}, \kbd{p1}. The garbage here consists of \kbd{p1} and
611
\kbd{p3}, which are separated by \kbd{p2}. But if we compute \kbd{p3}
612
\emph{before} \kbd{p2} then the garbage becomes connected, and we get the
613
following program with garbage collecting:
614
%
615
\bprog
616
ltop = avma; p1 = gsqr(x); p3 = gsqr(y);
617
lbot = avma; z = cgetg(3, t_VEC);
618
gel(z, 1) = gadd(p1,y);
619
gel(z, 2) = gadd(p3,x); z = gerepile(ltop,lbot,z);
620
@eprog\noindent Finishing by \kbd{z = gerepileupto(ltop, z)} would be ok as
621
well. Beware that
622
\bprog
623
ltop = avma; p1 = gadd(gsqr(x), y); p3 = gadd(gsqr(y), x);
624
z = cgetg(3, t_VEC);
625
gel(z, 1) = p1;
626
gel(z, 2) = p3; z = gerepileupto(ltop,z); /*@Ccom WRONG */
627
@eprog\noindent
628
is a disaster since \kbd{p1} and \kbd{p3} are created before
629
\kbd{z}, so the call to \kbd{gerepileupto} overwrites them, leaving
630
\kbd{gel(z, 1)} and \kbd{gel(z, 2)} pointing at random data! The following
631
does work:
632
\bprog
633
ltop = avma; p1 = gsqr(x); p3 = gsqr(y);
634
lbot = avma; z = mkvec2(gadd(p1,y), gadd(p3,x));
635
z = gerepile(ltop,lbot,z);
636
@eprog\noindent but is very subtly wrong in the sense that
637
\kbd{z = gerepileupto(ltop, z)} would \emph{not} work. The reason being
638
that \kbd{mkvec2} creates the root \kbd{z} of the vector \emph{after}
639
its arguments have been evaluated, creating the components of \kbd{z}
640
too early; \kbd{gerepile} does not care, but the created \kbd{z} is a time
641
bomb which will explode on any later \kbd{gerepileupto}.
642
On the other hand
643
\bprog
644
ltop = avma; z = cgetg(3, t_VEC);
645
gel(z, 1) = gadd(gsqr(x), y);
646
gel(z, 2) = gadd(gsqr(y), x); z = gerepileupto(ltop,z); /*@Ccom INEFFICIENT */
647
@eprog\noindent
648
leaves the results of \kbd{gsqr(x)} and \kbd{gsqr(y)} on the stack (and
649
lets \kbd{gerepileupto} update them for naught). Finally, the most elegant
650
and efficient version (with respect to time and memory use) is as follows
651
\bprog
652
z = cgetg(3, t_VEC);
653
ltop = avma; gel(z, 1) = gerepileupto(ltop, gadd(gsqr(x), y));
654
ltop = avma; gel(z, 2) = gerepileupto(ltop, gadd(gsqr(y), x));
655
@eprog\noindent
656
which avoids updating the container \kbd{z} and cleans up its components
657
individually, as soon as they are computed.
658
659
\misctitle{One last example} Let us compute the product of two complex
660
numbers $x$ and $y$, using the $3M$ method which requires 3 multiplications
661
instead of the obvious 4. Let $z = x*y$, and set $x = x_r + i*x_i$ and
662
similarly for $y$ and $z$. We compute $p_1 = x_r*y_r$, $p_2=x_i*y_i$,
663
$p_3=(x_r+x_i)*(y_r+y_i)$, and then we have $z_r=p_1-p_2$,
664
$z_i=p_3-(p_1+p_2)$. The program is as follows:
665
%
666
\bprog
667
ltop = avma;
668
p1 = gmul(gel(x,1), gel(y,1));
669
p2 = gmul(gel(x,2), gel(y,2));
670
p3 = gmul(gadd(gel(x,1), gel(x,2)), gadd(gel(y,1), gel(y,2)));
671
p4 = gadd(p1,p2);
672
lbot = avma; z = cgetg(3, t_COMPLEX);
673
gel(z, 1) = gsub(p1,p2);
674
gel(z, 2) = gsub(p3,p4); z = gerepile(ltop,lbot,z);
675
@eprog
676
677
\misctitle{Exercise} Write a function which multiplies a matrix by a column
678
vector. Hint: start with a \kbd{cgetg} of the result, and use \kbd{gerepile}
679
whenever a coefficient of the result vector is computed. You can look at the
680
answer in \kbd{src/basemath/RgV.c:RgM\_RgC\_mul()}.
681
682
\subsubsec{gerepileall}
683
684
Let us now see why we may need the \kbd{gerepileall} variants. Although it
685
is not an infrequent occurrence, we do not give a specific example but a
686
general one: suppose that we want to do a computation (usually inside a
687
larger function) producing more than one PARI object as a result, say two for
688
instance. Then even if we set up the work properly, before cleaning up we
689
have a stack which has the desired results \kbd{z1}, \kbd{z2} (say), and
690
then connected garbage from lbot to ltop. If we write
691
\bprog
692
z1 = gerepile(ltop, lbot, z1);
693
@eprog\noindent
694
then the stack is cleaned, the pointers fixed up, but we have lost the
695
address of \kbd{z2}. This is where we need the \idx{gerepileall}
696
function:
697
\bprog
698
gerepileall(ltop, 2, &z1, &z2)
699
@eprog
700
\noindent copies \kbd{z1} and \kbd{z2} to new locations, cleans the stack
701
from \kbd{ltop} to the old \kbd{avma}, and updates the pointers \kbd{z1} and
702
\kbd{z2}. Here we do not assume anything about the stack: the garbage can be
703
disconnected and \kbd{z1}, \kbd{z2} need not be at the bottom of the stack.
704
If all of these assumptions are in fact satisfied, then we can call
705
\kbd{gerepilemanysp} instead, which is usually faster since we do not need
706
the initial copy (on the other hand, it is less cache friendly).
707
708
A most important usage is ``random'' garbage collection during loops
709
whose size requirements we cannot (or do not bother to) control in advance:
710
\bprog
711
pari_sp av = avma;
712
GEN x, y;
713
while (...)
714
{
715
garbage(); x = anything();
716
garbage(); y = anything(); garbage();
717
if (gc_needed(av,1)) /*@Ccom memory is running low (half spent since entry) */
718
gerepileall(av, 2, &x, &y);
719
}
720
@eprog
721
\noindent Here we assume that only \kbd{x} and \kbd{y} are needed from one
722
iteration to the next. As it would be costly to call gerepile once for each
723
iteration, we only do it when it seems to have become necessary.
724
725
More precisely, the macro \tet{stack_lim}\kbd{(av,$n$)} denotes an address
726
where $2^{n-1} / (2^{n-1}+1)$ of the remaining stack space since reference
727
point \kbd{av} is exhausted ($1/2$ for $n=1$, $2/3$ for $n=2$). The test
728
\tet{gc_needed}\kbd{(av,$n$)} becomes true whenever \kbd{avma} drops below
729
that address.
730
731
\subsec{Comments}
732
733
First, \kbd{gerepile} has turned out to be a flexible and fast garbage
734
collector for number-theoretic computations, which compares favorably with
735
more sophisticated methods used in other systems. Our benchmarks indicate
736
that the price paid for using \kbd{gerepile} and \kbd{gerepile}-related
737
copies, when properly used, is usually less than 1\% of the total
738
running time, which is quite acceptable!
739
740
Second, it is of course harder on the programmer, and quite error-prone
741
if you do not stick to a consistent PARI programming style. If all seems
742
lost, just use \tet{gerepilecopy} (or \tet{gerepileall}) to fix up the stack
743
for you. You can always optimize later when you have sorted out exactly which
744
routines are crucial and what objects need to be preserved and their usual
745
sizes.
746
747
\smallskip If you followed us this far, congratulations, and rejoice: the
748
rest is much easier.
749
750
\section{Creation of PARI objects, assignments, conversions}
751
752
\subsec{Creation of PARI objects}\sidx{creation}
753
The basic function which creates a PARI object is
754
755
\fun{GEN}{cgetg}{long l, long t}
756
$l$ specifies the number of longwords to be allocated to the
757
object, and $t$ is the type of the object, in symbolic
758
form (see \secref{se:impl} for the list of these). The precise effect of
759
this function is as follows: it first creates on the PARI \emph{stack} a
760
chunk of memory of size \kbd{length} longwords, and saves the address of the
761
chunk which it will in the end return. If the stack has been used up, a
762
message to the effect that ``the PARI stack overflows'' is printed,
763
and an error raised. Otherwise, it sets the type and length of the PARI object.
764
In effect, it fills its first codeword (\kbd{z[0]}). Many PARI
765
objects also have a second codeword (types \typ{INT}, \typ{REAL},
766
\typ{PADIC}, \typ{POL}, and \typ{SER}). In case you want to produce one of
767
those from scratch, which should be exceedingly rare, \emph{it is your
768
responsibility to fill this second codeword}, either explicitly (using the
769
macros described in \secref{se:impl}), or implicitly using an assignment
770
statement (using \kbd{gaffect}).
771
772
Note that the length argument $l$ is predetermined for a number of types:
773
3 for types \typ{INTMOD}, \typ{FRAC}, \typ{COMPLEX}, \typ{POLMOD},
774
\typ{RFRAC}, 4 for type \typ{QUAD}, and 5 for type \typ{PADIC} and \typ{QFB}.
775
However for the sake of efficiency, \kbd{cgetg} does not check this: disasters
776
will occur if you give an incorrect length for those types.
777
778
\misctitle{Notes} 1) The main use of this function is create efficiently
779
a constant object, or to prepare for later assignments (see
780
\secref{se:assign}). Most of the time you will use \kbd{GEN} objects as they
781
are created and returned by PARI functions. In this case you do not need to
782
use \kbd{cgetg} to create space to hold them.
783
784
\noindent 2) For the creation of leaves, i.e.~\typ{INT} or \typ{REAL},
785
786
\fun{GEN}{cgeti}{long length}
787
788
\fun{GEN}{cgetr}{long length}
789
790
\noindent should be used instead of \kbd{cgetg(length, t\_INT)} and
791
\kbd{cgetg(length, t\_REAL)} respectively. Finally
792
793
\fun{GEN}{cgetc}{long prec}
794
795
\noindent creates a \typ{COMPLEX} whose real and imaginary part are
796
\typ{REAL}s allocated by \kbd{cgetr(prec)}.
797
798
\misctitle{Examples} 1) Both \kbd{z = cgeti(DEFAULTPREC)} and
799
\kbd{cgetg(DEFAULTPREC, t\_INT)} create a \typ{INT} whose ``precision'' is
800
\kbd{bit\_accuracy(DEFAULTPREC)} = 64. This means \kbd{z} can hold rational
801
integers of absolute value less than $2^{64}$. Note that in both cases, the
802
second codeword is \emph{not} filled. Of course we could use numerical
803
values, e.g.~\kbd{cgeti(4)}, but this would have different meanings on
804
different machines as \kbd{bit\_accuracy(4)} equals 64 on 32-bit machines,
805
but 128 on 64-bit machines.
806
807
\noindent 2) The following creates a \emph{complex number} whose real and
808
imaginary parts can hold real numbers of precision
809
$\kbd{bit\_accuracy(MEDDEFAULTPREC)} = 96\hbox{ bits:}$
810
%
811
\bprog
812
z = cgetg(3, t_COMPLEX);
813
gel(z, 1) = cgetr(MEDDEFAULTPREC);
814
gel(z, 2) = cgetr(MEDDEFAULTPREC);
815
@eprog\noindent
816
or simply \kbd{z = cgetc(MEDDEFAULTPREC)}.
817
818
\noindent 3) To create a matrix object for $4\times 3$ matrices:
819
%
820
\bprog
821
z = cgetg(4, t_MAT);
822
for(i=1; i<4; i++) gel(z, i) = cgetg(5, t_COL);
823
@eprog\noindent
824
or simply \kbd{z = zeromatcopy(4, 3)}, which further initializes all entries
825
to \kbd{gen\_0}.
826
827
These last two examples illustrate the fact that since PARI types are
828
recursive, all the branches of the tree must be created. The function
829
\teb{cgetg} creates only the ``root'', and other calls to \kbd{cgetg} must be
830
made to produce the whole tree. For matrices, a common mistake is to think
831
that \kbd{z = cgetg(4, t\_MAT)} (for example) creates the root of the
832
matrix: one needs also to create the column vectors of the matrix (obviously,
833
since we specified only one dimension in the first \kbd{cgetg}!). This is
834
because a matrix is really just a row vector of column vectors (hence a
835
priori not a basic type), but it has been given a special type number so that
836
operations with matrices become possible.
837
838
Finally, to facilitate input of constant objects when speed is not paramount,
839
there are four \tet{varargs} functions:
840
841
\fun{GEN}{mkintn}{long n, ...}
842
returns the nonnegative \typ{INT} whose development in base $2^{32}$
843
is given by the following $n$ 32bit-words (\kbd{unsigned int}).
844
\bprog
845
mkintn(3, a2, a1, a0);
846
@eprog
847
\noindent returns $a_2 2^{64} + a_1 2^{32} + a_0$.
848
849
\fun{GEN}{mkpoln}{long n, ...}
850
Returns the \typ{POL} whose $n$ coefficients (\kbd{GEN}) follow, in order of
851
decreasing degree.
852
\bprog
853
mkpoln(3, gen_1, gen_2, gen_0);
854
@eprog
855
\noindent returns the polynomial $X^2 + 2X$ (in variable $0$, use
856
\tet{setvarn} if you want other variable numbers). Beware that $n$ is the
857
number of coefficients, hence \emph{one more} than the degree.
858
859
\fun{GEN}{mkvecn}{long n, ...}
860
returns the \typ{VEC} whose $n$ coefficients (\kbd{GEN}) follow.
861
862
\fun{GEN}{mkcoln}{long n, ...}
863
returns the \typ{COL} whose $n$ coefficients (\kbd{GEN}) follow.
864
865
\misctitle{Warning} Contrary to the policy of general PARI functions, the
866
latter three functions do \emph{not} copy their arguments, nor do they produce
867
an object a priori suitable for \tet{gerepileupto}. For instance
868
\bprog
869
/*@Ccom gerepile-safe: components are universal objects */
870
z = mkvecn(3, gen_1, gen_0, gen_2);
871
872
/*@Ccom not OK for gerepileupto: stoi(3) creates component before root */
873
z = mkvecn(3, stoi(3), gen_0, gen_2);
874
875
/*@Ccom NO! First vector component \kbd{x} is destroyed */
876
x = gclone(gen_1);
877
z = mkvecn(3, x, gen_0, gen_2);
878
gunclone(x);
879
@eprog
880
881
\noindent The following function is also available as a special case of
882
\tet{mkintn}:
883
884
\fun{GEN}{uu32toi}{ulong a, ulong b}
885
886
Returns the \kbd{GEN} equal to $2^{32} a + b$, \emph{assuming} that
887
$a,b < 2^{32}$. This does not depend on \kbd{sizeof(long)}: the behavior is
888
as above on both $32$ and $64$-bit machines.
889
890
\subsec{Sizes}
891
892
\fun{long}{gsizeword}{GEN x} returns the total number of \B-bit words occupied
893
by the tree representing~\kbd{x}.
894
895
\fun{long}{gsizebyte}{GEN x} returns the total number of bytes occupied
896
by the tree representing~\kbd{x}, i.e.~\kbd{gsizeword(x)} multiplied by
897
\kbd{sizeof(long)}. This is normally useless since PARI functions use
898
a number of \emph{words} as input for lengths and precisions.
899
900
\subsec{Assignments}
901
Firstly, if \kbd{x} and \kbd{y} are both declared as \kbd{GEN} (i.e.~pointers
902
to something), the ordinary C assignment \kbd{y = x} makes perfect sense: we
903
are just moving a pointer around. However, physically modifying either
904
\kbd{x} or \kbd{y} (for instance, \kbd{x[1] = 0}) also changes the other
905
one, which is usually not desirable. \label{se:assign}
906
907
\misctitle{Very important note} Using the functions described in this
908
paragraph is inefficient and often awkward: one of the \tet{gerepile}
909
functions (see~\secref{se:garbage}) should be preferred. See the paragraph
910
end for one exception to this rule.
911
912
\noindent
913
The general PARI \idx{assignment} function is the function \teb{gaffect} with
914
the following syntax:
915
916
\fun{void}{gaffect}{GEN x, GEN y}
917
918
\noindent
919
Its effect is to assign the PARI object \kbd{x} into the \emph{preexisting}
920
object \kbd{y}. Both \kbd{x} and \kbd{y} must be \emph{scalar} types. For
921
convenience, vector or matrices of scalar types are also allowed.
922
923
This copies the whole structure of \kbd{x} into \kbd{y} so many conditions
924
must be met for the assignment to be possible. For instance it is allowed to
925
assign a \typ{INT} into a \typ{REAL}, but the converse is forbidden. For
926
that, you must use the truncation or rounding function of your choice,
927
e.g.\kbd{mpfloor}.
928
929
It can also happen that \kbd{y} is not large enough or does not have the proper
930
tree structure to receive the object \kbd{x}. For instance, let \kbd{y} the zero
931
integer with length equal to 2; then \kbd{y} is too small to accommodate any
932
nonzero \typ{INT}. In general common sense tells you what is possible,
933
keeping in mind the PARI philosophy which says that if it makes sense it is
934
valid. For instance, the assignment of an imprecise object into a precise one
935
does \emph{not} make sense. However, a change in precision of imprecise
936
objects is allowed, even if it \emph{increases} its accuracy: we complement
937
the ``mantissa'' with infinitely many $0$ digits in this case. (Mantissa
938
between quotes, because this is not restricted to \typ{REAL}s, it also
939
applies for $p$-adics for instance.)
940
941
All functions ending in ``\kbd{z}'' such as \teb{gaddz}
942
(see~\secref{se:low_level}) implicitly use this function. In fact what they
943
exactly do is record {\teb{avma}} (see~\secref{se:garbage}), perform the
944
required operation, \teb{gaffect} the result to the last operand, then
945
restore the initial \kbd{avma}.
946
947
You can assign ordinary C long integers into a PARI object (not necessarily
948
of type \typ{INT}) using
949
950
\fun{void}{gaffsg}{long s, GEN y}
951
952
\misctitle{Note} Due to the requirements mentioned above, it is usually
953
a bad idea to use \tet{gaffect} statements. There is one exception: for simple
954
objects (e.g.~leaves) whose size is controlled, they can be easier to use than
955
\kbd{gerepile}, and about as efficient.
956
957
\misctitle{Coercion} It is often useful to coerce an inexact object to a
958
given precision. For instance at the beginning of a routine where precision
959
can be kept to a minimum; otherwise the precision of the input is used in all
960
subsequent computations, which is inefficient if the latter is known to
961
thousands of digits. One may use the \kbd{gaffect} function for this, but it
962
is easier and more efficient to call
963
964
\fun{GEN}{gtofp}{GEN x, long prec} converts the complex number~\kbd{x}
965
(\typ{INT}, \typ{REAL}, \typ{FRAC}, \typ{QUAD} or \typ{COMPLEX}) to either
966
a \typ{REAL} or \typ{COMPLEX} whose components are \typ{REAL} of length
967
\kbd{prec}.
968
969
\subsec{Copy} It is also very useful to \idx{copy} a PARI object, not
970
just by moving around a pointer as in the \kbd{y = x} example, but by
971
creating a copy of the whole tree structure, without pre-allocating a
972
possibly complicated \kbd{y} to use with \kbd{gaffect}. The function which
973
does this is called \teb{gcopy}. Its syntax is:
974
975
\fun{GEN}{gcopy}{GEN x}
976
977
\noindent and the effect is to create a new copy of x on the PARI stack.
978
979
Sometimes, on the contrary, a quick copy of the skeleton of \kbd{x} is
980
enough, leaving pointers to the original data in \kbd{x} for the sake of
981
speed instead of making a full recursive copy. Use
982
\fun{GEN}{shallowcopy}{GEN x} for this. Note that the result is not suitable
983
for \tet{gerepileupto} !
984
985
Make sure at this point that you understand the difference between \kbd{y =
986
x}, \kbd{y = gcopy(x)}, \kbd{y = shallowcopy(x)} and \kbd{gaffect(x,y)}.
987
988
\subsec{Clones}\sidx{clone}\label{se:clone}
989
Sometimes, it is more efficient to create a \emph{persistent} copy of a PARI
990
object. This is not created on the stack but on the heap, hence unaffected by
991
\tet{gerepile} and friends. The function which does this is called
992
\teb{gclone}. Its syntax is:
993
994
\fun{GEN}{gclone}{GEN x}
995
996
A clone can be removed from the heap (thus destroyed) using
997
998
\fun{void}{gunclone}{GEN x}
999
1000
\noindent No PARI object should keep references to a clone which has been
1001
destroyed!
1002
1003
\subsec{Conversions}\sidx{conversions}
1004
The following functions convert C objects to PARI objects (creating them on
1005
the stack as usual):
1006
1007
\fun{GEN}{stoi}{long s}: C \kbd{long} integer (``small'') to \typ{INT}.
1008
1009
\fun{GEN}{dbltor}{double s}: C \kbd{double} to \typ{REAL}. The accuracy of
1010
the result is 19 decimal digits, i.e.~a type \typ{REAL} of length
1011
\kbd{DEFAULTPREC}, although on 32-bit machines only 16 of them are
1012
significant.
1013
1014
\noindent We also have the converse functions:
1015
1016
\fun{long}{itos}{GEN x}: \kbd{x} must be of type \typ{INT},
1017
1018
\fun{double}{rtodbl}{GEN x}: \kbd{x} must be of type \typ{REAL},
1019
1020
\noindent as well as the more general ones:
1021
1022
\fun{long}{gtolong}{GEN x},
1023
1024
\fun{double}{gtodouble}{GEN x}.
1025
1026
\section{Implementation of the PARI types}
1027
\label{se:impl}
1028
1029
\noindent
1030
We now go through each type and explain its implementation. Let \kbd{z} be a
1031
\kbd{GEN}, pointing at a PARI object. In the following paragraphs, we will
1032
constantly mix two points of view: on the one hand, \kbd{z} is treated as the
1033
C pointer it is, on the other, as PARI's handle on some mathematical entity,
1034
so we will shamelessly write $\kbd{z} \ne 0$ to indicate that the
1035
\emph{value} thus represented is nonzero (in which case the
1036
\emph{pointer}~\kbd{z} is certainly not \kbd{NULL}). We offer no apologies
1037
for this style. In fact, you had better feel comfortable juggling both views
1038
simultaneously in your mind if you want to write correct PARI programs.
1039
1040
Common to all the types is the first codeword \kbd{z[0]}, which we do not
1041
have to worry about since this is taken care of by \kbd{cgetg}. Its precise
1042
structure depends on the machine you are using, but it always contains the
1043
following data: the \emph{internal type number}\sidx{type number} attached
1044
to the symbolic type name, the \emph{length} of the root in longwords, and a
1045
technical bit which indicates whether the object is a clone or not (see
1046
\secref{se:clone}). This last one is used by \kbd{gp} for internal garbage
1047
collecting, you will not have to worry about it.
1048
1049
Some types have a second codeword, different for each type, which we will
1050
soon describe as we will shortly consider each of them in turn.
1051
1052
\noindent The first codeword is handled through the following \emph{macros}:
1053
1054
\fun{long}{typ}{GEN z} returns the type number of \kbd{z}.
1055
1056
\fun{void}{settyp}{GEN z, long n} sets the type number of \kbd{z} to
1057
\kbd{n} (you should not have to use this function if you use \kbd{cgetg}).
1058
1059
\fun{long}{lg}{GEN z} returns the length (in longwords) of the root of \kbd{z}.
1060
1061
\fun{long}{setlg}{GEN z, long l} sets the length of \kbd{z} to \kbd{l}; you
1062
should not have to use this function if you use \kbd{cgetg}.
1063
1064
\fun{void}{lg_increase}{GEN z} increase the length of \kbd{z} by 1; you
1065
should not have to use this function if you use \kbd{cgetg}.
1066
1067
\fun{long}{isclone}{GEN z} is \kbd{z} a clone?
1068
1069
\fun{void}{setisclone}{GEN z} sets the \emph{clone} bit.
1070
1071
\fun{void}{unsetisclone}{GEN z} clears the \emph{clone} bit.
1072
1073
\misctitle{Important remark} For the sake of efficiency, none of the
1074
codeword-handling macros check the types of their arguments even when there
1075
are stringent restrictions on their use. It is trivial to create invalid
1076
objects, or corrupt one of the ``universal constants'' (e.g. setting the sign
1077
of \kbd{gen\_0} to $1$), and they usually provide negligible savings.
1078
Use higher level functions whenever possible.
1079
1080
\misctitle{Remark} The clone bit is there so that \kbd{gunclone} can check
1081
it is deleting an object which was allocated by \kbd{gclone}. Miscellaneous
1082
vector entries are often cloned by \kbd{gp} so that a GP statement like
1083
\kbd{v[1] = x} does not involve copying the whole of \kbd{v}: the component
1084
\kbd{v[1]} is deleted if its clone bit is set, and is replaced by a clone of
1085
\kbd{x}. Don't set/unset yourself the clone bit unless you know what you are
1086
doing: in particular \emph{never} set the clone bit of a vector component
1087
when the said vector is scheduled to be uncloned. Hackish code may abuse the
1088
clone bit to tag objects for reasons unrelated to the above instead of using
1089
proper data structures. Don't do that.
1090
1091
\subsec{Type \typ{INT} (integer)}
1092
\sidx{integer}\kbdsidx{t_INT}this type has a second codeword \kbd{z[1]} which
1093
contains the following information:
1094
1095
the sign of \kbd{z}: coded as $1$, $0$ or $-1$ if $\kbd{z} > 0$, $\kbd{z} = 0$,
1096
$\kbd{z} < 0$ respectively.
1097
1098
the \emph{effective length} of \kbd{z}, i.e.~the total number of significant
1099
longwords. This means the following: apart from the integer 0, every integer
1100
is ``normalized'', meaning that the most significant mantissa longword is
1101
nonzero. However, the integer may have been created with a longer length.
1102
Hence the ``length'' which is in \kbd{z[0]} can be larger than the
1103
``effective length'' which is in \kbd{z[1]}.
1104
1105
\noindent This information is handled using the following macros:
1106
1107
\fun{long}{signe}{GEN z} returns the sign of \kbd{z}.
1108
1109
\fun{void}{setsigne}{GEN z, long s} sets the sign of \kbd{z} to \kbd{s}.
1110
1111
\fun{long}{lgefint}{GEN z} returns the \idx{effective length} of \kbd{z}.
1112
1113
\fun{void}{setlgefint}{GEN z, long l} sets the effective length
1114
of \kbd{z} to \kbd{l}.
1115
1116
The integer 0 can be recognized either by its sign being~0, or by its
1117
effective length being equal to~2. Now assume that $\kbd{z} \ne 0$, and let
1118
$$ |z| = \sum_{i = 0}^n z_i B^i,
1119
\quad\text{where}~z_n\ne 0~\text{and}~B = 2^{\kbd{BITS\_IN\_LONG}}.
1120
$$
1121
With these notations, $n$ is \kbd{lgefint(z) - 3}, and the mantissa of
1122
$\kbd{z}$ may be manipulated via the following interface:
1123
1124
\fun{GEN}{int_MSW}{GEN z} returns a pointer to the most significant word of
1125
\kbd{z}, $z_n$.
1126
1127
\fun{GEN}{int_LSW}{GEN z} returns a pointer to the least significant word of
1128
\kbd{z}, $z_0$.
1129
1130
\fun{GEN}{int_W}{GEN z, long i} returns the $i$-th significant word of
1131
\kbd{z}, $z_i$. Accessing the $i$-th significant word for $i > n$
1132
yields unpredictable results.
1133
1134
\fun{GEN}{int_W_lg}{GEN z, long i, long lz} returns the $i$-th significant
1135
word of \kbd{z}, $z_i$, assuming \kbd{lgefint(z)} is \kbd{lz} ($= n + 3$).
1136
Accessing the $i$-th significant word for $i > n$ yields unpredictable
1137
results.
1138
1139
\fun{GEN}{int_precW}{GEN z} returns the previous (less significant) word of
1140
\kbd{z}, $z_{i-1}$ assuming \kbd{z} points to $z_i$.
1141
1142
\fun{GEN}{int_nextW}{GEN z} returns the next (more significant) word of \kbd{z},
1143
$z_{i+1}$ assuming \kbd{z} points to $z_i$.
1144
1145
Unnormalized integers, such that $z_n$ is possibly $0$, are explicitly
1146
forbidden. To enforce this, one may write an arbitrary mantissa then call
1147
1148
\fun{void}{int_normalize}{GEN z, long known0}
1149
1150
\noindent normalizes in place a nonnegative integer (such that $z_n$ is
1151
possibly $0$), assuming at least the first \kbd{known0} words are zero.
1152
1153
\noindent For instance a binary \kbd{and} could be implemented in the
1154
following way:
1155
\bprog
1156
GEN AND(GEN x, GEN y) {
1157
long i, lx, ly, lout;
1158
long *xp, *yp, *outp; /* mantissa pointers */
1159
GEN out;
1160
1161
if (!signe(x) || !signe(y)) return gen_0;
1162
lx = lgefint(x); xp = int_LSW(x);
1163
ly = lgefint(y); yp = int_LSW(y); lout = min(lx,ly); /* > 2 */
1164
1165
out = cgeti(lout); out[1] = evalsigne(1) | evallgefint(lout);
1166
outp = int_LSW(out);
1167
for (i=2; i < lout; i++)
1168
{
1169
*outp = (*xp) & (*yp);
1170
outp = int_nextW(outp);
1171
xp = int_nextW(xp);
1172
yp = int_nextW(yp);
1173
}
1174
if ( !*int_MSW(out) ) out = int_normalize(out, 1);
1175
return out;
1176
}
1177
@eprog
1178
1179
\noindent This low-level interface is mandatory in order to write portable
1180
code since PARI can be compiled using various multiprecision kernels, for
1181
instance the native one or GNU MP, with incompatible internal structures
1182
(for one thing, the mantissa is oriented in different directions).
1183
1184
\subsec{Type \typ{REAL} (real number)}
1185
\kbdsidx{t_REAL}\sidx{real number}this type has a second codeword z[1] which
1186
also encodes its sign, obtained or set using the same functions as for a
1187
\typ{INT}, and a binary exponent. This exponent is handled using the
1188
following macros:
1189
1190
\fun{long}{expo}{GEN z} returns the exponent of \kbd{z}.
1191
This is defined even when \kbd{z} is equal to zero.
1192
1193
\fun{void}{setexpo}{GEN z, long e} sets the exponent of \kbd{z} to \kbd{e}.
1194
1195
\noindent Note the functions:
1196
1197
\fun{long}{gexpo}{GEN z} which tries to return an exponent for \kbd{z},
1198
even if \kbd{z} is not a real number.
1199
1200
\fun{long}{gsigne}{GEN z} which returns a sign for \kbd{z}, even when
1201
\kbd{z} is a real number of type \typ{INT}, \typ{FRAC} or \typ{REAL},
1202
an infinity (\typ{INFINITY}) or a \typ{QUAD} of positive discriminant.
1203
1204
The real zero is characterized by having its sign equal to 0. If \kbd{z} is
1205
not equal to~0, then it is represented as $2^e M$, where $e$ is the exponent,
1206
and $M\in [1, 2[$ is the mantissa of $z$, whose digits are stored in
1207
$\kbd{z[2]},\dots, \kbd{z[lg(z)-1]}$. For historical reasons, the \kbd{prec}
1208
parameter attached to floating point functions is measured in \B-bit words
1209
and is equal to the length of \kbd{x}: yes, this includes the two code words
1210
and depends on \kbd{sizeof(long)}. For clarity we advise to use
1211
\kbd{bit\_accuracy}, which computes the true length of the mantissa in bits,
1212
and convert between bits and \kbd{prec} using the \kbd{prec2nbits} and
1213
\kbd{nbits2prec} macros. But keep in mind that the accuracy of \typ{REAL}
1214
actually increases by increments of \B bits.
1215
1216
More precisely, let $m$ be the integer (\kbd{z[2]},\dots, \kbd{z[lg(z)-1]})
1217
in base \kbd{2\pow BITS\_IN\_LONG}; here, \kbd{z[2]} is the most significant
1218
longword and is normalized, i.e.~its most significant bit is~1. Then we have
1219
$M := m / 2^{\kbd{bit\_accuracy(lg(z))} - 1 - \kbd{expo}(z)}$.
1220
1221
\fun{GEN}{mantissa_real}{GEN z, long *e} returns the mantissa $m$ of $z$, and
1222
sets \kbd{*e} to the exponent $\kbd{bit\_accuracy(lg(z))}-1-\kbd{expo}(z)$,
1223
so that $z = m / 2^e$.
1224
1225
Thus, the real number $3.5$ to accuracy \kbd{bit\_accuracy(lg(z))} is
1226
represented as \kbd{z[0]} (encoding $\kbd{type} = \typ{REAL}$, \kbd{lg(z)}),
1227
\kbd{z[1]} (encoding $\kbd{sign} = 1$, $\kbd{expo} = 1$), $\kbd{z[2]} =
1228
\kbd{0xe0000000}$, $\kbd{z[3]} =\dots = \kbd{z[lg(z)-1]} = \kbd{0x0}$.
1229
1230
\subsec{Type \typ{INTMOD}}\kbdsidx{t_INTMOD}
1231
\kbd{z[1]} points to the modulus, and \kbd{z[2]} at the number representing
1232
the class \kbd{z}. Both are separate \kbd{GEN} objects, and both must be
1233
\typ{INT}s, satisfying the inequality $0 \le \kbd{z[2]} < \kbd{z[1]}$.
1234
1235
\subsec{Type \typ{FRAC} (rational number)}%
1236
\kbdsidx{t_FRAC}\sidx{rational number}
1237
\kbd{z[1]} points to the numerator $n$, and \kbd{z[2]} to the denominator
1238
$d$. Both must be of type \typ{INT} such that $n\neq 0$, $d > 0$ and
1239
$(n,d) = 1$.
1240
1241
\subsec{Type \typ{FFELT} (finite field element)}%
1242
\kbdsidx{t_FFELT}\sidx{finite field element} (Experimental)
1243
1244
Components of this type should normally not be accessed directly. Instead,
1245
finite field elements should be created using \kbd{ffgen}.
1246
1247
\noindent The second codeword \kbd{z[1]} determines the storage format of the
1248
element, among
1249
1250
\item \tet{t_FF_FpXQ}: \kbd{A=z[2]} and \kbd{T=z[3]} are \kbd{FpX},
1251
\kbd{p=z[4]} is a \typ{INT}, where $p$ is a prime number, $T$ is irreducible
1252
modulo $p$, and $\deg A < \deg T$.
1253
This represents the element $A\pmod{T}$ in $\F_p[X]/T$.
1254
1255
\item \tet{t_FF_Flxq}: \kbd{A=z[2]} and \kbd{T=z[3]} are \kbd{Flx},
1256
\kbd{l=z[4]} is a \typ{INT}, where $l$ is a prime number, $T$ is irreducible
1257
modulo $l$, and $\deg A < \deg T$ This represents the element $A\pmod{T}$ in
1258
$\F_l[X]/T$.
1259
1260
\item \tet{t_FF_F2xq}: \kbd{A=z[2]} and \kbd{T=z[3]} are \kbd{F2x},
1261
\kbd{l=z[4]} is the \typ{INT} $2$, $T$ is irreducible modulo $2$, and
1262
$\deg A < \deg T$. This represents the element $A\pmod{T}$ in $\F_2[X]/T$.
1263
1264
\subsec{Type \typ{COMPLEX} (complex number)}%
1265
\kbdsidx{t_COMPLEX}\sidx{complex number}
1266
\kbd{z[1]} points to the real part, and \kbd{z[2]} to the imaginary part.
1267
The components \kbd{z[1]} and \kbd{z[2]} must be of type
1268
\typ{INT}, \typ{REAL} or \typ{FRAC}. For historical reasons \typ{INTMOD}
1269
and \typ{PADIC} are also allowed (the latter for $p = 2$ or
1270
congruent to 3 mod 4 only), but one should rather use the more general
1271
\typ{POLMOD} construction.
1272
1273
\subsec{Type \typ{PADIC} ($p$-adic numbers)}%
1274
\sidx{p-adic number}\kbdsidx{t_PADIC} this type has a second codeword
1275
\kbd{z[1]} which contains the following information: the $p$-adic precision
1276
(the exponent of $p$ modulo which the $p$-adic unit corresponding to
1277
\kbd{z} is defined if \kbd{z} is not~0), i.e.~one less than the number of
1278
significant $p$-adic digits, and the exponent of \kbd{z}. This information
1279
can be handled using the following functions:
1280
1281
\fun{long}{precp}{GEN z} returns the $p$-adic precision of \kbd{z}. This is
1282
$0$ if $\kbd{z} = 0$.
1283
1284
\fun{void}{setprecp}{GEN z, long l} sets the $p$-adic precision of \kbd{z}
1285
to \kbd{l}.
1286
1287
\fun{long}{valp}{GEN z} returns the $p$-adic valuation of \kbd{z}
1288
(i.e. the exponent). This is defined even if \kbd{z} is equal to~0.
1289
1290
\fun{void}{setvalp}{GEN z, long e} sets the $p$-adic valuation of \kbd{z}
1291
to \kbd{e}.
1292
1293
In addition to this codeword, \kbd{z[2]} points to the prime $p$,
1294
\kbd{z[3]} points to $p^{\text{precp(z)}}$, and \kbd{z[4]} points to
1295
a\typ{INT} representing the $p$-adic unit attached to \kbd{z} modulo
1296
\kbd{z[3]} (and to zero if \kbd{z} is zero). To summarize, if $z\neq
1297
0$, we have the equality:
1298
$$ \kbd{z} = p^{\text{valp(z)}} * (\kbd{z[4]} + O(\kbd{z[3]})),\quad
1299
\text{where}\quad \kbd{z[3]} = O(p^{\text{precp(z)}}). $$
1300
1301
\subsec{Type \typ{QUAD} (quadratic number)}
1302
\sidx{quadratic number}\kbdsidx{t_QUAD}\kbd{z[1]} points to the canonical
1303
polynomial $P$ defining the quadratic field (as output by \tet{quadpoly}),
1304
\kbd{z[2]} to the ``real part'' and \kbd{z[3]} to the ``imaginary part''. The
1305
latter are of type \typ{INT}, \typ{FRAC}, \typ{INTMOD}, or \typ{PADIC} and
1306
are to be taken as the coefficients of \kbd{z} with respect to the canonical
1307
basis $(1,X)$ of $\Q[X]/(P(X))$. Exact complex numbers may be implemented as
1308
quadratics, but \typ{COMPLEX} is in general more versatile (\typ{REAL}
1309
components are allowed) and more efficient.
1310
1311
Operations involving a \typ{QUAD} and \typ{COMPLEX} are implemented by
1312
converting the \typ{QUAD} to a \typ{REAL} (or \typ{COMPLEX} with \typ{REAL}
1313
components) to the accuracy of the \typ{COMPLEX}. As a consequence,
1314
operations between \typ{QUAD} and \emph{exact} \typ{COMPLEX}s are not allowed.
1315
1316
\subsec{Type \typ{POLMOD} (polmod)}\kbdsidx{t_POLMOD}\sidx{polmod}
1317
as for \typ{INTMOD}s, \kbd{z[1]} points to the modulus, and \kbd{z[2]}
1318
to a polynomial representing the class of~\kbd{z}. Both must be of type
1319
\typ{POL} in the same variable, satisfying the inequality $\deg \kbd{z[2]}
1320
< \deg \kbd{z[1]}$. However, \kbd{z[2]} is allowed to be a simplification
1321
of such a polynomial, e.g.~a scalar. This is tricky considering the
1322
hierarchical structure of the variables; in particular, a polynomial in
1323
variable of \emph{lesser} priority (see \secref{se:vars}) than the
1324
modulus variable is valid, since it is considered as the constant term of
1325
a polynomial of degree 0 in the correct variable. On the other hand a
1326
variable of \emph{greater} priority is not acceptable.
1327
1328
\subsec{Type \typ{POL} (polynomial)}\kbdsidx{t_POL}\sidx{polynomial} this
1329
type has a second codeword. It contains a ``\emph{sign}'': 0 if the
1330
polynomial is equal to~0, and 1 if not (see however the important remark
1331
below) and a \emph{variable number} (e.g.~0 for $x$, 1 for $y$, etc\dots).
1332
1333
\noindent These data can be handled with the following macros: \teb{signe}
1334
and \teb{setsigne} as for \typ{INT} and \typ{REAL},
1335
1336
\fun{long}{varn}{GEN z} returns the variable number of the object \kbd{z},
1337
1338
\fun{void}{setvarn}{GEN z, long v} sets the variable number of \kbd{z} to
1339
\kbd{v}.
1340
1341
The variable numbers encode the relative priorities of variables, we will
1342
give more details in \secref{se:vars}. Note also the function
1343
\fun{long}{gvar}{GEN z} which tries to return a \idx{variable number} for
1344
\kbd{z}, even if \kbd{z} is not a polynomial or power series. The variable
1345
number of a scalar type is set by definition equal to \tet{NO_VARIABLE},
1346
which has lower priority than any other variable number.
1347
1348
The components \kbd{z[2]}, \kbd{z[3]},\dots \kbd{z[lg(z)-1]} point to the
1349
coefficients of the polynomial \emph{in ascending order}, with \kbd{z[2]}
1350
being the constant term and so on.
1351
1352
For a \typ{POL} of nonzero sign, \tet{degpol}, \tet{leading_coeff},
1353
\tet{constant_coeff}, return its degree, and a pointer to the leading,
1354
resp. constant, coefficient with respect to the main variable. Note that no
1355
copy is made on the PARI stack so the returned value is not safe for a basic
1356
\kbd{gerepile} call. Applied to any other type than \typ{POL}, the result is
1357
unspecified. Those three functions are still defined when the sign is $0$,
1358
see \secref{se:accessors} and \secref{se:polynomials}.
1359
1360
\fun{long}{degree}{GEN x} returns the degree of \kbd{x} with respect to its
1361
main variable even when \kbd{x} is not a polynomial (a rational function for
1362
instance). By convention, the degree of a zero polynomial is~$-1$.
1363
1364
\misctitle{Important remark}
1365
The leading coefficient of a \typ{POL} may be equal to zero:
1366
1367
\item it is not allowed to be an exact rational $0$, such as \tet{gen_0};
1368
1369
\item an exact nonrational $0$, like \kbd{Mod(0,2)}, is possible for
1370
constant polynomials, i.e. of length $3$ and no other coefficient: this
1371
carries information about the base ring for the polynomial;
1372
1373
\item an inexact $0$, like \kbd{0.E-38} or \kbd{O(3\pow 5)}, is always
1374
possible. Inexact zeroes do not correspond to an actual $0$, but to a
1375
very small coefficient according to some metric; we keep them to give
1376
information on how much cancellation occurred in previous computations.
1377
1378
A polynomial disobeying any of these rules is an invalid \emph{unnormalized}
1379
object. We advise \emph{not} to use low-level constructions to build a
1380
\typ{POL} coefficient by coefficient, such as
1381
\bprog
1382
GEN T = cgetg(4, t_POL);
1383
T[1] = evalvarn(0);
1384
gel(T, 2) = x;
1385
gel(T, 3) = y;
1386
@eprog\noindent But if you do and it is not clear whether the result will be
1387
normalized, call
1388
1389
\fun{GEN}{normalizepol}{GEN x} applied to an unnormalized \typ{POL}~\kbd{x}
1390
(with all coefficients correctly set except that \kbd{leading\_term(x)} might
1391
be zero), normalizes \kbd{x} correctly in place and returns~\kbd{x}. This
1392
functions sets \kbd{signe} (to $0$ or $1$) properly.
1393
1394
\misctitle{Caveat} A consequence of the remark above is that zero polynomials
1395
are characterized by the fact that their sign is~0. It is in general
1396
incorrect to check whether \kbd{lg(x)} is $2$ or \kbd{degpol(x)} $< 0$,
1397
although both tests are valid when the coefficient types are under control:
1398
for instance, when they are guaranteed to be \typ{INT}s or \typ{FRAC}s.
1399
The same remark applies to \typ{SER}s.
1400
1401
\subsec{Type \typ{SER} (power series)}
1402
\kbdsidx{t_SER}\sidx{power series}This type also has a second codeword, which
1403
encodes a ``\emph{sign}'', i.e.~0 if the power series is 0, and 1 if not, a
1404
\emph{variable number} as for polynomials, and an \emph{exponent}. This
1405
information can be handled with the following functions: \teb{signe},
1406
\teb{setsigne}, \teb{varn}, \teb{setvarn} as for polynomials, and \teb{valp},
1407
\teb{setvalp} for the exponent as for $p$-adic numbers. Beware: do \emph{not}
1408
use \teb{expo} and \teb{setexpo} on power series.
1409
1410
The coefficients \kbd{z[2]}, \kbd{z[3]},\dots \kbd{z[lg(z)-1]} point to
1411
the coefficients of \kbd{z} in ascending order. As for polynomials
1412
(see remark there), the sign of a \typ{SER} is $0$ if and only all
1413
its coefficients are equal to $0$. (The leading coefficient cannot be an
1414
integer $0$.) A series whose coefficients are integers equal to zero
1415
is represented as $O(x^n)$ (\kbd{zeroser}$(\var{vx},n)$).
1416
A series whose coefficients are exact zeroes, but not all of
1417
them integers (e.g. an \typ{INTMOD} such as \kbd{Mod(0,2)}) is
1418
represented as $z*x^{n-1} +O(x^n)$, where $z$ is the $0$ of the
1419
base ring, as per \tet{Rg_get_0}.
1420
1421
Note that the exponent of a power series can be negative, i.e.~we are
1422
then dealing with a Laurent series (with a finite number of negative
1423
terms).
1424
1425
\subsec{Type \typ{RFRAC} (rational function)}%
1426
\kbdsidx{t_RFRAC}\sidx{rational function} \kbd{z[1]} points to the
1427
numerator $n$, and \kbd{z[2]} on the denominator $d$. The denominator must be
1428
of type \typ{POL}, with variable of higher priority than the numerator. The
1429
numerator $n$ is not an exact $0$ and $(n,d) = 1$ (see \tet{gred_rfac2}).
1430
1431
\subsec{Type \typ{QFB} (binary quadratic form)}%
1432
\kbdsidx{t_QFB}\sidx{binary quadratic form} \kbd{z[1]}, \kbd{z[2]},
1433
\kbd{z[3]} point to the three coefficients of the form, and \kbd{z[4]} point to
1434
the form discriminant. All four are of type \typ{INT}.
1435
1436
\subsec{Type \typ{VEC} and \typ{COL} (vector)}%
1437
\kbdsidx{t_VEC}\kbdsidx{t_COL}\sidx{row vector}\sidx{column vector}
1438
\kbd{z[1]}, \kbd{z[2]},\dots \kbd{z[lg(z)-1]} point to the components of the
1439
vector.
1440
1441
\subsec{Type \typ{MAT} (matrix)}\kbdsidx{t_MAT}\sidx{matrix} \kbd{z[1]},
1442
\kbd{z[2]},\dots \kbd{z[lg(z)-1]} point to the column vectors of \kbd{z},
1443
i.e.~they must be of type \typ{COL} and of the same length.
1444
1445
\subsec{Type \typ{VECSMALL} (vector of small integers)}\kbdsidx{t_VECSMALL}
1446
\kbd{z[1]}, \kbd{z[2]},\dots \kbd{z[lg(z)-1]} are ordinary signed long
1447
integers. This type is used instead of a \typ{VEC} of \typ{INT}s for
1448
efficiency reasons, for instance to implement efficiently permutations,
1449
polynomial arithmetic and linear algebra over small finite fields, etc.
1450
1451
\subsec{Type \typ{STR} (character string)}%
1452
\kbdsidx{t_STR}\sidx{character string}
1453
1454
\fun{char *}{GSTR}{z} (= \kbd{(z+1)}) points to the first character of the
1455
(\kbd{NULL}-terminated) string.
1456
1457
\subsec{Type \typ{ERROR} (error context)}\kbdsidx{t_ERROR}
1458
This type holds error messages, as well as details about the error, as
1459
returned by the exception handling system. The second codeword \kbd{z[1]}
1460
contains the error type (an \kbd{int}, as passed to \tet{pari_err}). The
1461
subsequent words \kbd{z[2]},\dots \kbd{z[lg(z)-1]} are \kbd{GEN}s containing
1462
additional data, depending on the error type.\sidx{error context}
1463
1464
\subsec{Type \typ{CLOSURE} (closure)}\kbdsidx{t_CLOSURE}\sidx{closure}
1465
This type holds GP functions and closures, in compiled form. The internal detail
1466
of this type is subject to change each time the GP language evolves. Hence we
1467
do not describe it here and refer to the Developer's Guide. However
1468
functions to create or to evaluate \typ{CLOSURE}s are documented in
1469
\secref{se:closure}.
1470
1471
\fun{long}{closure_arity}{GEN C} returns the arity of the \typ{CLOSURE}.
1472
1473
\fun{long}{closure_is_variadic}{GEN C} returns $1$ if the closure \kbd{C} is
1474
variadic, $0$ else.
1475
1476
\subsec{Type \typ{INFINITY} (infinity)}\kbdsidx{t_INFINITY}\sidx{infinity}
1477
1478
This type has a single \typ{INT} component, which is either $1$ or $-1$,
1479
corresponding to $+\infty$ and $-\infty$ respectively.
1480
1481
\fun{GEN}{mkmoo}{} returns $-\infty$
1482
1483
\fun{GEN}{mkoo}{} returns $\infty$
1484
1485
\fun{long}{inf_get_sign}{GEN x} returns $1$ if $x$ is $+\infty$, and $-1$
1486
if $x$ is $-\infty$.
1487
1488
\subsec{Type \typ{LIST} (list)}\kbdsidx{t_LIST}\sidx{list}
1489
this type was introduced for specific \kbd{gp} use and is rather inefficient
1490
compared to a straightforward linked list implementation (it requires more
1491
memory, as well as many unnecessary copies). Hence we do not describe it
1492
here and refer to the Developer's Guide.
1493
1494
\misctitle{Implementation note} For the types including an exponent (or a
1495
valuation), we actually store a biased nonnegative exponent (bit-ORing the
1496
biased exponent to the codeword), obtained by adding a constant to the true
1497
exponent: either \kbd{HIGHEXPOBIT} (for \typ{REAL}) or \kbd{HIGHVALPBIT} (for
1498
\typ{PADIC} and \typ{SER}). Of course, this is encapsulated by the
1499
exponent/valuation-handling macros and needs not concern the library user.
1500
1501
\section{PARI variables}\label{se:vars}
1502
\subsec{Multivariate objects}
1503
\sidx{variable (priority)}
1504
1505
\noindent We now consider variables and formal computations. As we have seen
1506
in \secref{se:impl}, the codewords for types \typ{POL} and \typ{SER} encode a
1507
``variable number''. This is an integer, ranging from $0$ to \kbd{MAXVARN}.
1508
Relative priorities may be ascertained using
1509
1510
\fun{int}{varncmp}{long v, long w}
1511
1512
\noindent which is $>0$, $=0$, $<0$ whenever $v$ has lower, resp.~same,
1513
resp.~higher priority than $w$.
1514
1515
The way an object is considered in formal computations depends entirely on
1516
its ``principal variable number'' which is given by the function
1517
1518
\fun{long}{gvar}{GEN z}
1519
1520
\noindent which returns a \idx{variable number} for \kbd{z}, even if \kbd{z}
1521
is not a polynomial or power series. The variable number of a scalar type is
1522
set by definition equal to \tet{NO_VARIABLE} which has lower priority than any
1523
valid variable number. The variable number of a recursive type which is not a
1524
polynomial or power series is the variable number with highest priority among
1525
its components. But for polynomials and power series only the ``outermost''
1526
number counts (we directly access $\tet{varn}(x)$ in the codewords): the
1527
representation is not symmetrical at all.
1528
1529
Under \kbd{gp}, one needs not worry too much since the interpreter defines
1530
the variables as it sees them\footnote{*}{The first time a given identifier
1531
is read by the GP parser a new variable is created, and it is assigned a
1532
strictly lower priority than any variable in use at this point. On startup,
1533
before any user input has taken place, 'x' is defined in this way and has
1534
initially maximal priority (and variable number $0$).}
1535
%
1536
and do the right thing with the polynomials produced.
1537
1538
But in library mode, they are tricky objects if you intend to build
1539
polynomials yourself (and not just let PARI functions produce them, which is
1540
less efficient). For instance, it does not make sense to have a variable
1541
number occur in the components of a polynomial whose main variable has a
1542
lower priority, even though PARI cannot prevent you from doing it.
1543
1544
\subsec{Creating variables} A basic difficulty is to ``create'' a variable.
1545
Some initializations are needed before you can use a given integer $v$ as a
1546
variable number.
1547
1548
Initially, this is done for $0$ and $1$ (the variables \kbd{x} and
1549
\kbd{y} under \kbd{gp}), and $2,\dots,9$ (printed as \kbd{t2}, \dots
1550
\kbd{t9}), with decreasing priority.
1551
1552
\subsubsec{User variables}\sidx{variable (user)} When the program starts,
1553
\kbd{x} (number~$0$) and \kbd{y} (number~$1$) are the only available
1554
variables, numbers $2$ to $9$ (decreasing priority) are reserved for building
1555
polynomials with predictable priorities.
1556
1557
To define further ones, you may use
1558
1559
\fun{GEN}{varhigher}{const char *s}
1560
1561
\fun{GEN}{varlower}{const char *s}
1562
1563
to recover a monomial of degree $1$ in a new variable, which is guaranteed
1564
to have higer (resp.~lower) priority than all existing ones at the
1565
time of the function call. The variable is printed as $s$, but is not part
1566
of GP's interpreter: it is not a symbol bound to a value.
1567
1568
On the other hand
1569
1570
\fun{long}{fetch_user_var}{char *s}: inspects the user variable whose name is
1571
the string pointed to by \kbd{s}, creating it if needed, and returns its
1572
variable number.
1573
\bprog
1574
long v = fetch_user_var("y");
1575
GEN gy = pol_x(v);
1576
@eprog\noindent
1577
The function raises an exception if the name is already in use for an
1578
\tet{install}ed or built-in function, or an alias. This function
1579
is mostly useless since it returns a variable with unpredictable priority.
1580
Don't use it to create new variables.
1581
1582
\misctitle{Caveat} You can use \tet{gp_read_str}
1583
(see~\secref{se:gp_read_str}) to execute a GP command and create GP
1584
variables on the fly as needed:
1585
\bprog
1586
GEN gy = gp_read_str("'y"); /*@Ccom returns \kbd{pol\_x}($v$), for some $v$ */
1587
long v = varn(gy);
1588
@eprog\noindent
1589
But please note the quote \kbd{'y} in the above. Using \kbd{gp\_read\_str("y")}
1590
might work, but is dangerous, especially when programming functions to
1591
be used under \kbd{gp}. The latter reads the value of \kbd{y}, as
1592
\emph{currently} known by the \kbd{gp} interpreter, possibly creating it
1593
in the process. But if \kbd{y} has been modified by previous \kbd{gp}
1594
commands (e.g.~\kbd {y = 1}), then the value of \kbd{gy} is not what you
1595
expected it to be and corresponds instead to the current value of the
1596
\kbd{gp} variable (e.g.~\kbd{gen\_1}).
1597
1598
\fun{GEN}{fetch_var_value}{long v} returns a shallow copy of the current
1599
value of the variable numbered $v$. Returns \kbd{NULL} if that variable
1600
number is unknown to the interpreter, e.g. it is a user variable. Note
1601
that this may not be the same as \kbd{pol\_x(v)} if assignments have been
1602
performed in the interpreter.
1603
1604
\subsubsec{Temporary variables}\sidx{variable (temporary)}
1605
You can create temporary variables using
1606
1607
\fun{long}{fetch_var}{}
1608
returns a new variable with \emph{lower} priority than any variable currently
1609
in use.
1610
1611
\fun{long}{fetch_var_higher}{}
1612
returns a new variable with \emph{higher} priority than any variable
1613
currently in use.
1614
1615
\noindent
1616
After the statement \kbd{v = fetch\_var()}, you can use
1617
\kbd{pol\_1(v)} and \kbd{pol\_x(v)}. The variables created in this way have no
1618
identifier assigned to them though, and are printed as
1619
\kbd{t\text{number}}. You can assign a name to a temporary variable, after
1620
creating it, by calling the function
1621
1622
\fun{void}{name_var}{long n, char *s}
1623
1624
\noindent after which the output machinery will use the name \kbd{s} to
1625
represent the variable number~\kbd{n}. The GP parser will \emph{not}
1626
recognize it by that name, however, and calling this on a variable known
1627
to~\kbd{gp} raises an error. Temporary variables are meant to be used as free
1628
variables to build polynomials and power series, and you should never assign
1629
values or functions to them as you would do with variables under~\kbd{gp}.
1630
For that, you need a user variable.
1631
1632
All objects created by \kbd{fetch\_var} are on the heap and not on the stack,
1633
thus they are not subject to standard garbage collecting (they are not
1634
destroyed by a \kbd{gerepile} or \kbd{avma = ltop} statement). When you do
1635
not need a variable number anymore, you can delete it using
1636
1637
\fun{long}{delete_var}{}
1638
1639
\noindent which deletes the \emph{latest} temporary variable created and
1640
returns the variable number of the previous one (or simply returns 0 if none
1641
remain). Of course you should make sure that
1642
the deleted variable does not appear anywhere in the objects you use later
1643
on. Here is an example:
1644
1645
\bprog
1646
long first = fetch_var();
1647
long n1 = fetch_var();
1648
long n2 = fetch_var(); /*@Ccom prepare three variables for internal use */
1649
...
1650
/*@Ccom delete all variables before leaving */
1651
do { num = delete_var(); } while (num && num <= first);
1652
@eprog\noindent
1653
The (dangerous) statement
1654
1655
\bprog
1656
while (delete_var()) /*@Ccom empty */;
1657
@eprog\noindent
1658
removes all temporary variables in use.
1659
1660
\subsec{Comparing variables}
1661
1662
Let us go back to \kbd{varncmp}. There is an interesting corner case, when
1663
one of the compared variables (from \kbd{gvar}, say) is \kbd{NO\_VARIABLE}.
1664
In this case, \kbd{varncmp} declares it has lower priority than any other
1665
variable; of course, comparing \kbd{NO\_VARIABLE} with itself yields
1666
$0$ (same priority);
1667
1668
In addition to \kbd{varncmp} we have
1669
1670
\fun{long}{varnmax}{long v, long w} given two variable numbers (possibly
1671
\kbd{NO\_VARIABLE}), returns the variable with the highest priority.
1672
This function always returns a valid variable number unless it is comparing
1673
\kbd{NO\_VARIABLE} to itself.
1674
1675
\fun{long}{varnmin}{long x, long y} given two variable numbers (possibly
1676
\kbd{NO\_VARIABLE}), returns the variable with the lowest priority. Note
1677
that when comparing a true variable with \kbd{NO\_VARIABLE}, this function
1678
returns \kbd{NO\_VARIABLE}, which is not a valid variable number.
1679
1680
\section{Input and output}
1681
1682
\noindent
1683
Two important aspects have not yet been explained which are specific to
1684
library mode: input and output of PARI objects.
1685
1686
\subsec{Input}
1687
1688
\noindent
1689
For \idx{input}, PARI provides several powerful high level functions
1690
which enable you to input your objects as if you were under \kbd{gp}. In fact,
1691
it \emph{is} essentially the GP syntactical parser.
1692
1693
There are two similar functions available to parse a string:
1694
1695
\fun{GEN}{gp_read_str}{const char *s}\label{se:gp_read_str}
1696
1697
\fun{GEN}{gp_read_str_multiline}{const char *s, char *last}
1698
1699
\noindent
1700
Both functions read the whole string \kbd{s}. The function
1701
\kbd{gp\_read\_str} ignores newlines: it assumes that the input is one
1702
expression and returns the result of this expression.
1703
1704
The function \kbd{gp\_read\_str\_multiline} processes the text in the
1705
same way as the GP command \tet{read}: newlines are significant and can
1706
be used to separate expressions.
1707
The return value is that of the last nonempty expression evaluated.
1708
1709
In \kbd{gp\_read\_str\_multiline}, if \kbd{last} is not \kbd{NULL},
1710
then \kbd{*last} receives the last character from the \emph{filtered}
1711
input: this can be used to check if the last character was a semi-colon
1712
(to hide the output in interactive usage). If (and only if) the
1713
input contains no statements, then \kbd{*last} is set to \kbd{0}.
1714
1715
For both functions, \kbd{gp}'s metacommands \emph{are} recognized.
1716
1717
Two variants allow to specify a default precision while evaluating the
1718
string:
1719
1720
\fun{GEN}{gp_read_str_prec}{const char *s, long prec}
1721
As \kbd{gp\_read\_str}, but set the precision to \kbd{prec} words while evaluating $s$.
1722
1723
\fun{GEN}{gp_read_str_bitprec}{const char *s, long bitprec}
1724
As \kbd{gp\_read\_str}, but set the precision to \kbd{bitprec} bits while evaluating $s$.
1725
1726
\misctitle{Note} The obsolete form
1727
1728
\fun{GEN}{readseq}{char *t}
1729
1730
still exists for backward compatibility (assumes filtered input, without
1731
spaces or comments). Don't use it.
1732
1733
To read a \kbd{GEN} from a file, you can use the simpler interface
1734
1735
\fun{GEN}{gp_read_stream}{FILE *file}
1736
1737
\noindent which reads a character string of arbitrary length from the stream
1738
\kbd{file} (up to the first complete expression sequence), applies
1739
\kbd{gp\_read\_str} to it, and returns the resulting \kbd{GEN}. This way, you
1740
do not have to worry about allocating buffers to hold the string. To
1741
interactively input an expression, use \kbd{gp\_read\_stream(stdin)}.
1742
Return \kbd{NULL} when there are no more expressions to read (we reached
1743
EOF).
1744
1745
Finally, you can read in a whole file, as in GP's \tet{read} statement
1746
1747
\fun{GEN}{gp_read_file}{char *name}
1748
1749
\noindent As usual, the return value is that of the last nonempty expression
1750
evaluated. There is one technical exception: if \kbd{name} is a \emph{binary}
1751
file (from \tet{writebin}) containing more than one object, a \typ{VEC}
1752
containing them all is returned. This is because binary objects bypass the
1753
parser, hence reading them has no useful side effect.
1754
1755
\subsec{Output to screen or file, output to string}\sidx{output}
1756
1757
General output functions return nothing but print a character string as a
1758
side effect. Low level routines are available to write on PARI output stream
1759
\tet{pari_outfile} (\tet{stdout} by default):
1760
1761
\fun{void}{pari_putc}{char c}: write character \kbd{c} to the output stream.
1762
1763
\fun{void}{pari_puts}{char *s}: write \kbd{s} to the output stream.
1764
1765
\fun{void}{pari_flush}{}: flush output stream; most streams are buffered by
1766
default, this command makes sure that all characters output so are actually
1767
written.
1768
1769
\fun{void}{pari_printf}{const char *fmt, ...}: the most versatile such
1770
function. \kbd{fmt} is a character string similar to the one
1771
\tet{printf} uses. In there, \kbd{\%} characters have a special meaning, and
1772
describe how to print the remaining operands. In addition to the standard
1773
format types (see the GP function \tet{printf}), you can use the \emph{length
1774
modifier}~\kbd{P} (for PARI of course!) to specify that an argument is a
1775
\kbd{GEN}. For instance, the following are valid conversions for a \kbd{GEN}
1776
argument
1777
\bprog
1778
%Ps @com convert to \kbd{char*} (will print an arbitrary \kbd{GEN})
1779
%P.10s @com convert to \kbd{char*}, truncated to 10 chars
1780
%P.2f @com convert to floating point format with 2 decimals
1781
%P4d @com convert to integer, field width at least 4
1782
1783
pari_printf("x[%d] = %Ps is not invertible!\n", i, gel(x,i));
1784
@eprog\noindent
1785
Here \kbd{i} is an \kbd{int}, \kbd{x} a \kbd{GEN} which is not a leaf
1786
(presumably a vector, or a polynomial) and this would insert the value of its
1787
$i$-th \kbd{GEN} component: \kbd{gel(x,i)}.
1788
1789
\noindent Simple but useful variants to \kbd{pari\_printf} are
1790
1791
\fun{void}{output}{GEN x} prints \kbd{x} in raw format, followed by a
1792
newline and a buffer flush. This is more or less equivalent to
1793
\bprog
1794
pari_printf("%Ps\n", x);
1795
pari_flush();
1796
@eprog
1797
1798
\fun{void}{outmat}{GEN x} as above except if $x$ is a \typ{MAT}, in which
1799
case a multi-line display is used to display the matrix. This is prettier for
1800
small dimensions, but quickly becomes unreadable and cannot be pasted and
1801
reused for input. If all entries of $x$ are small integers, you may use the
1802
recursive features of \kbd{\%Pd} and obtain the same (or better) effect with
1803
\bprog
1804
pari_printf("%Pd\n", x);
1805
pari_flush();
1806
@eprog\noindent
1807
A variant like \kbd{"\%5Pd"} would improve alignment by imposing
1808
5 chars for each coefficient. Similarly if all entries are to be converted to
1809
floats, a format like \kbd{"\%5.1Pf"} could be useful.
1810
1811
1812
These functions write on (PARI's idea of) standard output, and must be used
1813
if you want your functions to interact nicely with \kbd{gp}. In most
1814
programs, this is not a concern and it is more flexible to write to an
1815
explicit \kbd{FILE*}, or to recover a character string:
1816
1817
\fun{void}{pari_fprintf}{FILE *file, const char *fmt, ...} writes the
1818
remaining arguments to stream \kbd{file} according to the format
1819
specification \kbd{fmt}.
1820
1821
\fun{char*}{pari_sprintf}{const char *fmt, ...} produces a string from the
1822
remaining arguments, according to the PARI format \kbd{fmt} (see \tet{printf}).
1823
This is the \kbd{libpari} equivalent of \tet{strprintf}, and returns a
1824
\kbd{malloc}'ed string, which must be freed by the caller. Note that contrary
1825
to the analogous \tet{sprintf} in the \kbd{libc} you do not provide a buffer
1826
(leading to all kinds of buffer overflow concerns); the function provided is
1827
actually closer to the GNU extension \kbd{asprintf}, although the latter has
1828
a different interface.
1829
1830
Simple variants of \tet{pari_sprintf} convert a \kbd{GEN} to a
1831
\kbd{malloc}'ed ASCII string, which you must still \kbd{free} after use:
1832
1833
\fun{char*}{GENtostr}{GEN x}, using the current default output format
1834
(\kbd{prettymat} by default).
1835
1836
\fun{char*}{GENtoTeXstr}{GEN x}, suitable for inclusion in a \TeX\ file.
1837
1838
1839
Note that we have \tet{va_list} analogs of the functions of \kbd{printf} type
1840
seen so far:
1841
1842
\fun{void}{pari_vprintf}{const char *fmt, va_list ap}
1843
1844
\fun{void}{pari_vfprintf}{FILE *file, const char *fmt, va_list ap}
1845
1846
\fun{char*}{pari_vsprintf}{const char *fmt, va_list ap}
1847
1848
\subsec{Errors}\sidx{error}\kbdsidx{e_MISC}
1849
1850
\noindent
1851
If you want your functions to issue error messages, you can use the general
1852
error handling routine \tet{pari_err}. The basic syntax is
1853
%
1854
\bprog
1855
pari_err(e_MISC, "error message");
1856
@eprog\noindent
1857
This prints the corresponding error message and exit the program (in
1858
library mode; go back to the \kbd{gp} prompt otherwise).\label{se:err} You can
1859
also use it in the more versatile guise
1860
\bprog
1861
pari_err(e_MISC, format, ...);
1862
@eprog\noindent
1863
where \kbd{format} describes the format to use to write the remaining
1864
operands, as in the \tet{pari_printf} function. For instance:
1865
\bprog
1866
pari_err(e_MISC, "x[%d] = %Ps is not invertible!", i, gel(x,i));
1867
@eprog\noindent
1868
The simple syntax seen above is just a special case with a constant format
1869
and no remaining arguments. The general syntax is
1870
1871
\fun{void}{pari_err}{numerr,...}
1872
1873
\noindent where \kbd{numerr} is a codeword which specifies the error class
1874
and what to do with the remaining arguments and what message to print.
1875
For instance, if $x$ is a \kbd{GEN} with internal type \typ{STR}, say,
1876
\kbd{pari\_err(e\_TYPE,"extgcd", $x$)} prints the message:
1877
\bprog
1878
*** incorrect type in extgcd (t_STR),
1879
@eprog\noindent See \secref{se:errors} for details. In the libpari code
1880
itself, the general-purpose \kbd{e\_MISC} is used sparingly: it is so
1881
flexible that the corresponding error contexts (\typ{ERROR}) become hard to
1882
use reliably. Other more rigid error types are generally more useful: for
1883
instance the error context attached to the \kbd{e\_TYPE} exception above is
1884
precisely documented and contains \kbd{"extgcd"} and $x$ (not only its type)
1885
as readily available components.
1886
1887
\subsec{Warnings}
1888
1889
\noindent To issue a warning, use
1890
1891
\fun{void}{pari_warn}{warnerr,...}
1892
In that case, of course, we do \emph{not} abort the computation, just print
1893
the requested message and go on. The basic example is
1894
%
1895
\bprog
1896
pari_warn(warner, "Strategy 1 failed. Trying strategy 2")
1897
@eprog\noindent
1898
which is the exact equivalent of \kbd{pari\_err(e\_MISC,...)} except that
1899
you certainly do not want to stop the program at this point, just inform the
1900
user that something important has occurred; in particular, this output would be
1901
suitably highlighted under \kbd{gp}, whereas a simple \kbd{printf} would not.
1902
1903
The valid \emph{warning} keywords are \tet{warner} (general), \tet{warnprec}
1904
(increasing precision), \tet{warnmem} (garbage collecting) and \tet{warnfile}
1905
(error in file operation), used as follows:
1906
\bprog
1907
pari_warn(warnprec, "bnfinit", newprec);
1908
pari_warn(warnmem, "bnfinit");
1909
pari_warn(warnfile, "close", "afile"); /* error when closing "afile" */
1910
@eprog
1911
1912
\subsec{Debugging output}\sidx{debugging}\sidx{format}\label{se:dbg_output}
1913
1914
For debugging output, you can use the standard output
1915
functions, \tet{output} and \tet{pari_printf} mainly. Corresponding to the
1916
\kbd{gp} metacommand \kbd{\b x}, you can also output the \idx{hexadecimal
1917
tree} attached to an object:
1918
1919
\fun{void}{dbgGEN}{GEN x, long nb = -1}, displays the recursive structure of
1920
\kbd{x}. If $\kbd{nb} = -1$, the full structure is printed, otherwise
1921
the leaves (nonrecursive components) are truncated to \kbd{nb} words.
1922
1923
\noindent The function \tet{output} is vital under debuggers, since none of
1924
them knows how to print PARI objects by default. Seasoned PARI developers
1925
add the following \kbd{gdb} macro to their \kbd{.gdbinit}:
1926
\bprog
1927
define oo
1928
call output((GEN)$arg0)
1929
end
1930
define xx
1931
call dbgGEN($arg0,-1)
1932
end
1933
@eprog\noindent
1934
Typing \kbd{i x} at a breakpoint in \kbd{gdb} then prints the value of the
1935
\kbd{GEN} \kbd{x} (provided the optimizer has not put it into a register, but
1936
it is rarely a good idea to debug optimized code).
1937
1938
\noindent
1939
The global variables \teb{DEBUGLEVEL} and \teb{DEBUGMEM} (corresponding
1940
to the default \teb{debug} and \teb{debugmem})
1941
are used throughout the PARI code to govern the amount of diagnostic and
1942
debugging output, depending on their values. You can use them to debug your
1943
own functions, especially if you \tet{install} the latter under \kbd{gp}.
1944
1945
\fun{void}{dbg_pari_heap}{void} print debugging statements about the PARI
1946
stack, heap, and number of variables used. Corresponds to \kbd{\bs s}
1947
under gp.
1948
1949
\subsec{Timers and timing output}
1950
1951
\noindent To handle timings in a reentrant way, PARI defines a dedicated data
1952
type, \tet{pari_timer}, together with the following methods:
1953
1954
\fun{void}{timer_start}{pari_timer *T} start (or reset) a timer.
1955
1956
\fun{long}{timer_delay}{pari_timer *T} returns the number of milliseconds
1957
elapsed since the timer was last reset. Resets the timer as a side effect.
1958
Assume $T$ was started by \kbd{timer\_start}.
1959
1960
\fun{long}{timer_get}{pari_timer *T} returns the number of milliseconds
1961
elapsed since the timer was last reset. Does \emph{not} reset the timer.
1962
Assume $T$ was started by \kbd{timer\_start}.
1963
1964
\fun{void}{walltimer_start}{pari_timer *T} start a timer, as if it
1965
had been started at the Unix epoch (see \kbd{getwalltime}).
1966
1967
\fun{long}{walltimer_delay}{pari_timer *T} returns the number of milliseconds
1968
elapsed since the timer was last checked. Assume $T$ was started by
1969
\kbd{walltimer\_start}.
1970
1971
\fun{long}{walltimer_get}{pari_timer *T} returns the number of milliseconds
1972
elapsed since the timer was last reset. Does \emph{not} reset the timer.
1973
Assume $T$ was started by \kbd{walltimer\_start}.
1974
1975
\fun{long}{timer_printf}{pari_timer *T, char *format,...} This diagnostics
1976
function is equivalent to the following code
1977
\bprog
1978
err_printf("Time ")
1979
... prints remaining arguments according to format ...
1980
err_printf(": %ld", timer_delay(T));
1981
@eprog\noindent Resets the timer as a side effect.
1982
1983
\noindent They are used as follows:
1984
\bprog
1985
pari_timer T;
1986
timer_start(&T); /* initialize timer */
1987
...
1988
printf("Total time: %ldms\n", timer_delay(&T));
1989
@eprog\noindent
1990
or
1991
\bprog
1992
pari_timer T;
1993
timer_start(&T);
1994
for (i = 1; i < 10; i++) {
1995
...
1996
timer_printf(&T, "for i = %ld (L[i] = %Ps)", i, gel(L,i));
1997
}
1998
@eprog
1999
2000
The following functions provided the same functionality, in a
2001
nonreentrant way, and are now deprecated.
2002
2003
\fun{long}{timer}{void}
2004
2005
\fun{long}{timer2}{void}
2006
2007
\fun{void}{msgtimer}{const char *format, ...}
2008
2009
The following function implements \kbd{gp}'s timer and should not be used in
2010
libpari programs:
2011
\fun{long}{gettime}{void} equivalent to \tet{timer_delay}$(T)$ attached to a
2012
private timer $T$.
2013
2014
\section{Iterators, Numerical integration, Sums, Products}
2015
\subsec{Iterators}
2016
Since it is easier to program directly simple loops in library mode, some GP
2017
iterators are mainly useful for GP programming. Here are the others:
2018
2019
\item \tet{fordiv} is a trivial iteration over a list produced by
2020
\tet{divisors}.
2021
2022
\item \tet{forell}, \tet{forqfvec} and \tet{forsubgroup} are currently not
2023
implemented as an iterator but as a procedure with callbacks.
2024
2025
\fun{void}{forell}{void *E, long fun(void*, GEN), GEN a, GEN b, long flag}
2026
goes through the same curves as \tet{forell(ell,a,b,,flag)}, calling
2027
\tet{fun(E, ell)} for each curve \kbd{ell}, stopping if \kbd{fun} returns a
2028
nonzero value.
2029
2030
\fun{void}{forqfvec}{void *E, long (*fun)(void *, GEN, GEN, double), GEN q, GEN b}:
2031
Evaluate \kbd{fun(E,U,v,m)} on all $v$ such that $q(U\*v)<b$, where $U$ is a
2032
\typ{MAT}, $v$ is a \typ{VECSMALL} and $m=q(v)$ is a C double. The function
2033
\kbd{fun} must return $0$, unless \kbd{forqfvec} should stop, in which case,
2034
it should return $1$.
2035
2036
\fun{void}{forqfvec1}{void *E, long (*fun)(void *, GEN), GEN q, GEN b}:
2037
Evaluate \kbd{fun(E,v)} on all $v$ such that $q(v)<b$, where $v$ is a
2038
\typ{COL}. The function \kbd{fun} must return $0$, unless \kbd{forqfvec}
2039
should stop, in which case, it should return $1$.
2040
2041
\fun{void}{forsubgroup}{void *E, long fun(void*, GEN), GEN G, GEN B}
2042
goes through the same subgroups as \tet{forsubgroup(H = G, B,)}, calling
2043
\tet{fun(E, H)} for each subgroup $H$, stopping if \kbd{fun} returns a
2044
nonzero value.
2045
2046
\item \tet{forprime} and \tet{forprimestep}, iterators over primes and
2047
primes in arithmetic progressions, for which we refer you to the
2048
next subsection.
2049
2050
\item \tet{forcomposite}, we provide an iterator over composite integers:
2051
2052
\fun{int}{forcomposite_init}{forcomposite_t *T, GEN a, GEN b} initialize an
2053
iterator $T$ over composite integers in $[a,b]$; over composites $\geq a$ if
2054
$b = \kbd{NULL}$. We must have $a\geq 0$. Return $0$ if the range is known to
2055
be empty from the start (as if $b < a$ or $b < 0$), and return $1$ otherwise.
2056
2057
\fun{GEN}{forcomposite_next}{forcomposite_t *T} returns the next composite in
2058
the range, assuming that $T$ was initialized by \tet{forcomposite_init}.
2059
2060
\item \tet{forvec}, for which we provide a convenient iterator. To
2061
initialize the analog of \kbd{forvec(X = v, ..., flag)}, call
2062
2063
\fun{int}{forvec_init}{forvec_t *T, GEN v, long flag}
2064
initialize an iterator $T$ over the vectors generated by
2065
\kbd{forvec(X = $v$,..., flag)}. This returns $0$ if this vector list is
2066
empty, and $1$ otherwise.
2067
2068
\fun{GEN}{forvec_next}{forvec_t *T} returns the next element in the
2069
\kbd{forvec} sequence, or \kbd{NULL} if we are done. The return value must be
2070
used immediately or copied since the next call to the iterator destroys it:
2071
the relevant vector is updated in place. The iterator works hard to not
2072
use up PARI stack, and is more efficient when all lower bounds in the
2073
initialization vector $v$ are integers. In that case, the cost is linear in
2074
the number of tuples enumerated, and you can expect to run over more than
2075
$10^9$ tuples per minute. If speed is critical and all integers involved
2076
would fit in $C$ \kbd{long}s, write a simple direct backtracking algorithm
2077
yourself.
2078
2079
\item \tet{forpart} is a variant of \kbd{forvec} which iterates over
2080
partitions. See the documentation of the \kbd{forpart} GP function for
2081
details. This function is available as a loop with callbacks:
2082
2083
\fun{void}{forpart}{void *data, long (*call)(void*,GEN), long k, GEN a, GEN n}
2084
2085
\noindent It is also available as an iterator:
2086
2087
\fun{void}{forpart_init}{forpart_t *T, long k, GEN a, GEN n} initializes an
2088
iterator over the partitions of $k$, with length restricted by $n$,
2089
and components restricted by $a$, either of which can be set to \kbd{NULL}
2090
to run without restriction.
2091
2092
\fun{GEN}{forpart_next}{forpart_t *T} returns the next partition, or
2093
\kbd{NULL} when all partitions have been exhausted.
2094
2095
\fun{GEN}{forpart_prev}{forpart_t *T} returns the previous partition, or
2096
\kbd{NULL} when all partitions have been exhausted.
2097
2098
In both cases, the partition must be used or copied before the next call
2099
since it is returned from a state array which will be modified in place.
2100
You may \emph{not} mix calls to \tet{forpart_next} and \tet{forpart_prev}:
2101
the first one called determines the ordering used to iterate over the
2102
partitions; you can not go back since the \tet{forpart_t} structure is used
2103
in incompatible ways.
2104
2105
\item \tet{forperm} to loop over permutations of $k$. See the documentation
2106
of the \kbd{forperm} GP function for details. This function is available as
2107
an iterator:
2108
2109
\fun{void}{forperm_init}{forperm_t *T, GEN k} initializes an iterator
2110
over the permutations of $k$ (\typ{INT}, \typ{VEC} or \typ{VECSMALL}).
2111
2112
\fun{GEN}{forperm_next}{forperm_t *T} returns the next permutation
2113
as a \typ{VECSMALL} or \kbd{NULL} whell all permutations have been
2114
exhausted. The permutation must be used or copied before the next call
2115
since it is returned from a state array which will be modified in place.
2116
2117
\item \tet{forsubset} to loop over subsets. See the documentation
2118
of the \kbd{forsubset} GP function for details. This function
2119
is available as two iterators:
2120
2121
\fun{void}{forallsubset_init}{forsubset_t *T, long n}
2122
2123
\fun{void}{forksubset_init}{forsubset_t *T, long n, long k}
2124
2125
\noindent It is also available in generic form:
2126
2127
\fun{void}{forsubset_init}{forsubset_t *T, GEN nk} where \kbd{nk} is either
2128
a \typ{INT} $n$ or a \typ{VEC} with two integral components $[n,k]$.
2129
2130
In all three cases, \fun{GEN}{forsubset_next}{forsubset_t *T} returns the
2131
next subset as a \typ{VECSMALL} or \kbd{NULL} when all subsets have been
2132
exhausted.
2133
2134
\subsec{Iterating over primes}\label{se:primeiter}
2135
2136
The library provides a high-level iterator, which stores its (private) data
2137
in a \kbd{struct} \tet{forprime_t} and runs over arbitrary ranges of primes,
2138
without ever overflowing.
2139
2140
The iterator has two flavors, one providing the successive primes as
2141
\kbd{ulong}s, the other as \kbd{GEN}. They are initialized as follows,
2142
where we expect to run over primes $\geq a$ and $\leq b$:
2143
2144
\fun{int}{u_forprime_init}{forprime_t *T, ulong a, ulong b} for the
2145
\kbd{ulong} variant, where $b = \kbd{ULONG\_MAX}$ means we will run through
2146
all primes representable in a \kbd{ulong} type.
2147
2148
\fun{int}{forprime_init}{forprime_t *T, GEN a, GEN b} for the \kbd{GEN}
2149
variant, where $b = \kbd{NULL}$ means $+\infty$.
2150
2151
\fun{int}{forprimestep_init}{forprime_t *T, GEN a, GEN b, GEN q} initialize
2152
an iterator $T$ over primes in an arithmetic progression, $p \geq a$ and $p
2153
\leq b$ (where $b = \kbd{NULL}$ means $+\infty$). The argument $q$ is
2154
either a \typ{INT} ($p \equiv a \pmod{q}$) or a \typ{INTMOD} \kbd{Mod(c,N)}
2155
and we restrict to that congruence class.
2156
2157
All variants return $1$ on success, and $0$ if the iterator would run over an
2158
empty interval (if $a > b$, for instance). They allocate the \tet{forprime_t}
2159
data structure on the PARI stack.
2160
2161
\noindent The successive primes are then obtained using
2162
2163
\fun{GEN}{forprime_next}{forprime_t *T}, returns \kbd{NULL} if no more primes
2164
are available in the interval and the next suitable prime as a \typ{INT}
2165
otherwise.
2166
2167
\fun{ulong}{u_forprime_next}{forprime_t *T}, returns $0$ if no more primes
2168
are available in the interval and fitting in an \kbd{ulong} and the next
2169
suitable prime otherwise.
2170
2171
These two functions leave alone the PARI stack, and write their state
2172
information in the preallocated \tet{forprime_t} struct. The typical usage is
2173
thus:
2174
\bprog
2175
forprime_t T;
2176
GEN p;
2177
pari_sp av = avma, av2;
2178
2179
forprime_init(&T, gen_2, stoi(1000));
2180
av2 = avma;
2181
while ( (p = forprime_next(&T)) )
2182
{
2183
...
2184
if ( prime_is_OK(p) ) break;
2185
avma = av2; /* delete garbage accumulated in this iteration */
2186
}
2187
avma = av; /* delete all */
2188
@eprog\noindent Of course, the final \kbd{avma = av} could be replaced
2189
by a \kbd{gerepile} call. Beware that swapping the
2190
\kbd{av2 = avma} and \tet{forprime_init} call would be incorrect: the
2191
first \kbd{avma = av2} would delete the \tet{forprime_t} structure!
2192
2193
\subsec{Parallel iterators}
2194
2195
Theses iterators loops over the values of a \typ{CLOSURE}
2196
taken at some data, where the evaluations are done in parallel.
2197
2198
\item \tet{parfor}. To initialize the analog of
2199
\kbd{parfor(i = a, b, ...)}, call
2200
2201
\fun{void}{parfor_init}{parfor_t *T, GEN a, GEN b, GEN code}
2202
initialize an iterator over the evaluation of \kbd{code} on the integers
2203
between $a$ and $b$.
2204
2205
\fun{GEN}{parfor_next}{parfor_t *T} returns a \typ{VEC} \kbd{[i,code(i)]} where
2206
$i$ is one of the integers and \kbd{code(i)} is the evaluation, \kbd{NULL} when
2207
all data have been exhausted. Once it happens, \kbd{parfor\_next} must not be
2208
called anymore with the same initialization.
2209
2210
\fun{void}{parfor_stop}{parfor_t *T} needs to be called when leaving the
2211
iterator before \kbd{parfor\_next} returned \kbd{NULL}.
2212
2213
The following returns an integer $1\leq i\leq N$ such that \kbd{fun(i)} is not
2214
zero, or \kbd{NULL}.
2215
\bprog
2216
GEN
2217
parfirst(GEN fun, GEN N)
2218
{
2219
parfor_t T;
2220
GEN e;
2221
parfor_init(&T, gen_1, N, fun);
2222
while ((e = parfor_next(&T)))
2223
{
2224
GEN i = gel(e,1), funi = gel(e,2);
2225
if (!gequal0(funi))
2226
{ /* found: stop the iterator and return the index */
2227
parfor_stop(&T);
2228
return i;
2229
}
2230
}
2231
return NULL; /* not found */
2232
}
2233
@eprog
2234
2235
\item \tet{parforeach}. To initialize the analog of
2236
\kbd{parforeach(V, X, ...)}, call
2237
2238
\fun{void}{parforeach_init}{parforeach_t *T, GEN V, GEN code}
2239
initialize an iterator over the evaluation of \kbd{code} on the components
2240
of $V$.
2241
2242
\fun{GEN}{parforeach_next}{parforeach_t *T} returns a \typ{VEC}
2243
\kbd{[V[i],code(V[i])]} where $V[i]$ is one of the components of $V$ and \kbd{code(V[i])} is
2244
the evaluation, \kbd{NULL} when all data have been exhausted. Once it happens,
2245
\kbd{parforprime\_next}
2246
must not be called anymore with the same initialization.
2247
2248
\fun{void}{parforeach_stop}{parforeach_t *T} needs to be called when leaving
2249
the iterator before \kbd{parforeach\_next} returned \kbd{NULL}.
2250
2251
\item \tet{parforprime}. To initialize the analog of
2252
\kbd{parforprime(p = a, b, ...)}, call
2253
2254
\fun{void}{parforprime_init}{parforprime_t *T, GEN a, GEN b, GEN code}
2255
initialize an iterator over the evaluation of \kbd{code} on the prime numbers
2256
between $a$ and $b$.
2257
2258
\item \tet{parforprimestep}. To initialize the analog of
2259
\kbd{parforprimestep(p = a, b, q, ...)}, call
2260
2261
\fun{void}{parforprimestep_init}{parforprime_t *T, GEN a, GEN b, GEN q, GEN code}
2262
initialize an iterator over the evaluation of \kbd{code} on the prime numbers
2263
between $a$ and $b$ in the congruence class defined by $q$.
2264
2265
\fun{GEN}{parforprime_next}{parforprime_t *T} returns a \typ{VEC}
2266
\kbd{[p,code(p)]} where $p$ is one of the prime numbers and \kbd{code(p)} is
2267
the evaluation, \kbd{NULL} when all data have been exhausted. Once it happens,
2268
\kbd{parforprime\_next}
2269
must not be called anymore with the same initialization.
2270
2271
\fun{void}{parforprime_stop}{parforprime_t *T} needs to be called when leaving
2272
the iterator before \kbd{parforprime\_next} returned \kbd{NULL}.
2273
2274
\item \tet{parforvec}. To initialize the analog of
2275
\kbd{parforvec(X = V, ..., flag)}, call
2276
2277
\fun{void}{parforvec_init}{parforvec_t *T, GEN V, GEN code, long flag}
2278
initialize an iterator over the evaluation of \kbd{code} on the vectors
2279
specified by \kbd{V} and \kbd{flag}, see \kbd{forvec} for detail.
2280
2281
\fun{GEN}{parforvec_next}{parforvec_t *T} returns a \typ{VEC} \kbd{[v,code(v)]}
2282
where $v$ is one of the vectors and \kbd{code(v)} is the evaluation, \kbd{NULL}
2283
when all data have been exhausted. Once it happens, \kbd{parforvec\_next}
2284
must not be called anymore with the same initialization.
2285
2286
\fun{void}{parforvec_stop}{parforvec_t *T} needs to be called when leaving the
2287
iterator before \kbd{parforvec\_next} returned \kbd{NULL}.
2288
2289
\subsec{Numerical analysis}
2290
2291
Numerical routines code a function (to be integrated, summed, zeroed, etc.)
2292
with two parameters named
2293
\bprog
2294
void *E;
2295
GEN (*eval)(void*, GEN)
2296
@eprog\noindent
2297
The second is meant to contain all auxiliary data needed by your function.
2298
The first is such that \kbd{eval(x, E)} returns your function evaluated at
2299
\kbd{x}. For instance, one may code the family of functions
2300
$f_t: x \to (x+t)^2$ via
2301
\bprog
2302
GEN fun(void *t, GEN x) { return gsqr(gadd(x, (GEN)t)); }
2303
@eprog\noindent
2304
One can then integrate $f_1$ between $a$ and $b$ with the call
2305
\bprog
2306
intnum((void*)stoi(1), &fun, a, b, NULL, prec);
2307
@eprog\noindent
2308
Since you can set \kbd{E} to a pointer to any \kbd{struct} (typecast to
2309
\kbd{void*}) the above mechanism handles arbitrary functions. For simple
2310
functions without extra parameters, you may set \kbd{E = NULL} and ignore
2311
that argument in your function definition.
2312
2313
\section{Catching exceptions}
2314
2315
\subsec{Basic use}
2316
2317
PARI provides a mechanism to trap exceptions generated via \kbd{pari\_err}
2318
using the \tet{pari_CATCH} construction. The basic usage is as follows
2319
\bprog
2320
pari_CATCH(err_code) {
2321
@com recovery branch
2322
}
2323
pari_TRY {
2324
@com main branch
2325
}
2326
pari_ENDCATCH
2327
@eprog\noindent This fragment executes the main branch, then the recovery
2328
branch \emph{if} exception \kbd{err\_code} is thrown, e.g. \kbd{e\_TYPE}.
2329
See \secref{se:errors} for the description of all error classes.
2330
The special error code \tet{CATCH_ALL} is available to catch all errors.
2331
2332
One can replace the \tet{pari_TRY} keyword by \tet{pari_RETRY}, in which case
2333
once the recovery branch is run, we run the main branch again, still catching
2334
the same exceptions.
2335
2336
\misctitle{Restrictions}
2337
2338
\item Such constructs can be nested without adverse effect, the innermost
2339
handler catching the exception.
2340
2341
\item It is \emph{valid} to leave either branch using \tet{pari_err}.
2342
2343
\item It is \emph{invalid} to use C flow control instructions (\kbd{break},
2344
\kbd{continue}, \kbd{return}) to directly leave either branch without seeing
2345
the \tet{pari_ENDCATCH} keyword. This would leave an invalid structure in the
2346
exception handler stack, and the next exception would crash.
2347
2348
\item In order to leave using \kbd{break}, \kbd{continue} or \kbd{return},
2349
one must precede the keyword by a call to
2350
2351
\fun{void}{pari_CATCH_reset}{} disable the current handler, allowing to leave
2352
without adverse effect.
2353
2354
\subsec{Advanced use}
2355
2356
In the recovery branch, the exception context can be examined via the
2357
following helper routines:
2358
2359
\fun{GEN}{pari_err_last}{} returns the exception context, as a \typ{ERROR}.
2360
The exception $E$ returned by \tet{pari_err_last} can be rethrown, using
2361
\bprog
2362
pari_err(0, E);
2363
@eprog
2364
2365
\fun{long}{err_get_num}{GEN E} returns the error symbolic name. E.g
2366
\kbd{e\_TYPE}.
2367
2368
\fun{GEN}{err_get_compo}{GEN E, long i} error $i$-th component, as documented
2369
in \secref{se:errors}.
2370
2371
\noindent For instance
2372
\bprog
2373
pari_CATCH(CATCH_ALL) { /* catch everything */
2374
GEN x, E = pari_err_last();
2375
long code = err_get_num(E);
2376
if (code != e_INV) pari_err(0, E); /* unexpected error, rethrow */
2377
x = err_get_compo(E, 2);
2378
/* e_INV has two components, 1: function name 2: noninvertible x */
2379
if (typ(x) != t_INTMOD) pari_err(0, E); /* unexpected type, rethrow */
2380
pari_CATCH_reset();
2381
return x; /* leave ! */
2382
@com @dots
2383
} pari_TRY {
2384
@com main branch
2385
}
2386
pari_ENDCATCH
2387
@eprog
2388
2389
\section{A complete program}
2390
\label{se:prog}
2391
2392
\noindent
2393
Now that the preliminaries are out of the way, the best way to learn how to
2394
use the library mode is to study a detailed example. We want to write a
2395
program which computes the gcd of two integers, together with the Bezout
2396
coefficients. We shall use the standard quadratic algorithm which is not
2397
optimal but is not too far from the one used in the PARI function
2398
\teb{bezout}.
2399
2400
Let $x,y$ two integers and initially
2401
$ \pmatrix{s_x & s_y \cr t_x & t_y } =
2402
\pmatrix{1 & 0 \cr 0 & 1}$, so that
2403
$$ \pmatrix{s_x & s_y \cr
2404
t_x & t_y }
2405
\pmatrix{x \cr y } =
2406
\pmatrix{x \cr y }.
2407
$$
2408
To apply the ordinary Euclidean algorithm to the right hand side,
2409
multiply the system from the left by
2410
$ \pmatrix{0 & 1 \cr 1 & -q }$,
2411
with $q = \kbd{floor}(x / y)$. Iterate until $y = 0$ in the right hand side,
2412
then the first line of the system reads
2413
$$ s_x x + s_y y = \gcd(x,y).$$
2414
In practice, there is no need to update $s_y$ and $t_y$ since
2415
$\gcd(x,y)$ and $s_x$ are enough to recover $s_y$. The following program
2416
is now straightforward. A couple of new functions appear in there, whose
2417
description can be found in the technical reference manual in Chapter 5,
2418
but whose meaning should be clear from their name and the context.
2419
2420
This program can be found in \kbd{examples/extgcd.c} together with a
2421
proper \kbd{Makefile}. You may ignore the first comment
2422
\bprog
2423
/*
2424
GP;install("extgcd", "GG&&", "gcdex", "./libextgcd.so");
2425
*/
2426
@eprog\noindent which instruments the program so that \kbd{gp2c-run extgcd.c}
2427
can import the \kbd{extgcd()} routine into an instance of the \kbd{gp}
2428
interpreter (under the name \kbd{gcdex}). See the \kbd{gp2c} manual for
2429
details.
2430
\newpage
2431
2432
\bprogfile{../examples/extgcd.c}
2433
2434
\noindent For simplicity, the inner loop does not include any garbage
2435
collection, hence memory use is quadratic in the size of the inputs instead
2436
of linear. Here is a better version of that loop:
2437
\bprog
2438
pari_sp av = avma;
2439
...
2440
while (!gequal0(b))
2441
{
2442
GEN r, q = dvmdii(a, b, &r), v = vx;
2443
2444
vx = subii(ux, mulii(q, vx));
2445
ux = v; a = b; b = r;
2446
if (gc_needed(av,1))
2447
gerepileall(av, 4, &a, &b, &ux, &vx);
2448
}
2449
@eprog
2450
\newpage
2451
2452