#
#   complex.rb - 
#   	$Release Version: 0.5 $
#   	$Revision: 1.3 $
#   	$Date: 1998/07/08 10:05:28 $
#   	by Keiju ISHITSUKA(SHL Japan Inc.)
#
# ----
#
# complex.rb implements the Complex class for complex numbers.  Additionally,
# some methods in other Numeric classes are redefined or added to allow greater
# interoperability with Complex numbers.
#
# Complex numbers can be created in the following manner:
# - <tt>Complex(a, b)</tt>
# - <tt>Complex.polar(radius, theta)</tt>
#   
# Additionally, note the following:
# - <tt>Complex::I</tt> (the mathematical constant <i>i</i>)
# - <tt>Numeric#im</tt> (e.g. <tt>5.im -> 0+5i</tt>)
#
# The following +Math+ module methods are redefined to handle Complex arguments.
# They will work as normal with non-Complex arguments.
#    sqrt exp cos sin tan log log10
#    cosh sinh tanh acos asin atan atan2 acosh asinh atanh
#


#
# Numeric is a built-in class on which Fixnum, Bignum, etc., are based.  Here
# some methods are added so that all number types can be treated to some extent
# as Complex numbers.
#
class Numeric
  #
  # Returns a Complex number <tt>(0,<i>self</i>)</tt>.
  #
  def im
    Complex(0, self)
  end
  
  #
  # The real part of a complex number, i.e. <i>self</i>.
  #
  def real
    self
  end
  
  #
  # The imaginary part of a complex number, i.e. 0.
  #
  def image
    0
  end
  alias imag image
  
  #
  # See Complex#arg.
  #
  def arg
    if self >= 0
      return 0
    else
      return Math::PI
    end
  end
  alias angle arg
  
  #
  # See Complex#polar.
  #
  def polar
    return abs, arg
  end
  
  #
  # See Complex#conjugate (short answer: returns <i>self</i>).
  #
  def conjugate
    self
  end
  alias conj conjugate
end


#
# Creates a Complex number.  +a+ and +b+ should be Numeric.  The result will be
# <tt>a+bi</tt>.
#
def Complex(a, b = 0)
  if b == 0 and (a.kind_of?(Complex) or defined? Complex::Unify)
    a
  else
    Complex.new( a.real-b.imag, a.imag+b.real )
  end
end

#
# The complex number class.  See complex.rb for an overview.
#
class Complex < Numeric
  @RCS_ID='-$Id: complex.rb,v 1.3 1998/07/08 10:05:28 keiju Exp keiju $-'

  undef step
  undef div, divmod
  undef floor, truncate, ceil, round

  def Complex.generic?(other) # :nodoc:
    other.kind_of?(Integer) or
    other.kind_of?(Float) or
    (defined?(Rational) and other.kind_of?(Rational))
  end

  #
  # Creates a +Complex+ number in terms of +r+ (radius) and +theta+ (angle).
  #
  def Complex.polar(r, theta)
    Complex(r*Math.cos(theta), r*Math.sin(theta))
  end

  #
  # Creates a +Complex+ number <tt>a</tt>+<tt>b</tt><i>i</i>.
  #
  def Complex.new!(a, b=0)
    new(a,b)
  end

  def initialize(a, b)
    raise TypeError, "non numeric 1st arg `#{a.inspect}'" if !a.kind_of? Numeric
    raise TypeError, "`#{a.inspect}' for 1st arg" if a.kind_of? Complex
    raise TypeError, "non numeric 2nd arg `#{b.inspect}'" if !b.kind_of? Numeric
    raise TypeError, "`#{b.inspect}' for 2nd arg" if b.kind_of? Complex
    @real = a
    @image = b
  end

  #
  # Addition with real or complex number.
  #
  def + (other)
    if other.kind_of?(Complex)
      re = @real + other.real
      im = @image + other.image
      Complex(re, im)
    elsif Complex.generic?(other)
      Complex(@real + other, @image)
    else
      x , y = other.coerce(self)
      x + y
    end
  end
  
  #
  # Subtraction with real or complex number.
  #
  def - (other)
    if other.kind_of?(Complex)
      re = @real - other.real
      im = @image - other.image
      Complex(re, im)
    elsif Complex.generic?(other)
      Complex(@real - other, @image)
    else
      x , y = other.coerce(self)
      x - y
    end
  end
  
  #
  # Multiplication with real or complex number.
  #
  def * (other)
    if other.kind_of?(Complex)
      re = @real*other.real - @image*other.image
      im = @real*other.image + @image*other.real
      Complex(re, im)
    elsif Complex.generic?(other)
      Complex(@real * other, @image * other)
    else
      x , y = other.coerce(self)
      x * y
    end
  end
  
  #
  # Division by real or complex number.
  #
  def / (other)
    if other.kind_of?(Complex)
      self*other.conjugate/other.abs2
    elsif Complex.generic?(other)
      Complex(@real/other, @image/other)
    else
      x, y = other.coerce(self)
      x/y
    end
  end
  
  def quo(other)
    Complex(@real.quo(1), @image.quo(1)) / other
  end

  #
  # Raise this complex number to the given (real or complex) power.
  #
  def ** (other)
    if other == 0
      return Complex(1)
    end
    if other.kind_of?(Complex)
      r, theta = polar
      ore = other.real
      oim = other.image
      nr = Math.exp!(ore*Math.log!(r) - oim * theta)
      ntheta = theta*ore + oim*Math.log!(r)
      Complex.polar(nr, ntheta)
    elsif other.kind_of?(Integer)
      if other > 0
	x = self
	z = x
	n = other - 1
	while n != 0
	  while (div, mod = n.divmod(2)
		 mod == 0)
	    x = Complex(x.real*x.real - x.image*x.image, 2*x.real*x.image)
	    n = div
	  end
	  z *= x
	  n -= 1
	end
	z
      else
	if defined? Rational
	  (Rational(1) / self) ** -other
	else
	  self ** Float(other)
	end
      end
    elsif Complex.generic?(other)
      r, theta = polar
      Complex.polar(r**other, theta*other)
    else
      x, y = other.coerce(self)
      x**y
    end
  end
  
  #
  # Remainder after division by a real or complex number.
  #
  def % (other)
    if other.kind_of?(Complex)
      Complex(@real % other.real, @image % other.image)
    elsif Complex.generic?(other)
      Complex(@real % other, @image % other)
    else
      x , y = other.coerce(self)
      x % y
    end
  end
  
