
/*
 * Computer Algebra Kit (c) 1993,98 by Comp.Alg.Objects.  All Rights Reserved.
 * $Id: matrix.m,v 1.8 1998/11/13 09:45:33 stes Exp $
 */

#include "cakit.h"

@implementation Matrix
- check
{
  int i, m, n;

  assert(rows && scalarZero && emptyVector);
  [rows check];
  [scalarZero check];
  [emptyVector check];

  assert([scalarZero isEqual:[emptyVector scalarZero]]);
  assert([emptyVector numScalars] == 0);

  m = [self numRows];
  n = [self numColumns];
  for (i = 0; i < m; i++) {
    assert([[[rows at:i] scalarZero] isEqual:scalarZero]);
    assert([[rows at:i] numScalars] == n);
  }

  return [super check];
}


+ diagonal:cltnOfScalars
{
  int i, n = [cltnOfScalars size];

  if (n == 0)
    [self error:"Collection Must Not Be Empty"];

  self = [self scalarZero:[[cltnOfScalars lastElement] zero] numRows:n numColumns:n];
  for (i = 0; i < n; i++)
    [self placeScalar:[[cltnOfScalars at:i] cheapCopy] at:i:i];

  assert([self isDiagonal]);
  return self;
}

+ circulant:cltnOfScalars
{
  int i, j, n = [cltnOfScalars size];

  if (n == 0)
    [self error:"Collection Must Not Be Empty"];

  self = [self scalarZero:[[cltnOfScalars lastElement] zero] numRows:n numColumns:n];
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
      [self placeScalar:[[cltnOfScalars at:(n - i + j) % n] cheapCopy] at:i:j];
    }
  }
  return self;
}

+ companion:cltnOfScalars
{
  id aZero;
  int i, n = [cltnOfScalars size];

  if (n == 0)
    [self error:"Collection Must Not Be Empty"];

  aZero = [[cltnOfScalars lastElement] zero];
  self = [self scalarZero:[aZero cheapCopy] numRows:n numColumns:n];
  for (i = 0; i < n - 1; i++)
    [self placeScalar:[aZero one] at:i + 1:i];
  for (i = 0; i < n; i++)
    [self placeScalar:[[cltnOfScalars at:i] negate] at:i:n - 1];
  [aZero free];
  return self;
}

+ hankel:rowScalars:colScalars
{
  id obj;
  int i, j, k, n, m;

  n = [colScalars size] + 1;
  m = [rowScalars size];
  self = [self scalarZero:[[rowScalars lastElement] zero] numRows:n numColumns:m];

  for (i = 0; i < n; i++) {
    for (j = 0; j < m; j++) {
      k = i + j;
      obj = (k < m) ? [rowScalars at : k]:[colScalars at:k - m];
      [self placeScalar:[obj cheapCopy] at:i:j];
    }
  }

  return self;
}

+ toeplitz:rowScalars:colScalars
{
  id obj;
  int i, j, k, n, m;

  n = [colScalars size];
  m = [rowScalars size] + 1;
  self = [self scalarZero:[[rowScalars lastElement] zero] numRows:n numColumns:m];

  for (i = 0; i < n; i++)
    for (j = 0; j < m; j++) {
      k = j - i + n - 1;
      obj = (k < n) ? [colScalars at : k]:[rowScalars at:k - n];
      [self placeScalar:[obj cheapCopy] at:i:j];
    }

  return self;
}

+ hilbert:(int)n scalarZero:aZero
{
  id obj;
  int i, j, k;

  self = [self scalarZero:aZero numRows:n numColumns:n];

  for (k = 1; k < (n + n); k++) {
    obj = [aZero intValue:k];
    obj = [obj inverseSelf];
    for (i = 0; i < n; i++)
      for (j = 0; j < n; j++) {
	if (i + j + 1 == k)
	  [self placeScalar:[obj cheapCopy] at:i:j];
      }
    [obj free];
  }

  return self;
}

+ hilbert:(int)n
{
  id aZero;

  aZero = [Integer new];
  aZero = [Fraction numerator:aZero denominator:[aZero one]];

  return [self hilbert:n scalarZero:aZero];
}


- _setUpScalarZero:aScalarZero numRows:(int)numRows numColumns:(int)numColumns
{
  if ([aScalarZero notZero])
    [self error:"ScalarZero Not Zero !"];

  scalarZero = aScalarZero;
  emptyVector = [Vector scalarZero:[scalarZero cheapCopy] numScalars:0];

  rows = [CACollection new:numRows];
  while (numRows--)
    [rows add:[emptyVector numScalars:numColumns]];
  return self;
}

