Commit | Line | Data |
66730be0 |
1 | # $RCSFile$ |
2 | # |
3 | # Complex numbers and associated mathematical functions |
0c721ce2 |
4 | # -- Raphael Manfredi, September 1996 |
5cd24f17 |
5 | # -- Jarkko Hietaniemi, March-April 1997 |
a0d0e21e |
6 | |
7 | require Exporter; |
5aabfad6 |
8 | package Math::Complex; |
a0d0e21e |
9 | |
0c721ce2 |
10 | use strict; |
11 | |
5aabfad6 |
12 | use vars qw($VERSION @ISA |
13 | @EXPORT %EXPORT_TAGS |
14 | $package $display |
5cd24f17 |
15 | $i $logn %logn); |
0c721ce2 |
16 | |
5aabfad6 |
17 | @ISA = qw(Exporter); |
18 | |
19 | $VERSION = 1.01; |
20 | |
21 | my @trig = qw( |
22 | pi |
5cd24f17 |
23 | sin cos tan |
5aabfad6 |
24 | csc cosec sec cot cotan |
25 | asin acos atan |
26 | acsc acosec asec acot acotan |
27 | sinh cosh tanh |
28 | csch cosech sech coth cotanh |
29 | asinh acosh atanh |
30 | acsch acosech asech acoth acotanh |
31 | ); |
32 | |
33 | @EXPORT = (qw( |
34 | i Re Im arg |
35 | sqrt exp log ln |
36 | log10 logn cbrt root |
37 | cplx cplxe |
38 | ), |
39 | @trig); |
40 | |
41 | %EXPORT_TAGS = ( |
42 | 'trig' => [@trig], |
66730be0 |
43 | ); |
a0d0e21e |
44 | |
a5f75d66 |
45 | use overload |
0c721ce2 |
46 | '+' => \&plus, |
47 | '-' => \&minus, |
48 | '*' => \&multiply, |
49 | '/' => \÷, |
66730be0 |
50 | '**' => \&power, |
51 | '<=>' => \&spaceship, |
52 | 'neg' => \&negate, |
0c721ce2 |
53 | '~' => \&conjugate, |
66730be0 |
54 | 'abs' => \&abs, |
55 | 'sqrt' => \&sqrt, |
56 | 'exp' => \&exp, |
57 | 'log' => \&log, |
58 | 'sin' => \&sin, |
59 | 'cos' => \&cos, |
0c721ce2 |
60 | 'tan' => \&tan, |
66730be0 |
61 | 'atan2' => \&atan2, |
62 | qw("" stringify); |
63 | |
64 | # |
65 | # Package globals |
66 | # |
67 | |
68 | $package = 'Math::Complex'; # Package name |
69 | $display = 'cartesian'; # Default display format |
70 | |
71 | # |
72 | # Object attributes (internal): |
73 | # cartesian [real, imaginary] -- cartesian form |
74 | # polar [rho, theta] -- polar form |
75 | # c_dirty cartesian form not up-to-date |
76 | # p_dirty polar form not up-to-date |
77 | # display display format (package's global when not set) |
78 | # |
79 | |
80 | # |
81 | # ->make |
82 | # |
83 | # Create a new complex number (cartesian form) |
84 | # |
85 | sub make { |
86 | my $self = bless {}, shift; |
87 | my ($re, $im) = @_; |
40da2db3 |
88 | $self->{'cartesian'} = [$re, $im]; |
66730be0 |
89 | $self->{c_dirty} = 0; |
90 | $self->{p_dirty} = 1; |
91 | return $self; |
92 | } |
93 | |
94 | # |
95 | # ->emake |
96 | # |
97 | # Create a new complex number (exponential form) |
98 | # |
99 | sub emake { |
100 | my $self = bless {}, shift; |
101 | my ($rho, $theta) = @_; |
102 | $theta += pi() if $rho < 0; |
40da2db3 |
103 | $self->{'polar'} = [abs($rho), $theta]; |
66730be0 |
104 | $self->{p_dirty} = 0; |
105 | $self->{c_dirty} = 1; |
106 | return $self; |
107 | } |
108 | |
109 | sub new { &make } # For backward compatibility only. |
110 | |
111 | # |
112 | # cplx |
113 | # |
114 | # Creates a complex number from a (re, im) tuple. |
115 | # This avoids the burden of writing Math::Complex->make(re, im). |
116 | # |
117 | sub cplx { |
118 | my ($re, $im) = @_; |
0c721ce2 |
119 | return $package->make($re, defined $im ? $im : 0); |
66730be0 |
120 | } |
121 | |
122 | # |
123 | # cplxe |
124 | # |
125 | # Creates a complex number from a (rho, theta) tuple. |
126 | # This avoids the burden of writing Math::Complex->emake(rho, theta). |
127 | # |
128 | sub cplxe { |
129 | my ($rho, $theta) = @_; |
0c721ce2 |
130 | return $package->emake($rho, defined $theta ? $theta : 0); |
66730be0 |
131 | } |
132 | |
133 | # |
134 | # pi |
135 | # |
136 | # The number defined as 2 * pi = 360 degrees |
137 | # |
5cd24f17 |
138 | |
139 | use constant pi => 4 * atan2(1, 1); |
140 | |
141 | # |
142 | # log2inv |
143 | # |
144 | # Used in log10(). |
145 | # |
146 | |
147 | use constant log10inv => 1 / log(10); |
66730be0 |
148 | |
149 | # |
150 | # i |
151 | # |
152 | # The number defined as i*i = -1; |
153 | # |
154 | sub i () { |
5cd24f17 |
155 | return $i if ($i); |
156 | $i = bless {}; |
40da2db3 |
157 | $i->{'cartesian'} = [0, 1]; |
5cd24f17 |
158 | $i->{'polar'} = [1, pi/2]; |
66730be0 |
159 | $i->{c_dirty} = 0; |
160 | $i->{p_dirty} = 0; |
161 | return $i; |
162 | } |
163 | |
164 | # |
165 | # Attribute access/set routines |
166 | # |
167 | |
0c721ce2 |
168 | sub cartesian {$_[0]->{c_dirty} ? |
169 | $_[0]->update_cartesian : $_[0]->{'cartesian'}} |
170 | sub polar {$_[0]->{p_dirty} ? |
171 | $_[0]->update_polar : $_[0]->{'polar'}} |
66730be0 |
172 | |
40da2db3 |
173 | sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] } |
174 | sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] } |
66730be0 |
175 | |
176 | # |
177 | # ->update_cartesian |
178 | # |
179 | # Recompute and return the cartesian form, given accurate polar form. |
180 | # |
181 | sub update_cartesian { |
182 | my $self = shift; |
40da2db3 |
183 | my ($r, $t) = @{$self->{'polar'}}; |
66730be0 |
184 | $self->{c_dirty} = 0; |
40da2db3 |
185 | return $self->{'cartesian'} = [$r * cos $t, $r * sin $t]; |
66730be0 |
186 | } |
187 | |
188 | # |
189 | # |
190 | # ->update_polar |
191 | # |
192 | # Recompute and return the polar form, given accurate cartesian form. |
193 | # |
194 | sub update_polar { |
195 | my $self = shift; |
40da2db3 |
196 | my ($x, $y) = @{$self->{'cartesian'}}; |
66730be0 |
197 | $self->{p_dirty} = 0; |
40da2db3 |
198 | return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0; |
199 | return $self->{'polar'} = [sqrt($x*$x + $y*$y), atan2($y, $x)]; |
66730be0 |
200 | } |
201 | |
202 | # |
203 | # (plus) |
204 | # |
205 | # Computes z1+z2. |
206 | # |
207 | sub plus { |
208 | my ($z1, $z2, $regular) = @_; |
209 | my ($re1, $im1) = @{$z1->cartesian}; |
0e505df1 |
210 | $z2 = cplx($z2) unless ref $z2; |
5cd24f17 |
211 | my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); |
66730be0 |
212 | unless (defined $regular) { |
213 | $z1->set_cartesian([$re1 + $re2, $im1 + $im2]); |
214 | return $z1; |
215 | } |
216 | return (ref $z1)->make($re1 + $re2, $im1 + $im2); |
217 | } |
218 | |
219 | # |
220 | # (minus) |
221 | # |
222 | # Computes z1-z2. |
223 | # |
224 | sub minus { |
225 | my ($z1, $z2, $inverted) = @_; |
226 | my ($re1, $im1) = @{$z1->cartesian}; |
0e505df1 |
227 | $z2 = cplx($z2) unless ref $z2; |
228 | my ($re2, $im2) = @{$z2->cartesian}; |
66730be0 |
229 | unless (defined $inverted) { |
230 | $z1->set_cartesian([$re1 - $re2, $im1 - $im2]); |
231 | return $z1; |
232 | } |
233 | return $inverted ? |
234 | (ref $z1)->make($re2 - $re1, $im2 - $im1) : |
235 | (ref $z1)->make($re1 - $re2, $im1 - $im2); |
0e505df1 |
236 | |
66730be0 |
237 | } |
238 | |
239 | # |
240 | # (multiply) |
241 | # |
242 | # Computes z1*z2. |
243 | # |
244 | sub multiply { |
245 | my ($z1, $z2, $regular) = @_; |
246 | my ($r1, $t1) = @{$z1->polar}; |
0e505df1 |
247 | $z2 = cplxe(abs($z2), $z2 >= 0 ? 0 : pi) unless ref $z2; |
248 | my ($r2, $t2) = @{$z2->polar}; |
66730be0 |
249 | unless (defined $regular) { |
250 | $z1->set_polar([$r1 * $r2, $t1 + $t2]); |
251 | return $z1; |
252 | } |
253 | return (ref $z1)->emake($r1 * $r2, $t1 + $t2); |
254 | } |
255 | |
256 | # |
0e505df1 |
257 | # _divbyzero |
0c721ce2 |
258 | # |
259 | # Die on division by zero. |
260 | # |
0e505df1 |
261 | sub _divbyzero { |
5cd24f17 |
262 | my $mess = "$_[0]: Division by zero.\n"; |
263 | |
264 | if (defined $_[1]) { |
265 | $mess .= "(Because in the definition of $_[0], the divisor "; |
266 | $mess .= "$_[1] " unless ($_[1] eq '0'); |
267 | $mess .= "is 0)\n"; |
268 | } |
269 | |
0c721ce2 |
270 | my @up = caller(1); |
5cd24f17 |
271 | |
272 | $mess .= "Died at $up[1] line $up[2].\n"; |
273 | |
274 | die $mess; |
0c721ce2 |
275 | } |
276 | |
277 | # |
66730be0 |
278 | # (divide) |
279 | # |
280 | # Computes z1/z2. |
281 | # |
282 | sub divide { |
283 | my ($z1, $z2, $inverted) = @_; |
284 | my ($r1, $t1) = @{$z1->polar}; |
0e505df1 |
285 | $z2 = cplxe(abs($z2), $z2 >= 0 ? 0 : pi) unless ref $z2; |
286 | my ($r2, $t2) = @{$z2->polar}; |
66730be0 |
287 | unless (defined $inverted) { |
0e505df1 |
288 | _divbyzero "$z1/0" if ($r2 == 0); |
66730be0 |
289 | $z1->set_polar([$r1 / $r2, $t1 - $t2]); |
290 | return $z1; |
291 | } |
0c721ce2 |
292 | if ($inverted) { |
0e505df1 |
293 | _divbyzero "$z2/0" if ($r1 == 0); |
0c721ce2 |
294 | return (ref $z1)->emake($r2 / $r1, $t2 - $t1); |
295 | } else { |
0e505df1 |
296 | _divbyzero "$z1/0" if ($r2 == 0); |
0c721ce2 |
297 | return (ref $z1)->emake($r1 / $r2, $t1 - $t2); |
298 | } |
66730be0 |
299 | } |
300 | |
301 | # |
0e505df1 |
302 | # _zerotozero |
303 | # |
304 | # Die on zero raised to the zeroth. |
305 | # |
306 | sub _zerotozero { |
307 | my $mess = "The zero raised to the zeroth power is not defined.\n"; |
308 | |
309 | my @up = caller(1); |
310 | |
311 | $mess .= "Died at $up[1] line $up[2].\n"; |
312 | |
313 | die $mess; |
314 | } |
315 | |
316 | # |
66730be0 |
317 | # (power) |
318 | # |
319 | # Computes z1**z2 = exp(z2 * log z1)). |
320 | # |
321 | sub power { |
322 | my ($z1, $z2, $inverted) = @_; |
ace5de91 |
323 | my $z1z = $z1 == 0; |
324 | my $z2z = $z2 == 0; |
325 | _zerotozero if ($z1z and $z2z); |
326 | if ($inverted) { |
327 | return 0 if ($z2z); |
328 | return 1 if ($z1z or $z2 == 1); |
329 | } else { |
330 | return 0 if ($z1z); |
331 | return 1 if ($z2z or $z1 == 1); |
332 | } |
0e505df1 |
333 | $z2 = cplx($z2) unless ref $z2; |
334 | unless (defined $inverted) { |
335 | my $z3 = exp($z2 * log $z1); |
336 | $z1->set_cartesian([@{$z3->cartesian}]); |
337 | return $z1; |
338 | } |
339 | return exp($z2 * log $z1) unless $inverted; |
340 | return exp($z1 * log $z2); |
66730be0 |
341 | } |
342 | |
343 | # |
344 | # (spaceship) |
345 | # |
346 | # Computes z1 <=> z2. |
347 | # Sorts on the real part first, then on the imaginary part. Thus 2-4i > 3+8i. |
348 | # |
349 | sub spaceship { |
350 | my ($z1, $z2, $inverted) = @_; |
5cd24f17 |
351 | my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0); |
352 | my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); |
66730be0 |
353 | my $sgn = $inverted ? -1 : 1; |
354 | return $sgn * ($re1 <=> $re2) if $re1 != $re2; |
355 | return $sgn * ($im1 <=> $im2); |
356 | } |
357 | |
358 | # |
359 | # (negate) |
360 | # |
361 | # Computes -z. |
362 | # |
363 | sub negate { |
364 | my ($z) = @_; |
365 | if ($z->{c_dirty}) { |
366 | my ($r, $t) = @{$z->polar}; |
367 | return (ref $z)->emake($r, pi + $t); |
368 | } |
369 | my ($re, $im) = @{$z->cartesian}; |
370 | return (ref $z)->make(-$re, -$im); |
371 | } |
372 | |
373 | # |
374 | # (conjugate) |
375 | # |
376 | # Compute complex's conjugate. |
377 | # |
378 | sub conjugate { |
379 | my ($z) = @_; |
380 | if ($z->{c_dirty}) { |
381 | my ($r, $t) = @{$z->polar}; |
382 | return (ref $z)->emake($r, -$t); |
383 | } |
384 | my ($re, $im) = @{$z->cartesian}; |
385 | return (ref $z)->make($re, -$im); |
386 | } |
387 | |
388 | # |
389 | # (abs) |
390 | # |
391 | # Compute complex's norm (rho). |
392 | # |
393 | sub abs { |
394 | my ($z) = @_; |
0c721ce2 |
395 | return abs($z) unless ref $z; |
66730be0 |
396 | my ($r, $t) = @{$z->polar}; |
397 | return abs($r); |
398 | } |
399 | |
400 | # |
401 | # arg |
402 | # |
403 | # Compute complex's argument (theta). |
404 | # |
405 | sub arg { |
406 | my ($z) = @_; |
0c721ce2 |
407 | return ($z < 0 ? pi : 0) unless ref $z; |
66730be0 |
408 | my ($r, $t) = @{$z->polar}; |
409 | return $t; |
410 | } |
411 | |
412 | # |
413 | # (sqrt) |
414 | # |
0c721ce2 |
415 | # Compute sqrt(z). |
66730be0 |
416 | # |
417 | sub sqrt { |
418 | my ($z) = @_; |
0c721ce2 |
419 | $z = cplx($z, 0) unless ref $z; |
66730be0 |
420 | my ($r, $t) = @{$z->polar}; |
421 | return (ref $z)->emake(sqrt($r), $t/2); |
422 | } |
423 | |
424 | # |
425 | # cbrt |
426 | # |
0c721ce2 |
427 | # Compute cbrt(z) (cubic root). |
66730be0 |
428 | # |
429 | sub cbrt { |
430 | my ($z) = @_; |
0c721ce2 |
431 | return cplx($z, 0) ** (1/3) unless ref $z; |
66730be0 |
432 | my ($r, $t) = @{$z->polar}; |
433 | return (ref $z)->emake($r**(1/3), $t/3); |
434 | } |
435 | |
436 | # |
0e505df1 |
437 | # _rootbad |
438 | # |
439 | # Die on bad root. |
440 | # |
441 | sub _rootbad { |
442 | my $mess = "Root $_[0] not defined, root must be positive integer.\n"; |
443 | |
444 | my @up = caller(1); |
445 | |
446 | $mess .= "Died at $up[1] line $up[2].\n"; |
447 | |
448 | die $mess; |
449 | } |
450 | |
451 | # |
66730be0 |
452 | # root |
453 | # |
454 | # Computes all nth root for z, returning an array whose size is n. |
455 | # `n' must be a positive integer. |
456 | # |
457 | # The roots are given by (for k = 0..n-1): |
458 | # |
459 | # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n)) |
460 | # |
461 | sub root { |
462 | my ($z, $n) = @_; |
0e505df1 |
463 | _rootbad($n) if ($n < 1 or int($n) != $n); |
66730be0 |
464 | my ($r, $t) = ref $z ? @{$z->polar} : (abs($z), $z >= 0 ? 0 : pi); |
465 | my @root; |
466 | my $k; |
467 | my $theta_inc = 2 * pi / $n; |
468 | my $rho = $r ** (1/$n); |
469 | my $theta; |
470 | my $complex = ref($z) || $package; |
471 | for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) { |
472 | push(@root, $complex->emake($rho, $theta)); |
a0d0e21e |
473 | } |
66730be0 |
474 | return @root; |
a0d0e21e |
475 | } |
476 | |
66730be0 |
477 | # |
478 | # Re |
479 | # |
480 | # Return Re(z). |
481 | # |
a0d0e21e |
482 | sub Re { |
66730be0 |
483 | my ($z) = @_; |
484 | return $z unless ref $z; |
485 | my ($re, $im) = @{$z->cartesian}; |
486 | return $re; |
a0d0e21e |
487 | } |
488 | |
66730be0 |
489 | # |
490 | # Im |
491 | # |
492 | # Return Im(z). |
493 | # |
a0d0e21e |
494 | sub Im { |
66730be0 |
495 | my ($z) = @_; |
496 | return 0 unless ref $z; |
497 | my ($re, $im) = @{$z->cartesian}; |
498 | return $im; |
a0d0e21e |
499 | } |
500 | |
66730be0 |
501 | # |
502 | # (exp) |
503 | # |
504 | # Computes exp(z). |
505 | # |
506 | sub exp { |
507 | my ($z) = @_; |
0c721ce2 |
508 | $z = cplx($z, 0) unless ref $z; |
66730be0 |
509 | my ($x, $y) = @{$z->cartesian}; |
510 | return (ref $z)->emake(exp($x), $y); |
511 | } |
512 | |
513 | # |
8c03c583 |
514 | # _logofzero |
515 | # |
516 | # Die on division by zero. |
517 | # |
518 | sub _logofzero { |
519 | my $mess = "$_[0]: Logarithm of zero.\n"; |
520 | |
521 | if (defined $_[1]) { |
522 | $mess .= "(Because in the definition of $_[0], the argument "; |
523 | $mess .= "$_[1] " unless ($_[1] eq '0'); |
524 | $mess .= "is 0)\n"; |
525 | } |
526 | |
527 | my @up = caller(1); |
528 | |
529 | $mess .= "Died at $up[1] line $up[2].\n"; |
530 | |
531 | die $mess; |
532 | } |
533 | |
534 | # |
66730be0 |
535 | # (log) |
536 | # |
537 | # Compute log(z). |
538 | # |
539 | sub log { |
540 | my ($z) = @_; |
0c721ce2 |
541 | $z = cplx($z, 0) unless ref $z; |
0c721ce2 |
542 | my ($x, $y) = @{$z->cartesian}; |
5cd24f17 |
543 | my ($r, $t) = @{$z->polar}; |
0c721ce2 |
544 | $t -= 2 * pi if ($t > pi() and $x < 0); |
545 | $t += 2 * pi if ($t < -pi() and $x < 0); |
66730be0 |
546 | return (ref $z)->make(log($r), $t); |
547 | } |
548 | |
549 | # |
0c721ce2 |
550 | # ln |
551 | # |
552 | # Alias for log(). |
553 | # |
554 | sub ln { Math::Complex::log(@_) } |
555 | |
556 | # |
66730be0 |
557 | # log10 |
558 | # |
559 | # Compute log10(z). |
560 | # |
5cd24f17 |
561 | |
66730be0 |
562 | sub log10 { |
563 | my ($z) = @_; |
5cd24f17 |
564 | |
565 | return log(cplx($z, 0)) * log10inv unless ref $z; |
66730be0 |
566 | my ($r, $t) = @{$z->polar}; |
5cd24f17 |
567 | return (ref $z)->make(log($r) * log10inv, $t * log10inv); |
66730be0 |
568 | } |
569 | |
570 | # |
571 | # logn |
572 | # |
573 | # Compute logn(z,n) = log(z) / log(n) |
574 | # |
575 | sub logn { |
576 | my ($z, $n) = @_; |
0c721ce2 |
577 | $z = cplx($z, 0) unless ref $z; |
66730be0 |
578 | my $logn = $logn{$n}; |
579 | $logn = $logn{$n} = log($n) unless defined $logn; # Cache log(n) |
0c721ce2 |
580 | return log($z) / $logn; |
66730be0 |
581 | } |
582 | |
583 | # |
584 | # (cos) |
585 | # |
586 | # Compute cos(z) = (exp(iz) + exp(-iz))/2. |
587 | # |
588 | sub cos { |
589 | my ($z) = @_; |
5cd24f17 |
590 | $z = cplx($z, 0) unless ref $z; |
66730be0 |
591 | my ($x, $y) = @{$z->cartesian}; |
592 | my $ey = exp($y); |
593 | my $ey_1 = 1 / $ey; |
0c721ce2 |
594 | return (ref $z)->make(cos($x) * ($ey + $ey_1)/2, |
595 | sin($x) * ($ey_1 - $ey)/2); |
66730be0 |
596 | } |
597 | |
598 | # |
599 | # (sin) |
600 | # |
601 | # Compute sin(z) = (exp(iz) - exp(-iz))/2. |
602 | # |
603 | sub sin { |
604 | my ($z) = @_; |
5cd24f17 |
605 | $z = cplx($z, 0) unless ref $z; |
66730be0 |
606 | my ($x, $y) = @{$z->cartesian}; |
607 | my $ey = exp($y); |
608 | my $ey_1 = 1 / $ey; |
0c721ce2 |
609 | return (ref $z)->make(sin($x) * ($ey + $ey_1)/2, |
610 | cos($x) * ($ey - $ey_1)/2); |
66730be0 |
611 | } |
612 | |
613 | # |
614 | # tan |
615 | # |
616 | # Compute tan(z) = sin(z) / cos(z). |
617 | # |
618 | sub tan { |
619 | my ($z) = @_; |
0c721ce2 |
620 | my $cz = cos($z); |
0e505df1 |
621 | _divbyzero "tan($z)", "cos($z)" if ($cz == 0); |
0c721ce2 |
622 | return sin($z) / $cz; |
66730be0 |
623 | } |
624 | |
625 | # |
0c721ce2 |
626 | # sec |
627 | # |
628 | # Computes the secant sec(z) = 1 / cos(z). |
629 | # |
630 | sub sec { |
631 | my ($z) = @_; |
632 | my $cz = cos($z); |
0e505df1 |
633 | _divbyzero "sec($z)", "cos($z)" if ($cz == 0); |
0c721ce2 |
634 | return 1 / $cz; |
635 | } |
636 | |
637 | # |
638 | # csc |
639 | # |
640 | # Computes the cosecant csc(z) = 1 / sin(z). |
641 | # |
642 | sub csc { |
643 | my ($z) = @_; |
644 | my $sz = sin($z); |
0e505df1 |
645 | _divbyzero "csc($z)", "sin($z)" if ($sz == 0); |
0c721ce2 |
646 | return 1 / $sz; |
647 | } |
648 | |
66730be0 |
649 | # |
0c721ce2 |
650 | # cosec |
66730be0 |
651 | # |
0c721ce2 |
652 | # Alias for csc(). |
653 | # |
654 | sub cosec { Math::Complex::csc(@_) } |
655 | |
656 | # |
657 | # cot |
658 | # |
659 | # Computes cot(z) = 1 / tan(z). |
660 | # |
661 | sub cot { |
66730be0 |
662 | my ($z) = @_; |
0c721ce2 |
663 | my $sz = sin($z); |
0e505df1 |
664 | _divbyzero "cot($z)", "sin($z)" if ($sz == 0); |
0c721ce2 |
665 | return cos($z) / $sz; |
66730be0 |
666 | } |
667 | |
668 | # |
0c721ce2 |
669 | # cotan |
670 | # |
671 | # Alias for cot(). |
672 | # |
673 | sub cotan { Math::Complex::cot(@_) } |
674 | |
675 | # |
66730be0 |
676 | # acos |
677 | # |
678 | # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)). |
679 | # |
680 | sub acos { |
681 | my ($z) = @_; |
0c721ce2 |
682 | $z = cplx($z, 0) unless ref $z; |
8c03c583 |
683 | my ($re, $im) = @{$z->cartesian}; |
684 | return atan2(sqrt(1 - $re * $re), $re) |
685 | if ($im == 0 and abs($re) <= 1.0); |
686 | my $acos = ~i * log($z + sqrt($z*$z - 1)); |
687 | if ($im == 0 || |
688 | (abs($re) < 1 && abs($im) < 1) || |
689 | (abs($re) > 1 && abs($im) > 1 |
690 | && !($re > 1 && $im > 1) |
691 | && !($re < -1 && $im < -1))) { |
692 | # this rule really, REALLY, must be simpler |
693 | return -$acos; |
694 | } |
695 | return $acos; |
66730be0 |
696 | } |
697 | |
698 | # |
699 | # asin |
700 | # |
701 | # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)). |
702 | # |
703 | sub asin { |
704 | my ($z) = @_; |
0c721ce2 |
705 | $z = cplx($z, 0) unless ref $z; |
8c03c583 |
706 | my ($re, $im) = @{$z->cartesian}; |
707 | return atan2($re, sqrt(1 - $re * $re)) |
708 | if ($im == 0 and abs($re) <= 1.0); |
0c721ce2 |
709 | return ~i * log(i * $z + sqrt(1 - $z*$z)); |
66730be0 |
710 | } |
711 | |
712 | # |
713 | # atan |
714 | # |
0c721ce2 |
715 | # Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)). |
66730be0 |
716 | # |
717 | sub atan { |
718 | my ($z) = @_; |
5cd24f17 |
719 | $z = cplx($z, 0) unless ref $z; |
8c03c583 |
720 | _divbyzero "atan(i)" if ( $z == i); |
721 | _divbyzero "atan(-i)" if (-$z == i); |
0c721ce2 |
722 | return i/2*log((i + $z) / (i - $z)); |
a0d0e21e |
723 | } |
724 | |
66730be0 |
725 | # |
0c721ce2 |
726 | # asec |
727 | # |
728 | # Computes the arc secant asec(z) = acos(1 / z). |
729 | # |
730 | sub asec { |
731 | my ($z) = @_; |
0e505df1 |
732 | _divbyzero "asec($z)", $z if ($z == 0); |
8c03c583 |
733 | $z = cplx($z, 0) unless ref $z; |
734 | my ($re, $im) = @{$z->cartesian}; |
735 | if ($im == 0 && abs($re) >= 1.0) { |
736 | my $ire = 1 / $re; |
737 | return atan2(sqrt(1 - $ire * $ire), $ire); |
738 | } |
739 | my $asec = acos(1 / $z); |
740 | return ~$asec if $re < 0 && $re > -1 && $im == 0; |
741 | return -$asec if $im && !($re > 0 && $im > 0) && !($re < 0 && $im < 0); |
742 | return $asec; |
0c721ce2 |
743 | } |
744 | |
745 | # |
5cd24f17 |
746 | # acsc |
0c721ce2 |
747 | # |
8c03c583 |
748 | # Computes the arc cosecant acsc(z) = asin(1 / z). |
0c721ce2 |
749 | # |
5cd24f17 |
750 | sub acsc { |
0c721ce2 |
751 | my ($z) = @_; |
0e505df1 |
752 | _divbyzero "acsc($z)", $z if ($z == 0); |
8c03c583 |
753 | $z = cplx($z, 0) unless ref $z; |
754 | my ($re, $im) = @{$z->cartesian}; |
755 | if ($im == 0 && abs($re) >= 1.0) { |
756 | my $ire = 1 / $re; |
757 | return atan2($ire, sqrt(1 - $ire * $ire)); |
758 | } |
759 | my $acsc = asin(1 / $z); |
760 | return ~$acsc if $re < 0 && $re > -1 && $im == 0; |
761 | return $acsc; |
0c721ce2 |
762 | } |
763 | |
764 | # |
5cd24f17 |
765 | # acosec |
66730be0 |
766 | # |
5cd24f17 |
767 | # Alias for acsc(). |
0c721ce2 |
768 | # |
5cd24f17 |
769 | sub acosec { Math::Complex::acsc(@_) } |
0c721ce2 |
770 | |
66730be0 |
771 | # |
0c721ce2 |
772 | # acot |
773 | # |
8c03c583 |
774 | # Computes the arc cotangent acot(z) = atan(1 / z) |
0c721ce2 |
775 | # |
776 | sub acot { |
66730be0 |
777 | my ($z) = @_; |
8c03c583 |
778 | _divbyzero "acot($z)" if ($z == 0); |
5cd24f17 |
779 | $z = cplx($z, 0) unless ref $z; |
8c03c583 |
780 | _divbyzero "acot(i)", if ( $z == i); |
781 | _divbyzero "acot(-i)" if (-$z == i); |
782 | return atan(1 / $z); |
66730be0 |
783 | } |
784 | |
785 | # |
0c721ce2 |
786 | # acotan |
787 | # |
788 | # Alias for acot(). |
789 | # |
790 | sub acotan { Math::Complex::acot(@_) } |
791 | |
792 | # |
66730be0 |
793 | # cosh |
794 | # |
795 | # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2. |
796 | # |
797 | sub cosh { |
798 | my ($z) = @_; |
0e505df1 |
799 | my $real; |
800 | unless (ref $z) { |
801 | $z = cplx($z, 0); |
802 | $real = 1; |
803 | } |
804 | my ($x, $y) = @{$z->cartesian}; |
66730be0 |
805 | my $ex = exp($x); |
806 | my $ex_1 = 1 / $ex; |
f4837644 |
807 | return cplx(0.5 * ($ex + $ex_1), 0) if $real; |
0c721ce2 |
808 | return (ref $z)->make(cos($y) * ($ex + $ex_1)/2, |
809 | sin($y) * ($ex - $ex_1)/2); |
66730be0 |
810 | } |
811 | |
812 | # |
813 | # sinh |
814 | # |
815 | # Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2. |
816 | # |
817 | sub sinh { |
818 | my ($z) = @_; |
0e505df1 |
819 | my $real; |
820 | unless (ref $z) { |
821 | $z = cplx($z, 0); |
822 | $real = 1; |
823 | } |
824 | my ($x, $y) = @{$z->cartesian}; |
66730be0 |
825 | my $ex = exp($x); |
826 | my $ex_1 = 1 / $ex; |
f4837644 |
827 | return cplx(0.5 * ($ex - $ex_1), 0) if $real; |
0c721ce2 |
828 | return (ref $z)->make(cos($y) * ($ex - $ex_1)/2, |
829 | sin($y) * ($ex + $ex_1)/2); |
66730be0 |
830 | } |
831 | |
832 | # |
833 | # tanh |
834 | # |
835 | # Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z). |
836 | # |
837 | sub tanh { |
838 | my ($z) = @_; |
0c721ce2 |
839 | my $cz = cosh($z); |
0e505df1 |
840 | _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0); |
0c721ce2 |
841 | return sinh($z) / $cz; |
66730be0 |
842 | } |
843 | |
844 | # |
0c721ce2 |
845 | # sech |
846 | # |
847 | # Computes the hyperbolic secant sech(z) = 1 / cosh(z). |
848 | # |
849 | sub sech { |
850 | my ($z) = @_; |
851 | my $cz = cosh($z); |
0e505df1 |
852 | _divbyzero "sech($z)", "cosh($z)" if ($cz == 0); |
0c721ce2 |
853 | return 1 / $cz; |
854 | } |
855 | |
856 | # |
857 | # csch |
858 | # |
859 | # Computes the hyperbolic cosecant csch(z) = 1 / sinh(z). |
66730be0 |
860 | # |
0c721ce2 |
861 | sub csch { |
862 | my ($z) = @_; |
863 | my $sz = sinh($z); |
0e505df1 |
864 | _divbyzero "csch($z)", "sinh($z)" if ($sz == 0); |
0c721ce2 |
865 | return 1 / $sz; |
866 | } |
867 | |
868 | # |
869 | # cosech |
870 | # |
871 | # Alias for csch(). |
872 | # |
873 | sub cosech { Math::Complex::csch(@_) } |
874 | |
66730be0 |
875 | # |
0c721ce2 |
876 | # coth |
877 | # |
878 | # Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z). |
879 | # |
880 | sub coth { |
66730be0 |
881 | my ($z) = @_; |
0c721ce2 |
882 | my $sz = sinh($z); |
0e505df1 |
883 | _divbyzero "coth($z)", "sinh($z)" if ($sz == 0); |
0c721ce2 |
884 | return cosh($z) / $sz; |
66730be0 |
885 | } |
886 | |
887 | # |
0c721ce2 |
888 | # cotanh |
889 | # |
890 | # Alias for coth(). |
891 | # |
892 | sub cotanh { Math::Complex::coth(@_) } |
893 | |
894 | # |
66730be0 |
895 | # acosh |
896 | # |
8c03c583 |
897 | # Computes the arc hyperbolic cosine acosh(z) = log(z +- sqrt(z*z-1)). |
66730be0 |
898 | # |
899 | sub acosh { |
900 | my ($z) = @_; |
5cd24f17 |
901 | $z = cplx($z, 0) unless ref $z; |
8c03c583 |
902 | my ($re, $im) = @{$z->cartesian}; |
903 | return log($re + sqrt(cplx($re*$re - 1, 0))) |
904 | if ($im == 0 && $re < 0); |
0c721ce2 |
905 | return log($z + sqrt($z*$z - 1)); |
66730be0 |
906 | } |
907 | |
908 | # |
909 | # asinh |
910 | # |
911 | # Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1)) |
912 | # |
913 | sub asinh { |
914 | my ($z) = @_; |
5cd24f17 |
915 | $z = cplx($z, 0) unless ref $z; |
0c721ce2 |
916 | return log($z + sqrt($z*$z + 1)); |
66730be0 |
917 | } |
918 | |
919 | # |
920 | # atanh |
921 | # |
922 | # Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)). |
923 | # |
924 | sub atanh { |
925 | my ($z) = @_; |
8c03c583 |
926 | _divbyzero 'atanh(1)', "1 - $z" if ($z == 1); |
927 | _logofzero 'atanh(-1)' if ($z == -1); |
5cd24f17 |
928 | $z = cplx($z, 0) unless ref $z; |
8c03c583 |
929 | my ($re, $im) = @{$z->cartesian}; |
930 | if ($im == 0 && $re > 1) { |
931 | return cplx(atanh(1 / $re), pi/2); |
932 | } |
933 | return log((1 + $z) / (1 - $z)) / 2; |
66730be0 |
934 | } |
935 | |
936 | # |
0c721ce2 |
937 | # asech |
938 | # |
939 | # Computes the hyperbolic arc secant asech(z) = acosh(1 / z). |
940 | # |
941 | sub asech { |
942 | my ($z) = @_; |
0e505df1 |
943 | _divbyzero 'asech(0)', $z if ($z == 0); |
8c03c583 |
944 | $z = cplx($z, 0) unless ref $z; |
945 | my ($re, $im) = @{$z->cartesian}; |
946 | if ($im == 0 && $re < 0) { |
947 | my $ire = 1 / $re; |
948 | return log($ire + sqrt(cplx($ire*$ire - 1, 0))); |
949 | } |
0c721ce2 |
950 | return acosh(1 / $z); |
951 | } |
952 | |
953 | # |
954 | # acsch |
66730be0 |
955 | # |
0c721ce2 |
956 | # Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z). |
66730be0 |
957 | # |
0c721ce2 |
958 | sub acsch { |
66730be0 |
959 | my ($z) = @_; |
0e505df1 |
960 | _divbyzero 'acsch(0)', $z if ($z == 0); |
0c721ce2 |
961 | return asinh(1 / $z); |
962 | } |
963 | |
964 | # |
965 | # acosech |
966 | # |
967 | # Alias for acosh(). |
968 | # |
969 | sub acosech { Math::Complex::acsch(@_) } |
970 | |
971 | # |
972 | # acoth |
973 | # |
974 | # Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)). |
975 | # |
976 | sub acoth { |
977 | my ($z) = @_; |
8c03c583 |
978 | _divbyzero 'acoth(1)', "$z - 1" if ($z == 1); |
979 | _logofzero 'acoth(-1)' if ($z == -1); |
5cd24f17 |
980 | $z = cplx($z, 0) unless ref $z; |
8c03c583 |
981 | my ($re, $im) = @{$z->cartesian}; |
982 | if ($im == 0 and abs($re) < 1) { |
983 | return cplx(acoth(1/$re) , pi/2); |
984 | } |
985 | return log((1 + $z) / ($z - 1)) / 2; |
66730be0 |
986 | } |
987 | |
988 | # |
0c721ce2 |
989 | # acotanh |
990 | # |
991 | # Alias for acot(). |
992 | # |
993 | sub acotanh { Math::Complex::acoth(@_) } |
994 | |
995 | # |
66730be0 |
996 | # (atan2) |
997 | # |
998 | # Compute atan(z1/z2). |
999 | # |
1000 | sub atan2 { |
1001 | my ($z1, $z2, $inverted) = @_; |
5cd24f17 |
1002 | my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0); |
1003 | my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); |
66730be0 |
1004 | my $tan; |
1005 | if (defined $inverted && $inverted) { # atan(z2/z1) |
1006 | return pi * ($re2 > 0 ? 1 : -1) if $re1 == 0 && $im1 == 0; |
1007 | $tan = $z2 / $z1; |
1008 | } else { |
1009 | return pi * ($re1 > 0 ? 1 : -1) if $re2 == 0 && $im2 == 0; |
1010 | $tan = $z1 / $z2; |
1011 | } |
1012 | return atan($tan); |
1013 | } |
1014 | |
1015 | # |
1016 | # display_format |
1017 | # ->display_format |
1018 | # |
1019 | # Set (fetch if no argument) display format for all complex numbers that |
1020 | # don't happen to have overrriden it via ->display_format |
1021 | # |
1022 | # When called as a method, this actually sets the display format for |
1023 | # the current object. |
1024 | # |
1025 | # Valid object formats are 'c' and 'p' for cartesian and polar. The first |
1026 | # letter is used actually, so the type can be fully spelled out for clarity. |
1027 | # |
1028 | sub display_format { |
1029 | my $self = shift; |
1030 | my $format = undef; |
1031 | |
1032 | if (ref $self) { # Called as a method |
1033 | $format = shift; |
0c721ce2 |
1034 | } else { # Regular procedure call |
66730be0 |
1035 | $format = $self; |
1036 | undef $self; |
1037 | } |
1038 | |
1039 | if (defined $self) { |
1040 | return defined $self->{display} ? $self->{display} : $display |
1041 | unless defined $format; |
1042 | return $self->{display} = $format; |
1043 | } |
1044 | |
1045 | return $display unless defined $format; |
1046 | return $display = $format; |
1047 | } |
1048 | |
1049 | # |
1050 | # (stringify) |
1051 | # |
1052 | # Show nicely formatted complex number under its cartesian or polar form, |
1053 | # depending on the current display format: |
1054 | # |
1055 | # . If a specific display format has been recorded for this object, use it. |
1056 | # . Otherwise, use the generic current default for all complex numbers, |
1057 | # which is a package global variable. |
1058 | # |
a0d0e21e |
1059 | sub stringify { |
66730be0 |
1060 | my ($z) = shift; |
1061 | my $format; |
1062 | |
1063 | $format = $display; |
1064 | $format = $z->{display} if defined $z->{display}; |
1065 | |
1066 | return $z->stringify_polar if $format =~ /^p/i; |
1067 | return $z->stringify_cartesian; |
1068 | } |
1069 | |
1070 | # |
1071 | # ->stringify_cartesian |
1072 | # |
1073 | # Stringify as a cartesian representation 'a+bi'. |
1074 | # |
1075 | sub stringify_cartesian { |
1076 | my $z = shift; |
1077 | my ($x, $y) = @{$z->cartesian}; |
1078 | my ($re, $im); |
1079 | |
55497cff |
1080 | $x = int($x + ($x < 0 ? -1 : 1) * 1e-14) |
1081 | if int(abs($x)) != int(abs($x) + 1e-14); |
1082 | $y = int($y + ($y < 0 ? -1 : 1) * 1e-14) |
1083 | if int(abs($y)) != int(abs($y) + 1e-14); |
1084 | |
66730be0 |
1085 | $re = "$x" if abs($x) >= 1e-14; |
1086 | if ($y == 1) { $im = 'i' } |
1087 | elsif ($y == -1) { $im = '-i' } |
40da2db3 |
1088 | elsif (abs($y) >= 1e-14) { $im = $y . "i" } |
66730be0 |
1089 | |
0c721ce2 |
1090 | my $str = ''; |
66730be0 |
1091 | $str = $re if defined $re; |
1092 | $str .= "+$im" if defined $im; |
1093 | $str =~ s/\+-/-/; |
1094 | $str =~ s/^\+//; |
1095 | $str = '0' unless $str; |
1096 | |
1097 | return $str; |
1098 | } |
1099 | |
1100 | # |
1101 | # ->stringify_polar |
1102 | # |
1103 | # Stringify as a polar representation '[r,t]'. |
1104 | # |
1105 | sub stringify_polar { |
1106 | my $z = shift; |
1107 | my ($r, $t) = @{$z->polar}; |
1108 | my $theta; |
0c721ce2 |
1109 | my $eps = 1e-14; |
66730be0 |
1110 | |
0c721ce2 |
1111 | return '[0,0]' if $r <= $eps; |
a0d0e21e |
1112 | |
66730be0 |
1113 | my $tpi = 2 * pi; |
1114 | my $nt = $t / $tpi; |
1115 | $nt = ($nt - int($nt)) * $tpi; |
1116 | $nt += $tpi if $nt < 0; # Range [0, 2pi] |
a0d0e21e |
1117 | |
0c721ce2 |
1118 | if (abs($nt) <= $eps) { $theta = 0 } |
1119 | elsif (abs(pi-$nt) <= $eps) { $theta = 'pi' } |
66730be0 |
1120 | |
55497cff |
1121 | if (defined $theta) { |
0c721ce2 |
1122 | $r = int($r + ($r < 0 ? -1 : 1) * $eps) |
1123 | if int(abs($r)) != int(abs($r) + $eps); |
1124 | $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps) |
1125 | if ($theta ne 'pi' and |
1126 | int(abs($theta)) != int(abs($theta) + $eps)); |
55497cff |
1127 | return "\[$r,$theta\]"; |
1128 | } |
66730be0 |
1129 | |
1130 | # |
1131 | # Okay, number is not a real. Try to identify pi/n and friends... |
1132 | # |
1133 | |
1134 | $nt -= $tpi if $nt > pi; |
1135 | my ($n, $k, $kpi); |
1136 | |
1137 | for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) { |
1138 | $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5); |
0c721ce2 |
1139 | if (abs($kpi/$n - $nt) <= $eps) { |
1140 | $theta = ($nt < 0 ? '-':''). |
1141 | ($k == 1 ? 'pi':"${k}pi").'/'.abs($n); |
66730be0 |
1142 | last; |
1143 | } |
1144 | } |
1145 | |
1146 | $theta = $nt unless defined $theta; |
1147 | |
0c721ce2 |
1148 | $r = int($r + ($r < 0 ? -1 : 1) * $eps) |
1149 | if int(abs($r)) != int(abs($r) + $eps); |
1150 | $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps) |
1151 | if ($theta !~ m(^-?\d*pi/\d+$) and |
1152 | int(abs($theta)) != int(abs($theta) + $eps)); |
55497cff |
1153 | |
66730be0 |
1154 | return "\[$r,$theta\]"; |
a0d0e21e |
1155 | } |
a5f75d66 |
1156 | |
1157 | 1; |
1158 | __END__ |
1159 | |
1160 | =head1 NAME |
1161 | |
66730be0 |
1162 | Math::Complex - complex numbers and associated mathematical functions |
a5f75d66 |
1163 | |
1164 | =head1 SYNOPSIS |
1165 | |
66730be0 |
1166 | use Math::Complex; |
5aabfad6 |
1167 | |
66730be0 |
1168 | $z = Math::Complex->make(5, 6); |
1169 | $t = 4 - 3*i + $z; |
1170 | $j = cplxe(1, 2*pi/3); |
a5f75d66 |
1171 | |
1172 | =head1 DESCRIPTION |
1173 | |
66730be0 |
1174 | This package lets you create and manipulate complex numbers. By default, |
1175 | I<Perl> limits itself to real numbers, but an extra C<use> statement brings |
1176 | full complex support, along with a full set of mathematical functions |
1177 | typically associated with and/or extended to complex numbers. |
1178 | |
1179 | If you wonder what complex numbers are, they were invented to be able to solve |
1180 | the following equation: |
1181 | |
1182 | x*x = -1 |
1183 | |
1184 | and by definition, the solution is noted I<i> (engineers use I<j> instead since |
1185 | I<i> usually denotes an intensity, but the name does not matter). The number |
1186 | I<i> is a pure I<imaginary> number. |
1187 | |
1188 | The arithmetics with pure imaginary numbers works just like you would expect |
1189 | it with real numbers... you just have to remember that |
1190 | |
1191 | i*i = -1 |
1192 | |
1193 | so you have: |
1194 | |
1195 | 5i + 7i = i * (5 + 7) = 12i |
1196 | 4i - 3i = i * (4 - 3) = i |
1197 | 4i * 2i = -8 |
1198 | 6i / 2i = 3 |
1199 | 1 / i = -i |
1200 | |
1201 | Complex numbers are numbers that have both a real part and an imaginary |
1202 | part, and are usually noted: |
1203 | |
1204 | a + bi |
1205 | |
1206 | where C<a> is the I<real> part and C<b> is the I<imaginary> part. The |
1207 | arithmetic with complex numbers is straightforward. You have to |
1208 | keep track of the real and the imaginary parts, but otherwise the |
1209 | rules used for real numbers just apply: |
1210 | |
1211 | (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i |
1212 | (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i |
1213 | |
1214 | A graphical representation of complex numbers is possible in a plane |
1215 | (also called the I<complex plane>, but it's really a 2D plane). |
1216 | The number |
1217 | |
1218 | z = a + bi |
1219 | |
1220 | is the point whose coordinates are (a, b). Actually, it would |
1221 | be the vector originating from (0, 0) to (a, b). It follows that the addition |
1222 | of two complex numbers is a vectorial addition. |
1223 | |
1224 | Since there is a bijection between a point in the 2D plane and a complex |
1225 | number (i.e. the mapping is unique and reciprocal), a complex number |
1226 | can also be uniquely identified with polar coordinates: |
1227 | |
1228 | [rho, theta] |
1229 | |
1230 | where C<rho> is the distance to the origin, and C<theta> the angle between |
1231 | the vector and the I<x> axis. There is a notation for this using the |
1232 | exponential form, which is: |
1233 | |
1234 | rho * exp(i * theta) |
1235 | |
1236 | where I<i> is the famous imaginary number introduced above. Conversion |
1237 | between this form and the cartesian form C<a + bi> is immediate: |
1238 | |
1239 | a = rho * cos(theta) |
1240 | b = rho * sin(theta) |
1241 | |
1242 | which is also expressed by this formula: |
1243 | |
1244 | z = rho * exp(i * theta) = rho * (cos theta + i * sin theta) |
1245 | |
1246 | In other words, it's the projection of the vector onto the I<x> and I<y> |
1247 | axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta> |
1248 | the I<argument> of the complex number. The I<norm> of C<z> will be |
1249 | noted C<abs(z)>. |
1250 | |
1251 | The polar notation (also known as the trigonometric |
1252 | representation) is much more handy for performing multiplications and |
1253 | divisions of complex numbers, whilst the cartesian notation is better |
1254 | suited for additions and substractions. Real numbers are on the I<x> |
1255 | axis, and therefore I<theta> is zero. |
1256 | |
1257 | All the common operations that can be performed on a real number have |
1258 | been defined to work on complex numbers as well, and are merely |
1259 | I<extensions> of the operations defined on real numbers. This means |
1260 | they keep their natural meaning when there is no imaginary part, provided |
1261 | the number is within their definition set. |
1262 | |
1263 | For instance, the C<sqrt> routine which computes the square root of |
1264 | its argument is only defined for positive real numbers and yields a |
1265 | positive real number (it is an application from B<R+> to B<R+>). |
1266 | If we allow it to return a complex number, then it can be extended to |
1267 | negative real numbers to become an application from B<R> to B<C> (the |
1268 | set of complex numbers): |
1269 | |
1270 | sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i |
1271 | |
1272 | It can also be extended to be an application from B<C> to B<C>, |
1273 | whilst its restriction to B<R> behaves as defined above by using |
1274 | the following definition: |
1275 | |
1276 | sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2) |
1277 | |
1278 | Indeed, a negative real number can be noted C<[x,pi]> |
1279 | (the modulus I<x> is always positive, so C<[x,pi]> is really C<-x>, a |
1280 | negative number) |
1281 | and the above definition states that |
1282 | |
1283 | sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i |
1284 | |
1285 | which is exactly what we had defined for negative real numbers above. |
a5f75d66 |
1286 | |
66730be0 |
1287 | All the common mathematical functions defined on real numbers that |
1288 | are extended to complex numbers share that same property of working |
1289 | I<as usual> when the imaginary part is zero (otherwise, it would not |
1290 | be called an extension, would it?). |
a5f75d66 |
1291 | |
66730be0 |
1292 | A I<new> operation possible on a complex number that is |
1293 | the identity for real numbers is called the I<conjugate>, and is noted |
1294 | with an horizontal bar above the number, or C<~z> here. |
a5f75d66 |
1295 | |
66730be0 |
1296 | z = a + bi |
1297 | ~z = a - bi |
a5f75d66 |
1298 | |
66730be0 |
1299 | Simple... Now look: |
a5f75d66 |
1300 | |
66730be0 |
1301 | z * ~z = (a + bi) * (a - bi) = a*a + b*b |
a5f75d66 |
1302 | |
66730be0 |
1303 | We saw that the norm of C<z> was noted C<abs(z)> and was defined as the |
1304 | distance to the origin, also known as: |
a5f75d66 |
1305 | |
66730be0 |
1306 | rho = abs(z) = sqrt(a*a + b*b) |
a5f75d66 |
1307 | |
66730be0 |
1308 | so |
1309 | |
1310 | z * ~z = abs(z) ** 2 |
1311 | |
1312 | If z is a pure real number (i.e. C<b == 0>), then the above yields: |
1313 | |
1314 | a * a = abs(a) ** 2 |
1315 | |
1316 | which is true (C<abs> has the regular meaning for real number, i.e. stands |
1317 | for the absolute value). This example explains why the norm of C<z> is |
1318 | noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet |
1319 | is the regular C<abs> we know when the complex number actually has no |
1320 | imaginary part... This justifies I<a posteriori> our use of the C<abs> |
1321 | notation for the norm. |
1322 | |
1323 | =head1 OPERATIONS |
1324 | |
1325 | Given the following notations: |
1326 | |
1327 | z1 = a + bi = r1 * exp(i * t1) |
1328 | z2 = c + di = r2 * exp(i * t2) |
1329 | z = <any complex or real number> |
1330 | |
1331 | the following (overloaded) operations are supported on complex numbers: |
1332 | |
1333 | z1 + z2 = (a + c) + i(b + d) |
1334 | z1 - z2 = (a - c) + i(b - d) |
1335 | z1 * z2 = (r1 * r2) * exp(i * (t1 + t2)) |
1336 | z1 / z2 = (r1 / r2) * exp(i * (t1 - t2)) |
1337 | z1 ** z2 = exp(z2 * log z1) |
1338 | ~z1 = a - bi |
1339 | abs(z1) = r1 = sqrt(a*a + b*b) |
1340 | sqrt(z1) = sqrt(r1) * exp(i * t1/2) |
1341 | exp(z1) = exp(a) * exp(i * b) |
1342 | log(z1) = log(r1) + i*t1 |
1343 | sin(z1) = 1/2i (exp(i * z1) - exp(-i * z1)) |
1344 | cos(z1) = 1/2 (exp(i * z1) + exp(-i * z1)) |
1345 | abs(z1) = r1 |
1346 | atan2(z1, z2) = atan(z1/z2) |
1347 | |
1348 | The following extra operations are supported on both real and complex |
1349 | numbers: |
1350 | |
1351 | Re(z) = a |
1352 | Im(z) = b |
1353 | arg(z) = t |
1354 | |
1355 | cbrt(z) = z ** (1/3) |
1356 | log10(z) = log(z) / log(10) |
1357 | logn(z, n) = log(z) / log(n) |
1358 | |
1359 | tan(z) = sin(z) / cos(z) |
0c721ce2 |
1360 | |
5aabfad6 |
1361 | csc(z) = 1 / sin(z) |
1362 | sec(z) = 1 / cos(z) |
0c721ce2 |
1363 | cot(z) = 1 / tan(z) |
66730be0 |
1364 | |
1365 | asin(z) = -i * log(i*z + sqrt(1-z*z)) |
1366 | acos(z) = -i * log(z + sqrt(z*z-1)) |
1367 | atan(z) = i/2 * log((i+z) / (i-z)) |
0c721ce2 |
1368 | |
5aabfad6 |
1369 | acsc(z) = asin(1 / z) |
1370 | asec(z) = acos(1 / z) |
8c03c583 |
1371 | acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i)) |
66730be0 |
1372 | |
1373 | sinh(z) = 1/2 (exp(z) - exp(-z)) |
1374 | cosh(z) = 1/2 (exp(z) + exp(-z)) |
0c721ce2 |
1375 | tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z)) |
1376 | |
5aabfad6 |
1377 | csch(z) = 1 / sinh(z) |
1378 | sech(z) = 1 / cosh(z) |
0c721ce2 |
1379 | coth(z) = 1 / tanh(z) |
66730be0 |
1380 | |
1381 | asinh(z) = log(z + sqrt(z*z+1)) |
1382 | acosh(z) = log(z + sqrt(z*z-1)) |
1383 | atanh(z) = 1/2 * log((1+z) / (1-z)) |
66730be0 |
1384 | |
5aabfad6 |
1385 | acsch(z) = asinh(1 / z) |
1386 | asech(z) = acosh(1 / z) |
0c721ce2 |
1387 | acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1)) |
1388 | |
1389 | I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>, I<coth>, |
1390 | I<acosech>, I<acotanh>, have aliases I<ln>, I<cosec>, I<cotan>, |
1391 | I<acosec>, I<acotan>, I<cosech>, I<cotanh>, I<acosech>, I<acotanh>, |
1392 | respectively. |
1393 | |
1394 | The I<root> function is available to compute all the I<n> |
66730be0 |
1395 | roots of some complex, where I<n> is a strictly positive integer. |
1396 | There are exactly I<n> such roots, returned as a list. Getting the |
1397 | number mathematicians call C<j> such that: |
1398 | |
1399 | 1 + j + j*j = 0; |
1400 | |
1401 | is a simple matter of writing: |
1402 | |
1403 | $j = ((root(1, 3))[1]; |
1404 | |
1405 | The I<k>th root for C<z = [r,t]> is given by: |
1406 | |
1407 | (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n) |
1408 | |
f4837644 |
1409 | The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In |
1410 | order to ensure its restriction to real numbers is conform to what you |
1411 | would expect, the comparison is run on the real part of the complex |
1412 | number first, and imaginary parts are compared only when the real |
1413 | parts match. |
66730be0 |
1414 | |
1415 | =head1 CREATION |
1416 | |
1417 | To create a complex number, use either: |
1418 | |
1419 | $z = Math::Complex->make(3, 4); |
1420 | $z = cplx(3, 4); |
1421 | |
1422 | if you know the cartesian form of the number, or |
1423 | |
1424 | $z = 3 + 4*i; |
1425 | |
1426 | if you like. To create a number using the trigonometric form, use either: |
1427 | |
1428 | $z = Math::Complex->emake(5, pi/3); |
1429 | $x = cplxe(5, pi/3); |
1430 | |
0c721ce2 |
1431 | instead. The first argument is the modulus, the second is the angle |
1432 | (in radians, the full circle is 2*pi). (Mnmemonic: C<e> is used as a |
1433 | notation for complex numbers in the trigonometric form). |
66730be0 |
1434 | |
1435 | It is possible to write: |
1436 | |
1437 | $x = cplxe(-3, pi/4); |
1438 | |
1439 | but that will be silently converted into C<[3,-3pi/4]>, since the modulus |
1440 | must be positive (it represents the distance to the origin in the complex |
1441 | plane). |
1442 | |
1443 | =head1 STRINGIFICATION |
1444 | |
1445 | When printed, a complex number is usually shown under its cartesian |
1446 | form I<a+bi>, but there are legitimate cases where the polar format |
1447 | I<[r,t]> is more appropriate. |
1448 | |
1449 | By calling the routine C<Math::Complex::display_format> and supplying either |
1450 | C<"polar"> or C<"cartesian">, you override the default display format, |
1451 | which is C<"cartesian">. Not supplying any argument returns the current |
1452 | setting. |
1453 | |
1454 | This default can be overridden on a per-number basis by calling the |
1455 | C<display_format> method instead. As before, not supplying any argument |
1456 | returns the current display format for this number. Otherwise whatever you |
1457 | specify will be the new display format for I<this> particular number. |
1458 | |
1459 | For instance: |
1460 | |
1461 | use Math::Complex; |
1462 | |
1463 | Math::Complex::display_format('polar'); |
1464 | $j = ((root(1, 3))[1]; |
1465 | print "j = $j\n"; # Prints "j = [1,2pi/3] |
1466 | $j->display_format('cartesian'); |
1467 | print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i" |
1468 | |
1469 | The polar format attempts to emphasize arguments like I<k*pi/n> |
1470 | (where I<n> is a positive integer and I<k> an integer within [-9,+9]). |
1471 | |
1472 | =head1 USAGE |
1473 | |
1474 | Thanks to overloading, the handling of arithmetics with complex numbers |
1475 | is simple and almost transparent. |
1476 | |
1477 | Here are some examples: |
1478 | |
1479 | use Math::Complex; |
1480 | |
1481 | $j = cplxe(1, 2*pi/3); # $j ** 3 == 1 |
1482 | print "j = $j, j**3 = ", $j ** 3, "\n"; |
1483 | print "1 + j + j**2 = ", 1 + $j + $j**2, "\n"; |
1484 | |
1485 | $z = -16 + 0*i; # Force it to be a complex |
1486 | print "sqrt($z) = ", sqrt($z), "\n"; |
1487 | |
1488 | $k = exp(i * 2*pi/3); |
1489 | print "$j - $k = ", $j - $k, "\n"; |
a5f75d66 |
1490 | |
5cd24f17 |
1491 | =head1 ERRORS DUE TO DIVISION BY ZERO |
5aabfad6 |
1492 | |
1493 | The division (/) and the following functions |
1494 | |
1495 | tan |
1496 | sec |
1497 | csc |
1498 | cot |
5cd24f17 |
1499 | asec |
1500 | acsc |
5aabfad6 |
1501 | atan |
1502 | acot |
1503 | tanh |
1504 | sech |
1505 | csch |
1506 | coth |
1507 | atanh |
1508 | asech |
1509 | acsch |
1510 | acoth |
1511 | |
1512 | cannot be computed for all arguments because that would mean dividing |
8c03c583 |
1513 | by zero or taking logarithm of zero. These situations cause fatal |
1514 | runtime errors looking like this |
5aabfad6 |
1515 | |
1516 | cot(0): Division by zero. |
5cd24f17 |
1517 | (Because in the definition of cot(0), the divisor sin(0) is 0) |
5aabfad6 |
1518 | Died at ... |
1519 | |
8c03c583 |
1520 | or |
1521 | |
1522 | atanh(-1): Logarithm of zero. |
1523 | Died at... |
1524 | |
1525 | For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>, |
1526 | C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the |
1527 | C<atanh>, C<acoth>, the argument cannot be C<1> (one). For the |
1528 | C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one). For the |
1529 | C<atan>, C<acot>, the argument cannot be C<i> (the imaginary unit). |
1530 | For the C<atan>, C<acoth>, the argument cannot be C<-i> (the negative |
1531 | imaginary unit). For the C<tan>, C<sec>, C<tanh>, C<sech>, the |
1532 | argument cannot be I<pi/2 + k * pi>, where I<k> is any integer. |
5cd24f17 |
1533 | |
a5f75d66 |
1534 | =head1 BUGS |
1535 | |
5cd24f17 |
1536 | Saying C<use Math::Complex;> exports many mathematical routines in the |
1537 | caller environment and even overrides some (C<sin>, C<cos>, C<sqrt>, |
1538 | C<log>, C<exp>). This is construed as a feature by the Authors, |
1539 | actually... ;-) |
66730be0 |
1540 | |
1541 | The code is not optimized for speed, although we try to use the cartesian |
1542 | form for addition-like operators and the trigonometric form for all |
1543 | multiplication-like operators. |
1544 | |
1545 | The arg() routine does not ensure the angle is within the range [-pi,+pi] |
1546 | (a side effect caused by multiplication and division using the trigonometric |
1547 | representation). |
a5f75d66 |
1548 | |
66730be0 |
1549 | All routines expect to be given real or complex numbers. Don't attempt to |
1550 | use BigFloat, since Perl has currently no rule to disambiguate a '+' |
1551 | operation (for instance) between two overloaded entities. |
a5f75d66 |
1552 | |
0c721ce2 |
1553 | =head1 AUTHORS |
a5f75d66 |
1554 | |
ace5de91 |
1555 | Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>> and |
1556 | Jarkko Hietaniemi <F<jhi@iki.fi>>. |
5cd24f17 |
1557 | |
1558 | =cut |
1559 | |
1560 | # eof |