Annotation of OpenXM_contrib/pari/doc/appb.tex, Revision 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
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>