- scalarZero:aScalarZero numRows:(int)numRows numColumns:(int)numColumns
{
  return [self notImplemented:_cmd];
}

+ scalarZero:aScalarZero numRows:(int)numRows numColumns:(int)numColumns
{
  return [[super new] _setUpScalarZero:aScalarZero numRows:numRows numColumns:numColumns];
}
- scalar:aScalar numRows:(int)numRows numColumns:(int)numColumns
{
  return [self notImplemented:_cmd];
}

+ scalar:aScalar numRows:(int)numRows numColumns:(int)numColumns
{
  int i = (numRows < numColumns) ? numRows : numColumns;

  self = [self scalarZero:[aScalar zero] numRows:numRows numColumns:numColumns];
  while (i--)
    [self placeScalar:[aScalar cheapCopy] at:i:i];
  /* we own aScalar now */
  [aScalar free];
  return self;
}
- copy
{
  self = [super copy];
  scalarZero = [scalarZero cheapCopy];
  emptyVector = [emptyVector cheapCopy];
  rows = [rows deepCopy];
  assert([self check]);
  return self;
}

- deepCopy
{
  self = [super deepCopy];
  scalarZero = [scalarZero deepCopy];
  emptyVector = [emptyVector deepCopy];
  rows = [rows deepCopy];
  assert([self check]);
  return self;
}

- clone
{
  self = [super clone];
  scalarZero = [scalarZero cheapCopy];
  emptyVector = [emptyVector cheapCopy];
  rows = nil;
  return self;
}

- rows:c
{
  self = [self clone];
  rows = c;
  assert([self check]);
  return self;
}

- scalarZero:aZero emptyVector:aVector rows:c
{
  self = [super clone];
  scalarZero = aZero;
  emptyVector = aVector;
  rows = c;
  assert([self check]);
  return self;
}

- numRows:(int)numRows numColumns:(int)numColumns
{
  id c = [rows capacity:numRows];

  while (numRows--)
    [c add:[emptyVector numScalars:numColumns]];
  return [self rows:c];
}

- emptyVector
{
  return [emptyVector cheapCopy];
}

- clear
{
  scalarZero = [scalarZero free];
  emptyVector = [emptyVector free];
  rows = [rows free];
  return [super clear];
}

- (BOOL)sameClass:b
{
  return isa == [b class] && [scalarZero isEqual:[b scalarZero]];
}

- (BOOL)inAdditiveSemiGroup		 {
  return YES;
}
- (BOOL)inAdditiveMonoid		 {
  return NO;
}
- (BOOL)inAdditiveGroup			 {
  return NO;
}

- (BOOL)inSemiGroup			 {
  return YES;
}
- (BOOL)inMonoid			 {
  return NO;
}
- (BOOL)inGroup				 {
  return NO;
}

- (BOOL)inEuclideanDomain		 {
  return NO;
}
- (BOOL)inIntegralDomain		 {
  return NO;
}
- (BOOL)inField				 {
  return NO;
}

- (BOOL)inOrderedSet			 {
  return NO;
}
- (int)characteristic			 {
  return [scalarZero characteristic];
}


- scalarZero
{
  return scalarZero;
}

- rows
{
  return rows;
}

+ classVectors
{
  return [self doesNotRecognize:_cmd];
}

- (int) numRows
{
  /* don't invoke check here as we call numRows from check */
  return [rows size];
}

- (int) numColumns
{
  /* don't invoke check here as we call numRows from check */
  id row = [rows lastElement];

  return (row) ? [row numScalars] : 0;
}

- checkDimensionEqual:b
{
  if ([self numRows] != [b numRows])
    [self error:"numRows not equal"];
  if ([self numColumns] != [b numColumns])
    [self error:"numColumns not equal"];
  return self;
}

- checkSquare
{
  if ([self numRows] != [self numColumns])
    [self error:"Matrix not Square"];
  return self;
}

- (BOOL) isEqual:aMatrix
{
  if (self == aMatrix) {
    return YES;
  } else {
    int m = [self numRows];

    if (m == [aMatrix numRows]) {
      while (m--)
	if ([[self rowAt:m] notEqual:[aMatrix rowAt:m]])
	  return NO;
      return YES;
    } else {
      return NO;
    }
  }
}

- (unsigned) hash
{
  return ([rows lastElement]) ? [[rows lastElement] hash] : 0;
}

- (BOOL) isDiagonal
{
  id ix;
  BOOL res;
  int i, j, m, n;

  m = [self numRows];
  n = [self numColumns];
  ix = [self eachSequence];
  res = YES;

  for (i = 0; i < m; i++)
    for (j = 0; j < n; j++)
      if (i != j)
	if ([[[ix at:i] at:j] notZero]) {
	  res = NO;
	  break;
	}
  ix = [ix free];
  return res;
}

