Pari internal

KIMURA Iwao at Toyama University 13th October 2000

Purpose of this talk

• To encourage writing programs in C, C++ with Pari library

• To explain Pari’s type system

• To explain how Pari memory management system is designed, is implemented and works

�

Purpose of this talk

�

� �

• Be happy with Pari.

About Pari-gp

Pari-gp is created by Henri Cohen on mid ’80s.

The current version of Pari-gp is maintained by C. Batut, K. Belabas, D. Bernardi, H. Cohen and M. Olivier.

Development process is a kind of Open Source model. Everybody can get entire source code, can submit bug fix, enhancements.

Capability

• multiple precision arithmetic, e.g. integers, rationales, finite fields, real numbers, padic numbers, polynomials of coefficients in such numbers, and algebraic extension (polynomials modulo polynomials) of those rings/fields

• elementary transcendental functions (sin, cos, exp, log,...)

• not so elementary transcendental functions (Riemann/Dedekind zeta, Weber’s function, polylog,...)

Capability (continued)

• Pari-gp has many functions to perform

– calculation about (absolute/relative) number fields. e.g. class numbers/groups, regulators/unit groups, discriminants, integral basis, class field theoretic functions...etc.

– calculation about elliptic curves over rationals/finite fields OK, let’s write nice programs in C with Pari library, but HOW?

You need...

• some knowledge about C

• good development environment, e.g. editor, C compiler (and make), debugger, cross referencer, source code beautifier...

You need...

� (We won’t enter

this topic) �

Emacs, gcc, gmake, global, GNU indent,...

You need...(continued)

• nice algorithm

• nice data structures

• Some Knowledge about PARI’s – Type system

– library functions (those parameters, return values...)

– Memory management

Tiny example 1

#include <stdio.h>
#include <pari.h>
int
main (void)
{
GEN i, j, k, x, y, z;
pari_init (500000, 2);
i = gun;
j = stoi (3);
k = gadd (i, j);
x = dbltor (.0);
y = dbltor (.1);
z = gadd (x, y);
printf ("1 + 3 = %s, 0.0 + 0.1 = %s\n",
        GENtostr (k), GENtostr (z));
return 0;
}

Tiny example 1 (result)

$ gcc -o tinyexample1 tinyexample1.c \
`pari-config --cflags` `pari-config --libs`
$ ./tinyexample1
1 + 3 = 4, 0.0 + 0.1 = 0.100000000000000005551

Pari’s type system

There is only one type GEN in Pari from C point of view.

But Pari has its own type system. This can be divided into two parts. One corresponds to numbers, for example, t INT, t INTMOD, t REAL, t COMPLEX, t PADIC,...

Another is considered as container. e.g. t VEC, t MAT, t POLY, t POLYMOD...

GEN internal

GEN is long ∗, or, array of long (abuse of notion).

If you declare z as GEN and initialize z correctly, z[0] contains some information of z (its type in Pari type system, length, clone bit).

Assignments

There are 3 kind of assignments in Pari. Let x, y be variables of type GEN and assume x represents some mathematical entity.

C level y = x; is legal. This is assignment of pointer.

gaffect() gaffect (x, y); is legal if y is properly initialized and typ(x) == typ (y).

gcopy() y = gcopy (x); is always legal. Copy of x is created on the Pari stack, y points it.

How one can construct and use numbers in Pari?

• There are universal constants in Pari. e.g. gzero, gun, gdeux, ghalf, geuler, gpi, gi,...

• Convert C numbers to Pari numbers. stoi() converts C long integers to Pari’s t INT, dbltor() converts C double numbers to t REAL,...

How one can use vectors, matrices,...?

There are two way to use such containers.

• To use flisseq(char ∗s), flisexpr(char ∗s). These are essentially gp. To make a matrix,

    GEN m = flisexpr ("[1,0;0,1]");

is legal at all.

• Another way is to construct vectors and matrices by hand. To do this, It is mandatory to know about Pari’s memory management system.

How one can know type of GEN object?

In C point of view, every object in Pari is of type GEN. But one wants to know, for example, is certain GEN object of type t INT? t REAL?...

There are functions by which we can inspect GEN object. typ(GEN z) returns its type. lg(GEN z) returns its length. For each Pari type, there are appropriate inspection functions.

Pari stack

Every program with Pari library must call pari init() in its early stage.

Pari’s memory management system is as follows. First of all, allocate huge chunk of memory as scratch-pad and place holder of temporary results. This is called Pari stack (Though it is allocated by malloc(3), so resides on heap).

�

�

� There are some functions which allocate some memory from Pari stack. One who wants to construct a matrix must, so,

• Create each entries in Pari stack

• Allocate memory which is nothing but matrix itself

• Insert each entries into the matrix

#include <stdio.h>
#include <pari.h>
int
main (void)
{
GEN v1, v2, mat1, mat2;
pari_init (500000, 2);
v1 = cgetg (3, t_COL);
v2 = cgetg (3, t_COL);
v1[1] = gzero;
v1[2] = gneg (gun);
v2[1] = gun;
v2[2] = gzero;
mat1 = cgetg (3, t_MAT);
mat1[1] = v1;
mat1[2] = v2;
mat2 = gmul (mat1, mat1);

printf ("mat1 = %s, mat2 = %s\n",
        GENtostr (mat1), GENtostr (mat2));
return 0;
}

