// smat.cc: implementation of sparse integer matrix class smat
//////////////////////////////////////////////////////////////////////////
//
// Copyright 1990-2012 John Cremona
// 
// This file is part of the eclib package.
// 
// eclib is free software; you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by the
// Free Software Foundation; either version 2 of the License, or (at your
// option) any later version.
// 
// eclib is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
// for more details.
// 
// You should have received a copy of the GNU General Public License
// along with eclib; if not, write to the Free Software Foundation,
// Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
// 
//////////////////////////////////////////////////////////////////////////
 
// Original version by Luiz Figueiredo
 
// ONLY to be included by smatrix.cc

void showrow(int*pos, scalar*val) // for debugging
{
  int d=pos[0]; 
  cout<<"[";
  int* posi=pos; posi++;
  scalar* vali=(scalar*)val;
  while(d--) cout<<"("<<(*posi++)<<","<<(*vali++)<<")";
  /*
  if(sizeof(scalar)==sizeof(int)) 
    {
      int* vali=(int*)val;
      while(d--) cout<<"("<<(*posi++)<<","<<(*vali++)<<")";
    }
  else 
    {
      long* vali=(long*)val;
      while(d--) cout<<"("<<(*posi++)<<","<<(*vali++)<<")";
    }
  */
  cout<<"]";
}

//#define DEBUG_MEM

// Definitions of member operators and functions:

smat::smat(int nr, int nc)
{
  nco = nc;
  nro = nr;
  col = new int * [nr];
  val = new scalar * [nr];
#ifdef DEBUG_MEM
  cout<<"Constructed an smat with (nr,nc)=("<<nr<<","<<nc<<"), with col="<<col<<", val="<<val<<endl;
#endif
  for( int i = 0; i < nr; i++ )
    {
      col[i] = new int [ 2 ];
      val[i] = new scalar [ 1 ];
      col[i][1] = col[i][0] = val[i][0] = 0;
    }
}

smat::smat(const smat& sm)
{
  nco = sm.nco;
  nro = sm.nro;
  col = new int * [nro];
  val = new scalar * [nro];
#ifdef DEBUG_MEM
  cout<<"Constructed an smat (copy constructor) with col="<<col<<", val="<<val<<endl;
#endif
  for( int i = 0; i < nro; i++ )
    { 
      int d = sm.col[i][0];
      col[i] = new int[ d+1 ];  
      val[i] = new scalar[ d ];  
      int *pos = col[i], *pi = sm.col[i];
      scalar *values = val[i], *vi = sm.val[i];
      *pos++ = *pi++;
      while (d--) { *values++ = *vi++; *pos++ = *pi++; }
    }
}

smat::smat(const mat& m)
{
  //  cout<<"Converting mat("<<m.nro<<"x"<<m.nco<<") to smat"<<endl;
  nco = m.nco;
  nro = m.nro;
  col = new int * [nro];
  val = new scalar * [nro];
#ifdef DEBUG_MEM
  cout<<"Constructed an smat (from a mat) with col="<<col<<", val="<<val<<endl;
#endif
  int i, j, k, l, p;
  for( i = 0; i < nro; i++ )
    {
      scalar *veci = m.entries + i*nco;
      for( j = 0, k = 0; j < nco; j++ ) if( *veci++ ) k++;
      col[i] = new int[ k+1 ];  
      val[i] = new scalar[ k ];  
      scalar *values = val[i]; int *pos = col[i];
      veci = m.entries + i*nco;
      *pos++ = k;
      for( l = 0, p = 1;  l < nco; l++, p++,veci++ ) 
	if( *veci ) { *values++ = *veci; *pos++ = p; }
    }
}

smat::~smat()
{
  for( int i = 0; i < nro; i++ ) { delete [] col[i]; delete [] val[i]; }
#ifdef DEBUG_MEM
  cout<<"Destroying an smat with col="<<col<<", val="<<val<<endl;
#endif
  delete [] col;
  delete [] val;
}

// member functions and operators

void smat::set_row( int i, int d, int* pos, scalar* values)
{
  if( col[i][0] != d ) {
    delete [] col[i]; delete [] val[i]; 
    col[i] = new int [d+1];
    val[i] = new scalar [d];
    col[i][0] = d;
  }
  for( int j = 0; j < d; j++ ) {
    col[i][j+1] = *pos++;
    val[i][j] = *values++;
  }
}