- (BOOL) isSymmetric
{
  id ix, jx;
  BOOL res;
  int i, j, m, n;

  m = [self numRows];
  n = [self numColumns];
  ix = [self eachSequence];
  jx = [self eachSequence];
  res = YES;

  for (i = 0; i < m; i++)
    for (j = 0; j < i; j++)
      if ([[[ix at:i] at:j] notEqual:[[jx at:j] at:i]]) {
	res = NO;
	break;
      }
  ix = [ix free];
  jx = [jx free];
  return res;
}

- (BOOL) isAntiSymmetric
{
  id ix, jx;
  BOOL res;
  int i, j, m, n;

  m = [self numRows];
  n = [self numColumns];
  ix = [self eachSequence];
  jx = [self eachSequence];
  res = YES;

  for (i = 0; i < m; i++)
    for (j = 0; j < i; j++)
      if ([[[ix at:i] at:j] notOpposite:[[jx at:j] at:i]]) {
	res = NO;
	break;
      }
  ix = [ix free];
  jx = [jx free];
  return res;
}


- insertRow:aVector
{
  id last;

  [self checkRefCountOne];
  last = [rows lastElement];
  if (last != nil && [last numScalars] != [aVector numScalars]) {
    [self error:"Rows Must Have Same Number of Scalars"];
  }
  [rows add:aVector];
  return [self invalidate];
}

- insertRow:aVector at:(int)i
{
  id last;

  [self checkRefCountOne];
  last = [rows lastElement];
  if (last != nil && [last numScalars] != [aVector numScalars]) {
    [self error:"Rows Must Have Same Number of Scalars"];
  }
  [rows at:i insert:aVector];
  return [self invalidate];
}

- insertColumn:aCollection
{
  int m, n;

  [self checkRefCountOne];

  m = [rows size];
  n = [aCollection size];

  if (m) {
    if (m != n)
      [self error:"Columns Must Have Same Number of Scalars"];
  } else {
    /* load enough empty vectors in the matrix */
    while (n--)
      [rows add:[self emptyVector]];
  }

  while (m--)
    [[rows at:m] insertScalar:[[aCollection at:m] cheapCopy]];
  [aCollection free];
  return [self invalidate];
}

- insertColumn:aCollection at:(int)i
{
  int m, n;

  [self checkRefCountOne];

  m = [rows size];
  n = [aCollection size];

  if (m) {
    if (m != n)
      [self error:"Columns Must Have Same Number of Scalars"];
  } else {
    /* load enough empty vectors in the matrix */
    while (n--)
      [rows add:[self emptyVector]];
  }

  while (m--)
    [[rows at:m] insertScalar:[[aCollection at:m] cheapCopy] at:i];
  [aCollection free];
  return [self invalidate];
}


- removeRow
{
  id r;

  [self checkRefCountOne];
  r = [rows removeLast];
  [self invalidate];
  return r;
}

- removeRowAt:(int)i
{
  id r;

  [self checkRefCountOne];
  r = [rows removeAt:i];
  [self invalidate];
  return r;
}

- removeColumn
{
  int j, m = [self numRows];

  if ([self numColumns]) {
    id column = [rows capacity:m];

    [self checkRefCountOne];
    for (j = 0; j < m; j++)
      [column add:[[rows at:j] removeScalar]];
    [self invalidate];
    return column;
  } else {
    return nil;
  }
}

- removeColumnAt:(int)i
{
  int j, m = [self numRows];

  if ([self numColumns]) {
    id column = [rows capacity:m];

    [self checkRefCountOne];
    for (j = 0; j < m; j++)
      [column add:[[rows at:j] removeScalarAt:i]];
    [self invalidate];
    return column;
  } else {
    return nil;
  }
}


- placeScalar:aScalar at:(int)i:(int)j
{
  [self checkRefCountOne];
  [[rows at:i] placeScalar:aScalar at:j];
  return [self invalidate];
}

- replaceScalarAt:(int)i:(int)j with:aScalar
{
  [self checkRefCountOne];
  [[rows at:i] replaceScalarAt:j with:aScalar];
  return [self invalidate];
}


- asNumerical
{
  id new;
  int i, m;
  id aZero, aVector;

  m = [self numRows];
  new = [rows capacity:m];
  for (i = 0; i < m; i++)
    [new add:[[rows at:i] asNumerical]];

  aZero = [scalarZero asNumerical];
  aVector = [emptyVector asNumerical];
  return [self scalarZero:aZero emptyVector:aVector rows:new];
}

