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     ! 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