void smat::setrow ( int i, const svec& v) // i counts from 1
{
  int j, d=v.entries.size();
  i--;
  if( col[i][0] != d ) {
    delete [] col[i]; delete [] val[i]; 
    col[i] = new int [d+1];
    val[i] = new scalar [d];
    col[i][0] = d;
  }
  map<int,scalar>::const_iterator vi;
  for(vi=v.entries.begin(), j=0;
      vi!=v.entries.end(); vi++, j++)
    {
      col[i][j+1] = vi->first;
      val[i][j] = vi->second;
    }
}

smat smat::select_rows(const vec& rows) const
  {
    int i,r, n=dim(rows);
    smat ans(n,nco);
    for(i=0; i<n; i++)
      {
	r=rows[i+1]-1;
	ans.set_row(i,col[r][0],col[r]+1,val[r]);
      }
    return ans;
  }

mat smat::as_mat( ) const
{
  //  cout<<"Converting smat to mat("<<nro<<"x"<<nco<<")"<<endl;
  mat ans( nro, nco ); 
  scalar *mi = ans.entries;
  for( int i = 0; i < nro; i++ )
    {
      int d = *col[i];
      scalar *values = val[i];
      int *posi = col[i] + 1;
      while( d-- )
	mi[ i*nco + (*posi++) - 1 ] = *values++;
    }
  return ans;
}

svec smat::row(int i) const // extract row i as an svec, i counts from 1
{
  i--;
  svec ans(nco);
  int d = *col[i];
  scalar *values = val[i];
  int *posi = col[i] + 1;
  while( d-- )
    ans.set( (*posi++),  (*values++));
  return ans;
}

scalar smat::elem(int i, int j)  const   /*returns (i,j) entry, 1 <= i <= nro
        				  * can only be used as a rvalue  */
{
 if( (0<i) && (i<=nro) && (0<j) && (j<=nco) )           
   {
     int d = *col[i-1];
     int *posi = col[i-1] + 1;
     scalar *veci = val[i-1]; 
     while( d-- )
       { 
	 if( j == *posi++ ) return *veci;
	 veci++;
       }
     return 0;
   }
 else 
   {
     cerr << "Bad indices in smat::operator ()\n"; 
     return 0;
   }
}


smat& smat::operator=(const smat& sm)
{
 if (this==&sm) return *this;
 nco = sm.nco;
 int i, n = sm.nro; 
 if (nro != n) // delete old space and replace with new.
   {            
     for( i = 0; i < nro; i++ ) { delete [] col[i]; delete [] val[i]; }
#ifdef DEBUG_MEM
  cout<<"in smat.=, destroying old smat with col="<<col<<", val="<<val<<endl;
#endif
     delete [] col;
     delete [] val;
     nro = n;
     col = new int * [nro]; 
     val = new scalar * [nro]; 
#ifdef DEBUG_MEM
  cout<<"in smat.=, creating new smat with col="<<col<<", val="<<val<<endl;
#endif
  for( i = 0; i < nro; i++ )
    {
      col[i] = new int [ 2 ];
      val[i] = new scalar [ 1 ];
      col[i][1] = col[i][0] = val[i][0] = 0;
    }
   }
 for( i = 0; i < nro; i++ )
   { 
     int d = *sm.col[i];
     if(d!=col[i][0])
       {
	 delete[]col[i]; delete[]val[i];
	 col[i] = new int [ d+1 ];
	 val[i] = new scalar [ d ];
	 col[i][0]=d;
       }
     scalar *values = val[i]; int *pos = col[i];
     scalar *vi = sm.val[i]; int *pi = sm.col[i];
     *pos++ = *pi++;
     while (d--) { *values++ = *vi++; *pos++ = *pi++; }
   }
  return *this;
}

