class Numeric
def im
Complex(0, self)
end
def real
self
end
def image
0
end
alias imag image
def arg
Math.atan2!(0, self)
end
alias angle arg
def polar
return abs, arg
end
def conjugate
self
end
alias conj conjugate
end
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
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) other.kind_of?(Integer) or
other.kind_of?(Float) or
(defined?(Rational) and other.kind_of?(Rational))
end
def Complex.polar(r, theta)
Complex(r*Math.cos(theta), r*Math.sin(theta))
end
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
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
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
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
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
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
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 abs
Math.hypot(@real, @image)
end
def abs2
@real*@real + @image*@image
end
def arg
Math.atan2!(@image, @real)
end
alias angle arg
def polar
return abs, arg
end
def conjugate
Complex(@real, -@image)
end
alias conj conjugate
def <=> (other)
self.abs <=> other.abs
end
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
def coerce(other)
if Complex.generic?(other)
return Complex.new!(other), self
else
super
end
end
def denominator
@real.denominator.lcm(@image.denominator)
end
def numerator
cd = denominator
Complex(@real.numerator*(cd/@real.denominator),
@image.numerator*(cd/@image.denominator))
end
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
def hash
@real.hash ^ @image.hash
end
def inspect
sprintf("Complex(%s, %s)", @real.inspect, @image.inspect)
end
I = Complex(0,1)
attr :real
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
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
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
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
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
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
def log(z)
if Complex.generic?(z) and z >= 0
log!(z)
else
r, theta = z.polar
Complex(log!(r.abs), theta)
end
end
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