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

#include "cakit.h"

@implementation Complex

+ real:aReal
{
  return [self real:aReal imaginary:[aReal zero]];
}

- _setUpReal:aReal imaginary:anotherReal
{
  real = aReal;
  imaginary = anotherReal;
  return self;
}

+ real:aReal imaginary:anotherReal
{
  return [[super new] _setUpReal:aReal imaginary:anotherReal];
}

- copy
{
  self = [super copy];
  real = [real cheapCopy];
  imaginary = [imaginary cheapCopy];
  return self;
}

- deepCopy
{
  self = [super deepCopy];
  real = [real deepCopy];
  imaginary = [imaginary deepCopy];
  return self;
}

- clear
{
  [real free];
  [imaginary free];
  return [super clear];
}

- real:aReal imaginary:anotherReal
{
  assert(aReal && anotherReal);
  self = [self clone];
  real = aReal;
  imaginary = anotherReal;
  return self;
}

- real:aReal
{
  return [self real:aReal imaginary:[imaginary zero]];
}

- imaginary:aReal
{
  return [self real:[real zero] imaginary:aReal];
}

- symbol:aSymbol
{
  if (strcmp([aSymbol str], "I") == 0) {
    return [self I];
  } else {
    return [self real:[real symbol:aSymbol]];
  }
}


- real
{
  return real;
}

- imaginary
{
  return imaginary;
}

- (BOOL) isReal
{
  return [imaginary isZero];
}

- (BOOL) notReal
{
  return [imaginary notZero];
}

- (unsigned) hash
{
  return ([real hash] << 16) + ([imaginary hash] & 0x0000ffff);
}

- (BOOL) isEqual:b
{
  if (self == b) {
    return YES;
  } else {
    return [real isEqual:[b real]] && [imaginary isEqual:[b imaginary]];
  }
}


- conjugate
{
  if ([self isReal]) {
    return [self cheapCopy];
  } else {
    return [self imaginary:[imaginary negate]];
  }
}

- conjugateSelf
{
  return idReplace(self, [self conjugate]);
}

- norm
{
  id norm;

  if ([self isReal]) {
    norm = [real square];
  } else {
    id t = [imaginary square];

    norm = [[real square] addSelf:t];
    [t free];
  }

  return norm;
}

- asReal
{
  assert([self isReal]);
  return [real cheapCopy];
}


- zero
{
  return [self real:[real zero] imaginary:[imaginary zero]];
}

- (BOOL) isZero
{
  return [real isZero] && [imaginary isZero];
}

- (BOOL) isOpposite:b
{
  return [real isOpposite:[b real]] && [imaginary isOpposite:[b imaginary]];
}

- negate
{
  if ([self isReal]) {
    return [self real:[real negate]];
  } else {
    return [self real:[real negate] imaginary:[imaginary negate]];
  }
}

- double
{
  if ([self isReal]) {
    return [self real:[real double]];
  } else {
    return [self real:[real double] imaginary:[imaginary double]];
  }
}

- add:b
{
  if (self == b) {
    return [self double];
  } else {
    id r = [real add:[b real]];
    id i = [imaginary add:[b imaginary]];

    return [self real:r imaginary:i];
  }
}

- subtract:b
{
  if (self == b) {
    return [self zero];
  } else {
    id r = [real subtract:[b real]];
    id i = [imaginary subtract:[b imaginary]];

    return [self real:r imaginary:i];
  }
}


- I
{
  return [self imaginary:[imaginary one]];
}

- one
{
  return [self real:[real one]];
}

- (BOOL) isOne
{
  if ([real isZero])
    return [imaginary isUnit];
  if ([imaginary isZero])
    return [real isUnit];
  return NO;
}

- (BOOL) isMinusOne
{
  if ([real isZero])
    return [imaginary isUnit];
  if ([imaginary isZero])
    return [real isUnit];
  return NO;
}

- (BOOL) isUnit
{
  if ([real isZero])
    return [imaginary isUnit];
  if ([imaginary isZero])
    return [real isUnit];
  return NO;
}

- square
{
  if ([self isZero])
    return [self zero];
  if ([self isOne] || [self isMinusOne])
    return [self one];

  if ([self isReal]) {
    return [self real:[real square]];
  } else {
    id r, i;

    r = [real square];
    i = [imaginary square];
    r = [r subtractSelf:i];
    [i free];
    i = [[real multiply:imaginary] doubleSelf];
    return [self real:r imaginary:i];
  }
}