smat& smat::operator+=(const smat& mat2)
{
  if (nro==mat2.nro)
    {
      for(int i = 0; i < nro; i++ )
	{
// 	  cout<<"Adding rows ";
// 	  showrow(col[i],val[i]); cout<<" and "; showrow(mat2.col[i],mat2.val[i]);
// 	  cout<<endl;
	  int d = *col[i], d2 = *mat2.col[i];
	  int *pos1 = col[i] + 1, *pos2 = mat2.col[i]+1;
 	  scalar *val1 = val[i], *val2 = mat2.val[i];
	  int *P = new int [ d + d2 + 1 ]; int* Pi=P+1;
	  scalar *V = new scalar [ d + d2 ]; scalar* Vi=V;
	  int k = 0;       /*k will be # of non-zero entries of sum*/
	  while( d && d2 )
	    { 
	      if( *pos1 < *pos2 ) 
		{  *Pi++ = *pos1++; *Vi++ = *val1++; d--; k++; }
	      else if(*pos2 < *pos1 )
		{  *Pi++ = *pos2++; *Vi++ = *val2++; d2--; k++; }
	      else {
// 		cout<<"two entries at position "<<(*pos1)<<", "<<(*val1)<<" and "<<(*val2)<<endl;
		*Pi = *pos1;
		scalar newval = (*val1++) + (*val2++);
// 		cout<<"sum = "<<newval<<endl;
		if( newval!=0 ) 
		  {
		    *Vi=newval;
// 		    cout<<"nonzero, putting "<<(*V)<<" into sum row in position "<<(*Pi)<<endl; 
		    Vi++; Pi++; k++;
		  }
// 		else cout<<"zero, omitting from sum row"<<endl;
		pos1++; pos2++; d--; d2--;
	      }
	    }
// 	  P[0] = k;
// 	  cout<<"intermediate result is ";
// 	  showrow(P,V);
// 	  cout<<endl;
	  while( d2-- )
	    { *Pi++ = *pos2++; *Vi++ = *val2++; k++; }
	  while( d -- )
	    { *Pi++ = *pos1++; *Vi++ = *val1++; k++; }
	  P[0] = k;
	  delete [] col[i]; col[i]=P;
	  delete [] val[i]; val[i]=V;
// 	  cout<<"Result is ";
// 	  showrow(col[i],val[i]);
// 	  cout<<endl;
	}
    }
  else cerr << "Incompatible smatrices in operator +=\n";
  return *this;
}

smat& smat::operator-=(const smat& mat2)
{
  if (nro==mat2.nro)
    {
      for(int i = 0; i < nro; i++ ) 
	{
	  int d = *col[i], d2 = *mat2.col[i];
	  int *pos1 = col[i] + 1, *pos2 = mat2.col[i]+1;
 	  scalar *val1 = val[i], *val2 = mat2.val[i];
	  int *P = new int [ d + d2 + 1 ]; int* Pi=P+1;
	  scalar *V = new scalar [ d + d2 ]; scalar* Vi=V;
	  int k = 0;       /*k will be # of non-zero entries of sum*/
	  while( d && d2 )
	    { 
	      if( *pos1 < *pos2 ) 
		{ *Pi++ = *pos1++; *Vi++ = *val1++; d--; k++; }
	      else if(*pos2 < *pos1 )
		{  *Pi++ = *pos2++; *Vi++ = -(*val2++); d2--; k++; }
	      else {
		*Pi = *pos1;
		scalar newval = (*val1++) - (*val2++);
		if( newval!=0 ) { *Vi++=newval; Pi++; k++; }
		pos1++; pos2++; d--; d2--;
	      }
	    }
	  while( d2-- )
	    { *Pi++ = *pos2++; *Vi++ = -(*val2++); k++; }
	  while( d -- )
	    { *Pi++ = *pos1++; *Vi++ = *val1++; k++; }
	  P[0] = k;
	  delete [] col[i]; col[i]=P;
	  delete [] val[i]; val[i]=V;
	}
    }
  else cerr << "Incompatible matrices in operator -=\n";
  return *this;
}

smat& smat::operator+= (const scalar& scal) // adds scalar*identity
{
  if(scal==0) return *this;
  int i, d, k;
  for(i = 0; i < nro; i++ )
    {
      d = *col[i];                      // length of old row
      int *pos1 = col[i] + 1;           // pointer to run along position vector
      scalar *val1 = val[i];            // pointer to run along value vector
      int *P = new int [ d + 2 ];       //  new position vector
      scalar *V = new scalar [ d + 1 ]; //  new value vector
      int* Pi=P+1;
      scalar* Vi=V;
      scalar newval;
      k = 0;           // k will be # of non-zero entries of new row
      while((d)&&(*pos1<(i+1)))  // just copy entries
	{ 
	  *Pi++ = *pos1++; *Vi++ = *val1++; k++; d--;
	}
      if(d&&(*pos1==(i+1)))     // add the scalar, see if it's zero
	{
	  newval = (*val1)+scal;
	  if( newval!=0) { *Vi++ = newval; *Pi++=*pos1; k++; }
	  pos1++; val1++; d--;
	}
      else // insert new entry
	{
	  *Vi++ = scal; *Pi++=(i+1); k++; 
	}
      while(d--) // copy remaining entries if necessary
	{
	  *Pi++ = *pos1++; *Vi++ = *val1++; k++; 
	}
      P[0] = k;
      delete [] col[i]; col[i]=P;
      delete [] val[i]; val[i]=V;
    }
  return *this;
}