- asModp:(unsigned short)p
{
  id new;
  int i, m;
  id aZero, aVector;

  m = [self numRows];
  new = [rows capacity:m];
  for (i = 0; i < m; i++)
    [new add:[[rows at:i] asModp:p]];

  aZero = [scalarZero asModp:p];
  aVector = [emptyVector asModp:p];
  return [self scalarZero:aZero emptyVector:aVector rows:new];
}

- commonDenominator
{
  id d, s = [self eachScalar];

  d = [s commonDenominator];
  [s free];
  return d;
}

- asIntegral
{
  int i, m;
  id aZero, aVector, new, seq;

  m = [self numRows];
  new = [rows capacity:m];
  for (i = 0; i < m; i++)
    [new add:[[rows at:i] asIntegral]];

  aZero = [scalarZero asIntegral];
  aVector = [emptyVector asIntegral];
  return [self scalarZero:aZero emptyVector:aVector rows:new];
}

- asFractional
{
  int i, m;
  id aZero, aVector, new;

  m = [self numRows];
  new = [rows capacity:m];
  for (i = 0; i < m; i++)
    [new add:[[rows at:i] asFractional]];

  aZero = [scalarZero asTotalFraction];
  aVector = [emptyVector asFractional];
  return [self scalarZero:aZero emptyVector:aVector rows:new];
}

- onCommonDenominator:(id *)denominator
{
  id n, d, e;

  d = [self commonDenominator];
  e = [scalarZero numerator:[d cheapCopy] denominator:[d one]];
  n = [self multiplyScalar:e];
  [e free];
  *denominator = d;
  return idReplace(n, [n asIntegral]);
}


- rowAt:(int)i
{
  return [rows at:i];
}

- eachRow
{
  return [rows eachElement];
}

- eachScalar
{
  return [CASequence over:[CAMatrixSequenceScalars content:[self cheapCopy]]];
}

- eachSequence
{
  return [CASequence over:[CAMatrixSequenceSequences content:[self cheapCopy]]];
}

- (float) floatValueAt:(int)i:(int)j
{
  return [[rows at:i] floatValueAt:j];
}


- zero
{
  id c;
  int i, m;

  c = [rows capacity:m = [self numRows]];
  for (i = 0; i < m; i++)
    [c add:[[self rowAt:i] zero]];
  return [self rows:c];
}

- negate
{
  id c;
  int i, m;

  c = [rows capacity:m = [self numRows]];
  for (i = 0; i < m; i++)
    [c add:[[self rowAt:i] negate]];
  return [self rows:c];
}

- double
{
  id c;
  int i, m;

  c = [rows capacity:m = [self numRows]];
  for (i = 0; i < m; i++)
    [c add:[[self rowAt:i] double]];
  return [self rows:c];
}

- add:b
{
  if (self == b) {
    return [self double];
  } else {
    id c;
    int i, m;

    [self checkDimensionEqual:b];
    c = [rows capacity:m = [self numRows]];
    for (i = 0; i < m; i++)
      [c add:[[self rowAt:i] add:[b rowAt:i]]];
    return [self rows:c];
  }
}

- subtract:b
{
  if (self == b) {
    return [self double];
  } else {
    id c;
    int i, m;

    [self checkDimensionEqual:b];
    c = [rows capacity:m = [self numRows]];
    for (i = 0; i < m; i++)
      [c add:[[self rowAt:i] subtract:[b rowAt:i]]];
    return [self rows:c];
  }
}

- addScalar:s
{
  if ([s isZero]) {
    return [self cheapCopy];
  } else {
    int i, m;
    id c;

    c = [rows capacity:m = [self numRows]];
    for (i = 0; i < m; i++)
      [c add:[[self rowAt:i] addScalar:s at:i]];
    return [self rows:c];
  }
}

- subtractScalar:s
{
  if ([s isZero]) {
    return [self cheapCopy];
  } else {
    int i, m;
    id c;

    c = [rows capacity:m = [self numRows]];
    for (i = 0; i < m; i++)
      [c add:[[self rowAt:i] subtractScalar:s at:i]];
    return [self rows:c];
  }
}


- one
{
  int m = [self numRows];
  int n = [self numColumns];

  self = [self zero];
  n = (n < m) ? n : m;
  while (n--)
    [self placeScalar:[scalarZero one] at:n:n];
  return self;
}

- square
{
  return [self multiply:self];
}

- multiply:b
{
  id c, bt;
  int i, j, m, n;

  assert([b check] && [self check]);

  m = [self numRows];
  n = [b numColumns];
  if ([self numColumns] != [b numRows]) {
    [self error:"-matrix multiply: numColumns!=numRows"];
  }
  bt = [b transpose];
  c = [self numRows:m numColumns:n];

  for (i = 0; i < m; i++) {
    id ri = [self rowAt:i];

    for (j = 0; j < n; j++) {
      [c placeScalar:[ri dotMultiply:[bt rowAt:j]] at:i:j];
    }
  }

  [bt free];
  return c;
}

