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