- multiplyReal:b
{
  if ([b isZero])
    return [self zero];
  if ([b isOne])
    return [self cheapCopy];
  if ([b isMinusOne])
    return [self negate];

  return [self real:[real multiply:b] imaginary:[imaginary multiply:b]];
}

- multiplySelfReal:b
{
  return idReplace(self, [self multiplyReal:b]);
}

- multiply:b
{
  if (self == b) {
    return [self square];
  } else {
    id ac, bd, bc, ad;

    if ([self isReal])
      return [b multiplyReal:real];
    if ([b isReal])
      return [self multiplyReal:[b real]];

    ac = [real multiply:[b real]];
    bd = [imaginary multiply:[b imaginary]];
    ad = [real multiply:[b imaginary]];
    bc = [imaginary multiply:[b real]];

    ac = [ac subtractSelf:bd];
    [bd free];
    ad = [ad addSelf:bc];
    [bc free];
    return [self real:ac imaginary:ad];
  }
}

- power:(int)n
{
  if ([self isReal]) {
    return [self real:[real power:n]];
  } else {
    id o, e;

    if (n == 0) {
      return [self one];
    }
    if (n == 1) {
      return [self cheapCopy];
    }
    if (n == 2) {
      return [self square];
    }
    o = [self one];
    e = [self cheapCopy];
    while (n > 1) {
      if (n & 1)
	o = [o multiplySelf:e];
      e = [e squareSelf];
      n >>= 1;
    }
    e = [e multiplySelf:o];
    o = [o free];

    return e;
  }
}


- divideReal:b
{
  id r, i;

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

  r = [real divide:b];
  if (!r)
    return nil;
  i = [imaginary divide:b];
  if (!i) {
    [r free];
    return nil;
  }
  return [self real:r imaginary:i];
}

- divideSelfReal:b
{
  return idReplace(self, [self divideReal:b]);
}

- divide:b
{
  if ([b isZero])
    return nil;
  if (self == b) {
    return [self one];
  } else {
    if ([b isReal]) {
      return [self divideReal:[b real]];
    } else {
      id bbar, norm, res;

      bbar = [b conjugate];
      norm = [b norm];
      /* this is not the same as multiplying by the inverse of b ! */
      res = [[self multiply:bbar] divideSelfReal:norm];
      [norm free];
      [bbar free];
      return res;
    }
  }
}

- inverse
{
  id inv, norm;

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

  norm = [self norm];
  inv = [[self conjugate] divideSelfReal:norm];
  [norm free];
  return inv;
}

- quotientReal:b
{
  id r, i;

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

  r = [real quotient:b];
  if (!r)
    return nil;
  i = [imaginary quotient:b];
  if (!i) {
    [r free];
    return nil;
  }
  return [self real:r imaginary:i];
}

- quotientSelfReal:b
{
  return idReplace(self, [self quotientReal:b]);
}

- quotient:b
{
  if ([b isZero])
    return nil;
  if (self == b) {
    return [self one];
  } else {
    if ([b isReal]) {
      return [self quotientReal:[b real]];
    } else {
      id bbar, norm, res;

      bbar = [b conjugate];
      norm = [b norm];
      res = [[self multiply:bbar] quotientSelfReal:norm];
      [norm free];
      [bbar free];
      return res;
    }
  }
}

- remainder:b quotient:(id *)q
{
  id rem, tmp, quo = [self quotient:b];

  tmp = [quo multiply:b];
  rem = [self subtract:tmp];
  tmp = [tmp free];
  if (q)
    *q = quo;
  else
    [quo free];
  return rem;
}


- gcd:b
{
  if ([self isZero])
    return [b cheapCopy];
  if ([b isZero])
    return [self cheapCopy];

  if ([self isReal] && [b isReal]) {
    id ra, rb, rc;

    ra = [self asReal];
    rb = [b asReal];
    rc = [ra gcd:rb];
    [ra free];
    [rb free];
    return [self real:rc];
  } else {
    id a = [self cheapCopy];

    b = [b cheapCopy];
    while ([b notZero]) {
      id c = [a remainderSelf:b];

      a = b;
      b = c;
    }
    [b free];
    return a;
  }
}


- printOn:(IOD)aFile
{
  if ([imaginary isZero]) {
    if ([real isZero]) {
      fprintf(aFile, "0");
    } else {
      [real printOn:aFile];
    }
  } else {
    if ([real notZero]) {
      [real printOn:aFile];
      fprintf(aFile, " + ");
    }
    fprintf(aFile, "I ");
    [imaginary printOn:aFile];
  }

  return self;
}

@end