- multiplyVector:aColumn
{
  return [aColumn multiplyLeftMatrix:self];
}


- multiplyScalar:b
{
  id c;
  int i, m;

  c = [rows capacity:m = [self numRows]];
  for (i = 0; i < m; i++)
    [c add:[[self rowAt:i] multiplyScalar:b]];
  return [self rows:c];
}


- divideScalar:b
{
  id c;
  int i, m;

  c = [rows capacity:m = [self numRows]];
  for (i = 0; i < m; i++) {
    id r = [[self rowAt:i] divideScalar:b];

    if (r)
      [c add:r];
    else {
      [c free];
      return nil;
    }
  }
  return [self rows:c];
}


- transpose
{
  id c, seq;
  int i, m = [self numRows];
  int j, n = [self numColumns];

  seq = [self eachSequence];
  c = [self numRows:n numColumns:m];

  for (i = 0; i < m; i++) {
    for (j = 0; j < n; j++) {
      [c placeScalar:[[[seq at:i] at:j] cheapCopy] at:j:i];
    }
  }

  [seq free];
  return c;
}
- transposeSelf
{
  return idReplace(self, [self transpose]);
}

- determinant
{
  id det;

  [self checkSquare];

  if ([scalarZero isFloatingPoint]) {
    det = [self error:"No Determinant Over Floating Points Yet !"];
  } else if ([scalarZero inFieldOfFractions]) {
    id n, d, m;

    m = [self onCommonDenominator:&d];
    d = [d powerSelf:[self numRows]];
    n = [m determinant];
    [m free];
    return [scalarZero numerator:n denominator:d];
  } else if ([scalarZero inField]) {
    det = [self determinantGauss];
  } else if ([scalarZero inIntegralDomain]) {
    det = [self determinantBareiss];
  } else {
    det = [self notImplemented:_cmd];
  }

  return det;
}

- determinantGauss
{
  int i, j, n, sign;
  id ri, rj, rn, p, q, det;

  /* create a scalar that will evolve into the determinant */

  det = [scalarZero one];
  rn = [rows deepCopy];
  n = [rn size];

  /* iterate over the rows */

  for (sign = +1, ri = nil, p = nil, i = 0; i < n; i++) {

    /* for the first row, search pivot in first column
     * for second row, search in second column etc
     */

    for (j = i; j < n; j++) {
      ri = [rn at:j];
      p = [ri removeScalar];
      if ([p isZero])
	p = [p free];
      else
	break;
    }

    /* quit if no pivot or swap pivot onto diagonal */

    if (j == n) {
      [det free];
      [rn free];
      return [scalarZero zero];
    }
    if (j != i) {
      sign = -sign;
      [rn swapElementsAt:i:j];
    }
    det = [det multiplySelf:p];
    p = [p inverseSelf];
    assert(p);

    /* reduction step - modify behind back of "rn" member at "j"
     * so don't forget to update the reference at position "j"
     */

    for (j++; j < n; j++) {
      rj = [rn at:j];
      q = [rj removeScalar];
      q = [q multiplySelf:p];
      rj = [rj subtractSelf:ri multiplyScalar:q];
      [rn at:j put:rj];
      q = [q free];
    }

    assert([rn check]);
    p = [p free];
  }

  [rn free];
  det = (sign > 0) ? det : [det negateSelf];
  return det;
}

- determinantBareiss
{
  int i, j, n, sign;
  id ri, rj, rn, p, q, det;

  /* create a scalar that will evolve into the determinant */

  det = [scalarZero one];
  rn = [rows deepCopy];
  n = [rn size];

  /* iterate over the rows */

  for (sign = +1, ri = nil, p = nil, i = 0; i < n; i++) {

    /* for the first row, search pivot in first column
     * for second row, search in second column etc
     */

    for (j = i; j < n; j++) {
      ri = [rn at:j];
      p = [ri removeScalar];
      if ([p isZero])
	p = [p free];
      else
	break;
    }

    /* quit if no pivot or swap pivot onto diagonal
     * unlike Gauss, have to modify rows with leading zero
     */

    if (j == n) {
      [det free];
      [rn free];
      return [scalarZero zero];
    }
    if (j != i) {
      int k;

      sign = -sign;
      for (k = i; k < j; k++) {
	id rk = [rn at:k];

	rk = [rk multiplySelfScalar:p divideScalar:det];
	[rn at:k put:rk];
      }
      [rn swapElementsAt:i:j];
    }
    /* reduction step - modify behind back of "rn" member at "j"
     * so don't forget to update the reference at position "j"
     * multiply by current pivot - divide by previous pivot
     */

    for (j++; j < n; j++) {
      rj = [rn at:j];
      q = [rj removeScalar];
      rj = [rj multiplySelfScalar:p subtract:ri multiplyScalar:q divideScalar:det];
      [rn at:j put:rj];
      q = [q free];
    }

    assert([rn check]);
    det = [det free];
    det = p;
  }

  [rn free];
  det = (sign > 0) ? det : [det negateSelf];
  return det;
}