void smat::sub_mod_p(const scalar& lambda, const scalar& p) 
// subtracts scalar*identity mod p
{
  this->operator-=(lambda);
  this->reduce_mod_p(p);
}

void smat::reduce_mod_p(const scalar& p)
{
  svec rowi;
  for(int i=1; i<=nro; i++)
    {
      rowi = row(i);
      rowi.reduce_mod_p(p);
      setrow(i,rowi);
    }
}

smat& smat::operator*=(scalar scal)
{
  if(scal==0) cerr<<"Attempt to multiply smat by 0\n"<<endl;
  int i, d; scalar *veci;
  for( i = 0; i < nro; i++)
    {
      d = *col[i];
      veci = val[i];
      while(d--) (*veci++) *= scal;
    }
  return *this;
}

smat& smat::mult_by_scalar_mod_p (scalar scal, const scalar& p)
{
  if(xmod(scal,p)==0) cerr<<"Attempt to multiply smat by 0\n"<<endl;
  int i, d; scalar *veci;
  for( i = 0; i < nro; i++)
    {
      d = *col[i];
      veci = val[i];
      while(d--) {(*veci) = xmodmul(*veci,scal,p); veci++;}
    }
  return *this;
}

smat& smat::operator/=(scalar scal)
{
  if(scal==0) cerr<<"Attempt to divide smat by 0\n"<<endl;
  int i, d; scalar *veci;
   for(  i = 0; i < nro; i++)
    {
      d = *col[i];
      veci = val[i];
      while(d--) (*veci++) /= scal;
    }
  return *this;
 }

mat smat::operator*( const mat& m )
{
  if( nco != m.nrows() )
    {
      cerr << "incompatible smat & mat in operator*\n";
      abort();
    }
  mat product( nro, m.ncols() );
  int i, j, d, t;
  scalar ans;
  for( i = 1; i <= nro; i++ ) 
    {
      d = col[i-1][0];
      for( j = 1; j <= m.ncols(); j++ ) 
	{
	  ans = 0;
	  for( t = 0; t < d; t++ ) ans += val[i-1][t]*m(col[i-1][t+1],j);
	  product(i,j) = ans;
	}
    }
  return product;
}

long smat::nullity(const scalar& lambda, scalar mod) // nullity of this-lambda*I
{
  smat sma(*this); sma-=lambda;  return sma.ncols()-sma.rank(mod);
}

// Definitions of non-member, friend operators and functions

svec operator* ( const smat& A, const svec& v )
{
  if( A.nco != dim(v) ) 
    { 
      cout << "incompatible smat*svec\n"; 
      cout << "Dimensions "<<dim(A)<<" and "<<dim(v)<<endl;
      abort();
    }
  int n = A.nro, j; scalar s;
  svec prod(n);
  for(j = 1; j<=n; j++)  
    {
      s = (A.row(j))*v;
      if(s) prod.entries[j]=s;
    }
  return prod;
}

vec operator*  (smat& m, const vec& v)
{
  int r = m.nrows(), c=m.ncols();
  if(c!=dim(v))
    {
      cout<<"Error in smat*vec:  wrong dimensions ("<<r<<"x"<<c<<")*"<<dim(v)<<endl;
      abort();
    }
  vec w(r);
  for(int i=1; i<=r; i++) w.set(i,m.row(i)*v);
  return w;
}

// (col) svec * smat

svec operator* ( const svec& v, const smat& A )
{
  if( v.d != A.nrows() ) 
    { 
      cout << "incompatible sizes in v*A\n"; 
      cout << "Dimensions "<<v.d<<" and "<<dim(A)<<endl;
      abort();
    }
  svec prod(A.ncols());
  map<int,scalar>::const_iterator vi;
  for(vi=v.entries.begin(); vi!=v.entries.end(); vi++)
    prod += (vi->second)*(A.row(vi->first));
  return prod;
}

#if(0)
svec mult_mod_p( const svec& v, const smat& A, const scalar& p  )
{
  if( v.d != A.nrows() ) 
    { 
      cout << "incompatible sizes in v*A\n"; 
      cout << "Dimensions "<<v.d<<" and "<<dim(A)<<endl;
      abort();
    }
  svec prod(A.ncols());
  map<int,scalar>::const_iterator vi;
  for(vi=v.entries.begin(); vi!=v.entries.end(); vi++)
    prod.add_scalar_times_mod_p(A.row(vi->first), vi->second,p);
  return prod;
}

