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

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: }