- solveVector:y
{
  [emptyVector checkSameClass:y];
  [self checkSquare];

  if ([y numScalars] != [self numRows]) {
    [self error:"Dimensions Not Compatible"];
  } else {
    id matrix;

    /* append column vector to matrix */

    matrix = [self copy];
    [matrix insertColumn:[y asCollection] at:[matrix numColumns]];

    if ([scalarZero inFieldOfFractions]) {
      id x = [matrix solveGaussDeterminant:NULL];

      [matrix free];
      return x;
    } else if ([scalarZero inField]) {
      id x = [matrix solveGaussDeterminant:NULL];

      [matrix free];
      return x;
    } else if ([scalarZero inIntegralDomain]) {
      id det, x;

      x = [matrix solveBareissDeterminant:&det];
      x = [x divideSelfScalar:det];
      [det free];
      [matrix free];
      return x;
    }
  }

  return [self notImplemented:_cmd];
}

- solveGaussDeterminant:(id *)detReference
{
  int i, j, n, sign;
  id ri, rj, p, q, rn, det, x;

  n = [self numRows];
  assert(n + 1 == [self numColumns]);
  rn = [rows deepCopy];

  det = (detReference) ? [scalarZero one] : nil;
  for (sign = (+1), ri = nil, p = nil, i = 0; i < n; i++) {

    for (j = i; j < n; j++) {
      ri = [rn at:j];
      p = [ri removeScalar];
      if ([p isZero])
	p = [p free];
      else
	break;
    }

    if (j == n) {
      [self error:"Matrix Singular"];
    }
    if (j != i) {
      sign = -sign;
      [rn swapElementsAt:i:j];
    }
    if (detReference)
      det = [det multiplySelf:p];
    p = [p inverseSelf];
    assert(p);
    ri = [ri multiplySelfScalar:p];
    [rn at:i put:ri];
    p = [p free];

    for (j++; j < n; j++) {
      rj = [rn at:j];
      q = [rj removeScalar];
      rj = [rj subtractSelf:ri multiplyScalar:q];
      [rn at:j put:rj];
      q = [q free];
    }

    assert([rn check]);
  }
  if (detReference)
    det = (sign > 0) ? det : [det negateSelf];

  x = [emptyVector capacity:n + 1];
  x = [x insertScalar:[scalarZero minusOne]];
  i = n;
  while (i--) {
    q = [x dotMultiply:[rn at:i]];
    [x insertScalar:[q negateSelf]];
  }

  rn = [rn free];

  if (detReference)
    *detReference = det;
  else
    [det free];
  return [x freeScalarAt:n];
}

- solveBareissDeterminant:(id *)detReference
{
  int i, j, n, sign;
  id ri, rj, rk, p, q, piv, rn, det, x;

  n = [self numRows];
  assert(n + 1 == [self numColumns]);
  rn = [rows deepCopy];

  det = [scalarZero one];
  piv = [rows capacity:n];

  for (sign = (+1), ri = nil, p = nil, i = 0; i < n; i++) {

    for (j = i; j < n; j++) {
      ri = [rn at:j];
      p = [ri removeScalar];
      if ([p isZero])
	p = [p free];
      else
	break;
    }

    if (j == n) {
      [self error:"Matrix Singular"];
    }
    if (j != i) {
      int k;

      for (k = i; k < j; k++) {
	rk = [[rn at:k] multiplySelfScalar:p divideScalar:det];
	[rn at:k put:rk];
      }
      sign = -sign;
      [rn swapElementsAt:i:j];
    }
    for (j++; j < n; j++) {
      rj = [rn at:j];
      q = [rj removeScalar];
      rj = [rj multiplySelfScalar:p subtract:ri multiplyScalar:q divideScalar:det];
      [rn at:j put:rj];
      q = [q free];
    }

    [piv add:[p cheapCopy]];
    assert([rn check]);
    det = [det free];
    det = p;
  }

  det = (sign > 0) ? det : [det negateSelf];

  x = [emptyVector capacity:n + 1];
  p = [piv removeLast];
  if (sign > 0) {
    assert([p isEqual:det]);
    x = [x insertScalar:[p negateSelf]];
    x = [x insertScalar:[[rn at:n - 1] removeScalar]];
  } else {
    assert([p isOpposite:det]);
    x = [x insertScalar:p];
    x = [x insertScalar:[[[rn at:n - 1] removeScalar] negateSelf]];
  }
  i = n - 1;
  while (i--) {
    q = [x dotMultiply:[rn at:i]];
    p = [piv removeLast];
    q = [[q divideSelf:p] negateSelf];
    [x insertScalar:q];
    p = [p free];
  }

  piv = [piv free];
  rn = [rn free];

  if (detReference)
    *detReference = det;
  else
    [det free];
  return [x freeScalarAt:n];
}