#else

svec mult_mod_p( const svec& v, const smat& A, const scalar& p  )
{
  if( v.d != A.nrows() )
    {
      cout << "incompatible sizes in v*A\n";
      cout << "Dimensions "<<v.d<<" and "<<dim(A)<<endl;
      abort();
    }
  vec prod(A.ncols());
  map<int,scalar>::const_iterator vi;
  for(vi=v.entries.begin(); vi!=v.entries.end(); vi++)
    {
      // prod.add_scalar_times_mod_p(A.row(vi->first), vi->second,p);
      int i = (vi->first)-1;        // the row of A to use (from 0)
      scalar c = vi->second;        // the coefficient tomultiply it by
      int d = A.col[i][0];          // #nonzero entries in this row
      int *posi = A.col[i] +1;      // pointer to array of columns
      scalar *values = A.val[i];    // pointer to array of values
      while (d--)
        prod.add_modp(*posi++,xmodmul(c,*values++,p),p);
    }
  return svec(prod);
}
#endif

svec mult_mod_p( const smat& A, const svec& v, const scalar& p  )
{
  if( v.d != A.ncols() )
    {
      cout << "incompatible sizes in A*v\n";
      cout << "Dimensions "<<dim(A)<<" and "<<v.d<<endl;
      abort();
    }
  svec w(A.nrows());
  int i;
  for(i=1; i<=A.nrows(); i++)
    w.set(i,dotmodp(A.row(i),v,p));
  return w;
}

vec mult_mod_p( const smat& A, const vec& v, const scalar& p  )
{
  if( dim(v) != A.ncols() )
    {
      cout << "incompatible sizes in A*v\n";
      cout << "Dimensions "<<dim(A)<<" and "<<dim(v)<<endl;
      abort();
    }
  vec w(A.nrows());
  int i;
  for(i=1; i<=A.nrows(); i++)
    w.set(i,dotmodp(A.row(i),v,p));
  return w;
}

smat operator* ( const smat& A, const smat& B )
{
  if( A.nco != B.nro ) { cerr << "incompatible smats in operator *\n"; abort();}
  int nro = A.nro, nco = B.nco;
  smat prod( nro, nco );

  for (int i=1; i<=nro; i++)
    prod.setrow(i, A.row(i)*B);
  return prod;
}

smat mult_mod_p ( const smat& A, const smat& B, const scalar& p )
{
  if( A.nco != B.nro ) { cerr << "incompatible smats in operator *\n"; abort();}
  int nro = A.nro, nco = B.nco;
  smat prod( nro, nco );

  for (int i=1; i<=nro; i++)
    prod.setrow(i, mult_mod_p(A.row(i),B,p));
  return prod;
}

smat transpose ( const smat& A )
{
  // 1. Count the number of entries in each column (as in operator*() below):
  int *colwts = new int[A.nco];
  int i, r;
  for(i=0; i<A.nco; i++) colwts[i]=0;
  for( r = 0; r <A.nro; r++ ) // counts # of elements in each col
    {
      int d = *A.col[r];
      int *p = A.col[r] + 1;
      while( d-- ) colwts[*p++ - 1]++;
    }
#if(0)
  cout<<"Column weights of A:\n";
  for(i=0; i<A.nco; i++) cout<<colwts[i]<<" ";
  cout<<endl;
#endif
  // 2. Make space for the output matrix:
  smat B(A.nco,A.nro);
  // Remove the default entries in B:
  for( int i = 0; i < B.nro; i++ ) { delete [] B.col[i]; delete [] B.val[i]; }
  // Replace with the correct sizes:
  for( i = 0; i < B.nro; i++ )
    {
      int d = colwts[i];
      B.col[i] = new int[ d+1 ];  
      B.val[i] = new scalar[ d ];  
      B.col[i][0] = d;
    }
  delete[]colwts;
  //3. Copy entries over.  aux[i] holds the number of entries so far
  //   put into row i of the transpose
  int * aux = new int [B.nro];
  int *a=aux;
  for( r = 0; r < B.nro; r++ ) *a++ = 0;
  for( r = 0; r < A.nro; r++ ) {
    int d = *A.col[r];
    //    cout<<"row "<<r<<" of A has "<<d<<" entries\n";
    scalar *v = A.val[r];
    int *p = A.col[r] + 1;
    while( d-- ) {
      int c = *p++ - 1;
      B.col[c][aux[c]+1] = r+1;
      B.val[c][aux[c]] = *v++;
      aux[c]++;
    }
#if(0)
    cout<<"After processing that row, aux = \n";
    for(i=0; i<A.nco; i++) cout<<aux[i]<<" ";
    cout<<endl;
#endif
  }
  delete[]aux;
  return B;
}

