Return to matexp.c CVS log

Up to [local] / OpenXM_contrib / pari / examples

Initial revision

#include "pari.h"

GEN
matexp(GEN x,long prec)
{
  long lx=lg(x),i,k,n, ltop = avma;
  GEN y,r,s,p1,p2;

  /* check that x is a square matrix */
  if (typ(x) != t_MAT) err(typeer,"matexp");
  if (lx == 1) return cgetg(1, t_MAT);
  if (lx != lg(x[1])) err(talker,"not a square matrix");

  /* convert x to real or complex of real and compute its L2 norm */
  s = gzero; r = cgetr(prec+1); affsr(1,r); x = gmul(r,x);
  for (i=1; i<lx; i++)
    s = gadd(s, gnorml2((GEN)x[i]));
  if (typ(s) == t_REAL) setlg(s,3);
  s = gsqrt(s,3); /* we do not need much precision on s */

  /* if s<1 we are happy */
  k = expo(s);
  if (k < 0) { n = 0; p1 = x; }
  else { n = k+1; p1 = gmul2n(x,-n); setexpo(s,-1); }

  /* initializations before the loop */
  y = gscalmat(r,lx-1); /* creates scalar matrix with r on diagonal */
  p2 = p1; r = s; k = 1;
  y = gadd(y,p2);

  /* now the main loop */
  while (expo(r) >= -BITS_IN_LONG*(prec-1))
  {
    k++; p2 = gdivgs(gmul(p2,p1),k);
    r = gdivgs(gmul(s,r),k); y = gadd(y,p2);
  }

  /* now square back n times if necessary */
  for (i=0; i<n; i++) y = gsqr(y);
  return gerepileupto(ltop,y);
}

int
main()
{
  long d, prec = 3;
  GEN x;

  /* take a stack of 10^6 bytes, no prime table */
  pari_init(1000000, 2);
  printf("precision of the computation in decimal digits:\n");
  d = itos(lisGEN(stdin));
  if (d > 0) prec = (long)(d*pariK1+3);

  printf("input your matrix in GP format:\n");
  x = matexp(lisGEN(stdin), prec);

  sor(x, 'g', d, 0);
  exit(0);
}