Tiny example 2 (continued)

$ ./tinyexample2
mat1 = [0, 1; -1, 0], mat2 = [-1, 0; 0, -1]

You will see some warnings when you compile this example. Ignore it.

Pari’s built-in functions

Most built-in functions takes some GEN’s and long as arguments. Those return values are, in most cases, GEN or long. This means, most built-in functions allocate some memory in Pari stack, and returns its address (== GEN).

The functions which return long don’t make any object on Pari stack.

Pari stack

(top−bot)/sizeof(long) is the 1st argument of pari init(). Avma stands for AVailable Memory Adress.

As computations continue, the Pari stack filled with temporary results and garbage. You must clean up Pari stack, otherwise error message,

Pari stack overflows!!!

will be raised (and all cumputation will be stopped).

GC, the easiest case

If you can abandon temporary results since certain points upto now, you are lucky.

/*  Some computation has done.
 *  You can abandon computation from now.
 */
    ltop = avma;
    somecomputation();
    makesgarbage();
    onparistack();
/*  now you want to collect garbage */
    avma = ltop;

Garbages are simply over-written.

GC, using cgiv(GEN).

If you can give back memory which is the most recently allocated, use cgiv(GEN).

GEN z = somefunction();
outbeaut (z);
/* now z is needless for you. */
cgiv (z);

But Life is not so simple!

Sometime your Pari stack is mixed up by garbage and some meaningful results...

If you can certify some interval of Pari stack is garbage, use gerepile(long ltop, long lbot, GEN z). gerepile() copies your important temporary data to somewhere, clean up Pari stack, and copies them back. gerepile(), before and after.

29

coding example

ltop = avma;
somefunc();
makesgarbage();
onmystack();
lbot = avma;
z = important();
z = gerepile (ltop, lbot, z);

Smart GC

ltop = avma;
limit = (avma + bot) / 2;
/* smart gc, taken from 4.4.2 of the Bible */
  if (avma < limit)
    {
      /* half of the stack is spent */
      lbot = avma;
      from = gerepile (ltop, lbot, gcopy (from));
      err (warner, "gerepile done! (avma = %ld)",
           avma);
    }

If you cannot predict when is the time for GC, you may be helped by coding in this way.

gerepile()’s family

gerepile() has some variants. gerepileupto (lbot,

z) is synonim to gerepile (top, lbot, z).

gerepilemany(), gerepilemanysp() are gerepiling function which update many significant objects.

Output your computation

char ∗ GENtostr (GEN z) stringify GEN object appropriately.

void outbeaut (GEN z) outputs GEN object to stdout with newline, buffer flush.

void output (GEN z) also outputs GEN object. This is useful when one traces a program using debugger.

Put your function into gp

You can put your function into gp, and use the function from gp script.

But this needs some knowledge about parser code (Instruction for gp parser what parameters the function takes and what kind of value are returned).

Please refer User’s guide, §4.9.2.

Tips

• To migrate your gp code to C program with Pari library, the easiest way is to use lisseq(char ∗ s) with gp program ‘s’, e.g. :

#include <stdio.h>
#include <pari.h>
int
main (void)
{
  GEN z;
  pari_init (800000, 5000);
  z = lisseq ("f=polsubcyclo(9,3);
  g=polcompositum(x^2+23,f);bnfclgp(g[1]);");
  outbeaut (z);
  return 0;
}

Tips (continued)

Returns following output (Class group of first layer of Z_{3}-extension of Q(√−23 ).

[9, [9], [
[6 0 0 4 0 2]
[0 2 0 0 0 0]
[0 0 2 0 0 0]
[0 0 0 1 0 0]
[0 0 0 0 1 0]
[0 0 0 0 0 1]
]]

Tips (Continued)

Use assert()! In C level, there is only one type GEN, and arguments of many library built-in functions are declared GEN. This means, compiler cannot warn meaningless function call:

  GEN a = stoi (10);        /* a is integer 10 */
  GEN b = cgetg (3, t_MAT); /* b is matrix */
  GEN c = gadd (a, b);      /* nonsense! */

Ofcourse, Pari screams error message but it is run-time. It is advised that, in your program, assertion like

  assert (typ (a) == t_INT && typ (b) == t_INT);

should be inserted, as a checkpoint.

Q & A

Q Does my gp program run faster if I rewrite it in C?

A It depends.

• If your gp program only calls a few builtin, No.

• If your gp program is very complicated, C version will be little bit faster. But beware Total Cost.

• If your gp program contains many loop constructs (for(), sum() and those variants), C version will be faster greatly.

To get the latest version of Pari-gp ftp ftp://megrez.math.u-bordeaux.fr/pub/pari cvs See CVS.txt in the source distribution Developers’ Mailing list: mailto://pari-dev@list.cr.yp.to

References

• H. Cohen, A course in computational algebraic number theory, Springer.

• —, Advanced topics in computational number theory, Springer.

• C. Batut, K. Belabas, D. Bernardi, H. Cohen, M. Olivier, User’s guide to Pari-gp.