Return to appb.tex CVS log

Up to [local] / OpenXM_contrib / pari / doc

Annotation of OpenXM_contrib/pari/doc/appb.tex, Revision 1.1.1.1

1.1       maekawa     1: \appendix{A Sample program and Makefile}
                      2:
                      3: We assume that you have installed the PARI library and include files as
                      4: explained in Appendix A or in the installation guide. If you chose
                      5: differently any of the directory names, change them accordingly in the
                      6: Makefiles.
                      7:
                      8: If the program example that we have given is in the file \kbd{matexp.c} (say
                      9: as the first of several matrix transcendental functions), then a sample
                     10: Makefile might look as follows. Note that the actual file
                     11: {\tt examples/Makefile} is much more elaborate and you should have a look at
                     12: it if you intend to use {\tt install()} on custom made functions, see
                     13: \secref{se:install}.
                     14:
                     15: \bprog
                     16: CC = cc
                     17: INCDIR = \includedir
                     18: LIBDIR = \libdir
                     19: CFLAGS = -O -I\$(INCDIR) -L\$(LIBDIR)
                     20: \h
                     21: all:\qquad matexp
                     22: \h
                     23: matexp:\qquad  matexp.c
                     24: \q\q \$(CC) \$(CFLAGS) -o matexp matexp.c -lpari -lm
                     25: \eprog
                     26:
                     27: \noindent We then give the listing of the program \kbd{examples/matexp.c}
                     28: seen in detail in \secref{se:prog}, with the slight modifications explained
                     29: at the end of that section.
                     30: %
                     31: \bprog
                     32: \#include <pari.h>
                     33: \h
                     34: GEN
                     35: matexp(GEN x, long prec)
                     36: \obr
                     37: \q  long lx=lg(x),i,k,n, ltop = avma;
                     38: \q  GEN y,r,s,p1,p2;
                     39: \h
                     40: \q /* {\rm check that x is a square matrix} */
                     41: \q if (typ(x) != t\_MAT) err(typeer,"matexp");
                     42: \q if (lx == 1) return cgetg(1, t\_MAT);
                     43: \q if (lx != lg(x[1])) err(talker,"not a square matrix");
                     44: \h
                     45: \q /* {\rm convert x to real or complex of real and compute its L2 norm} */
                     46: \q s = gzero; r = cgetr(prec+1); gaffsg(1,r); x = gmul(r,x);
                     47: \q for (i=1; i<lx; i++)
                     48: \q\q s = gadd(s, gnorml2((GEN)x[i]));
                     49: \q if (typ(s) == t\_REAL) setlg(s,3);
                     50: \q s = gsqrt(s,3); /* {\rm we do not need much precision on s} */
                     51: \h
                     52: \q /* {\rm if $\kbd{s}<1$ we are happy} */
                     53: \q k = expo(s);
                     54: \q if (k < 0) \obr\ n = 0; p1 = x; \cbr
                     55: \q else \obr\ n = k+1; p1 = gmul2n(x,-n); setexpo(s,-1); \cbr
                     56: \h
                     57: \q /* {\rm initializations before the loop} */
                     58: \q y = gscalmat(r,lx-1); /* {\rm creates scalar matrix with r on diagonal} */
                     59: \q p2 = p1; r = s; k = 1;
                     60: \q y = gadd(y,p2);
                     61: \h
                     62: \q /* {\rm now the main loop} */
                     63: \q while (expo(r) >= -BITS\_IN\_LONG*(prec-1))
                     64: \q\obr
                     65: \q\q k++; p2 = gdivgs(gmul(p2,p1),k);
                     66: \q\q r = gdivgs(gmul(s,r),k); y = gadd(y,p2);
                     67: \q\cbr
                     68: \h
                     69: \q /* {\rm now square back \kbd{n} times if necessary} */
                     70: \q for (i=0; i<n; i++) y = gsqr(y);
                     71: \q return gerepileupto(ltop,y);
                     72: \cbr
                     73: \h
                     74: int
                     75: main()
                     76: \obr
                     77: \q long d, prec = 3;
                     78: \q GEN x;
                     79: \h
                     80: \q /* {\rm take a stack of $10^6$ bytes, no prime table} */
                     81: \q pari\_init(1000000,2);
                     82: \q printf("precision of the computation in decimal digits:\bs n");
                     83: \q d = itos(lisGEN(stdin));
                     84: \q if (d > 0) prec = (long)(d*pariK1+3);
                     85: \h
                     86: \q printf("input your matrix in GP format:\bs n");
                     87: \q x = matexp(lisGEN(stdin), prec);
                     88: \h
                     89: \q sor(x, 'g', d, 0);
                     90: \q exit(0);
                     91: \cbr
                     92: \eprog\vfill\eject