Now that the preliminaries are out of the way, the best way to learn how to use the library mode is to work through a detailed nontrivial example of a main program. We will write a program which computes the exponential of a square matrix x. The complete listing is given in Appendix C, but each part of the program will be written and commented here. We will use an algorithm which is not optimal but is not far from the one used for the PARI function gexp (in fact in the function mpexp1). This consists in calculating the sum of the series:
We will now start writing our program. So as to be able to use it in other contexts, we will structure it in the following way: a main program which will do the input and output, and a function which we shall call matexp which does the real work. The main program is easy to write. It can be something like this:
#include <genpari.h> GEN matexp(); main() { GEN x,y; long prec,d; char s[512]; init(1000000,2); /* take a million bytes of memory for the stack */ printf("precision of the computation in decimal digits? "); scanf("%d",&d);if(d<0) prec=3;else prec=d*K1+3; printf("input your matrix in GP format: \n"); s[0]=0;while(!s[0]) gets(s);x=lisexpr(s); y=matexp(x,prec); sor(y,'g',d,0); }
The variable prec represents the length in longwords of the real numbers used. K1 is the constant (defined in gencom.h) equal to ln(10)/(32*ln(2)), which allows us to convert from a number of decimal digits to a number of longwords. The function lisexpr was mentioned earlier and avoids any trouble in the input. In fact, as was also mentioned, lisexpr can take any legal GP expression hence can do not only input but also computations. Note that we have used the construction s[0]=0;while(!s[0]) gets(s); to get our string, instead of scanf which would make trouble in this case.
Finally, sor is the general output routine. We have chosen to give d significant digits since this is what was asked. Note that there is a trick hidden here: if d was given to be negative, then the computation will be done in precision 3 (i.e. about 9.7 decimal digits) and in the function sor, giving a negative third argument outputs all the significant digits, hence nothing is wrong.
Now let us attack the main course, the function matexp:
GEN matexp(x,prec) GEN x; long prec; { GEN y,r,s,p1,p2; long tx=typ(x),lx=lg(x),i,k,n,lbot,ltop; /* check that x is a square matrix */ if(tx!=19) {printf("This expression is not a matrix \n");return(gzero);} if(lx!=lg(x[1])) {printf("This matrix is not square \n");return(gzero);} /* compute the L2 norm of x */ ltop=avma;s=gzero; for(i=1;i<lx;i++) s=gadd(s,gnorml2(x[i])); if(typ(s)==2) setlg(s,3); s=gsqrt(s,3); /* we do not need much precision on s */ /* if s<1, we are happy */ if(expo(s)<0) {n=0;p1=x;} else {n=expo(s)+1;p1=gmul2n(x,-n);s=gmul2n(s,-n);}
Before continuing, a few remarks are in order. First, after printing an error message we need to return a GEN value otherwise a fussy compiler will complain. Hence we return gzero, which in any case is an impossible result for an exponential.
Second, to compute the square of the L2-norm of x we just add the squares of the L2-norms of the column vectors which we obtain using the library function gnorml2. Had this function not existed, it would of course have been just as easy, but we would have had to make a double loop. Then, we take the square root of s, in precision 3. For that we use the function setlg which tells s to be in such a precision, only if s is of type real. Note that the matrix x is allowed to have complex entries, but the function gnorml2 guarantees that s is a nonnegative number of type real. If we had not known this fact, we would simply have added the instruction s = greal(s); just after the for loop.
Third, before starting this computation which will produce garbage on the stack, we have carefully kept the value of the stack pointer avma in ltop. Note that we are going to assume throughout that the garbage does not overflow the memory that we allocated on the stack. If it did, we would have two solutions. Either allocate more memory in the main program (for instance change 1000000 into 2000000), or do some gerepile along the way.
Fourth, note that we initialized the sum s to gzero, which is an exact zero. This is logical, but has some disadvantages: if all the entries of the matrix are integers (or rational numbers), the computation will be rather long, about twice as long as with real numbers of the same length. It would be better to initialize s to a real zero, using for instance the instructions:
s=cgetr(prec+1);gaffsg(0,s);.
However, this would not make much sense. In fact you should avoid making an assignment of an exact zero (essentially the integer zero) into a real number: which real zero is it going to give as a result? In fact a choice has been made, and it will give you the zero with exponent equal to -32 times the number of longwords in the mantissa, but this is not really satisfactory. Perhaps PARI should give an error message in that case, or at least a warning.
Fifth, we will want to look at the approximate size of real numbers, and the fastest way to do this is to look at their binary exponents. Hence we need to have s as a real number, and not as an integer or a rational number. This is done automatically when we take the square root.
Finally, note the use of the function . This function has the
following syntax:
GEN gmul2n(GEN x,long n);
The effect is simply to multiply x by 2n, where n can be positive or negative. This is much faster than gmul or gmulgs. Note that since n=expo(s)+1, the last gmul2n call could be replaced by setexpo(s,-1);.
There is another function with exactly the same syntax.
When n is nonnegative, the effects of these two functions is the same.
However when n is negative, gshift acts like a right shift of -n,
hence does not give the same effect on integers. The function gshift is
the PARI analogue of the C operators « and ».
We now come to the heart of the function. We have now a GEN p1 which points
to a certain matrix of which we want to take the exponential. We will want
to transform this matrix into a matrix with real (or complex of real) entries
before starting the computation. To do this, we simply multiply by the real
number 1 in precision prec+1 (to be safe). Then, to compute the series
we will keep three variables: a variable p2 which
at stage k will contain
/k!, a variable y which will contain
/i!, and the variable r which will contain
/k!. Note that we do not use Horner's rule. This is simply
because we are lazy and do not want to compute in advance the number of
terms that we need. We leave this modification (and many other improvements!)
to the reader. The program continues as follows:
/* initializations before the loop */
r=cgetr(prec+1);gaffsg(1,r);p1=gmul(r,p1);
y=gscalmat(r,lx-1); /* this creates the scalar matrix r
* */
p2=p1;r=s;k=1;
y=gadd(y,p2);
/* now the main loop */
while(expo(r) >= -32*(prec-1))
{k++;p2=gdivgs(gmul(p2,p1),k);r=gdivgs(gmul(s,r),k);y=gadd(y,p2);}
/* now square back n times if necessary */
if(!n){lbot=avma;y=gerepile(ltop,lbot,gcopy(y));}
else
{
for(i=0;i<n;i++){lbot=avma;y=gmul(y,y);}
y=gerepile(ltop,lbot,y);
}
return(y);
}
A few remarks once again. First note the use of the function
with the following syntax:
GEN gscalmat(GEN x, long l);
The effect of this function is to create the
×
scalar
matrix whose diagonal entries is the GEN x. Hence the length of
the matrix including the codeword will in fact be l+1.
There is a corresponding function
which takes a long as a first argument.
If we refer to what has been said above, the main loop is clear.
When we do the final squarings, according to the fundamental theorem on the use of gerepile we keep the value of avma in lbot just before the squaring, so that if it is the last one, lbot will indeed be the bottom address of the garbage pile, and gerepile will work. Note that it takes a completely negligible time to do this in each loop compared to a matrix squaring. However, when n is initially equal to 0, no squarings have to be done, and we have our final result ready but we lost the address of the bottom of the garbage pile. Hence we use the trick of copying y again on top of the stack. This is inefficient, but does the work. If we wanted to avoid this, the best thing to do would be to put the instruction lbot=avma just before both occurrences of the instruction y=gadd(p2,y);.
Remarks on the program: as such, the program should work most of the time if x is a square matrix with real or complex entries. Indeed, since essentially the first thing that we do is to multiply by the real number 1, the program should work for integer, real, rational, complex or quadratic entries. This is in accordance with the behavior of transcendental functions.
Furthermore, since this program is intended to be only an illustrative
example, it has been written a little sloppily. In particular many error checks
have been omitted, and the efficiency is far from optimal. An
evident improvement is to avoid the unnecessary gcopy by inserting a couple
of extra lbot=avma; instructions. Another improvement is to multiply
the matrix x by the real number 1 right at the beginning, speeding up the
computation of the L2-norm in many cases. These improvements are included
in the version given in Appendix C. Still
another improvement would come from a better choice of n. If the reader takes
a look at the implementation of the function in the file trans1.c
he can make himself the necessary changes. Finally, there exist
other algorithms of a different nature to compute the exponential of a matrix.
Remark: While writing this program, we have seen a few new functions
(the complete list of available functions is given by alphabetical order in the index). However, if you care to look at the file gencom.h, you will notice that
many more functions are defined. But in every case these missing functions
are particular cases of general functions. For example, we have the function
gneg, which takes the negative of a PARI object. However, there also
exist functions like negi (for type integer), negr
(for type real), and mpneg (for type integer or real).
These functions can of course be called by the user but we feel that the few microseconds lost in calling more general functions (in this case gneg) is compensated by the fact that one needs to remember a much smaller number of functions, and also because there is a hidden danger here: the type of the objects that you use, if they are themselves results of a previous computation, is not completely predetermined. For instance the multiplication of a type real by a type integer usually gives a result of type real, except when the integer is 0, in which case according to the PARI philosophy the result is the exact integer 0. Hence if afterwards you call a function which specifically needs a real type argument, you are in trouble.
If you really want to use these functions, their names are self-explanatory once you know that i stands for a PARI integer, r for a PARI real, mp for i or r, s for an ordinary 32-bit signed long, z (as a suffix) meaning as usual that the result is not created on the PARI stack but assigned to a preexisting GEN given as an extra argument.
For completeness, in Chapter 5 we have given a description of all these low-level functions.
Please note that in the present version the names of the functions are not always consistent. This will be changed. Hence anyone programming in PARI must be aware that the names of almost all functions that he uses will change eventually. If the need arises (i.e. if there really are people out there who delve into the innards of PARI), updated versions with no name changes will be released.