Return to matexp.c CVS log

Up to [local] / OpenXM_contrib / pari / examples

Annotation of OpenXM_contrib/pari/examples/matexp.c, Revision 1.1

1.1     ! maekawa     1: #include "pari.h"
        !             2:
        !             3: GEN
        !             4: matexp(GEN x,long prec)
        !             5: {
        !             6:   long lx=lg(x),i,k,n, ltop = avma;
        !             7:   GEN y,r,s,p1,p2;
        !             8:
        !             9:   /* check that x is a square matrix */
        !            10:   if (typ(x) != t_MAT) err(typeer,"matexp");
        !            11:   if (lx == 1) return cgetg(1, t_MAT);
        !            12:   if (lx != lg(x[1])) err(talker,"not a square matrix");
        !            13:
        !            14:   /* convert x to real or complex of real and compute its L2 norm */
        !            15:   s = gzero; r = cgetr(prec+1); affsr(1,r); x = gmul(r,x);
        !            16:   for (i=1; i<lx; i++)
        !            17:     s = gadd(s, gnorml2((GEN)x[i]));
        !            18:   if (typ(s) == t_REAL) setlg(s,3);
        !            19:   s = gsqrt(s,3); /* we do not need much precision on s */
        !            20:
        !            21:   /* if s<1 we are happy */
        !            22:   k = expo(s);
        !            23:   if (k < 0) { n = 0; p1 = x; }
        !            24:   else { n = k+1; p1 = gmul2n(x,-n); setexpo(s,-1); }
        !            25:
        !            26:   /* initializations before the loop */
        !            27:   y = gscalmat(r,lx-1); /* creates scalar matrix with r on diagonal */
        !            28:   p2 = p1; r = s; k = 1;
        !            29:   y = gadd(y,p2);
        !            30:
        !            31:   /* now the main loop */
        !            32:   while (expo(r) >= -BITS_IN_LONG*(prec-1))
        !            33:   {
        !            34:     k++; p2 = gdivgs(gmul(p2,p1),k);
        !            35:     r = gdivgs(gmul(s,r),k); y = gadd(y,p2);
        !            36:   }
        !            37:
        !            38:   /* now square back n times if necessary */
        !            39:   for (i=0; i<n; i++) y = gsqr(y);
        !            40:   return gerepileupto(ltop,y);
        !            41: }
        !            42:
        !            43: int
        !            44: main()
        !            45: {
        !            46:   long d, prec = 3;
        !            47:   GEN x;
        !            48:
        !            49:   /* take a stack of 10^6 bytes, no prime table */
        !            50:   pari_init(1000000, 2);
        !            51:   printf("precision of the computation in decimal digits:\n");
        !            52:   d = itos(lisGEN(stdin));
        !            53:   if (d > 0) prec = (long)(d*pariK1+3);
        !            54:
        !            55:   printf("input your matrix in GP format:\n");
        !            56:   x = matexp(lisGEN(stdin), prec);
        !            57:
        !            58:   sor(x, 'g', d, 0);
        !            59:   exit(0);
        !            60: }