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

#include "cakit.h"

@implementation Vector
- check
{
  assert(carrier && scalarZero);
  [carrier check];
  [scalarZero check];
  assert([scalarZero isZero]);
  assert([[carrier scalarZero] isEqual:scalarZero]);
  return [super check];
}


+ collection:aCltn
{
  if ([aCltn isEmpty]) {
    return [self error:"-collection: cltn must not be empty"];
  } else {
    id res;
    int n = [aCltn size];

    res = [self scalarZero:[[aCltn lastElement] zero] numScalars:n];
    while (n--)
      [res insertScalar:[[aCltn at:n] cheapCopy]];
    return res;
  }
}

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

- _setUpScalarZero:aScalarZero numScalars:(int)numScalars
{
  id classCarrier;

  scalarZero = [aScalarZero cheapCopy];
  classCarrier = [aScalarZero class_vector];
  carrier = [classCarrier scalarZero:aScalarZero numScalars:numScalars];

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

+ scalarZero:aScalarZero numScalars:(int)numScalars
{
  return [[super new] _setUpScalarZero:aScalarZero numScalars:numScalars];
}

- copy
{
  self = [super copy];
  carrier = [carrier copy];
  scalarZero = [scalarZero cheapCopy];
  return self;
}

- deepCopy
{
  self = [super deepCopy];
  carrier = [carrier deepCopy];
  scalarZero = [scalarZero deepCopy];
  return self;
}

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

- numScalars:(int)numScalars
{
  return [self over:[carrier numScalars:numScalars]];
}

- capacity:(int)aCapacity
{
  return [self over:[carrier capacity:aCapacity]];
}

- collection:aCltn
{
  int n = [aCltn size];

  /* the order of members in a vector/collection is opposite */

  self = [self capacity:n];
  while (n--)
    [self insertScalar:[[aCltn at:n] cheapCopy]];

  [aCltn free];
  return self;
}

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

- over:aCarrier
{
  self = [self clone];
  carrier = aCarrier;
  assert([self check]);
  return self;
}

- (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;
}

- (int) numScalars
{
  return [carrier numScalars];
}

- carrier
{
  return carrier;
}

- (unsigned) hash
{
  return [carrier hash];
}

- (BOOL) isEqual:b
{
  return (self == b) ? YES : [carrier isEqual:[b carrier]];
}

- (BOOL) notEqual:b
{
  return (self == b) ? NO : [carrier notEqual:[b carrier]];
}


- insertScalar:aScalar
{
  if (aScalar) {
    [self checkRefCountOne];
    [carrier insertScalar:aScalar];
    return [self invalidate];
  } else {
    return nil;
  }
}

- insertScalar:aScalar at:(int)i
{
  if (aScalar) {
    [self checkRefCountOne];
    [carrier insertScalar:aScalar at:i];
    return [self invalidate];
  } else {
    return self;
  }
}


- removeScalar
{
  id s;

  [self checkRefCountOne];
  s = [carrier removeScalar];
  [self invalidate];
  return s;
}

- removeScalarAt:(int)i
{
  id s;

  [self checkRefCountOne];
  s = [carrier removeScalarAt:i];
  [self invalidate];
  return s;
}

- freeScalarAt:(int)i
{
  [[self removeScalarAt:i] free];
  return self;
}


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

- replaceScalarAt:(int)i with:aScalar
{
  id s;

  [self checkRefCountOne];
  s = [carrier replaceScalarAt:i with:aScalar];
  [self invalidate];
  return s;
}

- asCollection
{
  id c, seq, scalar;

  c = [CACollection new];

  seq = [self eachScalar];
  while (scalar = [seq next])
    [c add:[scalar cheapCopy]];
  seq = [seq free];

  assert([c size] == [self numScalars]);
  return c;
}

- asNumerical
{
  int n;
  id v, seq;

  n = [self numScalars];
  v = [Vector scalarZero:[scalarZero asNumerical] numScalars:n];

  seq = [self eachScalar];
  while (n--)
    [v placeScalar:[[seq at:n] asNumerical] at:n];
  seq = [seq free];

  return v;
}

- asModp:(unsigned short)p
{
  int n;
  id v, seq;

  n = [self numScalars];
  v = [Vector scalarZero:[scalarZero asModp:p] numScalars:n];

  seq = [self eachScalar];
  while (n--) {
    id s = [[seq at:n] asModp:p];

    if (s == nil)
      [self error:"asModp: modulo failed"];
    [v placeScalar:s at:n];
  }
  seq = [seq free];

  return v;
}

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

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

- 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]);
}