- inverse
{
  [self checkSquare];
  return [[self one] divideSelf:self];
}

- divide:b
{
  if ([scalarZero inField])
    return [self divideGauss:b];
  return [self notImplemented:_cmd];
}

- divideGauss:b
{
  int i, j, n;
  id p, ri, rj, iri, irj, q, t, rn, irn, triangle;

  n = [self numRows];
  assert(n == [self numColumns] && n == [b numRows] && n == [b numColumns]);
  rn = [[b rows] deepCopy];
  irn = [rows deepCopy];
  triangle = [rows capacity:n];

  /*          - - -                 1 | |
   *     rn   - - -   ->  triangle    1 | 
   *          - - -                     1
   */

  for (ri = nil, p = nil, i = 0; i < n; i++) {

    t = [rn capacity:n];
    for (j = 0; j < i; j++) {
      rj = [rn at:j];
      [t add:[rj removeScalar]];
      assert([[rn at:j] checkRefCountOne]);
    }

    for (j = i; j < n; j++) {
      ri = [rn at:j];
      p = [ri removeScalar];
      if ([p isZero])
	p = [p free];
      else
	break;
    }

    if (j == n) {
      return nil;
    }
    if (j != i) {
      [rn swapElementsAt:i:j];
      [irn swapElementsAt:i:j];
    }
    p = [p inverseSelf];
    assert(p);
    iri = [irn at:i];
    ri = [ri multiplySelfScalar:p];
    iri = [iri multiplySelfScalar:p];
    [rn at:i put:ri];
    [irn at:i put:iri];
    p = [p free];

    for (j++; j < n; j++) {
      rj = [rn at:j];
      q = [rj removeScalar];
      rj = [rj subtractSelf:ri multiplyScalar:q];
      [rn at:j put:rj];
      irj = [irn at:j];
      irj = [irj subtractSelf:iri multiplyScalar:q];
      [irn at:j put:irj];
      q = [q free];
    }

    [triangle add:t];
    assert([rn check]);
    assert([irn check]);
  }

  while (i--) {
    t = [triangle at:i];
    iri = [irn at:i];
    for (j = i - 1; j >= 0; j--) {
      q = [t removeLast];
      irj = [irn at:j];
      irj = [irj subtractSelf:iri multiplyScalar:q];
      [irn at:j put:irj];
      q = [q free];
    }
  }

  triangle = [triangle free];
  rn = [rn free];
  return [self rows:irn];
}

- (int) rank
{
  BOOL *rowmark;
  id q, qj, ri, rk, rm;
  int i, j, k, m, n, rank;

  m = [self numRows];
  n = [self numColumns];
  rm = [rows deepCopy];
  qj = [rows capacity:m];

  rowmark = (BOOL *) OC_Malloc(m * sizeof(BOOL));

  for (i = 0; i < m; i++)
    rowmark[i] = NO;

  for (q = nil, ri = nil, rank = n, j = 0; j < n; j++) {

    for (i = 0; i < m; i++) {
      ri = [rm at:i];
      q = [ri removeScalar];
      if (rowmark[i] == NO && [q notZero])
	break;
      else
	[qj add:q];
    }

    if (i == m) {
      rank--;
    } else {
      rowmark[i] = YES;

      q = [q inverseSelf];
      assert(q);
      ri = [ri multiplySelfScalar:q];
      [rm at:i put:ri];
      q = [q free];

      for (k = i - 1; k >= 0; k--) {
	q = [qj removeLast];
	rk = [rm at:k];
	[rm at:k put:[rk subtractSelf:ri multiplyScalar:q]];
	q = [q free];
      }

      for (k = i + 1; k < m; k++) {
	rk = [rm at:k];
	q = [rk removeScalar];
	[rm at:k put:[rk subtractSelf:ri multiplyScalar:q]];
	q = [q free];
      }
    }

  }

  OC_Free(rowmark);
  qj = [qj free];
  rm = [rm free];
  return rank;
}

