home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 2 / RISC_DISC_2.iso / pd_share / program / code / pari / example / mattrans_c < prev    next >
Encoding:
Text File  |  1991-02-12  |  1.4 KB  |  66 lines
  1. #include <genpari.h>
  2.  
  3. GEN matexp();
  4.  
  5. main()
  6. {
  7.   GEN x,y;
  8.   long prec,d;
  9.   char s[512];
  10.  
  11.   init(1000000,2); /* take a million bytes of memory for the stack */
  12.   printf("precision of the computation in decimal digits? ");
  13.   scanf("%d",&d);if(d<0) prec=3;else prec=d*K1+3;
  14.   printf("input your matrix in GP format:\n");
  15.   s[0]=0;while(!s[0]) gets(s); x=lisexpr(s);
  16.   y=matexp(x,prec);
  17.   sor(y,'g',d,0,0);
  18. }
  19.  
  20. GEN matexp(x,prec)
  21.      GEN x;
  22.      long prec;
  23.  
  24. {
  25.   GEN y,r,s,p1,p2;
  26.   long tx=typ(x),lx=lg(x),i,k,n,lbot,ltop;
  27.  
  28. /* check that x is a square matrix */
  29.  
  30.   if(tx!=19) {printf("This expression is not a matrix\n");return(gzero);}
  31.   if(lx!=lg(x[1])) {printf("This matrix is not square\n");return(gzero);}
  32.  
  33. /* compute the L2 norm of x */
  34.  
  35.   ltop=avma;s=gzero;
  36.   r=cgetr(prec+1);gaffsg(1,r);p1=gmul(r,x);
  37.   for(i=1;i<lx;i++) s=gadd(s,gnorml2(p1[i]));
  38.   if(typ(s)==2) setlg(s,3);
  39.   s=gsqrt(s,3); /* we do not need much precision on s */
  40.  
  41. /* if s<1 we are happy */
  42.  
  43.   if(expo(s)<0) n=0;
  44.   else {n=expo(s)+1;p1=gmul2n(p1,-n);setexpo(s,-1);}
  45.  
  46. /* initializations before the loop */
  47.  
  48.   y=gscalmat(r,lx-1); /* this creates the scalar matrix r*Id */
  49.   p2=p1;r=s;k=1;
  50.   lbot=avma;y=gadd(y,p2);
  51.  
  52. /* now the main loop */
  53.  
  54.   while(expo(r) >= -32*(prec-1))
  55.     {
  56.       k++;p2=gdivgs(gmul(p2,p1),k);r=gdivgs(gmul(s,r),k);
  57.       lbot=avma;y=gadd(y,p2);
  58.     }
  59.  
  60. /* now square back n times if necessary */
  61.  
  62.   for(i=0;i<n;i++) {lbot=avma;y=gmul(y,y);}
  63.   y=gerepile(ltop,lbot,y);
  64.   return(y);
  65. }
  66.