home *** CD-ROM | disk | FTP | other *** search
/ Aminet 7 / Aminet 7 - August 1995.iso / Aminet / misc / math / gpamiga_1_38_3.lha / examples / mattrans.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-02  |  1.4 KB  |  64 lines
  1. #include "genpari.h"
  2.  
  3. GEN matexp(GEN,long);
  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=(long)(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);
  18.   return 0;
  19. }
  20.  
  21. GEN matexp(GEN x,long prec)
  22. {
  23.   GEN y,r,s,p1,p2;
  24.   long tx=typ(x),lx=lg(x),i,k,n,lbot,ltop;
  25.  
  26. /* check that x is a square matrix */
  27.  
  28.   if(tx!=19) {printf("This expression is not a matrix\n");return(gzero);}
  29.   if(lx!=lg((GEN)x[1])) {printf("Not a square matrix\n");return(gzero);}
  30.  
  31. /* compute the L2 norm of x */
  32.  
  33.   ltop=avma;s=gzero;
  34.   r=cgetr(prec+1);gaffsg(1,r);p1=gmul(r,x);
  35.   for(i=1;i<lx;i++) s=gadd(s,gnorml2((GEN)p1[i]));
  36.   if(typ(s)==2) setlg(s,3);
  37.   s=gsqrt(s,3); /* we do not need much precision on s */
  38.  
  39. /* if s<1 we are happy */
  40.  
  41.   if(expo(s)<0) n=0;
  42.   else {n=expo(s)+1;p1=gmul2n(p1,-n);setexpo(s,-1);}
  43.  
  44. /* initializations before the loop */
  45.  
  46.   y=gscalmat(r,lx-1); /* this creates the scalar matrix r*Id */
  47.   p2=p1;r=s;k=1;
  48.   lbot=avma;y=gadd(y,p2);
  49.  
  50. /* now the main loop */
  51.  
  52.   while(expo(r) >= -32*(prec-1))
  53.     {
  54.       k++;p2=gdivgs(gmul(p2,p1),k);r=gdivgs(gmul(s,r),k);
  55.       lbot=avma;y=gadd(y,p2);
  56.     }
  57.  
  58. /* now square back n times if necessary */
  59.  
  60.   for(i=0;i<n;i++) {lbot=avma;y=gmul(y,y);}
  61.   y=gerepile(ltop,lbot,y);
  62.   return(y);
  63. }
  64.