- asIntegral
{
  int n;
  id v, seq;

  n = [self numScalars];
  v = [Vector scalarZero:[scalarZero asIntegral] numScalars:n];

  seq = [self eachScalar];
  while (n--)
    [v placeScalar:[[seq at:n] asIntegral] at:n];
  seq = [seq free];

  return v;
}

- asFractional
{
  int n;
  id v, seq;

  n = [self numScalars];
  v = [Vector scalarZero:[scalarZero asTotalFraction] numScalars:n];

  seq = [self eachScalar];
  while (n--)
    [v placeScalar:[[seq at:n] asTotalFraction] at:n];
  seq = [seq free];

  return v;
}


- eachScalar
{
  return [carrier eachScalar];
}

- (float) floatValueAt:(int)i
{
  return [carrier floatValueAt:i];
}

- (int) intValueAt:(int)i
{
  return [carrier intValueAt:i];
}


- zero
{
  return [self numScalars:[self numScalars]];
}

- (BOOL) isZero
{
  return [carrier isZero];
}

- (BOOL) notZero
{
  return [carrier notZero];
}

- (BOOL) isOpposite:b
{
  [self checkSameClass:b];

  if (self == b) {
    return [self isZero] || [self isCharacteristicTwo];
  } else {
    return [carrier isOpposite:[b carrier]];
  }
}

- (BOOL) notOpposite:b
{
  if (self == b) {
    return [self notZero] && [self notCharacteristicTwo];
  } else {
    return [carrier notOpposite:[b carrier]];
  }
}

- negate
{
  if ([self isZero] || [self isCharacteristicTwo]) {
    return [self cheapCopy];
  } else {
    return [self over:[carrier negate]];
  }
}

- negateSelf
{
  if (refCount > 1)
    return idReplace(self, [self negate]);

  if ([self isZero] || [self isCharacteristicTwo]) {
    return self;
  } else {
    carrier = [carrier negateSelf];
    return [self invalidate];
  }
}

- _double:(int)v
{
  if ([self isZero] || [self isCharacteristicTwo]) {
    return [self zero];
  } else {
    return [self over:[carrier _double:v]];
  }
}

- _doubleSelf:(int)v
{
  if (refCount > 1)
    return idReplace(self, [self _double:v]);

  if ([self isZero] || [self isCharacteristicTwo]) {
    return idReplace(self, [self zero]);
  } else {
    carrier = [carrier _doubleSelf:v];
    return [self invalidate];
  }
}

- double
{
  if ([self isZero] || [self isCharacteristicTwo]) {
    return [self zero];
  } else {
    return [self over:[carrier double]];
  }
}

- doubleSelf
{
  if (refCount > 1)
    return idReplace(self, [self double]);

  if ([self isZero] || [self isCharacteristicTwo]) {
    return idReplace(self, [self zero]);
  } else {
    carrier = [carrier doubleSelf];
    return [self invalidate];
  }
}

- _add:(int)v:b:(int)w
{
  assert([self sameClass:b] && v * v == 1 && w * w == 1);

  if (self == b) {
    return (v == w) ? [self _double : v]:[self zero];
  } else {
    if ([self inAdditiveMonoid]) {
      if ([self isZero])
	return [b _negate:w];
      if ([b isZero])
	return [self _negate:v];
    }
    return [self over:[carrier _add:v:[b carrier]:w]];
  }
}

- _addSelf:(int)v:b:(int)w
{
  assert([self sameClass:b] && v * v == 1 && w * w == 1);
  if (refCount > 1)
    return idReplace(self, [self _add:v:b:w]);

  if (self == b) {
    return (v == w) ? [self _doubleSelf : v]:idReplace(self, [self zero]);
  } else {
    if ([self inAdditiveMonoid]) {
      if ([b isZero])
	return [self _negateSelf:v];
      if ([self isZero])
	return idReplace(self, [b _negate:w]);
    }
    carrier = [carrier _addSelf:v:[b carrier]:w];
    return [self invalidate];
  }
}

- add:b
{
  [self checkSameClass:b];

  if (self == b) {
    return [self double];
  } else {
    if ([self inAdditiveMonoid]) {
      if ([self isZero])
	return [b cheapCopy];
      if ([b isZero])
	return [self cheapCopy];
    }
    return [self over:[carrier add:[b carrier]]];
  }
}

