lib/complex.rb

   1 #
   2 #   complex.rb -
   3 #       $Release Version: 0.5 $
   4 #       $Revision: 1.3 $
   5 #       $Date: 1998/07/08 10:05:28 $
   6 #       by Keiju ISHITSUKA(SHL Japan Inc.)
   7 #
   8 # ----
   9 #
  10 # complex.rb implements the Complex class for complex numbers.  Additionally,
  11 # some methods in other Numeric classes are redefined or added to allow greater
  12 # interoperability with Complex numbers.
  13 #
  14 # Complex numbers can be created in the following manner:
  15 # - <tt>Complex(a, b)</tt>
  16 # - <tt>Complex.polar(radius, theta)</tt>
  17 #
  18 # Additionally, note the following:
  19 # - <tt>Complex::I</tt> (the mathematical constant <i>i</i>)
  20 # - <tt>Numeric#im</tt> (e.g. <tt>5.im -> 0+5i</tt>)
  21 #
  22 # The following +Math+ module methods are redefined to handle Complex arguments.
  23 # They will work as normal with non-Complex arguments.
  24 #    sqrt exp cos sin tan log log10
  25 #    cosh sinh tanh acos asin atan atan2 acosh asinh atanh
  26 #
  27
  28
  29 #
  30 # Numeric is a built-in class on which Fixnum, Bignum, etc., are based.  Here
  31 # some methods are added so that all number types can be treated to some extent
  32 # as Complex numbers.
  33 #
  34 class Numeric
  35   #
  36   # Returns a Complex number <tt>(0,<i>self</i>)</tt>.
  37   #
  38   def im
  39     Complex(0, self)
  40   end
  41
  42   #
  43   # The real part of a complex number, i.e. <i>self</i>.
  44   #
  45   def real
  46     self
  47   end
  48
  49   #
  50   # The imaginary part of a complex number, i.e. 0.
  51   #
  52   def image
  53     0
  54   end
  55   alias imag image
  56
  57   #
  58   # See Complex#arg.
  59   #
  60   def arg
  61     if self >= 0
  62       return 0
  63     else
  64       return Math::PI
  65     end
  66   end
  67   alias angle arg
  68
  69   #
  70   # See Complex#polar.
  71   #
  72   def polar
  73     return abs, arg
  74   end
  75
  76   #
  77   # See Complex#conjugate (short answer: returns <i>self</i>).
  78   #
  79   def conjugate
  80     self
  81   end
  82   alias conj conjugate
  83 end
  84
  85
  86 #
  87 # Creates a Complex number.  +a+ and +b+ should be Numeric.  The result will be
  88 # <tt>a+bi</tt>.
  89 #
  90 def Complex(a, b = 0)
  91   if b == 0 and (a.kind_of?(Complex) or defined? Complex::Unify)
  92     a
  93   else
  94     Complex.new( a.real-b.imag, a.imag+b.real )
  95   end
  96 end
  97
  98 #
  99 # The complex number class.  See complex.rb for an overview.
 100 #
 101 class Complex < Numeric
 102   @RCS_ID='-$Id: complex.rb,v 1.3 1998/07/08 10:05:28 keiju Exp keiju $-'
 103
 104   undef step
 105   undef div, divmod
 106   undef floor, truncate, ceil, round
 107
 108   def Complex.generic?(other) # :nodoc:
 109     other.kind_of?(Integer) or
 110     other.kind_of?(Float) or
 111     (defined?(Rational) and other.kind_of?(Rational))
 112   end
 113
 114   #
 115   # Creates a +Complex+ number in terms of +r+ (radius) and +theta+ (angle).
 116   #
 117   def Complex.polar(r, theta)
 118     Complex(r*Math.cos(theta), r*Math.sin(theta))
 119   end
 120
 121   #
 122   # Creates a +Complex+ number <tt>a</tt>+<tt>b</tt><i>i</i>.
 123   #
 124   def Complex.new!(a, b=0)
 125     new(a,b)
 126   end
 127
 128   def initialize(a, b)
 129     raise TypeError, "non numeric 1st arg `#{a.inspect}'" if !a.kind_of? Numeric
 130     raise TypeError, "`#{a.inspect}' for 1st arg" if a.kind_of? Complex
 131     raise TypeError, "non numeric 2nd arg `#{b.inspect}'" if !b.kind_of? Numeric
 132     raise TypeError, "`#{b.inspect}' for 2nd arg" if b.kind_of? Complex
 133     @real = a
 134     @image = b
 135   end
 136
 137   #
 138   # Addition with real or complex number.
 139   #
 140   def + (other)
 141     if other.kind_of?(Complex)
 142       re = @real + other.real
 143       im = @image + other.image
 144       Complex(re, im)
 145     elsif Complex.generic?(other)
 146       Complex(@real + other, @image)
 147     else
 148       x , y = other.coerce(self)
 149       x + y
 150     end
 151   end
 152
 153   #
 154   # Subtraction with real or complex number.
 155   #
 156   def - (other)
 157     if other.kind_of?(Complex)
 158       re = @real - other.real
 159       im = @image - other.image
 160       Complex(re, im)
 161     elsif Complex.generic?(other)
 162       Complex(@real - other, @image)
 163     else
 164       x , y = other.coerce(self)
 165       x - y
 166     end
 167   end
 168
 169   #
 170   # Multiplication with real or complex number.
 171   #
 172   def * (other)
 173     if other.kind_of?(Complex)
 174       re = @real*other.real - @image*other.image
 175       im = @real*other.image + @image*other.real
 176       Complex(re, im)
 177     elsif Complex.generic?(other)
 178       Complex(@real * other, @image * other)
 179     else
 180       x , y = other.coerce(self)
 181       x * y
 182     end
 183   end
 184
 185   #
 186   # Division by real or complex number.
 187   #
 188   def / (other)
 189     if other.kind_of?(Complex)
 190       self*other.conjugate/other.abs2
 191     elsif Complex.generic?(other)
 192       Complex(@real/other, @image/other)
 193     else
 194       x, y = other.coerce(self)
 195       x/y
 196     end
 197   end
 198
 199   def quo(other)
 200     Complex(@real.quo(1), @image.quo(1)) / other
 201   end
 202
 203   #
 204   # Raise this complex number to the given (real or complex) power.
 205   #
 206   def ** (other)
 207     if other == 0
 208       return Complex(1)
 209     end
 210     if other.kind_of?(Complex)
 211       r, theta = polar
 212       ore = other.real
 213       oim = other.image
 214       nr = Math.exp!(ore*Math.log!(r) - oim * theta)
 215       ntheta = theta*ore + oim*Math.log!(r)
 216       Complex.polar(nr, ntheta)
 217     elsif other.kind_of?(Integer)
 218       if other > 0
 219         x = self
 220         z = x
 221         n = other - 1
 222         while n != 0
 223           while (div, mod = n.divmod(2)
 224                  mod == 0)
 225             x = Complex(x.real*x.real - x.image*x.image, 2*x.real*x.image)
 226             n = div
 227           end
 228           z *= x
 229           n -= 1
 230         end
 231         z
 232       else
 233         if defined? Rational
 234           (Rational(1) / self) ** -other
 235         else
 236           self ** Float(other)
 237         end
 238       end
 239     elsif Complex.generic?(other)
 240       r, theta = polar
 241       Complex.polar(r**other, theta*other)
 242     else
 243       x, y = other.coerce(self)
 244       x**y
 245     end
 246   end
 247
 248   #
 249   # Remainder after division by a real or complex number.
 250   #
 251   def % (other)
 252     if other.kind_of?(Complex)
 253       Complex(@real % other.real, @image % other.image)
 254     elsif Complex.generic?(other)
 255       Complex(@real % other, @image % other)
 256     else
 257       x , y = other.coerce(self)
 258       x % y
 259     end
 260   end
 261
 262 #--
 263 #    def divmod(other)
 264 #      if other.kind_of?(Complex)
 265 #        rdiv, rmod = @real.divmod(other.real)
 266 #        idiv, imod = @image.divmod(other.image)
 267 #        return Complex(rdiv, idiv), Complex(rmod, rmod)
 268 #      elsif Complex.generic?(other)
 269 #        Complex(@real.divmod(other), @image.divmod(other))
 270 #      else
 271 #        x , y = other.coerce(self)
 272 #        x.divmod(y)
 273 #      end
 274 #    end
 275 #++
 276
 277   #
 278   # Absolute value (aka modulus): distance from the zero point on the complex
 279   # plane.
 280   #
 281   def abs
 282     Math.hypot(@real, @image)
 283   end
 284
 285   #
 286   # Square of the absolute value.
 287   #
 288   def abs2
 289     @real*@real + @image*@image
 290   end
 291
 292   #
 293   # Argument (angle from (1,0) on the complex plane).
 294   #
 295   def arg
 296     Math.atan2!(@image, @real)
 297   end
 298   alias angle arg
 299
 300   #
 301   # Returns the absolute value _and_ the argument.
 302   #
 303   def polar
 304     return abs, arg
 305   end
 306
 307   #
 308   # Complex conjugate (<tt>z + z.conjugate = 2 * z.real</tt>).
 309   #
 310   def conjugate
 311     Complex(@real, -@image)
 312   end
 313   alias conj conjugate
 314
 315   #
 316   # Compares the absolute values of the two numbers.
 317   #
 318   def <=> (other)
 319     self.abs <=> other.abs
 320   end
 321
 322   #
 323   # Test for numerical equality (<tt>a == a + 0<i>i</i></tt>).
 324   #
 325   def == (other)
 326     if other.kind_of?(Complex)
 327       @real == other.real and @image == other.image
 328     elsif Complex.generic?(other)
 329       @real == other and @image == 0
 330     else
 331       other == self
 332     end
 333   end
 334
 335   #
 336   # Attempts to coerce +other+ to a Complex number.
 337   #
 338   def coerce(other)
 339     if Complex.generic?(other)
 340       return Complex.new!(other), self
 341     else
 342       super
 343     end
 344   end
 345
 346   #
 347   # FIXME
 348   #
 349   def denominator
 350     @real.denominator.lcm(@image.denominator)
 351   end
 352
 353   #
 354   # FIXME
 355   #
 356   def numerator
 357     cd = denominator
 358     Complex(@real.numerator*(cd/@real.denominator),
 359             @image.numerator*(cd/@image.denominator))
 360   end
 361
 362   #
 363   # Standard string representation of the complex number.
 364   #
 365   def to_s
 366     if @real != 0
 367       if defined?(Rational) and @image.kind_of?(Rational) and @image.denominator != 1
 368         if @image >= 0
 369           @real.to_s+"+("+@image.to_s+")i"
 370         else
 371           @real.to_s+"-("+(-@image).to_s+")i"
 372         end
 373       else
 374         if @image >= 0
 375           @real.to_s+"+"+@image.to_s+"i"
 376         else
 377           @real.to_s+"-"+(-@image).to_s+"i"
 378         end
 379       end
 380     else
 381       if defined?(Rational) and @image.kind_of?(Rational) and @image.denominator != 1
 382         "("+@image.to_s+")i"
 383       else
 384         @image.to_s+"i"
 385       end
 386     end
 387   end
 388
 389   #
 390   # Returns a hash code for the complex number.
 391   #
 392   def hash
 393     @real.hash ^ @image.hash
 394   end
 395
 396   #
 397   # Returns "<tt>Complex(<i>real</i>, <i>image</i>)</tt>".
 398   #
 399   def inspect
 400     sprintf("Complex(%s, %s)", @real.inspect, @image.inspect)
 401   end
 402
 403
 404   #
 405   # +I+ is the imaginary number.  It exists at point (0,1) on the complex plane.
 406   #
 407   I = Complex(0,1)
 408
 409   # The real part of a complex number.
 410   attr :real
 411
 412   # The imaginary part of a complex number.
 413   attr :image
 414   alias imag image
 415
 416 end
 417
 418 class Integer
 419
 420   unless defined?(1.numerator)
 421     def numerator() self end
 422     def denominator() 1 end
 423
 424     def gcd(other)
 425       min = self.abs
 426       max = other.abs
 427       while min > 0
 428         tmp = min
 429         min = max % min
 430         max = tmp
 431       end
 432       max
 433     end
 434
 435     def lcm(other)
 436       if self.zero? or other.zero?
 437         0
 438       else
 439         (self.div(self.gcd(other)) * other).abs
 440       end
 441     end
 442
 443   end
 444
 445 end
 446
 447 module Math
 448   alias sqrt! sqrt
 449   alias exp! exp
 450   alias log! log
 451   alias log10! log10
 452   alias cos! cos
 453   alias sin! sin
 454   alias tan! tan
 455   alias cosh! cosh
 456   alias sinh! sinh
 457   alias tanh! tanh
 458   alias acos! acos
 459   alias asin! asin
 460   alias atan! atan
 461   alias atan2! atan2
 462   alias acosh! acosh
 463   alias asinh! asinh
 464   alias atanh! atanh
 465
 466   # Redefined to handle a Complex argument.
 467   def sqrt(z)
 468     if Complex.generic?(z)
 469       if z >= 0
 470         sqrt!(z)
 471       else
 472         Complex(0,sqrt!(-z))
 473       end
 474     else
 475       if z.image < 0
 476         sqrt(z.conjugate).conjugate
 477       else
 478         r = z.abs
 479         x = z.real
 480         Complex( sqrt!((r+x)/2), sqrt!((r-x)/2) )
 481       end
 482     end
 483   end
 484
 485   # Redefined to handle a Complex argument.
 486   def exp(z)
 487     if Complex.generic?(z)
 488       exp!(z)
 489     else
 490       Complex(exp!(z.real) * cos!(z.image), exp!(z.real) * sin!(z.image))
 491     end
 492   end
 493
 494   # Redefined to handle a Complex argument.
 495   def cos(z)
 496     if Complex.generic?(z)
 497       cos!(z)
 498     else
 499       Complex(cos!(z.real)*cosh!(z.image),
 500               -sin!(z.real)*sinh!(z.image))
 501     end
 502   end
 503
 504   # Redefined to handle a Complex argument.
 505   def sin(z)
 506     if Complex.generic?(z)
 507       sin!(z)
 508     else
 509       Complex(sin!(z.real)*cosh!(z.image),
 510               cos!(z.real)*sinh!(z.image))
 511     end
 512   end
 513
 514   # Redefined to handle a Complex argument.
 515   def tan(z)
 516     if Complex.generic?(z)
 517       tan!(z)
 518     else
 519       sin(z)/cos(z)
 520     end
 521   end
 522
 523   def sinh(z)
 524     if Complex.generic?(z)
 525       sinh!(z)
 526     else
 527       Complex( sinh!(z.real)*cos!(z.image), cosh!(z.real)*sin!(z.image) )
 528     end
 529   end
 530
 531   def cosh(z)
 532     if Complex.generic?(z)
 533       cosh!(z)
 534     else
 535       Complex( cosh!(z.real)*cos!(z.image), sinh!(z.real)*sin!(z.image) )
 536     end
 537   end
 538
 539   def tanh(z)
 540     if Complex.generic?(z)
 541       tanh!(z)
 542     else
 543       sinh(z)/cosh(z)
 544     end
 545   end
 546
 547   # Redefined to handle a Complex argument.
 548   def log(z)
 549     if Complex.generic?(z) and z >= 0
 550       log!(z)
 551     else
 552       r, theta = z.polar
 553       Complex(log!(r.abs), theta)
 554     end
 555   end
 556
 557   # Redefined to handle a Complex argument.
 558   def log10(z)
 559     if Complex.generic?(z)
 560       log10!(z)
 561     else
 562       log(z)/log!(10)
 563     end
 564   end
 565
 566   def acos(z)
 567     if Complex.generic?(z) and z >= -1 and z <= 1
 568       acos!(z)
 569     else
 570       -1.0.im * log( z + 1.0.im * sqrt(1.0-z*z) )
 571     end
 572   end
 573
 574   def asin(z)
 575     if Complex.generic?(z) and z >= -1 and z <= 1
 576       asin!(z)
 577     else
 578       -1.0.im * log( 1.0.im * z + sqrt(1.0-z*z) )
 579     end
 580   end
 581
 582   def atan(z)
 583     if Complex.generic?(z)
 584       atan!(z)
 585     else
 586       1.0.im * log( (1.0.im+z) / (1.0.im-z) ) / 2.0
 587     end
 588   end
 589
 590   def atan2(y,x)
 591     if Complex.generic?(y) and Complex.generic?(x)
 592       atan2!(y,x)
 593     else
 594       -1.0.im * log( (x+1.0.im*y) / sqrt(x*x+y*y) )
 595     end
 596   end
 597
 598   def acosh(z)
 599     if Complex.generic?(z) and z >= 1
 600       acosh!(z)
 601     else
 602       log( z + sqrt(z*z-1.0) )
 603     end
 604   end
 605
 606   def asinh(z)
 607     if Complex.generic?(z)
 608       asinh!(z)
 609     else
 610       log( z + sqrt(1.0+z*z) )
 611     end
 612   end
 613
 614   def atanh(z)
 615     if Complex.generic?(z) and z >= -1 and z <= 1
 616       atanh!(z)
 617     else
 618       log( (1.0+z) / (1.0-z) ) / 2.0
 619     end
 620   end
 621
 622   module_function :sqrt!
 623   module_function :sqrt
 624   module_function :exp!
 625   module_function :exp
 626   module_function :log!
 627   module_function :log
 628   module_function :log10!
 629   module_function :log10
 630   module_function :cosh!
 631   module_function :cosh
 632   module_function :cos!
 633   module_function :cos
 634   module_function :sinh!
 635   module_function :sinh
 636   module_function :sin!
 637   module_function :sin
 638   module_function :tan!
 639   module_function :tan
 640   module_function :tanh!
 641   module_function :tanh
 642   module_function :acos!
 643   module_function :acos
 644   module_function :asin!
 645   module_function :asin
 646   module_function :atan!
 647   module_function :atan
 648   module_function :atan2!
 649   module_function :atan2
 650   module_function :acosh!
 651   module_function :acosh
 652   module_function :asinh!
 653   module_function :asinh
 654   module_function :atanh!
 655   module_function :atanh
 656
 657 end
 658
 659 # Documentation comments:
 660 #  - source: original (researched from pickaxe)
 661 #  - a couple of fixme's
 662 #  - RDoc output for Bignum etc. is a bit short, with nothing but an
 663 #    (undocumented) alias.  No big deal.