int operator==(const smat& sm1, const smat& sm2)
{
   int nr = sm1.nro, i;
   int equal = ( nr == sm2.nro );
   if(!equal) return 0;

   for( i = 0; i < nr && equal; i++ )
     {
       int d1 = *sm1.col[i], d2 = *sm2.col[i];
       if( d1 != d2 ) {return 0;}
       
       scalar *sm1val = sm1.val[i], *sm2val= sm2.val[i];
       int *sm1pos= sm1.col[i] + 1;
       int *sm2pos = sm2.col[i] + 1;
       while (equal && d1--) equal = ((*sm1val++)==(*sm2val++));
       while (equal && d2--) equal = ((*sm1pos++)==(*sm2pos++));
     }
   return equal;
}

int eqmodp(const smat& sm1, const smat& sm2, const scalar& p)
{
   int nr = sm1.nro, i;
   int equal = ( nr == sm2.nro );
   if(!equal) return 0;

   for( i = 0; i < nr && equal; i++ )
     {
       int d1 = *sm1.col[i], d2 = *sm2.col[i];
       if( d1!=d2 ) {return 0;}
       
       scalar *sm1val = sm1.val[i], *sm2val= sm2.val[i];
       int *sm1pos= sm1.col[i] + 1;
       int *sm2pos = sm2.col[i] + 1;
       while (equal && d2--) equal = ((*sm1pos++)==(*sm2pos++));
       while (equal && d1--) equal = (xmod((*sm1val++)-(*sm2val++),p)==0);
     }
   return equal;
}

ostream& operator << (ostream& s, const smat& sm)
{
  for( int i = 0; i < sm.nro; i++ )
    {
      cout << "row[" << i+1 << "] =";
      int d = *sm.col[i]; 
      int *posi = sm.col[i] + 1; scalar *veci = sm.val[i];
      int n = d-1 > 0 ? d-1 : 0;
      s << "{ ";
      s << "values " << "[";
      if( d > 0 ) s << *veci++; 
      while ( n-- ) s << "," << (*veci++); 
      s << "]";
      s << "   positions: " << "[";
      if( d > 0 ) { s << *posi++; d = d-1; }
      while ( d-- ) { s << "," << (*posi++); }
      s << "]    }" << endl;
    }
  return s;
}


istream& operator >> (istream& s, smat& sm)
{
  int *pos = new int [ sm.nco ];
  scalar *values = new scalar [ sm.nco ];
  int r, k, count;
  for( r = 0; r < sm.nro; r++ )
    { 
      cout << "input row " << r+1 << endl;
      int *p = pos; scalar *v = values;
      s >> k;
      for( count = 0; k != 0; s >> k )
	{
	  *v++ = k;
	  s >> k;
	  if( k ) *p++ = k;
	  else { cerr << "enter zero as a value!!!\n"; abort(); }
	  count++;
	}
      
      delete [] sm.col[r];
      delete [] sm.val[r];
      sm.col[r] = new int [ count + 1 ];
      sm.val[r] = new scalar [ count ];
      
      sm.col[r][0] = count;
      p = pos;
      v = values;
      for( k = 0; k < count; k++ ) 
	{ sm.col[r][k+1] = *p++; sm.val[r][k] = *v++; }
    }
  delete [] pos;
  delete [] values;
  return s;
}


    
// Definition of non-friend functions

smat operator+(const smat& sm)
{
        return sm;
}

smat operator-(const smat& sm)
{
        return (-1)*sm;
}

smat operator+(const smat& sm1, const smat& sm2)
{
  smat ans(sm1); ans+=sm2;  return ans;
}

smat operator-(const smat& sm1, const smat& sm2) 
{
  smat ans(sm1); ans-=sm2;  return ans;
}

smat operator*(scalar scal, const smat& sm)
{
  smat ans(sm); ans*=scal;  return ans;
}

smat operator/(const smat& sm, scalar scal)
{
  smat ans(sm); ans/=scal;  return ans;
}

