requires: {}.
provides: {#Magnitude. #Number. #Integer. #SmallInteger. #Float. #Fraction. #Complex. #Infinity. #Epsilon}.

prototypes ensureDelegatedNamespace: #numerics.

numerics addPrototype: #Magnitude derivedFrom: {Cloneable}.
"A Magnitude is any object that compares linearly (<, =, >) with other like
objects."

numerics addPrototype: #Number derivedFrom: {Magnitude}.
"A Number is a Magnitude without dimension."

n@(Number traits) zero
"This method returns the zero element for the numeric type.
So, n * n zero = n zero."
[0].

n@(Number traits) unit
"This method returns the unit element for the numeric type.
So, n * n unit = n."
[1].

n@(Number traits) negated
"The argument's additive inverse generically."
[n zero - n].

n@(Number traits) isZero
"This method tests for zero generically."
[n = n zero].

n@(Number traits) isPositive
"Whether it's positive generically."
[n > n zero].

n@(Number traits) isNegative
"Whether it's negative generically."
[n < n zero].

n@(Number traits) sign
"Answers 1 for positive, -1 for negative, 0 for zero."
[
  n isZero
    ifTrue: [0]
    ifFalse: [n isPositive ifTrue: [1] ifFalse: [-1]]
].

n@(Number traits) squared
"The square of the number."
[n * n].

n@(Number traits) divideByZero
"The default handler for division by zero."
[
  [(DivideByZero of: n) signal] handleWith: DivideByZero ReturnInput
].

numerics addPrototype: #Integer derivedFrom: {Number}.

SmallInteger traits atSlotNamed: #parent0 put: Integer traits.
SmallInteger traits addImmutableSlot: #minimumValue valued: -16r40000000.
SmallInteger traits addImmutableSlot: #maximumValue valued: 16r3FFFFFFF.
SmallInteger traits addImmutableSlot: #level valued: 0.
"Level of the number type for coercion purposes."

i@(SmallInteger traits) highBit
[| n bit |
  i = SmallInteger minimumValue
    ifTrue: [^ (i bitNot highBit + 1)].
  bit: 0.
  n: i abs.
  n > 16rFFFF
    ifTrue: [bit: bit + 16. n: (n bitShift: -16)].
  n > 16rFF
    ifTrue: [bit: bit + 8. n: (n bitShift: -8)].
  n > 16rF
    ifTrue: [bit: bit + 4. n: (n bitShift: -4)].
  n > 2r11
    ifTrue: [bit: bit + 2. n: (n bitShift: -2)].
  n > 2r1
    ifTrue: [bit: bit + 1].
  bit
].

x@(SmallInteger traits) bitShiftOverflow: y@(SmallInteger traits)
[
  (x as: BigInteger) bitShift: y
].

x@(SmallInteger traits) addOverflow: y@(SmallInteger traits)
[
  (x as: BigInteger) + (y as: BigInteger)
].

x@(SmallInteger traits) subtractOverflow: y@(SmallInteger traits)
[
  (x as: BigInteger) - (y as: BigInteger)
].

x@(SmallInteger traits) multiplyOverflow: y@(SmallInteger traits)
[
  (x as: BigInteger) * (y as: BigInteger)
].

numerics addPrototype: #Float derivedFrom: {Number}.

n@(Float traits) zero
"This method returns the zero element for the numeric type.
So, n * n zero = n zero."
[0.0].

n@(Float traits) unit
"This method returns the unit element for the numeric type.
So, n * n unit = n."
[1.0].

x@(Float traits) isNaN
"Simple, endianness-independent test for Not-a-Number."
[
  x ~= x
].

n@(Number traits) isNaN
"Compatibility with Float."
[False].

SingleFloat traits atSlotNamed: #parent0 put: Float traits.
SingleFloat traits addImmutableSlot: #level valued: 3.
"Level of the number type for coercion purposes."

SingleFloat traits addImmutableSlot: #tolerance valued: 0.0001.

_@(Number traits) tolerance [0].

f@(SingleFloat traits) hash
[
  (((f byteAt: 0) 
    bitOr: ((f byteAt: 1) bitShift: 8))
    bitOr: ((f byteAt: 2) bitShift: 16))
    bitXor: (f byteAt: 3)
].

x@(SingleFloat traits) = y@(SingleFloat traits)
[
  (x byteAt: 0) = (y byteAt: 0)
    and: [(x byteAt: 1) = (y byteAt: 1)]
    and: [(x byteAt: 2) = (y byteAt: 2)]
    and: [(x byteAt: 3) = (y byteAt: 3)]
].

x@(Number traits) isCloseTo: y@(Number traits)
"Whether the two Floats are close to each other in value."
[
  x isZero ifTrue: [^ x abs < x tolerance].
  y isZero ifTrue: [^ y abs < y tolerance].
  x isNaN == y isNaN
    or: [(x - y) abs <= ((x abs max: y abs) * x tolerance)]
].

"FIX ME: Make these primitives or use #[] when available to avoid endianness oddities!"
f@(SingleFloat traits) significand
[
  ((f byteAt: 0) 
    bitOr: ((f byteAt: 1) bitShift: 8))
    bitOr: (((f byteAt: 2) bitAnd: 16r7F) bitShift: 16)
].

f@(SingleFloat traits) exponent
[
  ((f byteAt: 2) bitShift: -7)
    bitOr: (((f byteAt: 3) bitAnd: 16r7F) bitShift: 1) 
].
      
f@(SingleFloat traits) significandSize
[
  23
].

f@(SingleFloat traits) exponentSize
[
  8
].

f@(SingleFloat traits) bias
[
  127
].

f@(SingleFloat traits) normalizedWithSignificand: significand exponent: exponent
[| newF |
  newF: f clone.
  newF byteAt: 0 put: significand.
  newF byteAt: 1 put: (significand bitShift: -8).
  newF byteAt: 2 put: ((significand bitShift: -16) bitOr: (exponent bitShift: 7)).
  newF byteAt: 3 put: (exponent bitShift: -1).
  newF
].

f@(SingleFloat traits) denormalizedWithSignificand: significand
[
  newF: f clone.
  newF byteAt: 0 put: significand.
  newF byteAt: 1 put: (significand bitShift: -8).
  newF byteAt: 2 put: (significand bitShift: -16).
  newF
].

f@(Float traits) truncated
[| significand exponent n |
  significand: f significand.
  exponent: f exponent.
  exponent = (f bias * 2 + 1)
    ifTrue:
      [error: 'Can\'t coerce floating-point Inf or NaN to an integer.'].
  exponent isZero
    ifTrue:
      [
        significand isZero
          ifTrue: [^ 0].
        f isNegative
          ifTrue: [^ -1]
          ifFalse: [^ 0]
      ].
  exponent: exponent - f bias.
  significand: ((1 bitShift: f significandSize) bitOr: significand).
  n: (significand bitShift: exponent - f significandSize).
  f isNegative
    ifTrue: [n negated]
    ifFalse: [n]
].

f@(Float traits) as: _@(Integer traits)
[
  f truncated
].

i@(Integer traits) as: _@(Float traits)
[| f n length significand bits newF |
  i isZero
    ifTrue: [^ SingleFloat zero].
"  i highBit <= SingleFloat significandSize
    ifTrue: [f: SingleFloat]
    ifFalse: [f: DoubleFloat].
"
  f: SingleFloat.
  n: i abs.
  length: n highBit.
  length > f bias
    ifTrue:
      [error: 'Integer is too large to coerce to floating-point value.'].
  significand:
   ((n bitShift: f significandSize - length)
     bitAnd: (1 bitShift: f significandSize) - 1).
  newF: (f normalizedWithSignificand: significand exponent: length + f bias).
  i isNegative
   ifTrue: [newF negated]
   ifFalse: [newF]
].

"Constants"
numerics addSlot: #Pi valued: 3.141592654.
numerics addSlot: #E valued: 2.718281828.

n@(Number traits) sqrt
"The square root of the number."
[
  n isNegative ifTrue: [(n as: Complex) sqrt]
    ifFalse: [n raisedTo: 1 / 2]
].

n@(Integer traits) succ
"Answer the Integer's successor."
[n + 1].

n@(Integer traits) pred
"Answer the Integer's predecessor."
[n - 1].

n@(Integer traits) factorial
"The standard recursive definition."
[
  n <= 1
    ifTrue: [1]
    ifFalse: [n * ((n - 1) factorial)]
].

x@(Integer traits) gcd: y@(Integer traits)
"Euclid's algorithm for finding the greatest common denominator."
[| n m |
  n: x.
  m: y.
  [n isZero] whileFalse: [n: m \\ (m: n)].
  m abs
].

x@(Number traits) quo: y@(Number traits)
"Division rounded towards zero."
[
  (x / y) truncated
].

x@(Number traits) rem: y@(Number traits)
"Find the modulus (remainder) of a pair of numbers. Floats will return a non-
Integer difference."
[
  x - ((x quo: y) * y)
].

x@(Number traits) mod: y@(Number traits)
"An alias for rem:."
[
  x rem: y
].

x@(Number traits) isEven
"Answer whether the Number is even; modulo tells us the remainder of division."
[
  (x rem: 2) isZero
].

x@(Number traits) isOdd
[
  x isEven not
].

x@(Number traits) truncated
"The greatest lesser Integer value."
[
  x quo: 1
].

x@(Integer traits) truncated
"Integers have no fractional part."
[x].

i@(Integer traits) as: _@(Integer traits)
"By default integers merely coerce to themselves."
[i].

x@(Number traits) ln
[
  (x as: Float) ln
].

x@(Number traits) exp
[
  (x as: Float) exp
].

x@(Number traits) raisedTo: y@(Number traits)
[
  (x as: Float) raisedTo: (y as: Float)
].

x@(Number traits) ceiling
"The smallest greater Integer value."
[| trunc |
  trunc: x truncated.
  (x = trunc or: [x isNegative]) 
    ifTrue: [trunc] 
    ifFalse: [trunc + 1]
].

x@(Number traits) floor
"The largest lesser Integer value."
[| trunc |
  trunc: x truncated.
  (x = trunc or: [x isPositive])
    ifTrue: [trunc] 
    ifFalse: [trunc - 1]
].

x@(Integer traits) ceiling
"Integers have no fractional part."
[x].

x@(Integer traits) floor
"Integers have no fractional part."
[x].

i@(Integer traits) highBit
[(i ln / 2 ln) ceiling].

i@(Integer traits) lowBit
[
  i = SmallInteger minimumValue
    ifTrue: [i highBit]
    ifFalse: [(i bitXor: i - 1) highBit]
].

x@(Integer traits) << y@(Integer traits)
[
  x bitShift: y
].

x@(Integer traits) >> y@(Integer traits)
[
  x bitShift: y negated
].

x@(Integer traits) rounded
[x].

x@(Float traits) rounded
"Answer the integer nearest to x."
[(x + 0.5) floor].

x@(Number traits) roundTo: q@(Number traits)
"Answer the nearest number which is a multiple of q."
[(x / q) rounded * q].

x@(Number traits) fractionPart
[
  x - x truncated
].

x@(Number traits) \\ y@(Number traits)
[
  x - (x // y * y)
].

x@(Number traits) reciprocal
"Return the multiplicative inverse."
[
  x isZero
    ifTrue: [x divideByZero]
    ifFalse: [1 / x]
].

x@(Number traits) max: y@(Number traits)
"Return the greater argument."
[
  x > y
    ifTrue: [x]
    ifFalse: [y]
].

x@(Number traits) min: y@(Number traits)
"Return the lesser argument."
[
  x < y
    ifTrue: [x]
    ifFalse: [y]
].

x@(Number traits) coerceTo: y@(Number traits)
"Convert x to the same type as y, only in order of ascending level."
[
  x level < y level ifTrue: [x as: y] ifFalse: [x]
].

"Generic arithmetic methods. These work by catching any arithmetic
operation with dissimilar types and converting the less generic type to
the more generic type. They then dispatch the converted arguments to the
method which handles that operation for matched types.
E.g. Float + Complex and Complex + Float will both dispatch to
Complex + Complex"
x@(Number traits) + y@(Number traits) [(x coerceTo: y) + (y coerceTo: x)].
x@(Number traits) - y@(Number traits) [(x coerceTo: y) - (y coerceTo: x)].
x@(Number traits) * y@(Number traits) [(x coerceTo: y) * (y coerceTo: x)].
x@(Number traits) / y@(Number traits) [(x coerceTo: y) / (y coerceTo: x)].
x@(Number traits) < y@(Number traits) [(x coerceTo: y) < (y coerceTo: x)].
x@(Number traits) = y@(Number traits) [(x coerceTo: y) = (y coerceTo: x)].

x@(Integer traits) bitAnd: y@(Integer traits) [(x coerceTo: y) bitAnd: (y coerceTo: x)].
x@(Integer traits) bitOr: y@(Integer traits) [(x coerceTo: y) bitOr: (y coerceTo: x)].
x@(Integer traits) bitXor: y@(Integer traits) [(x coerceTo: y) bitXor: (y coerceTo: x)].

x@(Integer traits) bitNot
[
  x negated - 1
].

x@(Integer traits) bitShift: y@(Integer traits)
[
  (x * (2 raisedTo: y)) truncated
].

x@(Integer traits) / y@(Integer traits)
"Create and reduce a new Fraction for the pair."
[
  y isZero
    ifTrue: [^ x divideByZero].
  (x rem: y) isZero
    ifTrue: [x quo: y]
    ifFalse: [(Fraction newFor: x over: y) reduced]
].

x@(Integer traits) // y@(Integer traits)
"Integer division."
[| q |
  y isZero ifTrue: [^ x divideByZero].
  q: (x quo: y).
  (x isNegative ~= y isNegative and: [q * y ~= x])
    ifTrue: [q - 1]
    ifFalse: [q]
].

x@(Integer traits) lcm: y@(Integer traits)
"Find the least common multiple of the pair."
[
  x // (x gcd: y) * y
].

n@(Integer traits) taken: r@(Integer traits)
"Calculates the number of combinations of N elements taken R at a time.
It's zero outside of Pascal's triangle."
[| num denom |
  (r isNegative or: [r > n]) ifTrue: [^ 0].
  num: 1.
  n above: (r max: n - r)
    do: [| :factor | num: num * factor].
  denom: 1.
  1 to: (r min: n - r)
    do: [| :factor | denom: denom * factor].
  num // denom
].

n@(Integer traits) hashMultiply
"A copying of Squeak's hash-multiplier. TODO: To Be Replaced."
[| low28 low14 |
  low28: (n bitAnd: 268435455).
  low14: (low28 bitAnd: 16383).
  9741 * low14 + ((9741 * (low28 bitShift: -14) + (101 * low14)
		     bitAnd: 16383) * 16384) bitAnd: 268435455
].

x@(Number traits) raisedTo: y@(Integer traits)
[
  y isZero ifTrue: [^ x unit].
  (x isZero or: [y = 1]) ifTrue: [^ x].
  y isPositive
    ifTrue: [(x * x raisedTo: y // 2) * (x raisedTo: y \\ 2)]
    ifFalse: [(x raisedTo: y abs) reciprocal]
].

numerics addPrototype: #Fraction derivedFrom: {Number}.
Fraction addSlot: #numerator.
Fraction addSlot: #denominator.
Fraction traits addImmutableSlot: #level valued: 2.

"Fractions obviously represent reduced ratios of Integers. They are not unique
objects per value, however, and they are lazily reduced."

x@(Fraction traits) numerator: numer@(Integer traits) denominator: denom@(Integer traits)
"Initialize the Fraction."
[
  denom isNegative
    ifTrue: [
      x numerator: numer negated.  
      x denominator: denom negated.
      ^ x].
  x numerator: numer.
  x denominator: denom.
  x
].

x@(Fraction traits) newFor: numer over: denom
[
  x clone numerator: numer denominator: denom
].

x@(Fraction traits) printOn: s
"Print out a representation of the Fraction."
[
  x numerator printOn: s.
  s ; ' / '.
  x denominator printOn: s.
  x
].

x@(Fraction traits) truncated
"The nearest lesser Integer to the Fraction."
[
  x numerator quo: x denominator
].

x@(Fraction traits) * y@(Fraction traits)
[| xn xd yn yd d1 d2 |
  d1: ((xn: x numerator) gcd: (yd: y denominator)).
  d2: ((xd: x denominator) gcd: (yn: y numerator)).
  ((d2 = xd and: [d1 = yd])
     ifTrue: [xn // d1 * yn // d2]
     ifFalse: [x newFor: xn // d1 * (yn // d2)
		    over: xd // d2 * (yd // d1)])
].

x@(Fraction traits) + y@(Fraction traits)
[| n d d1 d2 |
  d: (x denominator gcd: y denominator).
  n: (x numerator * (d1: y denominator // d) + (y numerator * (d2: x denominator // d))).
  d1: d1 * d2.
  n: n // (d2: (n gcd: d)).
  (d: d1 * (d // d2)) = 1 ifTrue: [^ n].
  x newFor: n over: d
].

x@(Fraction traits) - y@(Fraction traits)
[x + y negated].

x@(Fraction traits) reduced
[| gcd numer denom |
  x numerator isZero ifTrue: [^ 0].
  gcd: (x numerator abs gcd: x denominator).
  numer: x numerator // gcd.
  denom: x denominator // gcd.
  denom = 1 ifTrue: [^ numer].
  x newFor: numer over: denom
].

x@(Fraction traits) negated
[
  x newFor: x numerator negated over: x denominator
].

x@(Fraction traits) reciprocal
[| numer |
  numer: x numerator.
  numer isZero ifTrue: [^ x denominator divideByZero].
  numer = 1 ifTrue: [^ x denominator].
  numer = -1 ifTrue: [^ x denominator negated].
  x newFor: x denominator over: numer
].

x@(Fraction traits) < y@(Fraction traits)
[
  (x numerator * y denominator) < (y numerator * x denominator)
].

x@(Fraction traits) = y@(Fraction traits)
[
  x numerator = y numerator and: [x denominator = y denominator]
].

x@(Fraction traits) / y@(Fraction traits)
[x * y reciprocal].

x@(Fraction traits) as: y@(Float traits)
[(x numerator as: y) / (x denominator as: y)].

x@(Fraction traits) as: _@(Integer traits)
[x truncated].

x@(Integer traits) as: y@(Fraction traits)
[y newFor: x over: 1].

f@(Fraction traits) raisedTo: n@(Number traits)
"Raise a Fraction to a power n"
[
  n isZero ifTrue: [^ f unit].
  n isPositive
    ifTrue: [(f numerator raisedTo: n) / (f denominator raisedTo: n)]
    ifFalse: [(f raisedTo: n negated) reciprocal]
].

x@(Fraction traits) rounded
"Answer the integer nearest to x."
[(x + (1 / 2)) floor].

x@(Magnitude traits) < y@(Magnitude traits)  [overrideThis].
x@(Magnitude traits) <= y@(Magnitude traits) [(y < x) not].
x@(Magnitude traits) >= y@(Magnitude traits) [(x < y) not].
x@(Magnitude traits) > y@(Magnitude traits)  [y < x].

x@(Magnitude traits) between: min@(Magnitude traits) and: max@(Magnitude traits)
"Return whether the first argument is between the other two."
[
  x >= min and: [x <= max]
].

x@(Magnitude traits) min: y@(Magnitude traits)
"Return the lesser argument."
[
  x < y
    ifTrue: [x]
    ifFalse: [y]
].

x@(Magnitude traits) max: y@(Magnitude traits)
"Return the greater argument."
[
  x > y
    ifTrue: [x]
    ifFalse: [y]
].

x@(Magnitude traits) min: y@(Magnitude traits) max: z@(Magnitude traits)
[
  (x min: y) max: z
].

n@(Magnitude traits) abs
"The absolute value."
[
  n isNegative
    ifTrue: [n negated]
    ifFalse: [n]
].

numerics addSlot: #LargeUnbounded valued: Magnitude clone.
"A quantity which is unbounded but still a Number: as-large-as-you-like. It's
definitely not an infinity."

x@(Magnitude traits) < _@LargeUnbounded
[True].

_@LargeUnbounded < _@LargeUnbounded
[False].

numerics addSlot: #PositiveInfinity valued: Magnitude clone.

x@(Magnitude traits) < pi@PositiveInfinity
[
  pi ~= x
].

pi@PositiveInfinity * x
[
  x isZero ifTrue: [^ x].
  x isNegative
    ifTrue: [pi negated]
    ifFalse: [pi]
].

_@PositiveInfinity printOn: s
[
  s ; 'positive infinity'
].

numerics addSlot: #NegativeInfinity valued: Magnitude clone.

ni@NegativeInfinity < x@(Magnitude traits)
[
  ni ~= x
].

ni@NegativeInfinity * x
[
  x isZero ifTrue: [^ x].
  x isNegative
    ifTrue: [ni negated]
    ifFalse: [ni]
].

_@NegativeInfinity printOn: s
[
  s ; 'negative infinity'
].

_@PositiveInfinity negated
[NegativeInfinity].

_@NegativeInfinity negated
[PositiveInfinity].

_ - pi@PositiveInfinity
[
  pi negated
].

_ - ni@NegativeInfinity
[
  ni negated
].

numerics addSlot: #PositiveEpsilon valued: Magnitude clone.
"The reciprocal of PositiveInfinity."

pe@PositiveEpsilon = x
[pe == x].

_@PositiveEpsilon isZero
[False].

_@PositiveEpsilon isPositive
[True].

_@(Magnitude traits) / _@PositiveInfinity
[PositiveEpsilon].

_@(Magnitude traits) / _@PositiveEpsilon
[PositiveInfinity].

_@(Magnitude traits) * pe@PositiveEpsilon
[pe].

x@(Magnitude traits) + _@PositiveEpsilon
[x].

_@PositiveEpsilon as: x@(Magnitude traits)
[0.0 as: x].

numerics addSlot: #NegativeEpsilon valued: Magnitude clone.
"The reciprocal of NegativeInfinity."

ne@NegativeEpsilon = x
[ne == x].

_@NegativeEpsilon isZero
[False].

_@NegativeEpsilon isNegative
[True].

_@(Magnitude traits) / _@NegativeInfinity
[NegativeEpsilon].

_@(Magnitude traits) / _@NegativeEpsilon
[NegativeInfinity].

_@(Magnitude traits) * ne@NegativeEpsilon
[ne].

x@(Magnitude traits) + _@NegativeEpsilon
[x].

_@NegativeEpsilon as: x@(Magnitude traits)
[
  -0.0 as: x
].

_@PositiveEpsilon negated
[NegativeEpsilon].

_@NegativeEpsilon negated
[PositiveEpsilon].

numerics addPrototype: #Complex derivedFrom: {Number}.
"Represents a complex number, having a real and imaginary part, which should
not be altered per instance."
Complex addSlot: #real valued: 0.
"The real component."
Complex addSlot: #imaginary valued: 0.
"The imaginary component."
Complex traits addSlot: #level valued: 4.

c1@(Complex traits) = c2@(Complex traits)
[
  c1 real = c2 real and: [c1 imaginary = c2 imaginary]
].

c@(Complex traits) = n@(Magnitude traits)
[
  c imaginary isZero and: [c real = n]
].

c@(Complex traits) hash
[
  c real hash bitXor: c imaginary hash
].

c@(Complex traits) real: real imaginary: imaginary
[| newC |
  newC: c clone.
  newC real: real.
  newC imaginary: imaginary.
  newC
].

r@(Magnitude traits) plusImaginary: i@(Magnitude traits)
"A convenient method for creating Complex quantities."
[Complex real: r imaginary: i].

x@(Magnitude traits) as: c@(Complex traits)
"Creates the corresponding complex number for the value, using a corresponding
zero coordinate."
[
  c real: x imaginary: x zero
].

c@(Complex traits) as: x@(Magnitude traits)
"Converts to a number if possible, but returns itself otherwise."
[
  c imaginary isZero
    ifTrue: [c real]
    ifFalse: [c]
].

c@(Complex traits) unit
[
  c real unit as: c
].

c@(Complex traits) zero
[
  c real zero as: c
].

c@(Complex traits) isZero
[
  c real isZero and: [c imaginary isZero]
].

c@(Complex traits) sign
"Returns a complex number of unit dimensions in the same quadrant."
[
  c real: c real sign imaginary: c imaginary sign
].

c@(Complex traits) radius
[
  (c real squared + c imaginary squared) sqrt
].

c@(Complex traits) abs
"The magnitude of the complex as a vector."
[
  c radius
].

c1@(Complex traits) < c2@(Complex traits)
"This is highly dubious, since Complex arithmetic and relationships are
totally different than linearly-ordered Magnitudes. This is here solely to
satisfy Magnitude protocol, embedding Complex into Real numbers."
[
  deprecated.
  c1 real < c2 real
].

numerics addSlot: #I valued: (Complex real: 0 imaginary: -1).
"The square root of -1."

c@(Complex traits) * _@I
"The complex number multiplied by -1 sqrt."
[
  c real: c imaginary negated imaginary: c real
].

c@(Complex traits) theta
"The angle between the complex coordinates and the positive real-axis."
[
  c real isZero
     ifTrue:
       [c imaginary isPositive
	  ifTrue: [Pi / 2]
	  ifFalse: [Pi * 3 / 2]]
     ifFalse:
       [| tan result |
	  tan: c imaginary / c real.
	  result: tan arcTan.
	  c real isPositive
	    ifTrue: [theta]
	    ifFalse: [Pi + result]]
].

c@(Complex traits) conjugated
"Returns the complex conjugate."
[
  c real: c real imaginary: c imaginary negated
].

c1@(Complex traits) + c2@(Complex traits)
[
  c1 real: c1 real + c2 real
     imaginary: c1 imaginary + c2 imaginary
].

c1@(Complex traits) - c2@(Complex traits)
[
  c1 real: c1 real - c2 real
     imaginary: c1 imaginary - c2 imaginary
].

c1@(Complex traits) * c2@(Complex traits)
[
  c1 real: c2 real * c1 real - (c1 imaginary * c2 imaginary)
     imaginary: c1 real * c2 imaginary + (c1 imaginary * c2 real)
].

c1@(Complex traits) / c2@(Complex traits)
[| denom |
  denom: c2 real squared + c2 imaginary squared.
  c1 real: c1 real * c2 real + (c1 imaginary * c2 imaginary) / denom
     imaginary: c2 real * c1 imaginary - (c2 imaginary * c1 real) / denom
].

c@(Complex traits) sqrt
[| r |
  (c imaginary isZero and: [c real isNegative not])
     ifTrue: [^ c real sqrt].
  c real: c imaginary / 2 / (r: (c radius - c real / 2) sqrt)
    imaginary: r
].

c@(Complex traits) exp
[| zero im |
  zero: ((c real zero + c imaginary zero) as: Real).
  im: (zero as: c imaginary).
  (zero as: c real) exp * (im cos + (im sin * i))
].

c@(Complex traits) ln
[
  c real: c radius ln imaginary: c theta
].

c@(Complex traits) cos
[
  (c * I) exp + (c * I) negated exp / 2
].

c@(Complex traits) sin
[
  (c * I) exp - (c * I) negated exp / 2
].

c@(Complex traits) tan
[
  c sin / c cos
].

c@(Complex traits) raisedTo: n
[
  n isZero ifTrue: [^ 1].
  n isPositive
    ifTrue: [(n * c ln) exp]
    ifFalse: [(c raisedTo: n negated) reciprocal]
].

c@(Complex traits) printOn: s
[
  s ; '('.
  c real printOn: s.
  s ; ' + '.
  c imaginary printOn: s.
  s ; ' i)'.
  s
].

"Yuck! This should be in dimensioned.slate but is here for now"
"TODO: Move this"
n@(Number traits) degreesToRadians
[n * Pi / 180].

n@(Number traits) radiansToDegrees
[n * 180 / Pi].

"A portable way to represent floating point numbers without loss of precision
http://research.microsoft.com/~hollasch/cgindex/coding/portfloat.html
Would be useful in #storeOn: #loadFrom: for floats"

n@(Float traits) storeOn: s
"Stores a float in portable encoding"
[].

n@(Float traits) loadFrom: s
"Loads a float in portable encoding"
[].