- addSelf:b
{
  assert([self sameClass:b]);
  if (refCount > 1)
    return idReplace(self, [self add:b]);

  if (self == b) {
    return [self doubleSelf];
  } else {
    if ([self inAdditiveMonoid]) {
      if ([b isZero])
	return self;
      if ([self isZero])
	return idReplace(self, [b cheapCopy]);
    }
    carrier = [carrier addSelf:[b carrier]];
    return [self invalidate];
  }
}

- subtract:b
{
  [self checkSameClass:b];

  if (self == b) {
    return [self zero];
  } else {
    if ([self inAdditiveMonoid]) {
      if ([self isZero])
	return [b negate];
      if ([b isZero])
	return [self cheapCopy];
    }
    return [self over:[carrier subtract:[b carrier]]];
  }
}

- subtractSelf:b
{
  assert([self sameClass:b]);
  if (refCount > 1)
    return idReplace(self, [self subtract:b]);

  if (self == b) {
    return idReplace(self, [self zero]);
  } else {
    if ([self inAdditiveMonoid]) {
      if ([b isZero])
	return self;
      if ([self isZero])
	return idReplace(self, [b negate]);
    }
    carrier = [carrier subtractSelf:[b carrier]];
    return [self invalidate];
  }
}

- addScalar:s at:(int)i
{
  if (0 <= i && i < [self numScalars]) {
    id sum, aSequence = [self eachScalar];

    sum = [[aSequence at:i] add:s];
    aSequence = [aSequence free];

    self = [self copy];
    [self placeScalar:sum at:i];
    return self;
  } else {
    return [self error:"-addScalar:at: index out of bounds"];
  }
}
- subtractScalar:s at:(int)i
{
  if (0 <= i && i < [self numScalars]) {
    id sum, aSequence = [self eachScalar];

    sum = [[aSequence at:i] subtract:s];
    aSequence = [aSequence free];

    self = [self copy];
    [self placeScalar:sum at:i];
    return self;
  } else {
    return [self error:"-addScalar:at: index out of bounds"];
  }
}


- multiplyScalar:s
{
  assert([scalarZero sameClass:s]);

  if ([s isZero])
    return [self zero];
  if ([s isOne])
    return [self cheapCopy];
  if ([s isMinusOne])
    return [self negate];

  return [self over:[carrier multiplyScalar:s]];
}

- multiplySelfScalar:s
{
  assert([scalarZero sameClass:s]);
  if (refCount > 1)
    return idReplace(self, [self multiplyScalar:s]);

  if ([s isZero])
    return idReplace(self, [self zero]);
  if ([s isOne])
    return self;
  if ([s isMinusOne])
    return [self negateSelf];

  carrier = [carrier multiplySelfScalar:s];
  return [self invalidate];
}

- divideScalar:s
{
  id tmp;

  assert([scalarZero sameClass:s]);

  if ([s isZero])
    return nil;
  if ([s isOne])
    return [self cheapCopy];
  if ([s isMinusOne])
    return [self negate];

  if (tmp = [carrier divideScalar:s]) {
    return [self over:tmp];
  } else {
    return nil;
  }
}

- divideSelfScalar:s
{
  id tmp;

  assert([scalarZero sameClass:s]);
  if (refCount > 1)
    return idReplace(self, [self divideScalar:s]);

  if ([s isZero])
    return nil;
  if ([s isOne])
    return self;
  if ([s isMinusOne])
    return [self negateSelf];

  if (tmp = [carrier divideSelfScalar:s]) {
    carrier = tmp;
    return [self invalidate];
  } else {
    return nil;
  }
}

- _add:(int)v:B multiplyScalar:b:(int)w
{
  assert([self sameClass:B] && [scalarZero sameClass:b] && v * v == 1 && w * w == 1);

  if ([b isZero])
    return [self _negate:v];
  if ([b isOne])
    return [self _add:v:B:+w];
  if ([b isMinusOne])
    return [self _add:v:B:-w];

  return [self over:[carrier _add:v:[B carrier] multiplyScalar:b:w]];
}

- _addSelf:(int)v:B multiplyScalar:b:(int)w
{
  if (refCount > 1)
    return idReplace(self, [self _add:v:B multiplyScalar:b:w]);
  assert([self sameClass:B] && [scalarZero sameClass:b] && v * v == 1 && w * w == 1);

  if ([b isZero])
    return [self _negateSelf:v];
  if ([b isOne])
    return [self _addSelf:v:B:+w];
  if ([b isMinusOne])
    return [self _addSelf:v:B:-w];

  carrier = [carrier _addSelf:v:[B carrier] multiplyScalar:b:w];
  return [self invalidate];
}