int operator!=(const smat& sm1, const smat& sm2)
{
        return !(sm1==sm2);
}

int get_population(const smat& m )
{
  int r,d,count = 0;
  for( r = 0; r < m.nro; r++ ) 
    {
      d = *(m.col[r]);
      if(d==0) continue;
      int *pos = m.col[r] + 1;
      while( d-- ) { count += ( *pos++ != 0 );}
    }
  return count;
}

smat sidmat(scalar n)  // identity matrix
{
  smat I(n,n); // creates enough space
  for( int i = 0; i < n; i++ )
    {
      I.col[i][0] = 1;   // one entry in this row
      I.col[i][1] = i+1; // ...it's in column i+1
      I.val[i][0] = 1;   // ...its value is 1
    }
  return I;
}

int liftmat(const smat& mm, scalar pr, smat& m, scalar& dd, int trace)
{
  scalar modulus=pr,n,d; long nr,nc; dd=1;
  int succ=0,success=1;
  float lim=floor(sqrt(pr/2.0));
  m = mm;
  if(trace)
    {
      cout << "Lifting mod-p smat;  smat mod "<<pr<<" is:\n";
      cout << m.as_mat();
      cout << "Now lifting back to Q.\n";
      cout << "lim = " << lim << "\n";
    }
  for(nr=0; nr<m.nro; nr++)
    for(nc=0; nc<m.col[nr][0]; nc++)
      {
	succ = modrat(m.val[nr][nc],modulus,lim,n,d);
	success = success && succ;
	dd=lcm(d,dd);
      }
  if(!success)
    {
      //cout << "Problems encountered with modrat lifting of smat." << endl;
      return 0;
    }
  dd=abs(dd);
  if(trace) cout << "Common denominator = " << dd << "\n";
  for(nr=0; nr<m.nro; nr++)
    for(nc=0; nc<m.col[nr][0]; nc++)
      m.val[nr][nc] = mod(xmodmul(dd,(m.val[nr][nc]),pr),pr);
  if(trace) cout << "Lifted smat = " << m.as_mat() << "\n";
  return 1;
}

//#define DEBUG_CHINESE

int liftmats_chinese(const smat& m1, scalar pr1, const smat& m2, scalar pr2,
                     smat& m, scalar& dd)
{
  long modulus=(long)pr1*(long)pr2,n,d,mij;
  long nr,nc,u,v;
  float lim=floor(sqrt(modulus/2.0));

  dd = bezout(pr1,pr2,u,v); //==1
  if (dd!=1) return 0;

  // First time through: compute CRTs, common denominator and success flag
  m = m1; // NB We assume that m1 and m2 have nonzero entries in the same places
  for(nr=0; nr<m1.nro; nr++)
    for(nc=0; nc<m1.col[nr][0]; nc++)
      {
        mij = mod(v*m1.val[nr][nc],pr1)*pr2 + mod(u*m2.val[nr][nc],pr2)*pr1;
        mij = mod(mij,modulus);
#ifdef DEBUG_CHINESE
        if (((mij-m1.val[nr][nc])%pr1)||((mij-m2.val[nr][nc])%pr2))
          {
            cout<< "bad CRT(["<<m1.val[nr][nc]<<","<<m2.val[nr][nc]<<"],["<<pr1<<","<<pr2<<"]) = "<<mij<<endl;
          }
#endif
        m.val[nr][nc] = mij;
	if (modrat(mij,modulus,lim,n,d))
          dd=lcm(d,dd);
        else
          {
#ifdef DEBUG_CHINESE
            cout<<"CRT("<<m1.val[nr][nc]<<","<<m2.val[nr][nc]<<")="<<mij<<" (mod "<<modulus<<") fails to lift (lim="<<lim<<")\n";
            cout << "Problems encountered in chinese lifting of smat modulo "<<pr1<<" and "<<pr2<< endl;
#endif
            return 0;
          }
      }
  dd=abs(dd);
#ifdef DEBUG_CHINESE
  cout << "Common denominator = " << dd << "\n";
#endif
  // Second time through: rescale
  for(nr=0; nr<m.nro; nr++)
    for(nc=0; nc<m.col[nr][0]; nc++)
      {
        m.val[nr][nc] = mod(xmodmul((dd/d),(long)m.val[nr][nc],modulus),modulus);
      }
  return 1;
}