#--
#    def divmod(other)
#      if other.kind_of?(Complex)
#        rdiv, rmod = @real.divmod(other.real)
#        idiv, imod = @image.divmod(other.image)
#        return Complex(rdiv, idiv), Complex(rmod, rmod)
#      elsif Complex.generic?(other)
#        Complex(@real.divmod(other), @image.divmod(other))
#      else
#        x , y = other.coerce(self)
#        x.divmod(y)
#      end
#    end
#++
  
  #
  # Absolute value (aka modulus): distance from the zero point on the complex
  # plane.
  #
  def abs
    Math.hypot(@real, @image)
  end
  
  #
  # Square of the absolute value.
  #
  def abs2
    @real*@real + @image*@image
  end
  
  #
  # Argument (angle from (1,0) on the complex plane).
  #
  def arg
    Math.atan2!(@image, @real)
  end
  alias angle arg
  
  #
  # Returns the absolute value _and_ the argument.
  #
  def polar
    return abs, arg
  end
  
  #
  # Complex conjugate (<tt>z + z.conjugate = 2 * z.real</tt>).
  #
  def conjugate
    Complex(@real, -@image)
  end
  alias conj conjugate
  
  #
  # Compares the absolute values of the two numbers.
  #
  def <=> (other)
    self.abs <=> other.abs
  end
  
  #
  # Test for numerical equality (<tt>a == a + 0<i>i</i></tt>).
  #
  def == (other)
    if other.kind_of?(Complex)
      @real == other.real and @image == other.image
    elsif Complex.generic?(other)
      @real == other and @image == 0
    else
      other == self
    end
  end

  #
  # Attempts to coerce +other+ to a Complex number.
  #
  def coerce(other)
    if Complex.generic?(other)
      return Complex.new!(other), self
    else
      super
    end
  end

  #
  # FIXME
  #
  def denominator
    @real.denominator.lcm(@image.denominator)
  end
  
  #
  # FIXME
  #
  def numerator
    cd = denominator
    Complex(@real.numerator*(cd/@real.denominator),
	    @image.numerator*(cd/@image.denominator))
  end
  
  #
  # Standard string representation of the complex number.
  #
  def to_s
    if @real != 0
      if defined?(Rational) and @image.kind_of?(Rational) and @image.denominator != 1
	if @image >= 0
	  @real.to_s+"+("+@image.to_s+")i"
	else
	  @real.to_s+"-("+(-@image).to_s+")i"
	end
      else
	if @image >= 0
	  @real.to_s+"+"+@image.to_s+"i"
	else
	  @real.to_s+"-"+(-@image).to_s+"i"
	end
      end
    else
      if defined?(Rational) and @image.kind_of?(Rational) and @image.denominator != 1
	"("+@image.to_s+")i"
      else
	@image.to_s+"i"
      end
    end
  end
  
  #
  # Returns a hash code for the complex number.
  #
  def hash
    @real.hash ^ @image.hash
  end
  
  #
  # Returns "<tt>Complex(<i>real</i>, <i>image</i>)</tt>".
  #
  def inspect
    sprintf("Complex(%s, %s)", @real.inspect, @image.inspect)
  end

  
  #
  # +I+ is the imaginary number.  It exists at point (0,1) on the complex plane.
  #
  I = Complex(0,1)
  
  # The real part of a complex number.
  attr :real

  # The imaginary part of a complex number.
  attr :image
  alias imag image
  
end

class Integer

  unless defined?(1.numerator)
    def numerator() self end
    def denominator() 1 end

    def gcd(other)
      min = self.abs
      max = other.abs
      while min > 0
        tmp = min
        min = max % min
        max = tmp
      end
      max
    end

    def lcm(other)
      if self.zero? or other.zero?
        0
      else
        (self.div(self.gcd(other)) * other).abs
      end
    end

  end

