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