Return to matexp.c CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / examples

Annotation of OpenXM_contrib/pari-2.2/examples/matexp.c, Revision 1.1.1.1

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