// Possible FLINT_LEVEL values are as follows.
//
// 0: no FLINT support (or a version <2.3)
// 1: support for 64-bit nmod_mat (standard from version 2.3)
// 2: support for 32-bit hmod_mat (non-standard, from version 2.3)
//
// The configure script should have detected whether the functions
// nmod_mat_rref and/or hmod_mat_rref are available to be used here.
//
#if FLINT_LEVEL!=0
//#define TRACE_FLINT_RREF
#include "flint/fmpz.h"
#if (SCALAR_OPTION==1)&&(FLINT_LEVEL==2)
#include "flint/hmod_mat.h"
#undef uscalar
#undef mod_mat
#undef mod_mat_init
#undef mod_mat_clear
#undef mod_mat_entry
#undef mod_mat_nrows
#undef mod_mat_ncols
#undef mod_mat_rref
#undef mod_mat_mul
#define uscalar hlimb_t // unsigned int
#define mod_mat hmod_mat_t // uses unsigned ints
#define mod_mat_init hmod_mat_init
#define mod_mat_clear hmod_mat_clear
#define mod_mat_entry hmod_mat_entry
#define mod_mat_nrows hmod_mat_nrows
#define mod_mat_ncols hmod_mat_ncols
#define mod_mat_rref hmod_mat_rref
#define mod_mat_mul hmod_mat_mul
#else
#include "flint/nmod_mat.h"
#undef uscalar
#undef mod_mat
#undef mod_mat_init
#undef mod_mat_clear
#undef mod_mat_entry
#undef mod_mat_nrows
#undef mod_mat_ncols
#undef mod_mat_rref
#undef mod_mat_mul
#define uscalar mp_limb_t // unsigned long
#define mod_mat nmod_mat_t // uses unsigned longs
#define mod_mat_init nmod_mat_init
#define mod_mat_clear nmod_mat_clear
#define mod_mat_entry nmod_mat_entry
#define mod_mat_nrows nmod_mat_nrows
#define mod_mat_ncols nmod_mat_ncols
#define mod_mat_rref nmod_mat_rref
#define mod_mat_mul nmod_mat_mul
#endif
#include "flint/profiler.h"
#include "eclib/timer.h"

// FLINT has two types for modular matrices: standard in FLINT-2.3 has
// nmod_mat_t with entries of type mp_limb_t (unsigned long);
// non-standard (in an optional branch) is hmod_mat_t, with entries
// hlimb_t (unsigned int).  We use the former when scalar=long and the
// latter when scalar=int, provided that the optional functions are
// present, which should have been determined by the configure script.
// The unsigned scalar types are #define'd as uscalar.

void mod_mat_from_smat(mod_mat& A, const smat& M, scalar pr)
{
  long nr=M.nrows(), nc=M.ncols();
  long i, j;

  // copy of the modulus for FLINT
  uscalar mod = (uscalar)pr;

  // create flint matrix copy of M:
  mod_mat_init(A, nr, nc, mod);
  for(i=0; i<nr; i++)
    for(j=0; j<nc; j++)
      mod_mat_entry(A,i,j) = (uscalar)posmod(M.elem(i+1,j+1),pr);
}

smat smat_from_mod_mat(const mod_mat& A, const scalar& p) //scalar just to fix return type
{
  long nr=mod_mat_nrows(A), nc=mod_mat_ncols(A);

  // create matrix copy of A:
  smat M(nr, nc);
  long i, j;
  for(i=0; i<nr; i++)
    {
      svec rowi(nc);
      for(j=0; j<nc; j++)
        rowi.set(j+1, mod_mat_entry(A,i,j));
      M.setrow(i+1,rowi);
    }
  return M;
}

smat mult_mod_p_flint ( const smat& A, const smat& B, const scalar& pr )
{
  if( A.ncols() != B.nrows() )
    {
      cerr << "incompatible smats in operator *\n";
      abort();
    }
  mod_mat AA, BB, CC;
  mod_mat_from_smat(AA,A,pr);
  mod_mat_from_smat(BB,B,pr);
  mod_mat_init(CC, A.nrows(), B.ncols(), pr);
  // timer T;
  // T.start();
  // mod_mat_mul(CC,AA,BB);
  // T.stop();
  // cout<<"mult_mod_p_flint time (size "<<dim(A)<<"x"<<dim(B)<<"): ";
  // T.show();
  smat C = smat_from_mod_mat(CC, pr);
  mod_mat_clear(AA);
  mod_mat_clear(BB);
  mod_mat_clear(CC);
  return C;
}

#endif