- (int) nullity
{
  return [self numColumns] - [self rank];
}

- rowReducedFormRank:(int *)rankReference indexRank:(int**)indexRankReference
{
  BOOL *rowmark;
  id q, qj, ri, rk, rm, rrf;
  int i, j, k, m, n, rank, *indexrank;

  m = [self numRows];
  n = [self numColumns];
  rm = [rows deepCopy];

  rrf = [rows capacity:n];
  rowmark = (BOOL *) OC_Malloc(sizeof(BOOL) * m);

  for (i = 0; i < m; i++)
    rowmark[i] = NO;
  indexrank = (int *) OC_Malloc(sizeof(int) * n);

  for (q = nil, ri = nil, rank = n, j = 0; j < n; j++) {

    qj = [rows capacity:m];

    for (i = 0; i < m; i++) {
      ri = [rm at:i];
      q = [ri removeScalar];
      if (rowmark[i] == NO && [q notZero]) {
	break;
      } else {
	[qj add:q];
      }
    }

    if (i == m) {
      indexrank[j] = -1;
      rank--;
    } else {
      indexrank[j] = i;
      rowmark[i] = YES;

      q = [q inverseSelf];
      assert(q);
      ri = [ri multiplySelfScalar:q];
      [rm at:i put:ri];
      q = [q free];
      [qj add:[scalarZero one]];

      for (k = 0; k < i; k++) {
	q = [qj at:k put:[scalarZero zero]];
	rk = [[rm at:k] subtractSelf:ri multiplyScalar:q];
	[rm at:k put:rk];
	q = [q free];
      }

      for (k = i + 1; k < m; k++) {
	rk = [rm at:k];
	q = [rk removeScalar];
	rk = [rk subtractSelf:ri multiplyScalar:q];
	[rm at:k put:rk];
	q = [q free];
	[qj add:[scalarZero zero]];
      }
    }

    [rrf add:qj];
  }

  OC_Free(rowmark);
  *rankReference = rank;
  *indexRankReference = indexrank;
  return rrf;
}

- kernel
{
  id q, c, ker, rrf;
  int i, ix, j, k, n, rank, nullity, *indexrank;

  rrf = [self rowReducedFormRank:&rank indexRank:&indexrank];
  n = [self numColumns];
  nullity = n - rank;

  ker = [rows capacity:nullity];

  for (j = k = 1; j < nullity; j++, k++) {
    c = [emptyVector numScalars:n];
    while (indexrank[k] >= 0)
      k++;
    for (i = 0; i < k; i++) {
      ix = indexrank[i];
      if (ix < 0) {
	q = [scalarZero zero];
      } else {
	q = [[[rrf at:i] at:ix] negate];
      }
      [c placeScalar:q at:i];
    }
    [ker add:c];
  }

  OC_Free(indexrank);
  rrf = [rrf free];
  return ker;
}

- image
{
  id qk, c, image, rrf;
  int i, j, k, m, rank, *indexrank;

  rrf = [self rowReducedFormRank:&rank indexRank:&indexrank];
  m = [self numRows];

  image = [rows capacity:rank];

  for (j = k = 1; j < rank; j++, k++) {
    while (indexrank[k] < 0)
      k++;
    qk = [rrf at:k];
    c = [emptyVector numScalars:m];
    i = m;
    while (i--)
      [c placeScalar:[[qk at:i] cheapCopy] at:i];
    [image add:c];
  }

  OC_Free(indexrank);
  rrf = [rrf free];
  return image;
}


- trace
{
  int n;
  id seq, trace;

  [self checkSquare];

  trace = [scalarZero zero];
  n = [self numRows];

  seq = [self eachSequence];
  while (n--)
    trace = [trace addSelf:[[seq at:n] at:n]];
  seq = [seq free];

  return trace;
}

- adjoint
{
  int k, n;
  id adj, c;

  n = [self numRows];
  k = [scalarZero characteristic];
  if (0 < k && k < n)
    [self error:"-adjoint: Characteristic to Small For This Implementation"];

  c = [self trace];
  c = idReplace(c, [c negate]);
  adj = [self addScalar:c];
  c = [c free];

  for (k = 2; k < n; k++) {
    assert([adj check]);
    adj = idReplace(adj, [self multiply:adj]);
    c = [adj trace];
    c = idReplace(c, [c divideIntValue:-k]);
    adj = idReplace(adj, [adj addScalar:c]);
    c = [c free];
  }

  return adj;
}

- adjointSelf
{
  return idReplace(self, [self adjoint]);
}


- printOn:(IOD)aFile
{
  [rows printOn:aFile];
  return self;
}

@end

