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
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>