\title{Appendix C : A complete program}

We give here the listing of the program seen in detail in section 4.3,
with the slight modifications explained at the end of that section.

{\tt \obeylines\parskip=0pt plus 1pt
\hbox{}
\#include <genpari.h>
\hbox{}
GEN matexp();
\hbox{}
main()
$\{$
\quad GEN x,y;
\quad long prec,d;
\quad char s[512];
\hbox{}
\quad init(1000000,2); /* take a million bytes of memory for the stack */
\quad printf("precision of the computation in decimal digits? ");
\quad scanf("\%d",\&d);if(d<0) prec=3;else prec=d*K1+3;
\quad printf("input your matrix in GP format:$\backslash$n");
\quad s[0]=0;while(!s[0]) gets(s);x=lisexpr(s);
\quad y=matexp(x,prec);
\quad sor(y,'g',d,0);
$\}$
\hbox{}
GEN matexp(x,prec)
\qquad GEN x;
\qquad long prec;
\hbox{}
$\{$
\quad GEN y,r,s,p1,p2;
\quad long tx=typ(x),lx=lg(x),i,k,n,lbot,ltop;
\hbox{}
/* check that x is a square matrix */
\hbox{}
\quad if(tx!=19) $\{$printf("This expression is not a matrix$\backslash$n");return(gzero);$\}$
\quad if(lx!=lg(x[1])) $\{$printf("This matrix is not square$\backslash$n");return(gzero);$\}$
\hbox{}
/* convert x to real or complex of real and compute its L2 norm */
\hbox{}
\quad ltop=avma;s=gzero;
\quad r=cgetr(prec+1);gaffsg(1,r);p1=gmul(r,x);
\quad for(i=1;i<lx;i++) s=gadd(s,gnorml2(p1[i]));
\quad if(typ(s)==2) setlg(s,3);
\quad s=gsqrt(s,3); /* we do not need much precision on s */
\hbox{}
/* if s<1, we are happy */
\hbox{}
\quad if(expo(s)<0) n=0;
\quad else $\{$n=expo(s)+1;p1=gmul2n(p1,-n);setexpo(s,-1);$\}$
\hbox{}
/* initialisations before the loop */
\hbox{}
\quad y=gscalmat(r,lx-1); /* this creates the scalar matrix r$*{\tt I_{\text {lx-1}}}$ */
\quad p2=p1;r=s;k=1;
\quad lbot=avma;y=gadd(y,p2);
\hbox{}
/* now the main loop */
\hbox{}
\quad while(expo(r) >= -32*(prec-1))
\quad $\{$
\qquad k++;p2=gdivgs(gmul(p2,p1),k);r=gdivgs(gmul(s,r),k);
\qquad lbot=avma;y=gadd(y,p2);
\quad $\}$
\hbox{}
/* now square back n times if necessary */
\hbox{}
\quad for(i=0;i<n;i++) $\{$lbot=avma;y=gmul(y,y);$\}$
\quad y=gerepile(ltop,lbot,y);
\quad return(y);
$\}$
}

\vfill\eject