source: trunk/essentials/dev-lang/perl/lib/Math/Complex.pm@ 3951

Last change on this file since 3951 was 3181, checked in by bird, 19 years ago

perl 5.8.8
File size: 44.8 KB
Line 
1#
2# Complex numbers and associated mathematical functions
3# -- Raphael Manfredi Since Sep 1996
4# -- Jarkko Hietaniemi Since Mar 1997
5# -- Daniel S. Lewart Since Sep 1997
6#
7
8package Math::Complex;
9
10use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $Inf);
11
12$VERSION = 1.35;
13
14BEGIN {
15 unless ($^O eq 'unicosmk') {
16 my $e = $!;
17 # We do want an arithmetic overflow, Inf INF inf Infinity:.
18 undef $Inf unless eval <<'EOE' and $Inf =~ /^inf(?:inity)?$/i;
19 local $SIG{FPE} = sub {die};
20 my $t = CORE::exp 30;
21 $Inf = CORE::exp $t;
22EOE
23 if (!defined $Inf) { # Try a different method
24 undef $Inf unless eval <<'EOE' and $Inf =~ /^inf(?:inity)?$/i;
25 local $SIG{FPE} = sub {die};
26 my $t = 1;
27 $Inf = $t + "1e99999999999999999999999999999999";
28EOE
29 }
30 $! = $e; # Clear ERANGE.
31 }
32 $Inf = "Inf" if !defined $Inf || !($Inf > 0); # Desperation.
33}
34
35use strict;
36
37my $i;
38my %LOGN;
39
40# Regular expression for floating point numbers.
41# These days we could use Scalar::Util::lln(), I guess.
42my $gre = qr'\s*([\+\-]?(?:(?:(?:\d+(?:_\d+)*(?:\.\d*(?:_\d+)*)?|\.\d+(?:_\d+)*)(?:[eE][\+\-]?\d+(?:_\d+)*)?))|inf)'i;
43
44require Exporter;
45
46@ISA = qw(Exporter);
47
48my @trig = qw(
49 pi
50 tan
51 csc cosec sec cot cotan
52 asin acos atan
53 acsc acosec asec acot acotan
54 sinh cosh tanh
55 csch cosech sech coth cotanh
56 asinh acosh atanh
57 acsch acosech asech acoth acotanh
58 );
59
60@EXPORT = (qw(
61 i Re Im rho theta arg
62 sqrt log ln
63 log10 logn cbrt root
64 cplx cplxe
65 atan2
66 ),
67 @trig);
68
69@EXPORT_OK = qw(decplx);
70
71%EXPORT_TAGS = (
72 'trig' => [@trig],
73);
74
75use overload
76 '+' => \&plus,
77 '-' => \&minus,
78 '*' => \&multiply,
79 '/' => \&divide,
80 '**' => \&power,
81 '==' => \&numeq,
82 '<=>' => \&spaceship,
83 'neg' => \&negate,
84 '~' => \&conjugate,
85 'abs' => \&abs,
86 'sqrt' => \&sqrt,
87 'exp' => \&exp,
88 'log' => \&log,
89 'sin' => \&sin,
90 'cos' => \&cos,
91 'tan' => \&tan,
92 'atan2' => \&atan2,
93 qw("" stringify);
94
95#
96# Package "privates"
97#
98
99my %DISPLAY_FORMAT = ('style' => 'cartesian',
100 'polar_pretty_print' => 1);
101my $eps = 1e-14; # Epsilon
102
103#
104# Object attributes (internal):
105# cartesian [real, imaginary] -- cartesian form
106# polar [rho, theta] -- polar form
107# c_dirty cartesian form not up-to-date
108# p_dirty polar form not up-to-date
109# display display format (package's global when not set)
110#
111
112# Die on bad *make() arguments.
113
114sub _cannot_make {
115 die "@{[(caller(1))[3]]}: Cannot take $_[0] of '$_[1]'.\n";
116}
117
118sub _make {
119 my $arg = shift;
120 my ($p, $q);
121
122 if ($arg =~ /^$gre$/) {
123 ($p, $q) = ($1, 0);
124 } elsif ($arg =~ /^(?:$gre)?$gre\s*i\s*$/) {
125 ($p, $q) = ($1 || 0, $2);
126 } elsif ($arg =~ /^\s*\(\s*$gre\s*(?:,\s*$gre\s*)?\)\s*$/) {
127 ($p, $q) = ($1, $2 || 0);
128 }
129
130 if (defined $p) {
131 $p =~ s/^\+//;
132 $p =~ s/^(-?)inf$/"${1}9**9**9"/e;
133 $q =~ s/^\+//;
134 $q =~ s/^(-?)inf$/"${1}9**9**9"/e;
135 }
136
137 return ($p, $q);
138}
139
140sub _emake {
141 my $arg = shift;
142 my ($p, $q);
143
144 if ($arg =~ /^\s*\[\s*$gre\s*(?:,\s*$gre\s*)?\]\s*$/) {
145 ($p, $q) = ($1, $2 || 0);
146 } elsif ($arg =~ m!^\s*\[\s*$gre\s*(?:,\s*([-+]?\d*\s*)?pi(?:/\s*(\d+))?\s*)?\]\s*$!) {
147 ($p, $q) = ($1, ($2 eq '-' ? -1 : ($2 || 1)) * pi() / ($3 || 1));
148 } elsif ($arg =~ /^\s*\[\s*$gre\s*\]\s*$/) {
149 ($p, $q) = ($1, 0);
150 } elsif ($arg =~ /^\s*$gre\s*$/) {
151 ($p, $q) = ($1, 0);
152 }
153
154 if (defined $p) {
155 $p =~ s/^\+//;
156 $q =~ s/^\+//;
157 $p =~ s/^(-?)inf$/"${1}9**9**9"/e;
158 $q =~ s/^(-?)inf$/"${1}9**9**9"/e;
159 }
160
161 return ($p, $q);
162}
163
164#
165# ->make
166#
167# Create a new complex number (cartesian form)
168#
169sub make {
170 my $self = bless {}, shift;
171 my ($re, $im);
172 if (@_ == 0) {
173 ($re, $im) = (0, 0);
174 } elsif (@_ == 1) {
175 return (ref $self)->emake($_[0])
176 if ($_[0] =~ /^\s*\[/);
177 ($re, $im) = _make($_[0]);
178 } elsif (@_ == 2) {
179 ($re, $im) = @_;
180 }
181 if (defined $re) {
182 _cannot_make("real part", $re) unless $re =~ /^$gre$/;
183 }
184 $im ||= 0;
185 _cannot_make("imaginary part", $im) unless $im =~ /^$gre$/;
186 $self->set_cartesian([$re, $im ]);
187 $self->display_format('cartesian');
188
189 return $self;
190}
191
192#
193# ->emake
194#
195# Create a new complex number (exponential form)
196#
197sub emake {
198 my $self = bless {}, shift;
199 my ($rho, $theta);
200 if (@_ == 0) {
201 ($rho, $theta) = (0, 0);
202 } elsif (@_ == 1) {
203 return (ref $self)->make($_[0])
204 if ($_[0] =~ /^\s*\(/ || $_[0] =~ /i\s*$/);
205 ($rho, $theta) = _emake($_[0]);
206 } elsif (@_ == 2) {
207 ($rho, $theta) = @_;
208 }
209 if (defined $rho && defined $theta) {
210 if ($rho < 0) {
211 $rho = -$rho;
212 $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
213 }
214 }
215 if (defined $rho) {
216 _cannot_make("rho", $rho) unless $rho =~ /^$gre$/;
217 }
218 $theta ||= 0;
219 _cannot_make("theta", $theta) unless $theta =~ /^$gre$/;
220 $self->set_polar([$rho, $theta]);
221 $self->display_format('polar');
222
223 return $self;
224}
225
226sub new { &make } # For backward compatibility only.
227
228#
229# cplx
230#
231# Creates a complex number from a (re, im) tuple.
232# This avoids the burden of writing Math::Complex->make(re, im).
233#
234sub cplx {
235 return __PACKAGE__->make(@_);
236}
237
238#
239# cplxe
240#
241# Creates a complex number from a (rho, theta) tuple.
242# This avoids the burden of writing Math::Complex->emake(rho, theta).
243#
244sub cplxe {
245 return __PACKAGE__->emake(@_);
246}
247
248#
249# pi
250#
251# The number defined as pi = 180 degrees
252#
253sub pi () { 4 * CORE::atan2(1, 1) }
254
255#
256# pit2
257#
258# The full circle
259#
260sub pit2 () { 2 * pi }
261
262#
263# pip2
264#
265# The quarter circle
266#
267sub pip2 () { pi / 2 }
268
269#
270# deg1
271#
272# One degree in radians, used in stringify_polar.
273#
274
275sub deg1 () { pi / 180 }
276
277#
278# uplog10
279#
280# Used in log10().
281#
282sub uplog10 () { 1 / CORE::log(10) }
283
284#
285# i
286#
287# The number defined as i*i = -1;
288#
289sub i () {
290 return $i if ($i);
291 $i = bless {};
292 $i->{'cartesian'} = [0, 1];
293 $i->{'polar'} = [1, pip2];
294 $i->{c_dirty} = 0;
295 $i->{p_dirty} = 0;
296 return $i;
297}
298
299#
300# ip2
301#
302# Half of i.
303#
304sub ip2 () { i / 2 }
305
306#
307# Attribute access/set routines
308#
309
310sub cartesian {$_[0]->{c_dirty} ?
311 $_[0]->update_cartesian : $_[0]->{'cartesian'}}
312sub polar {$_[0]->{p_dirty} ?
313 $_[0]->update_polar : $_[0]->{'polar'}}
314
315sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{c_dirty} = 0;
316 $_[0]->{'cartesian'} = $_[1] }
317sub set_polar { $_[0]->{c_dirty}++; $_[0]->{p_dirty} = 0;
318 $_[0]->{'polar'} = $_[1] }
319
320#
321# ->update_cartesian
322#
323# Recompute and return the cartesian form, given accurate polar form.
324#
325sub update_cartesian {
326 my $self = shift;
327 my ($r, $t) = @{$self->{'polar'}};
328 $self->{c_dirty} = 0;
329 return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)];
330}
331
332#
333#
334# ->update_polar
335#
336# Recompute and return the polar form, given accurate cartesian form.
337#
338sub update_polar {
339 my $self = shift;
340 my ($x, $y) = @{$self->{'cartesian'}};
341 $self->{p_dirty} = 0;
342 return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
343 return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y),
344 CORE::atan2($y, $x)];
345}
346
347#
348# (plus)
349#
350# Computes z1+z2.
351#
352sub plus {
353 my ($z1, $z2, $regular) = @_;
354 my ($re1, $im1) = @{$z1->cartesian};
355 $z2 = cplx($z2) unless ref $z2;
356 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
357 unless (defined $regular) {
358 $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
359 return $z1;
360 }
361 return (ref $z1)->make($re1 + $re2, $im1 + $im2);
362}
363
364#
365# (minus)
366#
367# Computes z1-z2.
368#
369sub minus {
370 my ($z1, $z2, $inverted) = @_;
371 my ($re1, $im1) = @{$z1->cartesian};
372 $z2 = cplx($z2) unless ref $z2;
373 my ($re2, $im2) = @{$z2->cartesian};
374 unless (defined $inverted) {
375 $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
376 return $z1;
377 }
378 return $inverted ?
379 (ref $z1)->make($re2 - $re1, $im2 - $im1) :
380 (ref $z1)->make($re1 - $re2, $im1 - $im2);
381
382}
383
384#
385# (multiply)
386#
387# Computes z1*z2.
388#
389sub multiply {
390 my ($z1, $z2, $regular) = @_;
391 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
392 # if both polar better use polar to avoid rounding errors
393 my ($r1, $t1) = @{$z1->polar};
394 my ($r2, $t2) = @{$z2->polar};
395 my $t = $t1 + $t2;
396 if ($t > pi()) { $t -= pit2 }
397 elsif ($t <= -pi()) { $t += pit2 }
398 unless (defined $regular) {
399 $z1->set_polar([$r1 * $r2, $t]);
400 return $z1;
401 }
402 return (ref $z1)->emake($r1 * $r2, $t);
403 } else {
404 my ($x1, $y1) = @{$z1->cartesian};
405 if (ref $z2) {
406 my ($x2, $y2) = @{$z2->cartesian};
407 return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2);
408 } else {
409 return (ref $z1)->make($x1*$z2, $y1*$z2);
410 }
411 }
412}
413
414#
415# _divbyzero
416#
417# Die on division by zero.
418#
419sub _divbyzero {
420 my $mess = "$_[0]: Division by zero.\n";
421
422 if (defined $_[1]) {
423 $mess .= "(Because in the definition of $_[0], the divisor ";
424 $mess .= "$_[1] " unless ("$_[1]" eq '0');
425 $mess .= "is 0)\n";
426 }
427
428 my @up = caller(1);
429
430 $mess .= "Died at $up[1] line $up[2].\n";
431
432 die $mess;
433}
434
435#
436# (divide)
437#
438# Computes z1/z2.
439#
440sub divide {
441 my ($z1, $z2, $inverted) = @_;
442 if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
443 # if both polar better use polar to avoid rounding errors
444 my ($r1, $t1) = @{$z1->polar};
445 my ($r2, $t2) = @{$z2->polar};
446 my $t;
447 if ($inverted) {
448 _divbyzero "$z2/0" if ($r1 == 0);
449 $t = $t2 - $t1;
450 if ($t > pi()) { $t -= pit2 }
451 elsif ($t <= -pi()) { $t += pit2 }
452 return (ref $z1)->emake($r2 / $r1, $t);
453 } else {
454 _divbyzero "$z1/0" if ($r2 == 0);
455 $t = $t1 - $t2;
456 if ($t > pi()) { $t -= pit2 }
457 elsif ($t <= -pi()) { $t += pit2 }
458 return (ref $z1)->emake($r1 / $r2, $t);
459 }
460 } else {
461 my ($d, $x2, $y2);
462 if ($inverted) {
463 ($x2, $y2) = @{$z1->cartesian};
464 $d = $x2*$x2 + $y2*$y2;
465 _divbyzero "$z2/0" if $d == 0;
466 return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d);
467 } else {
468 my ($x1, $y1) = @{$z1->cartesian};
469 if (ref $z2) {
470 ($x2, $y2) = @{$z2->cartesian};
471 $d = $x2*$x2 + $y2*$y2;
472 _divbyzero "$z1/0" if $d == 0;
473 my $u = ($x1*$x2 + $y1*$y2)/$d;
474 my $v = ($y1*$x2 - $x1*$y2)/$d;
475 return (ref $z1)->make($u, $v);
476 } else {
477 _divbyzero "$z1/0" if $z2 == 0;
478 return (ref $z1)->make($x1/$z2, $y1/$z2);
479 }
480 }
481 }
482}
483
484#
485# (power)
486#
487# Computes z1**z2 = exp(z2 * log z1)).
488#
489sub power {
490 my ($z1, $z2, $inverted) = @_;
491 if ($inverted) {
492 return 1 if $z1 == 0 || $z2 == 1;
493 return 0 if $z2 == 0 && Re($z1) > 0;
494 } else {
495 return 1 if $z2 == 0 || $z1 == 1;
496 return 0 if $z1 == 0 && Re($z2) > 0;
497 }
498 my $w = $inverted ? &exp($z1 * &log($z2))
499 : &exp($z2 * &log($z1));
500 # If both arguments cartesian, return cartesian, else polar.
501 return $z1->{c_dirty} == 0 &&
502 (not ref $z2 or $z2->{c_dirty} == 0) ?
503 cplx(@{$w->cartesian}) : $w;
504}
505
506#
507# (spaceship)
508#
509# Computes z1 <=> z2.
510# Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.
511#
512sub spaceship {
513 my ($z1, $z2, $inverted) = @_;
514 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
515 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
516 my $sgn = $inverted ? -1 : 1;
517 return $sgn * ($re1 <=> $re2) if $re1 != $re2;
518 return $sgn * ($im1 <=> $im2);
519}
520
521#
522# (numeq)
523#
524# Computes z1 == z2.
525#
526# (Required in addition to spaceship() because of NaNs.)
527sub numeq {
528 my ($z1, $z2, $inverted) = @_;
529 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
530 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
531 return $re1 == $re2 && $im1 == $im2 ? 1 : 0;
532}
533
534#
535# (negate)
536#
537# Computes -z.
538#
539sub negate {
540 my ($z) = @_;
541 if ($z->{c_dirty}) {
542 my ($r, $t) = @{$z->polar};
543 $t = ($t <= 0) ? $t + pi : $t - pi;
544 return (ref $z)->emake($r, $t);
545 }
546 my ($re, $im) = @{$z->cartesian};
547 return (ref $z)->make(-$re, -$im);
548}
549
550#
551# (conjugate)
552#
553# Compute complex's conjugate.
554#
555sub conjugate {
556 my ($z) = @_;
557 if ($z->{c_dirty}) {
558 my ($r, $t) = @{$z->polar};
559 return (ref $z)->emake($r, -$t);
560 }
561 my ($re, $im) = @{$z->cartesian};
562 return (ref $z)->make($re, -$im);
563}
564
565#
566# (abs)
567#
568# Compute or set complex's norm (rho).
569#
570sub abs {
571 my ($z, $rho) = @_;
572 unless (ref $z) {
573 if (@_ == 2) {
574 $_[0] = $_[1];
575 } else {
576 return CORE::abs($z);
577 }
578 }
579 if (defined $rho) {
580 $z->{'polar'} = [ $rho, ${$z->polar}[1] ];
581 $z->{p_dirty} = 0;
582 $z->{c_dirty} = 1;
583 return $rho;
584 } else {
585 return ${$z->polar}[0];
586 }
587}
588
589sub _theta {
590 my $theta = $_[0];
591
592 if ($$theta > pi()) { $$theta -= pit2 }
593 elsif ($$theta <= -pi()) { $$theta += pit2 }
594}
595
596#
597# arg
598#
599# Compute or set complex's argument (theta).
600#
601sub arg {
602 my ($z, $theta) = @_;
603 return $z unless ref $z;
604 if (defined $theta) {
605 _theta(\$theta);
606 $z->{'polar'} = [ ${$z->polar}[0], $theta ];
607 $z->{p_dirty} = 0;
608 $z->{c_dirty} = 1;
609 } else {
610 $theta = ${$z->polar}[1];
611 _theta(\$theta);
612 }
613 return $theta;
614}
615
616#
617# (sqrt)
618#
619# Compute sqrt(z).
620#
621# It is quite tempting to use wantarray here so that in list context
622# sqrt() would return the two solutions. This, however, would
623# break things like
624#
625# print "sqrt(z) = ", sqrt($z), "\n";
626#
627# The two values would be printed side by side without no intervening
628# whitespace, quite confusing.
629# Therefore if you want the two solutions use the root().
630#
631sub sqrt {
632 my ($z) = @_;
633 my ($re, $im) = ref $z ? @{$z->cartesian} : ($z, 0);
634 return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re)
635 if $im == 0;
636 my ($r, $t) = @{$z->polar};
637 return (ref $z)->emake(CORE::sqrt($r), $t/2);
638}
639
640#
641# cbrt
642#
643# Compute cbrt(z) (cubic root).
644#
645# Why are we not returning three values? The same answer as for sqrt().
646#
647sub cbrt {
648 my ($z) = @_;
649 return $z < 0 ?
650 -CORE::exp(CORE::log(-$z)/3) :
651 ($z > 0 ? CORE::exp(CORE::log($z)/3): 0)
652 unless ref $z;
653 my ($r, $t) = @{$z->polar};
654 return 0 if $r == 0;
655 return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);
656}
657
658#
659# _rootbad
660#
661# Die on bad root.
662#
663sub _rootbad {
664 my $mess = "Root '$_[0]' illegal, root rank must be positive integer.\n";
665
666 my @up = caller(1);
667
668 $mess .= "Died at $up[1] line $up[2].\n";
669
670 die $mess;
671}
672
673#
674# root
675#
676# Computes all nth root for z, returning an array whose size is n.
677# `n' must be a positive integer.
678#
679# The roots are given by (for k = 0..n-1):
680#
681# z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
682#
683sub root {
684 my ($z, $n, $k) = @_;
685 _rootbad($n) if ($n < 1 or int($n) != $n);
686 my ($r, $t) = ref $z ?
687 @{$z->polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
688 my $theta_inc = pit2 / $n;
689 my $rho = $r ** (1/$n);
690 my $cartesian = ref $z && $z->{c_dirty} == 0;
691 if (@_ == 2) {
692 my @root;
693 for (my $i = 0, my $theta = $t / $n;
694 $i < $n;
695 $i++, $theta += $theta_inc) {
696 my $w = cplxe($rho, $theta);
697 # Yes, $cartesian is loop invariant.
698 push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
699 }
700 return @root;
701 } elsif (@_ == 3) {
702 my $w = cplxe($rho, $t / $n + $k * $theta_inc);
703 return $cartesian ? cplx(@{$w->cartesian}) : $w;
704 }
705}
706
707#
708# Re
709#
710# Return or set Re(z).
711#
712sub Re {
713 my ($z, $Re) = @_;
714 return $z unless ref $z;
715 if (defined $Re) {
716 $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ];
717 $z->{c_dirty} = 0;
718 $z->{p_dirty} = 1;
719 } else {
720 return ${$z->cartesian}[0];
721 }
722}
723
724#
725# Im
726#
727# Return or set Im(z).
728#
729sub Im {
730 my ($z, $Im) = @_;
731 return 0 unless ref $z;
732 if (defined $Im) {
733 $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ];
734 $z->{c_dirty} = 0;
735 $z->{p_dirty} = 1;
736 } else {
737 return ${$z->cartesian}[1];
738 }
739}
740
741#
742# rho
743#
744# Return or set rho(w).
745#
746sub rho {
747 Math::Complex::abs(@_);
748}
749
750#
751# theta
752#
753# Return or set theta(w).
754#
755sub theta {
756 Math::Complex::arg(@_);
757}
758
759#
760# (exp)
761#
762# Computes exp(z).
763#
764sub exp {
765 my ($z) = @_;
766 my ($x, $y) = @{$z->cartesian};
767 return (ref $z)->emake(CORE::exp($x), $y);
768}
769
770#
771# _logofzero
772#
773# Die on logarithm of zero.
774#
775sub _logofzero {
776 my $mess = "$_[0]: Logarithm of zero.\n";
777
778 if (defined $_[1]) {
779 $mess .= "(Because in the definition of $_[0], the argument ";
780 $mess .= "$_[1] " unless ($_[1] eq '0');
781 $mess .= "is 0)\n";
782 }
783
784 my @up = caller(1);
785
786 $mess .= "Died at $up[1] line $up[2].\n";
787
788 die $mess;
789}
790
791#
792# (log)
793#
794# Compute log(z).
795#
796sub log {
797 my ($z) = @_;
798 unless (ref $z) {
799 _logofzero("log") if $z == 0;
800 return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
801 }
802 my ($r, $t) = @{$z->polar};
803 _logofzero("log") if $r == 0;
804 if ($t > pi()) { $t -= pit2 }
805 elsif ($t <= -pi()) { $t += pit2 }
806 return (ref $z)->make(CORE::log($r), $t);
807}
808
809#
810# ln
811#
812# Alias for log().
813#
814sub ln { Math::Complex::log(@_) }
815
816#
817# log10
818#
819# Compute log10(z).
820#
821
822sub log10 {
823 return Math::Complex::log($_[0]) * uplog10;
824}
825
826#
827# logn
828#
829# Compute logn(z,n) = log(z) / log(n)
830#
831sub logn {
832 my ($z, $n) = @_;
833 $z = cplx($z, 0) unless ref $z;
834 my $logn = $LOGN{$n};
835 $logn = $LOGN{$n} = CORE::log($n) unless defined $logn; # Cache log(n)
836 return &log($z) / $logn;
837}
838
839#
840# (cos)
841#
842# Compute cos(z) = (exp(iz) + exp(-iz))/2.
843#
844sub cos {
845 my ($z) = @_;
846 return CORE::cos($z) unless ref $z;
847 my ($x, $y) = @{$z->cartesian};
848 my $ey = CORE::exp($y);
849 my $sx = CORE::sin($x);
850 my $cx = CORE::cos($x);
851 my $ey_1 = $ey ? 1 / $ey : $Inf;
852 return (ref $z)->make($cx * ($ey + $ey_1)/2,
853 $sx * ($ey_1 - $ey)/2);
854}
855
856#
857# (sin)
858#
859# Compute sin(z) = (exp(iz) - exp(-iz))/2.
860#
861sub sin {
862 my ($z) = @_;
863 return CORE::sin($z) unless ref $z;
864 my ($x, $y) = @{$z->cartesian};
865 my $ey = CORE::exp($y);
866 my $sx = CORE::sin($x);
867 my $cx = CORE::cos($x);
868 my $ey_1 = $ey ? 1 / $ey : $Inf;
869 return (ref $z)->make($sx * ($ey + $ey_1)/2,
870 $cx * ($ey - $ey_1)/2);
871}
872
873#
874# tan
875#
876# Compute tan(z) = sin(z) / cos(z).
877#
878sub tan {
879 my ($z) = @_;
880 my $cz = &cos($z);
881 _divbyzero "tan($z)", "cos($z)" if $cz == 0;
882 return &sin($z) / $cz;
883}
884
885#
886# sec
887#
888# Computes the secant sec(z) = 1 / cos(z).
889#
890sub sec {
891 my ($z) = @_;
892 my $cz = &cos($z);
893 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
894 return 1 / $cz;
895}
896
897#
898# csc
899#
900# Computes the cosecant csc(z) = 1 / sin(z).
901#
902sub csc {
903 my ($z) = @_;
904 my $sz = &sin($z);
905 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
906 return 1 / $sz;
907}
908
909#
910# cosec
911#
912# Alias for csc().
913#
914sub cosec { Math::Complex::csc(@_) }
915
916#
917# cot
918#
919# Computes cot(z) = cos(z) / sin(z).
920#
921sub cot {
922 my ($z) = @_;
923 my $sz = &sin($z);
924 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
925 return &cos($z) / $sz;
926}
927
928#
929# cotan
930#
931# Alias for cot().
932#
933sub cotan { Math::Complex::cot(@_) }
934
935#
936# acos
937#
938# Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
939#
940sub acos {
941 my $z = $_[0];
942 return CORE::atan2(CORE::sqrt(1-$z*$z), $z)
943 if (! ref $z) && CORE::abs($z) <= 1;
944 $z = cplx($z, 0) unless ref $z;
945 my ($x, $y) = @{$z->cartesian};
946 return 0 if $x == 1 && $y == 0;
947 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
948 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
949 my $alpha = ($t1 + $t2)/2;
950 my $beta = ($t1 - $t2)/2;
951 $alpha = 1 if $alpha < 1;
952 if ($beta > 1) { $beta = 1 }
953 elsif ($beta < -1) { $beta = -1 }
954 my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
955 my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
956 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
957 return (ref $z)->make($u, $v);
958}
959
960#
961# asin
962#
963# Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
964#
965sub asin {
966 my $z = $_[0];
967 return CORE::atan2($z, CORE::sqrt(1-$z*$z))
968 if (! ref $z) && CORE::abs($z) <= 1;
969 $z = cplx($z, 0) unless ref $z;
970 my ($x, $y) = @{$z->cartesian};
971 return 0 if $x == 0 && $y == 0;
972 my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
973 my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
974 my $alpha = ($t1 + $t2)/2;
975 my $beta = ($t1 - $t2)/2;
976 $alpha = 1 if $alpha < 1;
977 if ($beta > 1) { $beta = 1 }
978 elsif ($beta < -1) { $beta = -1 }
979 my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
980 my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
981 $v = -$v if $y > 0 || ($y == 0 && $x < -1);
982 return (ref $z)->make($u, $v);
983}
984
985#
986# atan
987#
988# Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
989#
990sub atan {
991 my ($z) = @_;
992 return CORE::atan2($z, 1) unless ref $z;
993 my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
994 return 0 if $x == 0 && $y == 0;
995 _divbyzero "atan(i)" if ( $z == i);
996 _logofzero "atan(-i)" if (-$z == i); # -i is a bad file test...
997 my $log = &log((i + $z) / (i - $z));
998 return ip2 * $log;
999}
1000
1001#
1002# asec
1003#
1004# Computes the arc secant asec(z) = acos(1 / z).
1005#
1006sub asec {
1007 my ($z) = @_;
1008 _divbyzero "asec($z)", $z if ($z == 0);
1009 return acos(1 / $z);
1010}
1011
1012#
1013# acsc
1014#
1015# Computes the arc cosecant acsc(z) = asin(1 / z).
1016#
1017sub acsc {
1018 my ($z) = @_;
1019 _divbyzero "acsc($z)", $z if ($z == 0);
1020 return asin(1 / $z);
1021}
1022
1023#
1024# acosec
1025#
1026# Alias for acsc().
1027#
1028sub acosec { Math::Complex::acsc(@_) }
1029
1030#
1031# acot
1032#
1033# Computes the arc cotangent acot(z) = atan(1 / z)
1034#
1035sub acot {
1036 my ($z) = @_;
1037 _divbyzero "acot(0)" if $z == 0;
1038 return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z)
1039 unless ref $z;
1040 _divbyzero "acot(i)" if ($z - i == 0);
1041 _logofzero "acot(-i)" if ($z + i == 0);
1042 return atan(1 / $z);
1043}
1044
1045#
1046# acotan
1047#
1048# Alias for acot().
1049#
1050sub acotan { Math::Complex::acot(@_) }
1051
1052#
1053# cosh
1054#
1055# Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
1056#
1057sub cosh {
1058 my ($z) = @_;
1059 my $ex;
1060 unless (ref $z) {
1061 $ex = CORE::exp($z);
1062 return $ex ? ($ex + 1/$ex)/2 : $Inf;
1063 }
1064 my ($x, $y) = @{$z->cartesian};
1065 $ex = CORE::exp($x);
1066 my $ex_1 = $ex ? 1 / $ex : $Inf;
1067 return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
1068 CORE::sin($y) * ($ex - $ex_1)/2);
1069}
1070
1071#
1072# sinh
1073#
1074# Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
1075#
1076sub sinh {
1077 my ($z) = @_;
1078 my $ex;
1079 unless (ref $z) {
1080 return 0 if $z == 0;
1081 $ex = CORE::exp($z);
1082 return $ex ? ($ex - 1/$ex)/2 : "-$Inf";
1083 }
1084 my ($x, $y) = @{$z->cartesian};
1085 my $cy = CORE::cos($y);
1086 my $sy = CORE::sin($y);
1087 $ex = CORE::exp($x);
1088 my $ex_1 = $ex ? 1 / $ex : $Inf;
1089 return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
1090 CORE::sin($y) * ($ex + $ex_1)/2);
1091}
1092
1093#
1094# tanh
1095#
1096# Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
1097#
1098sub tanh {
1099 my ($z) = @_;
1100 my $cz = cosh($z);
1101 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
1102 return sinh($z) / $cz;
1103}
1104
1105#
1106# sech
1107#
1108# Computes the hyperbolic secant sech(z) = 1 / cosh(z).
1109#
1110sub sech {
1111 my ($z) = @_;
1112 my $cz = cosh($z);
1113 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
1114 return 1 / $cz;
1115}
1116
1117#
1118# csch
1119#
1120# Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
1121#
1122sub csch {
1123 my ($z) = @_;
1124 my $sz = sinh($z);
1125 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
1126 return 1 / $sz;
1127}
1128
1129#
1130# cosech
1131#
1132# Alias for csch().
1133#
1134sub cosech { Math::Complex::csch(@_) }
1135
1136#
1137# coth
1138#
1139# Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
1140#
1141sub coth {
1142 my ($z) = @_;
1143 my $sz = sinh($z);
1144 _divbyzero "coth($z)", "sinh($z)" if $sz == 0;
1145 return cosh($z) / $sz;
1146}
1147
1148#
1149# cotanh
1150#
1151# Alias for coth().
1152#
1153sub cotanh { Math::Complex::coth(@_) }
1154
1155#
1156# acosh
1157#
1158# Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
1159#
1160sub acosh {
1161 my ($z) = @_;
1162 unless (ref $z) {
1163 $z = cplx($z, 0);
1164 }
1165 my ($re, $im) = @{$z->cartesian};
1166 if ($im == 0) {
1167 return CORE::log($re + CORE::sqrt($re*$re - 1))
1168 if $re >= 1;
1169 return cplx(0, CORE::atan2(CORE::sqrt(1 - $re*$re), $re))
1170 if CORE::abs($re) < 1;
1171 }
1172 my $t = &sqrt($z * $z - 1) + $z;
1173 # Try Taylor if looking bad (this usually means that
1174 # $z was large negative, therefore the sqrt is really
1175 # close to abs(z), summing that with z...)
1176 $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1177 if $t == 0;
1178 my $u = &log($t);
1179 $u->Im(-$u->Im) if $re < 0 && $im == 0;
1180 return $re < 0 ? -$u : $u;
1181}
1182
1183#
1184# asinh
1185#
1186# Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z+1))
1187#
1188sub asinh {
1189 my ($z) = @_;
1190 unless (ref $z) {
1191 my $t = $z + CORE::sqrt($z*$z + 1);
1192 return CORE::log($t) if $t;
1193 }
1194 my $t = &sqrt($z * $z + 1) + $z;
1195 # Try Taylor if looking bad (this usually means that
1196 # $z was large negative, therefore the sqrt is really
1197 # close to abs(z), summing that with z...)
1198 $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
1199 if $t == 0;
1200 return &log($t);
1201}
1202
1203#
1204# atanh
1205#
1206# Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
1207#
1208sub atanh {
1209 my ($z) = @_;
1210 unless (ref $z) {
1211 return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
1212 $z = cplx($z, 0);
1213 }
1214 _divbyzero 'atanh(1)', "1 - $z" if (1 - $z == 0);
1215 _logofzero 'atanh(-1)' if (1 + $z == 0);
1216 return 0.5 * &log((1 + $z) / (1 - $z));
1217}
1218
1219#
1220# asech
1221#
1222# Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
1223#
1224sub asech {
1225 my ($z) = @_;
1226 _divbyzero 'asech(0)', "$z" if ($z == 0);
1227 return acosh(1 / $z);
1228}
1229
1230#
1231# acsch
1232#
1233# Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
1234#
1235sub acsch {
1236 my ($z) = @_;
1237 _divbyzero 'acsch(0)', $z if ($z == 0);
1238 return asinh(1 / $z);
1239}
1240
1241#
1242# acosech
1243#
1244# Alias for acosh().
1245#
1246sub acosech { Math::Complex::acsch(@_) }
1247
1248#
1249# acoth
1250#
1251# Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
1252#
1253sub acoth {
1254 my ($z) = @_;
1255 _divbyzero 'acoth(0)' if ($z == 0);
1256 unless (ref $z) {
1257 return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
1258 $z = cplx($z, 0);
1259 }
1260 _divbyzero 'acoth(1)', "$z - 1" if ($z - 1 == 0);
1261 _logofzero 'acoth(-1)', "1 + $z" if (1 + $z == 0);
1262 return &log((1 + $z) / ($z - 1)) / 2;
1263}
1264
1265#
1266# acotanh
1267#
1268# Alias for acot().
1269#
1270sub acotanh { Math::Complex::acoth(@_) }
1271
1272#
1273# (atan2)
1274#
1275# Compute atan(z1/z2), minding the right quadrant.
1276#
1277sub atan2 {
1278 my ($z1, $z2, $inverted) = @_;
1279 my ($re1, $im1, $re2, $im2);
1280 if ($inverted) {
1281 ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1282 ($re2, $im2) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
1283 } else {
1284 ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
1285 ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
1286 }
1287 if ($im1 || $im2) {
1288 # In MATLAB the imaginary parts are ignored.
1289 # warn "atan2: Imaginary parts ignored";
1290 # http://documents.wolfram.com/mathematica/functions/ArcTan
1291 # NOTE: Mathematica ArcTan[x,y] while atan2(y,x)
1292 my $s = $z1 * $z1 + $z2 * $z2;
1293 _divbyzero("atan2") if $s == 0;
1294 my $i = &i;
1295 my $r = $z2 + $z1 * $i;
1296 return -$i * &log($r / &sqrt( $s ));
1297 }
1298 return CORE::atan2($re1, $re2);
1299}
1300
1301#
1302# display_format
1303# ->display_format
1304#
1305# Set (get if no argument) the display format for all complex numbers that
1306# don't happen to have overridden it via ->display_format
1307#
1308# When called as an object method, this actually sets the display format for
1309# the current object.
1310#
1311# Valid object formats are 'c' and 'p' for cartesian and polar. The first
1312# letter is used actually, so the type can be fully spelled out for clarity.
1313#
1314sub display_format {
1315 my $self = shift;
1316 my %display_format = %DISPLAY_FORMAT;
1317
1318 if (ref $self) { # Called as an object method
1319 if (exists $self->{display_format}) {
1320 my %obj = %{$self->{display_format}};
1321 @display_format{keys %obj} = values %obj;
1322 }
1323 }
1324 if (@_ == 1) {
1325 $display_format{style} = shift;
1326 } else {
1327 my %new = @_;
1328 @display_format{keys %new} = values %new;
1329 }
1330
1331 if (ref $self) { # Called as an object method
1332 $self->{display_format} = { %display_format };
1333 return
1334 wantarray ?
1335 %{$self->{display_format}} :
1336 $self->{display_format}->{style};
1337 }
1338
1339 # Called as a class method
1340 %DISPLAY_FORMAT = %display_format;
1341 return
1342 wantarray ?
1343 %DISPLAY_FORMAT :
1344 $DISPLAY_FORMAT{style};
1345}
1346
1347#
1348# (stringify)
1349#
1350# Show nicely formatted complex number under its cartesian or polar form,
1351# depending on the current display format:
1352#
1353# . If a specific display format has been recorded for this object, use it.
1354# . Otherwise, use the generic current default for all complex numbers,
1355# which is a package global variable.
1356#
1357sub stringify {
1358 my ($z) = shift;
1359
1360 my $style = $z->display_format;
1361
1362 $style = $DISPLAY_FORMAT{style} unless defined $style;
1363
1364 return $z->stringify_polar if $style =~ /^p/i;
1365 return $z->stringify_cartesian;
1366}
1367
1368#
1369# ->stringify_cartesian
1370#
1371# Stringify as a cartesian representation 'a+bi'.
1372#
1373sub stringify_cartesian {
1374 my $z = shift;
1375 my ($x, $y) = @{$z->cartesian};
1376 my ($re, $im);
1377
1378 my %format = $z->display_format;
1379 my $format = $format{format};
1380
1381 if ($x) {
1382 if ($x =~ /^NaN[QS]?$/i) {
1383 $re = $x;
1384 } else {
1385 if ($x =~ /^-?$Inf$/oi) {
1386 $re = $x;
1387 } else {
1388 $re = defined $format ? sprintf($format, $x) : $x;
1389 }
1390 }
1391 } else {
1392 undef $re;
1393 }
1394
1395 if ($y) {
1396 if ($y =~ /^(NaN[QS]?)$/i) {
1397 $im = $y;
1398 } else {
1399 if ($y =~ /^-?$Inf$/oi) {
1400 $im = $y;
1401 } else {
1402 $im =
1403 defined $format ?
1404 sprintf($format, $y) :
1405 ($y == 1 ? "" : ($y == -1 ? "-" : $y));
1406 }
1407 }
1408 $im .= "i";
1409 } else {
1410 undef $im;
1411 }
1412
1413 my $str = $re;
1414
1415 if (defined $im) {
1416 if ($y < 0) {
1417 $str .= $im;
1418 } elsif ($y > 0 || $im =~ /^NaN[QS]?i$/i) {
1419 $str .= "+" if defined $re;
1420 $str .= $im;
1421 }
1422 } elsif (!defined $re) {
1423 $str = "0";
1424 }
1425
1426 return $str;
1427}
1428
1429
1430#
1431# ->stringify_polar
1432#
1433# Stringify as a polar representation '[r,t]'.
1434#
1435sub stringify_polar {
1436 my $z = shift;
1437 my ($r, $t) = @{$z->polar};
1438 my $theta;
1439
1440 my %format = $z->display_format;
1441 my $format = $format{format};
1442
1443 if ($t =~ /^NaN[QS]?$/i || $t =~ /^-?$Inf$/oi) {
1444 $theta = $t;
1445 } elsif ($t == pi) {
1446 $theta = "pi";
1447 } elsif ($r == 0 || $t == 0) {
1448 $theta = defined $format ? sprintf($format, $t) : $t;
1449 }
1450
1451 return "[$r,$theta]" if defined $theta;
1452
1453 #
1454 # Try to identify pi/n and friends.
1455 #
1456
1457 $t -= int(CORE::abs($t) / pit2) * pit2;
1458
1459 if ($format{polar_pretty_print} && $t) {
1460 my ($a, $b);
1461 for $a (2..9) {
1462 $b = $t * $a / pi;
1463 if ($b =~ /^-?\d+$/) {
1464 $b = $b < 0 ? "-" : "" if CORE::abs($b) == 1;
1465 $theta = "${b}pi/$a";
1466 last;
1467 }
1468 }
1469 }
1470
1471 if (defined $format) {
1472 $r = sprintf($format, $r);
1473 $theta = sprintf($format, $theta) unless defined $theta;
1474 } else {
1475 $theta = $t unless defined $theta;
1476 }
1477
1478 return "[$r,$theta]";
1479}
1480
14811;
1482__END__
1483
1484=pod
1485
1486=head1 NAME
1487
1488Math::Complex - complex numbers and associated mathematical functions
1489
1490=head1 SYNOPSIS
1491
1492 use Math::Complex;
1493
1494 $z = Math::Complex->make(5, 6);
1495 $t = 4 - 3*i + $z;
1496 $j = cplxe(1, 2*pi/3);
1497
1498=head1 DESCRIPTION
1499
1500This package lets you create and manipulate complex numbers. By default,
1501I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1502full complex support, along with a full set of mathematical functions
1503typically associated with and/or extended to complex numbers.
1504
1505If you wonder what complex numbers are, they were invented to be able to solve
1506the following equation:
1507
1508 x*x = -1
1509
1510and by definition, the solution is noted I<i> (engineers use I<j> instead since
1511I<i> usually denotes an intensity, but the name does not matter). The number
1512I<i> is a pure I<imaginary> number.
1513
1514The arithmetics with pure imaginary numbers works just like you would expect
1515it with real numbers... you just have to remember that
1516
1517 i*i = -1
1518
1519so you have:
1520
1521 5i + 7i = i * (5 + 7) = 12i
1522 4i - 3i = i * (4 - 3) = i
1523 4i * 2i = -8
1524 6i / 2i = 3
1525 1 / i = -i
1526
1527Complex numbers are numbers that have both a real part and an imaginary
1528part, and are usually noted:
1529
1530 a + bi
1531
1532where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1533arithmetic with complex numbers is straightforward. You have to
1534keep track of the real and the imaginary parts, but otherwise the
1535rules used for real numbers just apply:
1536
1537 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1538 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1539
1540A graphical representation of complex numbers is possible in a plane
1541(also called the I<complex plane>, but it's really a 2D plane).
1542The number
1543
1544 z = a + bi
1545
1546is the point whose coordinates are (a, b). Actually, it would
1547be the vector originating from (0, 0) to (a, b). It follows that the addition
1548of two complex numbers is a vectorial addition.
1549
1550Since there is a bijection between a point in the 2D plane and a complex
1551number (i.e. the mapping is unique and reciprocal), a complex number
1552can also be uniquely identified with polar coordinates:
1553
1554 [rho, theta]
1555
1556where C<rho> is the distance to the origin, and C<theta> the angle between
1557the vector and the I<x> axis. There is a notation for this using the
1558exponential form, which is:
1559
1560 rho * exp(i * theta)
1561
1562where I<i> is the famous imaginary number introduced above. Conversion
1563between this form and the cartesian form C<a + bi> is immediate:
1564
1565 a = rho * cos(theta)
1566 b = rho * sin(theta)
1567
1568which is also expressed by this formula:
1569
1570 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1571
1572In other words, it's the projection of the vector onto the I<x> and I<y>
1573axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1574the I<argument> of the complex number. The I<norm> of C<z> will be
1575noted C<abs(z)>.
1576
1577The polar notation (also known as the trigonometric
1578representation) is much more handy for performing multiplications and
1579divisions of complex numbers, whilst the cartesian notation is better
1580suited for additions and subtractions. Real numbers are on the I<x>
1581axis, and therefore I<theta> is zero or I<pi>.
1582
1583All the common operations that can be performed on a real number have
1584been defined to work on complex numbers as well, and are merely
1585I<extensions> of the operations defined on real numbers. This means
1586they keep their natural meaning when there is no imaginary part, provided
1587the number is within their definition set.
1588
1589For instance, the C<sqrt> routine which computes the square root of
1590its argument is only defined for non-negative real numbers and yields a
1591non-negative real number (it is an application from B<R+> to B<R+>).
1592If we allow it to return a complex number, then it can be extended to
1593negative real numbers to become an application from B<R> to B<C> (the
1594set of complex numbers):
1595
1596 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1597
1598It can also be extended to be an application from B<C> to B<C>,
1599whilst its restriction to B<R> behaves as defined above by using
1600the following definition:
1601
1602 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1603
1604Indeed, a negative real number can be noted C<[x,pi]> (the modulus
1605I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative
1606number) and the above definition states that
1607
1608 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1609
1610which is exactly what we had defined for negative real numbers above.
1611The C<sqrt> returns only one of the solutions: if you want the both,
1612use the C<root> function.
1613
1614All the common mathematical functions defined on real numbers that
1615are extended to complex numbers share that same property of working
1616I<as usual> when the imaginary part is zero (otherwise, it would not
1617be called an extension, would it?).
1618
1619A I<new> operation possible on a complex number that is
1620the identity for real numbers is called the I<conjugate>, and is noted
1621with a horizontal bar above the number, or C<~z> here.
1622
1623 z = a + bi
1624 ~z = a - bi
1625
1626Simple... Now look:
1627
1628 z * ~z = (a + bi) * (a - bi) = a*a + b*b
1629
1630We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1631distance to the origin, also known as:
1632
1633 rho = abs(z) = sqrt(a*a + b*b)
1634
1635so
1636
1637 z * ~z = abs(z) ** 2
1638
1639If z is a pure real number (i.e. C<b == 0>), then the above yields:
1640
1641 a * a = abs(a) ** 2
1642
1643which is true (C<abs> has the regular meaning for real number, i.e. stands
1644for the absolute value). This example explains why the norm of C<z> is
1645noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1646is the regular C<abs> we know when the complex number actually has no
1647imaginary part... This justifies I<a posteriori> our use of the C<abs>
1648notation for the norm.
1649
1650=head1 OPERATIONS
1651
1652Given the following notations:
1653
1654 z1 = a + bi = r1 * exp(i * t1)
1655 z2 = c + di = r2 * exp(i * t2)
1656 z = <any complex or real number>
1657
1658the following (overloaded) operations are supported on complex numbers:
1659
1660 z1 + z2 = (a + c) + i(b + d)
1661 z1 - z2 = (a - c) + i(b - d)
1662 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1663 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1664 z1 ** z2 = exp(z2 * log z1)
1665 ~z = a - bi
1666 abs(z) = r1 = sqrt(a*a + b*b)
1667 sqrt(z) = sqrt(r1) * exp(i * t/2)
1668 exp(z) = exp(a) * exp(i * b)
1669 log(z) = log(r1) + i*t
1670 sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
1671 cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
1672 atan2(y, x) = atan(y / x) # Minding the right quadrant, note the order.
1673
1674The definition used for complex arguments of atan2() is
1675
1676 -i log((x + iy)/sqrt(x*x+y*y))
1677
1678The following extra operations are supported on both real and complex
1679numbers:
1680
1681 Re(z) = a
1682 Im(z) = b
1683 arg(z) = t
1684 abs(z) = r
1685
1686 cbrt(z) = z ** (1/3)
1687 log10(z) = log(z) / log(10)
1688 logn(z, n) = log(z) / log(n)
1689
1690 tan(z) = sin(z) / cos(z)
1691
1692 csc(z) = 1 / sin(z)
1693 sec(z) = 1 / cos(z)
1694 cot(z) = 1 / tan(z)
1695
1696 asin(z) = -i * log(i*z + sqrt(1-z*z))
1697 acos(z) = -i * log(z + i*sqrt(1-z*z))
1698 atan(z) = i/2 * log((i+z) / (i-z))
1699
1700 acsc(z) = asin(1 / z)
1701 asec(z) = acos(1 / z)
1702 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
1703
1704 sinh(z) = 1/2 (exp(z) - exp(-z))
1705 cosh(z) = 1/2 (exp(z) + exp(-z))
1706 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1707
1708 csch(z) = 1 / sinh(z)
1709 sech(z) = 1 / cosh(z)
1710 coth(z) = 1 / tanh(z)
1711
1712 asinh(z) = log(z + sqrt(z*z+1))
1713 acosh(z) = log(z + sqrt(z*z-1))
1714 atanh(z) = 1/2 * log((1+z) / (1-z))
1715
1716 acsch(z) = asinh(1 / z)
1717 asech(z) = acosh(1 / z)
1718 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1719
1720I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>,
1721I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>,
1722I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>,
1723I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>,
1724C<rho>, and C<theta> can be used also as mutators. The C<cbrt>
1725returns only one of the solutions: if you want all three, use the
1726C<root> function.
1727
1728The I<root> function is available to compute all the I<n>
1729roots of some complex, where I<n> is a strictly positive integer.
1730There are exactly I<n> such roots, returned as a list. Getting the
1731number mathematicians call C<j> such that:
1732
1733 1 + j + j*j = 0;
1734
1735is a simple matter of writing:
1736
1737 $j = ((root(1, 3))[1];
1738
1739The I<k>th root for C<z = [r,t]> is given by:
1740
1741 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1742
1743You can return the I<k>th root directly by C<root(z, n, k)>,
1744indexing starting from I<zero> and ending at I<n - 1>.
1745
1746The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1747order to ensure its restriction to real numbers is conform to what you
1748would expect, the comparison is run on the real part of the complex
1749number first, and imaginary parts are compared only when the real
1750parts match.
1751
1752=head1 CREATION
1753
1754To create a complex number, use either:
1755
1756 $z = Math::Complex->make(3, 4);
1757 $z = cplx(3, 4);
1758
1759if you know the cartesian form of the number, or
1760
1761 $z = 3 + 4*i;
1762
1763if you like. To create a number using the polar form, use either:
1764
1765 $z = Math::Complex->emake(5, pi/3);
1766 $x = cplxe(5, pi/3);
1767
1768instead. The first argument is the modulus, the second is the angle
1769(in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a
1770notation for complex numbers in the polar form).
1771
1772It is possible to write:
1773
1774 $x = cplxe(-3, pi/4);
1775
1776but that will be silently converted into C<[3,-3pi/4]>, since the
1777modulus must be non-negative (it represents the distance to the origin
1778in the complex plane).
1779
1780It is also possible to have a complex number as either argument of the
1781C<make>, C<emake>, C<cplx>, and C<cplxe>: the appropriate component of
1782the argument will be used.
1783
1784 $z1 = cplx(-2, 1);
1785 $z2 = cplx($z1, 4);
1786
1787The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also
1788understand a single (string) argument of the forms
1789
1790 2-3i
1791 -3i
1792 [2,3]
1793 [2,-3pi/4]
1794 [2]
1795
1796in which case the appropriate cartesian and exponential components
1797will be parsed from the string and used to create new complex numbers.
1798The imaginary component and the theta, respectively, will default to zero.
1799
1800The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also
1801understand the case of no arguments: this means plain zero or (0, 0).
1802
1803=head1 DISPLAYING
1804
1805When printed, a complex number is usually shown under its cartesian
1806style I<a+bi>, but there are legitimate cases where the polar style
1807I<[r,t]> is more appropriate. The process of converting the complex
1808number into a string that can be displayed is known as I<stringification>.
1809
1810By calling the class method C<Math::Complex::display_format> and
1811supplying either C<"polar"> or C<"cartesian"> as an argument, you
1812override the default display style, which is C<"cartesian">. Not
1813supplying any argument returns the current settings.
1814
1815This default can be overridden on a per-number basis by calling the
1816C<display_format> method instead. As before, not supplying any argument
1817returns the current display style for this number. Otherwise whatever you
1818specify will be the new display style for I<this> particular number.
1819
1820For instance:
1821
1822 use Math::Complex;
1823
1824 Math::Complex::display_format('polar');
1825 $j = (root(1, 3))[1];
1826 print "j = $j\n"; # Prints "j = [1,2pi/3]"
1827 $j->display_format('cartesian');
1828 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1829
1830The polar style attempts to emphasize arguments like I<k*pi/n>
1831(where I<n> is a positive integer and I<k> an integer within [-9, +9]),
1832this is called I<polar pretty-printing>.
1833
1834For the reverse of stringifying, see the C<make> and C<emake>.
1835
1836=head2 CHANGED IN PERL 5.6
1837
1838The C<display_format> class method and the corresponding
1839C<display_format> object method can now be called using
1840a parameter hash instead of just a one parameter.
1841
1842The old display format style, which can have values C<"cartesian"> or
1843C<"polar">, can be changed using the C<"style"> parameter.
1844
1845 $j->display_format(style => "polar");
1846
1847The one parameter calling convention also still works.
1848
1849 $j->display_format("polar");
1850
1851There are two new display parameters.
1852
1853The first one is C<"format">, which is a sprintf()-style format string
1854to be used for both numeric parts of the complex number(s). The is
1855somewhat system-dependent but most often it corresponds to C<"%.15g">.
1856You can revert to the default by setting the C<format> to C<undef>.
1857
1858 # the $j from the above example
1859
1860 $j->display_format('format' => '%.5f');
1861 print "j = $j\n"; # Prints "j = -0.50000+0.86603i"
1862 $j->display_format('format' => undef);
1863 print "j = $j\n"; # Prints "j = -0.5+0.86603i"
1864
1865Notice that this affects also the return values of the
1866C<display_format> methods: in list context the whole parameter hash
1867will be returned, as opposed to only the style parameter value.
1868This is a potential incompatibility with earlier versions if you
1869have been calling the C<display_format> method in list context.
1870
1871The second new display parameter is C<"polar_pretty_print">, which can
1872be set to true or false, the default being true. See the previous
1873section for what this means.
1874
1875=head1 USAGE
1876
1877Thanks to overloading, the handling of arithmetics with complex numbers
1878is simple and almost transparent.
1879
1880Here are some examples:
1881
1882 use Math::Complex;
1883
1884 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1885 print "j = $j, j**3 = ", $j ** 3, "\n";
1886 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1887
1888 $z = -16 + 0*i; # Force it to be a complex
1889 print "sqrt($z) = ", sqrt($z), "\n";
1890
1891 $k = exp(i * 2*pi/3);
1892 print "$j - $k = ", $j - $k, "\n";
1893
1894 $z->Re(3); # Re, Im, arg, abs,
1895 $j->arg(2); # (the last two aka rho, theta)
1896 # can be used also as mutators.
1897
1898=head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1899
1900The division (/) and the following functions
1901
1902 log ln log10 logn
1903 tan sec csc cot
1904 atan asec acsc acot
1905 tanh sech csch coth
1906 atanh asech acsch acoth
1907
1908cannot be computed for all arguments because that would mean dividing
1909by zero or taking logarithm of zero. These situations cause fatal
1910runtime errors looking like this
1911
1912 cot(0): Division by zero.
1913 (Because in the definition of cot(0), the divisor sin(0) is 0)
1914 Died at ...
1915
1916or
1917
1918 atanh(-1): Logarithm of zero.
1919 Died at...
1920
1921For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1922C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the
1923logarithmic functions and the C<atanh>, C<acoth>, the argument cannot
1924be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be
1925C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be
1926C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
1927cannot be C<-i> (the negative imaginary unit). For the C<tan>,
1928C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
1929is any integer. atan2(0, 0) is undefined, and if the complex arguments
1930are used for atan2(), a division by zero will happen if z1**2+z2**2 == 0.
1931
1932Note that because we are operating on approximations of real numbers,
1933these errors can happen when merely `too close' to the singularities
1934listed above.
1935
1936=head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
1937
1938The C<make> and C<emake> accept both real and complex arguments.
1939When they cannot recognize the arguments they will die with error
1940messages like the following
1941
1942 Math::Complex::make: Cannot take real part of ...
1943 Math::Complex::make: Cannot take real part of ...
1944 Math::Complex::emake: Cannot take rho of ...
1945 Math::Complex::emake: Cannot take theta of ...
1946
1947=head1 BUGS
1948
1949Saying C<use Math::Complex;> exports many mathematical routines in the
1950caller environment and even overrides some (C<sqrt>, C<log>, C<atan2>).
1951This is construed as a feature by the Authors, actually... ;-)
1952
1953All routines expect to be given real or complex numbers. Don't attempt to
1954use BigFloat, since Perl has currently no rule to disambiguate a '+'
1955operation (for instance) between two overloaded entities.
1956
1957In Cray UNICOS there is some strange numerical instability that results
1958in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware.
1959The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
1960Whatever it is, it does not manifest itself anywhere else where Perl runs.
1961
1962=head1 AUTHORS
1963
1964Daniel S. Lewart <F<[email protected]>>
1965
1966Original authors Raphael Manfredi <F<[email protected]>> and
1967Jarkko Hietaniemi <F<[email protected]>>
1968
1969=cut
1970
19711;
1972
1973# eof
Note: See TracBrowser for help on using the repository browser.