diff options
Diffstat (limited to 'lib/complex.rb')
| -rw-r--r-- | lib/complex.rb | 490 |
1 files changed, 490 insertions, 0 deletions
diff --git a/lib/complex.rb b/lib/complex.rb new file mode 100644 index 0000000000..aa5d219d2f --- /dev/null +++ b/lib/complex.rb @@ -0,0 +1,490 @@ +# +# complex.rb - +# $Release Version: 0.5 $ +# $Revision: 1.1 $ +# $Date: 1996/11/11 04:25:19 $ +# by Keiju ISHITSUKA(SHL Japan Inc.) +# +# -- +# Usage: +# class Complex < Numeric +# +# Complex(x, y) --> x + yi +# y.im --> 0 + yi +# +# Complex::polar +# +# Complex::+ +# Complex::- +# Complex::* +# Complex::/ +# Complex::** +# Complex::% +# Complex::divmod +# Complex::abs +# Complex::abs2 +# Complex::arg +# Complex::polar +# Complex::conjugate +# Complex::<=> +# Complex::== +# Complex::to_i +# Complex::to_f +# Complex::to_r +# Complex::to_s +# +# Complex::I +# +# Numeric::im +# +# Math.sqrt +# Math.exp +# Math.cos +# Math.sin +# Math.tan +# Math.log +# Math.log10 +# Math.atan2 +# +# + +def Complex(a, b = 0) + if a.kind_of?(Complex) and b == 0 + a + elsif b == 0 and defined? Complex::Unify + a + else + Complex.new(a, b) + end +end + +class Complex < Numeric + + 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 initialize(a, b = 0) + @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 = a.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 = a.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 = a.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 = a.coerce(self) + x / y + end + 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.power!(other), theta * other) + else + x , y = a.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 = a.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, rdiv) + elsif Complex.generic?(other) + Complex(@real.divmod(other), @image.divmod(other)) + else + x , y = a.coerce(self) + x.divmod(y) + end + end + + def abs + Math.sqrt!((@real*@real + @image*@image).to_f) + end + + def abs2 + @real*@real + @image*@image + end + + def arg + Math.atan2(@image.to_f, @real.to_f) + end + + def polar + return abs, arg + end + + def conjugate + Complex(@real, -@image) + end + + 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 + x , y = a.coerce(self) + x == y + end + end + + def coerce(other) + if Complex.generic?(other) + return Complex.new(other), self + else + super + end + end + + def to_i + Complex(@real.to_i, @image.to_i) + end + + def to_f + Complex(@real.to_f, @image.to_f) + end + + def to_r + Complex(@real.to_r, @image.to_r) + 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 ^ @image + end + + I = Complex(0,1) + + attr :real + attr :image + +end + +class Numeric + def im + Complex(0, self) + end + + def real + self + end + + def image + 0 + end + + def arg + if self >= 0 + return 0 + else + return Math.atan2(1,1)*4 + end + end + + def polar + return abs, arg + end + + def conjugate + self + end +end + +class Fixnum + if not defined? Rational + alias power! ** + end + + def ** (other) + if self < 0 + Complex.new(self) ** other + else + if defined? Rational + if other >= 0 + self.power!(other) + else + Rational.new!(self,1)**other + end + else + self.power!(other) + end + end + end +end + +class Bignum + if not defined? Rational + alias power! ** + end +end + +class Float + alias power! ** +end + +module Math + alias sqrt! sqrt + alias exp! exp + alias cos! cos + alias sin! sin + alias tan! tan + alias log! log + alias log10! log10 + alias atan2! atan2 + + def sqrt(z) + if Complex.generic?(z) + if z >= 0 + sqrt!(z) + else + Complex(0,sqrt!(-z)) + end + else + z**Rational(1,2) + 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 cosh!(x) + (exp!(x) + exp!(-x))/2.0 + end + + def sinh!(x) + (exp!(x) - exp!(-x))/2.0 + 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 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 atan2(x, y) + if Complex.generic?(x) and Complex.generic?(y) + atan2!(x, y) + else + fail "Not yet implemented." + end + end + + def atanh!(x) + log((1.0 + x.to_f) / ( 1.0 - x.to_f)) / 2.0 + end + + def atan(z) + if Complex.generic?(z) + atan2!(z, 1) + elsif z.image == 0 + atan2(z.real,1) + else + a = z.real + b = z.image + + c = (a*a + b*b - 1.0) + d = (a*a + b*b + 1.0) + + Complex(atan2!((c + sqrt(c*c + 4.0*a*a)), 2.0*a), + atanh!((-d + sqrt(d*d - 4.0*b*b))/(2.0*b))) + end + end + + module_function :sqrt + module_function :sqrt! + module_function :exp! + module_function :exp + module_function :cosh! + module_function :cos! + module_function :cos + module_function :sinh! + module_function :sin! + module_function :sin + module_function :tan! + module_function :tan + module_function :log! + module_function :log + module_function :log10! + module_function :log + module_function :atan2! + module_function :atan2 +# module_function :atan! + module_function :atan + module_function :atanh! + +end + + |