Re: 2.0.7 bug?

[Gerhard Niklasch <nikl@mathematik.tu-muenchen.de>]

Following up further on
> Message-Id: <199803310158.DAA01161@pchelwig1.mathematik.tu-muenchen.de>
> In response to:
> > Message-Id: <199803302350.PAA28338@wish-bone.berkeley.edu>
> > Date: 	Tue, 31 Mar 1998 01:50:21 +0200
> > From: Roland Dreier <dreier@math.berkeley.edu>
> > 
> > ? matsnf(M)
> >   ***   segmentation fault: bug in GP (please report).

> If you play around with it a bit more you'll note that it happens if
> and only if the matrix is singular.  (I hope!)

Well.  smithall() calls hnfmod(), which is actually allhnfmod().  Twice.
If the input matrix is singular, it calls it with a zero modulus, in
which case allhnfmod() simply calls hnf(), which near the end has a short
piece of code `remove null columns'.

So smithall(), which assumes throughout that the matrix is square,
suddenly finds its toy a few columns short.  Oops.

Somewhat reluctantly (as it _might_ break things elsewhere) I've kludged
my way around this by (optionally) disabling the null-columns removal.
[Other functions, in particular hnfperm(), depend on the removal;  it
should only be disabled for the particular calling sequence out of
smithall().]

This exposed _another_, unrelated, buglet (see end of patch) -- smithall()
now ran to completion, but returned something which triggered segfaults
when GP tried to simplify() it or, with simplification turned off, to
gclone() it into the history array!

This is most assuredly an Unofficial Patch -- try at your own risk,
but do try it and give it a good knock around.

Enjoy, Gerhard

PS: One gratuitous speeling errror gratuitously fixed en route. :^)


bash$ diff -u src/basemath/base1.c.orig src/basemath/base1.c
--- src/basemath/base1.c.orig   Tue Mar 31 13:50:44 1998
+++ src/basemath/base1.c        Tue Mar 31 15:57:09 1998
@@ -219,7 +219,7 @@
 }
 
 GEN
-hnf(GEN x)
+hnf0(GEN x, long remove)       /* remove: throw away lin.dep.columns, GN */
 {
   long av0 = avma, s,li,co,av,tetpil,i,j,k,def,ldef,lim;
   GEN p1,u,v,d,denx,a,b;
@@ -268,18 +268,27 @@
       tetpil=avma; x=gerepile(av,tetpil,gcopy(x));
     }
   }
-  /* remove null columns */
-  for (i=1,j=1; j<co; j++)
-    if (!gcmp0((GEN)x[j])) x[i++] = x[j];
-  setlg(x,i); tetpil=avma; 
+  if (remove)
+  {                            /* remove null columns */
+    for (i=1,j=1; j<co; j++)
+      if (!gcmp0((GEN)x[j])) x[i++] = x[j];
+    setlg(x,i);
+  }
+  tetpil=avma; 
   x = denx? gdiv(x,denx): gcopy(x);
   return gerepile(av0,tetpil,x);
 }
 
+GEN
+hnf(GEN x)
+{ return hnf0(x,1); }          /* default is to remove null columns */
+
 #define cmod(x,u,us2) \
   {GEN a=modii((GEN)x,u); if (cmpii(a,us2)>0) a=subii(a,u); x=(long)a;}
 
-/* dm = multiple if diag element (usually detint(x)) */
+/* dm = multiple of diag element (usually detint(x)) */
+/* flag&1: don't/do append dm*matid.  flag&2: don't remove zero columns,
+   needed by smithall() --GN */
 static GEN
 allhnfmod(GEN x,GEN dm,long flag)
 {
@@ -287,7 +296,7 @@
   GEN a,b,w,p1,p2,d,u,v,dms2;
 
   if (typ(dm)!=t_INT) err(typeer,"allhnfmod");
-  if (!signe(dm)) return hnf(x);
+  if (!signe(dm)) return hnf0(x, !(flag&2));
   if (typ(x)!=t_MAT) err(typeer,"allhnfmod");
   if (DEBUGLEVEL>6) fprintferr("entering hnfmod");
 
@@ -296,7 +305,7 @@
   dms2=shifti(dm,-1); li=lg(x[1]);
   av=avma;
 
-  if (flag) 
+  if (flag&1)
   {
     p1 = cgetg(co,t_MAT); for (i=1; i<co; i++) p1[i]=x[i]; 
     if (li > co) err(talker,"nb lines > nb columns in hnfmod");
@@ -1075,7 +1084,7 @@
   if (ishnfall(x)) { if (all) { ml=idmat(n); mr=idmat(n); } }
   else
   {
-    p1=hnfmod(x,mdet);
+    p1=allhnfmod(x,mdet,3);    /* retain zero cols, don't append mdet*matid */
     if (all) { ml=idmat(n); mr=gauss(x,p1); }
     x=p1;
   }
@@ -1093,7 +1102,7 @@
     for (j=1; j<=n; j++) { p3[j]=ml[p2[j]]; p4[j]=mr[p2[j]]; }
     ml=p3; mr=p4;
   }
-  p1=hnfmod(x,mdet);
+  p1=allhnfmod(x,mdet,3);
   if (all) mr=gmul(mr,gauss(x,p1));
   x=p1; mun = negi(gun);
   
@@ -1223,7 +1232,7 @@
   tetpil=avma; z=cgetg(n+1,t_VEC); j=n;
   for (k=n; k; k--)
     if (signe(gcoeff(x,k,k))) z[j--]=labsi(gcoeff(x,k,k));
-  for (   ; k; k--) z[j--]=zero;
+  for (   ; j; j--) z[j]=zero;
   return gerepile(av,tetpil,z);
 }