Re: bug report/techn. question: incorrect type in ...

[Karim BELABAS <Karim.Belabas@math.u-psud.fr>]

[Maximilian Pitschi:]
> I'm trying to do some matix computations, e.g., matdet(...), with
> identical matrices M1 and M2, which depend on a parameter p.  See log
> below. To me, the matrices differ only in the way they were defined.

The following patch should get rid of _some_ of the problems linked to bad
guesses in the linear algebra package: the ones connected with polynomials
or rational functions (in particular your example with matdet/matsolve now
works ok).

[Note: I also changed slightly the maximal pivot strategy, comparing
exponents, instead of the numbers themselves. Hence this breaks the
`linear' bench, producing a different (in fact, slightly more precise)
result for 1/(1. * mathilbert(7))]

Cheers, Karim.

*** src/basemath/alglin1.c.orig	Thu Jul 23 15:40:23 1998
--- src/basemath/alglin1.c	Thu Jul 23 17:52:19 1998
***************
*** 608,613 ****
--- 608,648 ----
  /*                       Solve A*X=B (Gauss pivot)                 */
  /*                                                                 */
  /*******************************************************************/
+ 
+ /* Assume x is a non-empty matrix. Return 0 if maximal pivot should not be
+  * used, and the matrix precision (min real precision of coeffs) otherwise.
+  */
+ static long 
+ matprec(GEN x)
+ {
+   long tx,i,j,l, k = VERYBIGINT, lx = lg(x), ly = lg(x[1]);
+   GEN p1;
+   for (i=1; i<lx; i++)
+     for (j=1; j<ly; j++) 
+     {
+       p1 = gmael(x,i,j); tx = typ(p1);
+       if (!is_scalar_t(tx)) return 0;
+       l = precision(p1); if (l && l<k) k = l;
+     }
+   return (k==VERYBIGINT)? 0: k;
+ }
+ 
+ /* As above, returning 1 if the precision would be non-zero, 0 otherwise */
+ static long
+ use_maximal_pivot(GEN x)
+ {
+   long tx,i,j, lx = lg(x), ly = lg(x[1]);
+   GEN p1;
+   for (i=1; i<lx; i++)
+     for (j=1; j<ly; j++) 
+     {
+       p1 = gmael(x,i,j); tx = typ(p1);
+       if (!is_scalar_t(tx)) return 0;
+       if (precision(p1)) return 1;
+     }
+   return 0;
+ }
+ 
  static GEN
  check_b(GEN b, long nbli)
  {
***************
*** 645,651 ****
  GEN
  gauss(GEN a, GEN b)
  {
!   long inexact,ismat,nbli,nbco,i,j,k,av,av1,tetpil,lim;
    GEN p,m,u;
    /* nbli: nb lines of b = nb columns of a */
    /* nbco: nb columns of b (if matrix) */
--- 680,686 ----
  GEN
  gauss(GEN a, GEN b)
  {
!   long inexact,ismat,nbli,nbco,i,j,k,av,tetpil,lim;
    GEN p,m,u;
    /* nbli: nb lines of b = nb columns of a */
    /* nbco: nb columns of b (if matrix) */
***************
*** 664,706 ****
    a = dummycopy(a);
    b = check_b(b,nbli);
    nbco = lg(b)-1;
!   inexact = isinexactreal(a);
    ismat   = (typ(b)==t_MAT);
    if(DEBUGLEVEL>4)
      fprintferr("Entering gauss with inexact=%ld ismat=%ld\n",inexact,ismat);
  
    for (i=1; i<nbli; i++)
    {
-     long exchange;
- 
      /* k is the line where we find the pivot */
      p=gcoeff(a,i,i); k=i;
      if (inexact) /* maximal pivot */
      {
!       GEN p1, p2;
!       av1 = avma;
!       p2 = gabs(p,DEFAULTPREC);
        for (j=i+1; j<=nbli; j++)
        {
!         p1 = gabs(gcoeff(a,j,i),DEFAULTPREC);
!         if (gcmp(p1,p2)>0) { p2=p1; k=j; }
        }
!       if (gcmp0(p2)) err(matinv1);
!       exchange = (k > i);
!       avma = av1;
      }
!     else /* first non-zero pivot */
      {
!       exchange = gcmp0(p);
!       if (exchange)
!       {
!         do k++; while (k<=nbli && gcmp0(gcoeff(a,k,i)));
!         if (k>nbli) err(matinv1);
!       }
      }
  
!     /* exchange==1 if k<>i, we exchange the lines s.t. k=i */
!     if (exchange)
      {
        for (j=i; j<=nbli; j++)
        {
--- 699,731 ----
    a = dummycopy(a);
    b = check_b(b,nbli);
    nbco = lg(b)-1;
!   inexact = use_maximal_pivot(a);
    ismat   = (typ(b)==t_MAT);
    if(DEBUGLEVEL>4)
      fprintferr("Entering gauss with inexact=%ld ismat=%ld\n",inexact,ismat);
  
    for (i=1; i<nbli; i++)
    {
      /* k is the line where we find the pivot */
      p=gcoeff(a,i,i); k=i;
      if (inexact) /* maximal pivot */
      {
!       long e, ex = gexpo(p);
        for (j=i+1; j<=nbli; j++)
        {
!         e = gexpo(gcoeff(a,j,i));
!         if (e > ex) { ex=e; k=j; }
        }
!       if (gcmp0(gcoeff(a,k,i))) err(matinv1);
      }
!     else if (gcmp0(p)) /* first non-zero pivot */
      {
!       do k++; while (k<=nbli && gcmp0(gcoeff(a,k,i)));
!       if (k>nbli) err(matinv1);
      }
  
!     /* if (k!=i), exchange the lines s.t. k = i */
!     if (k != i)
      {
        for (j=i; j<=nbli; j++)
        {
***************
*** 1043,1057 ****
  static void
  gauss_get_prec(GEN x, long prec)
  {
!   long pr;
  
!   gauss_is_zero = &gcmp0;
!   if (!prec)
!   {
!     if (!isinexactreal(x)) return;
!     prec = DEFAULTPREC;
!   }
!   pr = gprecision(x); if (!pr) return;
    if (pr > prec) prec = pr;
  
    gauss_ex = BITS_IN_LONG - bit_accuracy(prec);
--- 1068,1076 ----
  static void
  gauss_get_prec(GEN x, long prec)
  {
!   long pr = matprec(x); 
  
!   if (!pr) { gauss_is_zero = &gcmp0; return; }
    if (pr > prec) prec = pr;
  
    gauss_ex = BITS_IN_LONG - bit_accuracy(prec);
***************
*** 1489,1501 ****
  static GEN
  det_simple_gauss(GEN a, long inexact)
  {
!   long nbco,i,j,k,av,av1,s,exchange;
    GEN x,p,m;
  
-   if (typ(a)!=t_MAT) err(mattype1,"det2");
-   nbco=lg(a)-1; if (!nbco) return gun;
-   if (nbco != lg(a[1])-1) err(mattype1,"det2");
- 
    av=avma; s=1; x=gun; a=dummycopy(a);
    for (i=1; i<nbco; i++)
    {
--- 1508,1516 ----
  static GEN
  det_simple_gauss(GEN a, long inexact)
  {
!   long i,j,k,av,av1,s, nbco = lg(a)-1;
    GEN x,p,m;
  
    av=avma; s=1; x=gun; a=dummycopy(a);
    for (i=1; i<nbco; i++)
    {
***************
*** 1512,1529 ****
        }
        p1 = gcoeff(a,i,k);
        if (gcmp0(p1)) { av1=avma; return gerepile(av,av1,gcopy(p1)); }
-       exchange = (k > i);
      }
!     else
      {
!       exchange = gcmp0(p);
!       if (exchange)
!       {
!         do k++; while(k<=nbco && gcmp0(gcoeff(a,i,k)));
!         if (k>nbco) { avma=av; return gzero; }
!       }
      }
!     if (exchange)
      {
        j=a[k]; a[k]=a[i]; a[i]=j;
        s = -s; p = gcoeff(a,i,i); 
--- 1527,1539 ----
        }
        p1 = gcoeff(a,i,k);
        if (gcmp0(p1)) { av1=avma; return gerepile(av,av1,gcopy(p1)); }
      }
!     else if (gcmp0(p))
      {
!       do k++; while(k<=nbco && gcmp0(gcoeff(a,i,k)));
!       if (k>nbco) { avma=av; return gzero; }
      }
!     if (k != i)
      {
        j=a[k]; a[k]=a[i]; a[i]=j;
        s = -s; p = gcoeff(a,i,i); 
***************
*** 1546,1554 ****
  }
  
  GEN
! det2(GEN x)
  {
!   return det_simple_gauss(x,isinexactreal(x));
  }
  
  /* determinant in a ring A: all computations are done within A
--- 1556,1568 ----
  }
  
  GEN
! det2(GEN a)
  {
!   long nbco = lg(a)-1;
!   if (typ(a)!=t_MAT) err(mattype1,"det2");
!   if (!nbco) return gun;
!   if (nbco != lg(a[1])-1) err(mattype1,"det2");
!   return det_simple_gauss(a,use_maximal_pivot(a));
  }
  
  /* determinant in a ring A: all computations are done within A
***************
*** 1557,1569 ****
  GEN
  det(GEN a)
  {
!   long nbco,i,j,k,av,av1,s;
    GEN p1,p,m,pprec;
  
-   if (isinexactreal(a)) return det_simple_gauss(a,1);
    if (typ(a)!=t_MAT) err(mattype1,"det");
!   nbco=lg(a)-1; if (!nbco) return gun;
    if (nbco != lg(a[1])-1) err(mattype1,"det");
  
    av=avma; a=dummycopy(a);
    pprec=gun; s=1;
--- 1571,1583 ----
  GEN
  det(GEN a)
  {
!   long nbco = lg(a)-1,i,j,k,av,av1,s;
    GEN p1,p,m,pprec;
  
    if (typ(a)!=t_MAT) err(mattype1,"det");
!   if (!nbco) return gun;
    if (nbco != lg(a[1])-1) err(mattype1,"det");
+   if (use_maximal_pivot(a)) return det_simple_gauss(a,1);
  
    av=avma; a=dummycopy(a);
    pprec=gun; s=1;
--
Karim Belabas                    email: Karim.Belabas@math.u-psud.fr
Dep. de Mathematiques, Bat. 425
Universite Paris-Sud             Tel: (00 33) 1 69 15 57 48
F-91405 Orsay (France)           Fax: (00 33) 1 69 15 60 19
--
PARI/GP Home Page: http://pari.home.ml.org