@@ -621,7 +621,7 @@ bigint egr_index(const vector<Point>& Plist, int real_too)
621621 imagematrix.push_back (im);
622622 vector<int > CG=CGS.ComponentGroup (*pi);
623623 for (unsigned int ni=0 ; ni<CG.size (); ni++, n++)
624- moduli.push_back (CG[ni]);
624+ moduli.push_back (CG[ni]);
625625 }
626626 mat m (Plist.size (),n);
627627 unsigned int j=0 , j1, j2, i;
@@ -806,7 +806,7 @@ bigint comp_map_image(const vector<int> moduli, const mat& image)
806806#endif
807807 if (modulus==1 ) continue ;
808808#ifdef DEBUG_INDEX
809- cout<<" Column = " <<m.col (j)<<endl;
809+ cout<<" Column " <<j<< " = " <<m.col (j)<<endl;
810810#endif
811811 for (i=1 ; (i<=npts); i++) m (i,j)=m (i,j)%modulus;
812812 long g=0 ,gm;
@@ -830,7 +830,7 @@ bigint comp_map_image(const vector<int> moduli, const mat& image)
830830#endif
831831 if (modulus==1 ) continue ;
832832#ifdef DEBUG_INDEX
833- cout<<" Column = " <<m.col (j)<<endl;
833+ cout<<" Column " <<j<< " = " <<m.col (j)<<endl;
834834#endif
835835 long colmin=modulus, imin=0 , r, q;
836836 while (abs (colmin)>1 )
@@ -843,7 +843,7 @@ bigint comp_map_image(const vector<int> moduli, const mat& image)
843843 if (abs (m (i,j))<abs (colmin)) colmin=m (i,j);
844844 for (i=1 ; i<=npts; i++) if (m (i,j)==colmin) {imin=i; break ;}
845845#ifdef DEBUG_INDEX
846- cout<<" colmin = " <<colmin<<" at imin=" <<imin<<endl;
846+ cout<<" colmin = " <<colmin<<" at row imin=" <<imin<<endl;
847847#endif
848848 for (i=1 ; i<=npts; i++)
849849 if (ndivides (colmin,m (i,j)))
@@ -856,31 +856,36 @@ bigint comp_map_image(const vector<int> moduli, const mat& image)
856856 for (jj=1 ; jj<=np; jj++) m (i,jj)-=q*m (imin,jj);
857857 }
858858#ifdef DEBUG_INDEX
859- cout<<" Column = " <<m.col (j)<<endl;
859+ cout<<" Column " <<j<< " = " <<m.col (j)<<endl;
860860#endif
861- }
861+ }
862862 // now |colmin|=1 and we can clear
863863#ifdef DEBUG_INDEX
864864 cout<<" colmin = " <<colmin<<" , clearing..." <<endl;
865- cout<<" Column = " <<m.col (j)<<endl;
865+ cout<<" Column " <<j<< " = " <<m.col (j)<<endl;
866866#endif
867- for (i=1 ; i<=npts; i++) if (m (i,j)==colmin) {imin=i; break ;}
867+ for (i=1 ; i<=npts; i++)
868+ if (m (i,j)==colmin)
869+ {
870+ imin=i; break ;
871+ }
868872 for (i=1 ; i<=npts; i++)
869873 if (i!=imin)
870- for (jj= 1 ; jj<=np; jj++)
871- {
872- if (colmin== 1 )
873- m (i,jj)-= m (i,j)* m (imin,jj);
874- else
875- m (i,jj)+= m (i,j)* m (imin,jj);
876- }
874+ {
875+ long piv = m (i,j)*colmin;
876+ for (jj= 1 ; jj<=np; jj++ )
877+ {
878+ m (i,jj) = ( m (i,jj)-piv* m (imin,jj)) % moduli[jj- 1 ];
879+ }
880+ }
877881 // now update matrix and ans:
878882#ifdef DEBUG_INDEX
883+ cout<<" Now matrix = " <<m<<endl;
879884 cout<<" Column = " <<m.col (j)<<endl;
880885 cout<<" multiplying ans and row " <<imin<<" by " <<modulus<<endl;
881886#endif
882887 ans*=modulus;
883- for (jj=1 ; jj<=np; jj++)
888+ for (jj=1 ; jj<=np; jj++)
884889 m (imin,jj)=(modulus*m (imin,jj))%moduli[jj-1 ];
885890#ifdef DEBUG_INDEX
886891 cout<<" Now matrix = " <<m<<endl;
0 commit comments