Karim BELABAS on Thu, 23 Jul 1998 18:20:57 +0200 (MET DST)

[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

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

```[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
```