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: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>