ALPHALINUX5/util/ALPHALINUX5/lib/ruby/1.6/complex.rb

   1 #
   2 #   complex.rb -
   3 #       $Release Version: 0.5 $
   4 #       $Revision: 1.1.1.1 $
   5 #       $Date: 2002/02/01 06:35:39 $
   6 #       by Keiju ISHITSUKA(SHL Japan Inc.)
   7 #
   8 # --
   9 #   Usage:
  10 #      class Complex < Numeric
  11 #
  12 #   Complex(x, y) --> x + yi
  13 #   y.im          --> 0 + yi
  14 #
  15 #   Complex::polar
  16 #
  17 #   Complex::+
  18 #   Complex::-
  19 #   Complex::*
  20 #   Complex::/
  21 #   Complex::**
  22 #   Complex::%
  23 #   Complex::divmod -- obsolete
  24 #   Complex::abs
  25 #   Complex::abs2
  26 #   Complex::arg
  27 #   Complex::polar
  28 #   Complex::conjugate
  29 #   Complex::<=>
  30 #   Complex::==
  31 #   Complex::to_i
  32 #   Complex::to_f
  33 #   Complex::to_r
  34 #   Complex::to_s
  35 #
  36 #   Complex::I
  37 #
  38 #   Numeric::im
  39 #
  40 #   Math.sqrt
  41 #   Math.exp
  42 #   Math.cos
  43 #   Math.sin
  44 #   Math.tan
  45 #   Math.log
  46 #   Math.log10
  47 #   Math.atan2
  48 #
  49 #
  50
  51 def Complex(a, b = 0)
  52   if a.kind_of?(Complex) and b == 0
  53     a
  54   elsif b.kind_of?(Complex)
  55     if a.kind_of?(Complex)
  56       Complex(a.real-b.image, a.image + b.real)
  57     else
  58       Complex(a-b.image, b.real)
  59     end
  60   elsif b == 0 and defined? Complex::Unify
  61     a
  62   else
  63     Complex.new(a, b)
  64   end
  65 end
  66
  67 class Complex < Numeric
  68   @RCS_ID='-$Id: complex.rb,v 1.1.1.1 2002/02/01 06:35:39 tacyas Exp $-'
  69
  70   def Complex.generic?(other)
  71     other.kind_of?(Integer) or
  72     other.kind_of?(Float) or
  73     (defined?(Rational) and other.kind_of?(Rational))
  74   end
  75
  76   def Complex.polar(r, theta)
  77     Complex(r*Math.cos(theta), r*Math.sin(theta))
  78   end
  79
  80   def initialize(a, b = 0)
  81     raise "non numeric 1st arg `#{a.inspect}'" if !a.kind_of? Numeric
  82     raise "non numeric 2nd arg `#{b.inspect}'" if !b.kind_of? Numeric
  83     @real = a
  84     @image = b
  85   end
  86
  87   def + (other)
  88     if other.kind_of?(Complex)
  89       re = @real + other.real
  90       im = @image + other.image
  91       Complex(re, im)
  92     elsif Complex.generic?(other)
  93       Complex(@real + other, @image)
  94     else
  95       x , y = other.coerce(self)
  96       x + y
  97     end
  98   end
  99
 100   def - (other)
 101     if other.kind_of?(Complex)
 102       re = @real - other.real
 103       im = @image - other.image
 104       Complex(re, im)
 105     elsif Complex.generic?(other)
 106       Complex(@real - other, @image)
 107     else
 108       x , y = other.coerce(self)
 109       x - y
 110     end
 111   end
 112
 113   def * (other)
 114     if other.kind_of?(Complex)
 115       re = @real*other.real - @image*other.image
 116       im = @real*other.image + @image*other.real
 117       Complex(re, im)
 118     elsif Complex.generic?(other)
 119       Complex(@real * other, @image * other)
 120     else
 121       x , y = other.coerce(self)
 122       x * y
 123     end
 124   end
 125
 126   def / (other)
 127     if other.kind_of?(Complex)
 128       self*other.conjugate/other.abs2
 129     elsif Complex.generic?(other)
 130       Complex(@real/other, @image/other)
 131     else
 132       x, y = other.coerce(self)
 133       x/y
 134     end
 135   end
 136
 137   def ** (other)
 138     if other == 0
 139       return Complex(1)
 140     end
 141     if other.kind_of?(Complex)
 142       r, theta = polar
 143       ore = other.real
 144       oim = other.image
 145       nr = Math.exp!(ore*Math.log!(r) - oim * theta)
 146       ntheta = theta*ore + oim*Math.log!(r)
 147       Complex.polar(nr, ntheta)
 148     elsif other.kind_of?(Integer)
 149       if other > 0
 150         x = self
 151         z = x
 152         n = other - 1
 153         while n != 0
 154           while (div, mod = n.divmod(2)
 155                  mod == 0)
 156             x = Complex(x.real*x.real - x.image*x.image, 2*x.real*x.image)
 157             n = div
 158           end
 159           z *= x
 160           n -= 1
 161         end
 162         z
 163       else
 164         if defined? Rational
 165           (Rational(1) / self) ** -other
 166         else
 167           self ** Float(other)
 168         end
 169       end
 170     elsif Complex.generic?(other)
 171       r, theta = polar
 172       Complex.polar(r.power!(other), theta * other)
 173     else
 174       x, y = other.coerce(self)
 175       x/y
 176     end
 177   end
 178
 179   def % (other)
 180     if other.kind_of?(Complex)
 181       Complex(@real % other.real, @image % other.image)
 182     elsif Complex.generic?(other)
 183       Complex(@real % other, @image % other)
 184     else
 185       x , y = other.coerce(self)
 186       x % y
 187     end
 188   end
 189
 190 #    def divmod(other)
 191 #      if other.kind_of?(Complex)
 192 #        rdiv, rmod = @real.divmod(other.real)
 193 #        idiv, imod = @image.divmod(other.image)
 194 #        return Complex(rdiv, idiv), Complex(rmod, rmod)
 195 #      elsif Complex.generic?(other)
 196 #        Complex(@real.divmod(other), @image.divmod(other))
 197 #      else
 198 #        x , y = other.coerce(self)
 199 #        x.divmod(y)
 200 #      end
 201 #    end
 202
 203   def abs
 204     Math.sqrt!((@real*@real + @image*@image).to_f)
 205   end
 206
 207   def abs2
 208     @real*@real + @image*@image
 209   end
 210
 211   def arg
 212     Math.atan2(@image.to_f, @real.to_f)
 213   end
 214
 215   def polar
 216     return abs, arg
 217   end
 218
 219   def conjugate
 220     Complex(@real, -@image)
 221   end
 222
 223   def <=> (other)
 224     self.abs <=> other.abs
 225   end
 226
 227   def == (other)
 228     if other.kind_of?(Complex)
 229       @real == other.real and @image == other.image
 230     elsif Complex.generic?(other)
 231       @real == other and @image == 0
 232     else
 233       x , y = other.coerce(self)
 234       x == y
 235     end
 236   end
 237
 238   def coerce(other)
 239     if Complex.generic?(other)
 240       return Complex.new(other), self
 241     else
 242       super
 243     end
 244   end
 245
 246   def to_i
 247     Complex(@real.to_i, @image.to_i)
 248   end
 249
 250   def to_f
 251     Complex(@real.to_f, @image.to_f)
 252   end
 253
 254   def to_r
 255     Complex(@real.to_r, @image.to_r)
 256   end
 257
 258   def denominator
 259     @real.denominator.lcm(@image.denominator)
 260   end
 261
 262   def numerator
 263     cd = denominator
 264     Complex(@real.numerator*(cd/@real.denominator),
 265             @image.numerator*(cd/@image.denominator))
 266   end
 267
 268   def to_s
 269     if @real != 0
 270       if defined?(Rational) and @image.kind_of?(Rational) and @image.denominator != 1
 271         if @image >= 0
 272           @real.to_s+"+("+@image.to_s+")i"
 273         else
 274           @real.to_s+"-("+(-@image).to_s+")i"
 275         end
 276       else
 277         if @image >= 0
 278           @real.to_s+"+"+@image.to_s+"i"
 279         else
 280           @real.to_s+"-"+(-@image).to_s+"i"
 281         end
 282       end
 283     else
 284       if defined?(Rational) and @image.kind_of?(Rational) and @image.denominator != 1
 285         "("+@image.to_s+")i"
 286       else
 287         @image.to_s+"i"
 288       end
 289     end
 290   end
 291
 292   def hash
 293     @real ^ @image
 294   end
 295
 296   def inspect
 297     sprintf("Complex(%s, %s)", @real.inspect, @image.inspect)
 298   end
 299
 300
 301   I = Complex(0,1)
 302
 303   attr :real
 304   attr :image
 305
 306 end
 307
 308 class Numeric
 309   def im
 310     Complex(0, self)
 311   end
 312
 313   def real
 314     self
 315   end
 316
 317   def image
 318     0
 319   end
 320
 321   def arg
 322     if self >= 0
 323       return 0
 324     else
 325       return Math.atan2(1,1)*4
 326     end
 327   end
 328
 329   def polar
 330     return abs, arg
 331   end
 332
 333   def conjugate
 334     self
 335   end
 336 end
 337
 338 class Fixnum
 339   if not defined? Rational
 340     alias power! **
 341   end
 342
 343   def ** (other)
 344     if self < 0
 345       Complex.new(self) ** other
 346     else
 347       if defined? Rational
 348         if other >= 0
 349           self.power!(other)
 350         else
 351           Rational.new!(self,1)**other
 352         end
 353       else
 354         self.power!(other)
 355       end
 356     end
 357   end
 358 end
 359
 360 class Bignum
 361   if not defined? Rational
 362     alias power! **
 363   end
 364 end
 365
 366 class Float
 367   alias power! **
 368 end
 369
 370 module Math
 371   alias sqrt! sqrt
 372   alias exp! exp
 373   alias cos! cos
 374   alias sin! sin
 375   alias tan! tan
 376   alias log! log
 377   alias log10! log10
 378   alias atan2! atan2
 379
 380   def sqrt(z)
 381     if Complex.generic?(z)
 382       if z >= 0
 383         sqrt!(z)
 384       else
 385         Complex(0,sqrt!(-z))
 386       end
 387     else
 388       z**Rational(1,2)
 389     end
 390   end
 391
 392   def exp(z)
 393     if Complex.generic?(z)
 394       exp!(z)
 395     else
 396       Complex(exp!(z.real) * cos!(z.image), exp!(z.real) * sin!(z.image))
 397     end
 398   end
 399
 400   def cosh!(x)
 401     (exp!(x) + exp!(-x))/2.0
 402   end
 403
 404   def sinh!(x)
 405     (exp!(x) - exp!(-x))/2.0
 406   end
 407
 408   def cos(z)
 409     if Complex.generic?(z)
 410       cos!(z)
 411     else
 412       Complex(cos!(z.real)*cosh!(z.image),
 413               -sin!(z.real)*sinh!(z.image))
 414     end
 415   end
 416
 417   def sin(z)
 418     if Complex.generic?(z)
 419       sin!(z)
 420     else
 421       Complex(sin!(z.real)*cosh!(z.image),
 422               cos!(z.real)*sinh!(z.image))
 423     end
 424   end
 425
 426   def tan(z)
 427     if Complex.generic?(z)
 428       tan!(z)
 429     else
 430       sin(z)/cos(z)
 431     end
 432   end
 433
 434   def log(z)
 435     if Complex.generic?(z) and z >= 0
 436       log!(z)
 437     else
 438       r, theta = z.polar
 439       Complex(log!(r.abs), theta)
 440     end
 441   end
 442
 443   def log10(z)
 444     if Complex.generic?(z)
 445       log10!(z)
 446     else
 447       log(z)/log!(10)
 448     end
 449   end
 450
 451   def atan2(x, y)
 452     if Complex.generic?(x) and Complex.generic?(y)
 453       atan2!(x, y)
 454     else
 455       fail "Not yet implemented."
 456     end
 457   end
 458
 459   def atanh!(x)
 460     log((1.0 + x.to_f) / ( 1.0 - x.to_f)) / 2.0
 461   end
 462
 463   def atan(z)
 464     if Complex.generic?(z)
 465       atan2!(z, 1)
 466     elsif z.image == 0
 467       atan2(z.real,1)
 468     else
 469       a = z.real
 470       b = z.image
 471
 472       c = (a*a + b*b - 1.0)
 473       d = (a*a + b*b + 1.0)
 474
 475       Complex(atan2!((c + sqrt(c*c + 4.0*a*a)), 2.0*a),
 476               atanh!((-d + sqrt(d*d - 4.0*b*b))/(2.0*b)))
 477     end
 478   end
 479
 480   module_function :sqrt
 481   module_function :sqrt!
 482   module_function :exp!
 483   module_function :exp
 484   module_function :cosh!
 485   module_function :cos!
 486   module_function :cos
 487   module_function :sinh!
 488   module_function :sin!
 489   module_function :sin
 490   module_function :tan!
 491   module_function :tan
 492   module_function :log!
 493   module_function :log
 494   module_function :log10!
 495   module_function :log
 496   module_function :atan2!
 497   module_function :atan2
 498 #  module_function :atan!
 499   module_function :atan
 500   module_function :atanh!
 501
 502 end