- _multiplyScalar:a:(int)v add:B:(int)w
{
  assert([scalarZero sameClass:a] && [self sameClass:B] && v * v == 1 && w * w == 1);

  if ([a isZero])
    return [B _negate:+w];
  if ([a isOne])
    return [self _add:+v:B:+w];
  if ([a isMinusOne])
    return [self _add:-v:B:+w];

  return [self over:[carrier _multiplyScalar:a:v add:[B carrier]:w]];
}

- _multiplySelfScalar:a:(int)v add:B:(int)w
{
  if (refCount > 1)
    return idReplace(self, [self _multiplyScalar:a:v add:B:w]);
  assert([scalarZero sameClass:a] && [self sameClass:B] && v * v == 1 && w * w == 1);

  if ([a isZero])
    return idReplace(self, [B _negate:+w]);
  if ([a isOne])
    return [self _addSelf:+v:B:+w];
  if ([a isMinusOne])
    return [self _addSelf:-v:B:+w];

  carrier = [carrier _multiplySelfScalar:a:v add:[B carrier]:w];
  return [self invalidate];
}

- _multiplyScalar:a:(int)v add:B multiplyScalar:b:(int)w
{
  assert([scalarZero sameClass:a] && [self sameClass:B] && [scalarZero sameClass:b] && v * v == 1 && w * w == 1);

  if ([a isZero])
    return [B _multiplyScalar:b:+w];
  if ([b isZero])
    return [self _multiplyScalar:a:+v];
  if ([a isOne])
    return [self _add:+v:B multiplyScalar:b:+w];
  if ([b isOne])
    return [self _multiplyScalar:a:+v add:B:+w];
  if ([a isMinusOne])
    return [self _add:-v:B multiplyScalar:b:+w];
  if ([b isMinusOne])
    return [self _multiplyScalar:a:+v add:B:-w];

  return [self over:[carrier _multiplyScalar:a:v add:[B carrier] multiplyScalar:b:w]];
}

- _multiplySelfScalar:a:(int)v add:B multiplyScalar:b:(int)w
{
  if (refCount > 1)
    return idReplace(self, [self _multiplyScalar:a:v add:B multiplyScalar:b:w]);
  assert([scalarZero sameClass:a] && [self sameClass:B] && [scalarZero sameClass:b] && v * v == 1 && w * w == 1);

  if ([a isZero])
    return idReplace(self, [B _multiplyScalar:b:+w]);
  if ([b isZero])
    return [self _multiplySelfScalar:a:+v];
  if ([a isOne])
    return [self _addSelf:+v:B multiplyScalar:b:+w];
  if ([b isOne])
    return [self _multiplySelfScalar:a:+v add:B:+w];
  if ([a isMinusOne])
    return [self _addSelf:-v:B multiplyScalar:b:+w];
  if ([b isMinusOne])
    return [self _multiplySelfScalar:a:+v add:B:-w];

  carrier = [carrier _multiplySelfScalar:a:v add:[B carrier] multiplyScalar:b:w];
  return [self invalidate];
}

- dotSquare
{
  return [carrier dotSquare];
}

- dotMultiply:aVector
{
  return [carrier dotMultiply:[aVector carrier]];
}

- multiplyLeftMatrix:aMatrix
{
  id v;
  int n;

  n = [aMatrix numRows];
  v = [self capacity:n];
  while (n--)
    [v insertScalar:[[aMatrix rowAt:n] dotMultiply:self]];
  return v;
}

- multiplySelfLeftMatrix:aMatrix
{
  return idReplace(self, [self multiplyLeftMatrix:aMatrix]);
}

- divideLeftMatrix:aMatrix
{
  return [aMatrix solveVector:self];
}
- divideSelfLeftMatrix:aMatrix
{
  return idReplace(self, [self divideLeftMatrix:aMatrix]);
}


- printOn:(IOD)aFile
{
  if ([self numScalars]) {
    id scalar, sequence = [self eachScalar];

    fprintf(aFile, "{");
    [[sequence firstElement] printOn:aFile];
    while (scalar = [sequence next]) {
      fprintf(aFile, ",");
      [scalar printOn:aFile];
    }
    [sequence free];
    fprintf(aFile, "}");
  } else {
    fprintf(aFile, "{}");
  }

  return self;
}

@end

