home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / math / gp / appc.tex < prev    next >
Encoding:
Text File  |  1993-07-28  |  2.0 KB  |  76 lines
  1. \title{Appendix C : A complete program}
  2.  
  3. We give here the listing of the program seen in detail in section 4.3,
  4. with the slight modifications explained at the end of that section.
  5.  
  6. {\tt \obeylines\parskip=0pt plus 1pt
  7. \hbox{}
  8. \#include <genpari.h>
  9. \hbox{}
  10. GEN matexp();
  11. \hbox{}
  12. main()
  13. $\{$
  14. \quad GEN x,y;
  15. \quad long prec,d;
  16. \quad char s[512];
  17. \hbox{}
  18. \quad init(1000000,2); /* take a million bytes of memory for the stack */
  19. \quad printf("precision of the computation in decimal digits? ");
  20. \quad scanf("\%d",\&d);if(d<0) prec=3;else prec=d*K1+3;
  21. \quad printf("input your matrix in GP format:$\backslash$n");
  22. \quad s[0]=0;while(!s[0]) gets(s);x=lisexpr(s);
  23. \quad y=matexp(x,prec);
  24. \quad sor(y,'g',d,0);
  25. $\}$
  26. \hbox{}
  27. GEN matexp(x,prec)
  28. \qquad GEN x;
  29. \qquad long prec;
  30. \hbox{}
  31. $\{$
  32. \quad GEN y,r,s,p1,p2;
  33. \quad long tx=typ(x),lx=lg(x),i,k,n,lbot,ltop;
  34. \hbox{}
  35. /* check that x is a square matrix */
  36. \hbox{}
  37. \quad if(tx!=19) $\{$printf("This expression is not a matrix$\backslash$n");return(gzero);$\}$
  38. \quad if(lx!=lg(x[1])) $\{$printf("This matrix is not square$\backslash$n");return(gzero);$\}$
  39. \hbox{}
  40. /* convert x to real or complex of real and compute its L2 norm */
  41. \hbox{}
  42. \quad ltop=avma;s=gzero;
  43. \quad r=cgetr(prec+1);gaffsg(1,r);p1=gmul(r,x);
  44. \quad for(i=1;i<lx;i++) s=gadd(s,gnorml2(p1[i]));
  45. \quad if(typ(s)==2) setlg(s,3);
  46. \quad s=gsqrt(s,3); /* we do not need much precision on s */
  47. \hbox{}
  48. /* if s<1, we are happy */
  49. \hbox{}
  50. \quad if(expo(s)<0) n=0;
  51. \quad else $\{$n=expo(s)+1;p1=gmul2n(p1,-n);setexpo(s,-1);$\}$
  52. \hbox{}
  53. /* initialisations before the loop */
  54. \hbox{}
  55. \quad y=gscalmat(r,lx-1); /* this creates the scalar matrix r$*{\tt I_{\text {lx-1}}}$ */
  56. \quad p2=p1;r=s;k=1;
  57. \quad lbot=avma;y=gadd(y,p2);
  58. \hbox{}
  59. /* now the main loop */
  60. \hbox{}
  61. \quad while(expo(r) >= -32*(prec-1))
  62. \quad $\{$
  63. \qquad k++;p2=gdivgs(gmul(p2,p1),k);r=gdivgs(gmul(s,r),k);
  64. \qquad lbot=avma;y=gadd(y,p2);
  65. \quad $\}$
  66. \hbox{}
  67. /* now square back n times if necessary */
  68. \hbox{}
  69. \quad for(i=0;i<n;i++) $\{$lbot=avma;y=gmul(y,y);$\}$
  70. \quad y=gerepile(ltop,lbot,y);
  71. \quad return(y);
  72. $\}$
  73. }
  74.  
  75. \vfill\eject
  76.