end

module Math
  alias sqrt! sqrt
  alias exp! exp
  alias log! log
  alias log10! log10
  alias cos! cos
  alias sin! sin
  alias tan! tan
  alias cosh! cosh
  alias sinh! sinh
  alias tanh! tanh
  alias acos! acos
  alias asin! asin
  alias atan! atan
  alias atan2! atan2
  alias acosh! acosh
  alias asinh! asinh
  alias atanh! atanh  

  # Redefined to handle a Complex argument.
  def sqrt(z)
    if Complex.generic?(z)
      if z >= 0
	sqrt!(z)
      else
	Complex(0,sqrt!(-z))
      end
    else
      if z.image < 0
	sqrt(z.conjugate).conjugate
      else
	r = z.abs
	x = z.real
	Complex( sqrt!((r+x)/2), sqrt!((r-x)/2) )
      end
    end
  end
  
  # Redefined to handle a Complex argument.
  def exp(z)
    if Complex.generic?(z)
      exp!(z)
    else
      Complex(exp!(z.real) * cos!(z.image), exp!(z.real) * sin!(z.image))
    end
  end
  
  # Redefined to handle a Complex argument.
  def cos(z)
    if Complex.generic?(z)
      cos!(z)
    else
      Complex(cos!(z.real)*cosh!(z.image),
	      -sin!(z.real)*sinh!(z.image))
    end
  end
    
  # Redefined to handle a Complex argument.
  def sin(z)
    if Complex.generic?(z)
      sin!(z)
    else
      Complex(sin!(z.real)*cosh!(z.image),
	      cos!(z.real)*sinh!(z.image))
    end
  end
  
  # Redefined to handle a Complex argument.
  def tan(z)
    if Complex.generic?(z)
      tan!(z)
    else
      sin(z)/cos(z)
    end
  end

  def sinh(z)
    if Complex.generic?(z)
      sinh!(z)
    else
      Complex( sinh!(z.real)*cos!(z.image), cosh!(z.real)*sin!(z.image) )
    end
  end

  def cosh(z)
    if Complex.generic?(z)
      cosh!(z)
    else
      Complex( cosh!(z.real)*cos!(z.image), sinh!(z.real)*sin!(z.image) )
    end
  end

  def tanh(z)
    if Complex.generic?(z)
      tanh!(z)
    else
      sinh(z)/cosh(z)
    end
  end
  
  # Redefined to handle a Complex argument.
  def log(z)
    if Complex.generic?(z) and z >= 0
      log!(z)
    else
      r, theta = z.polar
      Complex(log!(r.abs), theta)
    end
  end
  
  # Redefined to handle a Complex argument.
  def log10(z)
    if Complex.generic?(z)
      log10!(z)
    else
      log(z)/log!(10)
    end
  end

  def acos(z)
    if Complex.generic?(z) and z >= -1 and z <= 1
      acos!(z)
    else
      -1.0.im * log( z + 1.0.im * sqrt(1.0-z*z) )
    end
  end

  def asin(z)
    if Complex.generic?(z) and z >= -1 and z <= 1
      asin!(z)
    else
      -1.0.im * log( 1.0.im * z + sqrt(1.0-z*z) )
    end
  end

  def atan(z)
    if Complex.generic?(z)
      atan!(z)
    else
      1.0.im * log( (1.0.im+z) / (1.0.im-z) ) / 2.0
    end
  end

  def atan2(y,x)
    if Complex.generic?(y) and Complex.generic?(x)
      atan2!(y,x)
    else
      -1.0.im * log( (x+1.0.im*y) / sqrt(x*x+y*y) )
    end
  end

  def acosh(z)
    if Complex.generic?(z) and z >= 1
      acosh!(z)
    else
      log( z + sqrt(z*z-1.0) )
    end
  end

  def asinh(z)
    if Complex.generic?(z)
      asinh!(z)
    else
      log( z + sqrt(1.0+z*z) )
    end
  end

  def atanh(z)
    if Complex.generic?(z) and z >= -1 and z <= 1
      atanh!(z)
    else
      log( (1.0+z) / (1.0-z) ) / 2.0
    end
  end

  module_function :sqrt!
  module_function :sqrt
  module_function :exp!
  module_function :exp
  module_function :log!
  module_function :log
  module_function :log10!
  module_function :log10
  module_function :cosh!
  module_function :cosh
  module_function :cos!
  module_function :cos
  module_function :sinh!
  module_function :sinh
  module_function :sin!
  module_function :sin
  module_function :tan!
  module_function :tan
  module_function :tanh!
  module_function :tanh
  module_function :acos!
  module_function :acos
  module_function :asin!
  module_function :asin
  module_function :atan!
  module_function :atan
  module_function :atan2!
  module_function :atan2
  module_function :acosh!
  module_function :acosh
  module_function :asinh!
  module_function :asinh
  module_function :atanh!
  module_function :atanh
  
end

# Documentation comments:
#  - source: original (researched from pickaxe)
#  - a couple of fixme's
#  - RDoc output for Bignum etc. is a bit short, with nothing but an
#    (undocumented) alias.  No big deal.