Annotation of OpenXM_contrib/pari-2.2/examples/matexp.c, Revision 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: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>