// Modular matrix implementation file
// Implementor:  Michael J. Fischer
// Date: February 6, 2005

#include <lnv3/lnv3.h>
#include <nttl/gcd.h>
#include <nttl/inverse.h>
#include "error.h"
#include "matrix.h"

#define DEBUG(x) 

using namespace std;

//===============================================================

// Construct n x m matrix
Matrix::Matrix(int n, int m) : nrows_(n), ncols_(m) {
  allocate();
}

//---------------------------------------------------------------
// Copy constructor from submatrix
Matrix::Matrix(const Matrix& A, int  r1, int  r2, int  c1, int  c2) 
  : nrows_(r2-r1), ncols_(c2-c1) 
{
  allocate();
  for (int i=r1; i<r2; i++) {
    for (int j=c1; j<c2; j++) {
      rowhead_[i-r1][j-c1] = A[i][j];
    }
  }
}

//---------------------------------------------------------------
// Allocate matrix
void Matrix::allocate() 
{
  rowhead_ = new (int*)[nrows_];
  for (int i=0; i<nrows_; i++) {
    rowhead_[i] = new int[ncols_];
  }
}  

//---------------------------------------------------------------
// Destructor
Matrix::~Matrix() 
{
  for (int i=0; i<nrows_; i++) {
    delete [] rowhead_[i];
  }
  delete [] rowhead_;
}

//---------------------------------------------------------------
// Reduce modulo n
void Matrix::mod( const Matrix& A,  int n )
{
  for (int i=0; i<nrows(); i++) {
    for (int j=0; j<ncols(); j++) {
      rowhead_[i][j] = A[i][j] % n;
    }
  }
}

//---------------------------------------------------------------
// Chinese-remaindering
void Matrix::chinese( const Matrix& A1, int n1,
	      const Matrix& A2, int n2 )
{
  int n = n1*n2;
  int N1 = n2;
  int N2 = n1;
  int c1 = N1*Inverse( N1, n1 );
  int c2 = N2*Inverse( N2, n2 );

  for (int i=0; i<nrows(); i++) {
    for (int j=0; j<ncols(); j++) {
      rowhead_[i][j] = (c1*A1[i][j] + c2*A2[i][j]) % n;
    }
  }
}

//---------------------------------------------------------------
// Swap rows i and j if different
void Matrix::swapRows(int i, int j)
{
  if (i != j) {
    int* temp = rowhead_[i];
    rowhead_[i] = rowhead_[j];
    rowhead_[j] = temp;
  }
}

//---------------------------------------------------------------
// Multiply row i by c mod n, beginning in column col
void Matrix::scaleRow(int i, int c, int n, int col)
{
  for (int j=col; j<ncols(); j++) {
    rowhead_[i][j] = (rowhead_[i][j]*c) % n;
  }
}

//---------------------------------------------------------------
// Add row k * c  mod n to row i, beginning in column col
void Matrix::addMultToRow(int k, int c, int n, int i, int col)
{
  for (int j=col; j<ncols(); j++) {
    rowhead_[i][j] = (rowhead_[i][j] + rowhead_[k][j]*c) % n;
  }
}


//---------------------------------------------------------------
// Modular Gaussian elimination
// Performs row operations on matrix until no more are possible.
// Returns the number successfully completed.
// I.e., if k is returned, then the upper left k x k submatrix is
// the identity, the rest of the first k columns are zero, and
// column k+1 has no suitable pivot element below row k.
int Matrix::gauss(int modulus)
{
  int n = nrows();
  int m = ncols();
  int i,k;
  int v;			// inverse of pivot element

  // Perform Gaussian elimination

  for (k=0; k<m; k++) {

  // Find pivot element in column k relatively prime to n
    for (i=k; i<n && GCD(rowhead_[i][k], modulus)!=1; i++);
    if (i==n) break;		// no pivot element found

    DEBUG( cout << "Pivot element is " << k << endl );

    // swap rows i and k
    swapRows(i, k);

    // set diagonal element in pivot row to 1
    v = Inverse(rowhead_[k][k], modulus);
    scaleRow(k, v, modulus, k);

    // clear entire pivot column
    for (i=0; i<n; i++) {
      if (i!=k) {
	addMultToRow(k, modulus - rowhead_[i][k], modulus, i, k);
      }
    }
    DEBUG( cout << "After stage " << k << ":" << endl );
    DEBUG( cout << *this << endl );
  }

  return k;
}      

//---------------------------------------------------------------
// Modular Matrix inverse modulo 26
void Matrix::inverse26( const Matrix& A ) {
  int n = nrows();
  if ( n != ncols() ) throw Fatal("Attempt to invert non-square matrix");
 
  // Allocate n x 2n matrix for K
  int m = 2*n;
  Matrix K(n, m);

  // Copy A into left half of K
  for (int i=0; i<n; i++) {
    for (int j=0; j<n; j++) {
      K[i][j] = A[i][j];
    }
  }

  // Initialize right half of K to the identity
  for (int i=0; i<n; i++) {
    for (int j=n; j<m; j++) {
      K[i][j] = (i==(j-n)) ? 1 : 0;
    }
  }

  DEBUG( cout << "Matrix K before gauss():" << endl );
  DEBUG( cout << K << endl );

  // Do Chinese remaindering
  Matrix K1(n, m);
  Matrix K2(n, m);
  K1.mod(K, 2);
  K2.mod(K, 13);
  int rank1 = K1.gauss(2);
  int rank2 = K2.gauss(13);

  DEBUG( cout << "rank1 = " << rank1 << ", rank2 = " << rank2 << endl );

  if ( rank1 != n || rank2 != n ) {
    throw Fatal("Attempt to invert singular matrix");
  }

  K.chinese( K1, 2, K2, 13 );

  DEBUG( cout << "Matrix K after chinese():" << endl );
  DEBUG( cout << K << endl );

  // Copy answer back from right half of K
  for (int i=0; i<n; i++) {
    for (int j=0; j<n; j++) {
      rowhead_[i][j] = K[i][j+n];
    }
  }
}

//---------------------------------------------------------------
// Matrix print
void Matrix::print(ostream& out ) const
{
  for (int i=0; i<nrows(); i++) {
    for (int j=0; j<ncols(); j++) {
      out << setw(3) << rowhead_[i][j];
    }
    out << endl;
  }
}

//===============================================================
#if 0
int main()
{
  // read number of rows n
  int n;
  cin >> n;

  // Allocate n x n matrix for K
  Matrix K(n, n);

  // Read input into K
  for (int i=0; i<n; i++) {
    for (int j=0; j<n; j++) {
      cin >> K[i][j];
    }
  }

  cout << "K:" << endl << K << endl;

  // Compute K inverse
  Matrix KINV(n, n);
  KINV.inverse26(K);
  cout << "K inverse:" << endl << KINV << endl;

  return 0;
}
#endif