From 406c4314505509f74d73c689552c5c24293d963e Mon Sep 17 00:00:00 2001
From: John Cremona <john.cremona@gmail.com>
Date: Mon, 23 Mar 2015 10:10:52 +0000
Subject: [PATCH] fix saturation bug; bump version to 20150323, library to
 1.1.0

---
 configure.ac   |  4 ++--
 libsrc/egr.cc  | 37 +++++++++++++++++++++----------------
 tests/tsat3.cc |  9 +++++++--
 3 files changed, 30 insertions(+), 20 deletions(-)

diff --git a/configure.ac b/configure.ac
index 211cfa47..98887847 100644
--- a/configure.ac
+++ b/configure.ac
@@ -3,7 +3,7 @@
 # Process this file with autoconf to produce a configure script.
 
 AC_PREREQ([2.65])
-AC_INIT([eclib], [20150228], [john.cremona@gmail.com])
+AC_INIT([eclib], [20150323], [john.cremona@gmail.com])
 AM_INIT_AUTOMAKE([-Wall])
 AC_MSG_NOTICE([Configuring eclib...])
 AC_CONFIG_SRCDIR([libsrc])
@@ -24,7 +24,7 @@ LT_INIT
 # 3. If interfaces were removed (breaks backward compatibility): increment
 # current, and set both revision and age to zero.
 LT_CURRENT=1
-LT_REVISION=0
+LT_REVISION=1
 LT_AGE=0
 AC_SUBST(LT_CURRENT)
 AC_SUBST(LT_REVISION)
diff --git a/libsrc/egr.cc b/libsrc/egr.cc
index 9d383fa0..3a5abfd7 100644
--- a/libsrc/egr.cc
+++ b/libsrc/egr.cc
@@ -621,7 +621,7 @@ bigint egr_index(const vector<Point>& Plist, int real_too)
       imagematrix.push_back(im);
       vector<int> CG=CGS.ComponentGroup(*pi);
       for(unsigned int ni=0; ni<CG.size(); ni++, n++)
-	moduli.push_back(CG[ni]);	
+	moduli.push_back(CG[ni]);
     }
   mat m(Plist.size(),n);
   unsigned int j=0, j1, j2, i;
@@ -806,7 +806,7 @@ bigint comp_map_image(const vector<int> moduli, const mat& image)
 #endif
       if(modulus==1) continue;
 #ifdef DEBUG_INDEX
-      cout<<"Column = "<<m.col(j)<<endl;
+      cout<<"Column "<<j<<" = "<<m.col(j)<<endl;
 #endif
       for(i=1; (i<=npts); i++) m(i,j)=m(i,j)%modulus;
       long g=0,gm;
@@ -830,7 +830,7 @@ bigint comp_map_image(const vector<int> moduli, const mat& image)
 #endif
       if(modulus==1) continue;
 #ifdef DEBUG_INDEX
-      cout<<"Column = "<<m.col(j)<<endl;
+      cout<<"Column "<<j<<" = "<<m.col(j)<<endl;
 #endif
       long colmin=modulus, imin=0, r, q;
       while(abs(colmin)>1)
@@ -843,7 +843,7 @@ bigint comp_map_image(const vector<int> moduli, const mat& image)
 	      if(abs(m(i,j))<abs(colmin)) colmin=m(i,j);
 	  for(i=1; i<=npts; i++) if(m(i,j)==colmin) {imin=i; break;}
 #ifdef DEBUG_INDEX
-	  cout<<"colmin = "<<colmin<<" at imin="<<imin<<endl;
+	  cout<<"colmin = "<<colmin<<" at row imin="<<imin<<endl;
 #endif
 	  for(i=1; i<=npts; i++)
 	    if(ndivides(colmin,m(i,j)))
@@ -856,31 +856,36 @@ bigint comp_map_image(const vector<int> moduli, const mat& image)
 		for(jj=1; jj<=np; jj++) m(i,jj)-=q*m(imin,jj);
 	      }
 #ifdef DEBUG_INDEX
-	  cout<<"Column = "<<m.col(j)<<endl;
+	  cout<<"Column "<<j<<" = "<<m.col(j)<<endl;
 #endif
-	}     
+	}
       // now |colmin|=1 and we can clear
 #ifdef DEBUG_INDEX
       cout<<"colmin = "<<colmin<<", clearing..."<<endl;
-      cout<<"Column = "<<m.col(j)<<endl;
+      cout<<"Column "<<j<<" = "<<m.col(j)<<endl;
 #endif
-      for(i=1; i<=npts; i++) if(m(i,j)==colmin) {imin=i; break;}
+      for(i=1; i<=npts; i++)
+        if(m(i,j)==colmin)
+          {
+            imin=i; break;
+          }
       for(i=1; i<=npts; i++)
 	if(i!=imin)
-	  for(jj=1; jj<=np; jj++) 
-	    {
-	      if(colmin==1)
-		m(i,jj)-=m(i,j)*m(imin,jj);
-	      else
-		m(i,jj)+=m(i,j)*m(imin,jj);
-	    }
+          {
+            long piv = m(i,j)*colmin;
+            for(jj=1; jj<=np; jj++)
+              {
+		m(i,jj) = (m(i,jj)-piv*m(imin,jj)) % moduli[jj-1];
+              }
+          }
       // now update matrix and ans:
 #ifdef DEBUG_INDEX
+      cout<<"Now matrix = "<<m<<endl;
       cout<<"Column = "<<m.col(j)<<endl;
       cout<<"multiplying ans and row "<<imin<<" by "<<modulus<<endl;
 #endif
       ans*=modulus;
-      for(jj=1; jj<=np; jj++) 
+      for(jj=1; jj<=np; jj++)
 	m(imin,jj)=(modulus*m(imin,jj))%moduli[jj-1];
 #ifdef DEBUG_INDEX
       cout<<"Now matrix = "<<m<<endl;
diff --git a/tests/tsat3.cc b/tests/tsat3.cc
index 4945563c..8616bcce 100644
--- a/tests/tsat3.cc
+++ b/tests/tsat3.cc
@@ -125,9 +125,14 @@ int main()
     }
   cout<<"Input torsion points: "<<tpoints<<endl;
 
-  bigint index;
+  bigint index = BIGINT(1);
   vector<long> unsatprimes;
-  int success = saturate_points(C, points, index, unsatprimes, SAT_BND, use_egr, (verbose));
+  int success = 1;
+
+  if (npts)
+    {
+      success = saturate_points(C, points, index, unsatprimes, SAT_BND, use_egr, (verbose));
+    }
 
   if(success)
     {