This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
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 #
8 package Math::Complex;
10 use vars qw(\$VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS \$Inf);
12 \$VERSION = 1.42;
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 }
45 use strict;
47 my \$i;
48 my %LOGN;
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;
54 require Exporter;
56 @ISA = qw(Exporter);
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              );
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);
79 my @pi = qw(pi pi2 pi4 pip2 pip4 Inf);
81 @EXPORT_OK = @pi;
83 %EXPORT_TAGS = (
84     'trig' => [@trig],
85     'pi' => [@pi],
86 );
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;
108 #
109 # Package "privates"
110 #
112 my %DISPLAY_FORMAT = ('style' => 'cartesian',
113                       'polar_pretty_print' => 1);
114 my \$eps            = 1e-14;             # Epsilon
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 #
125 # Die on bad *make() arguments.
127 sub _cannot_make {
128     die "@{[(caller(1))]}: Cannot take \$_ of '\$_'.\n";
129 }
131 sub _make {
132     my \$arg = shift;
133     my (\$p, \$q);
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     }
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     }
150     return (\$p, \$q);
151 }
153 sub _emake {
154     my \$arg = shift;
155     my (\$p, \$q);
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     }
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     }
174     return (\$p, \$q);
175 }
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(\$_)
189             if (\$_ =~ /^\s*\[/);
190         (\$re, \$im) = _make(\$_);
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');
202     return \$self;
203 }
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(\$_)
217             if (\$_ =~ /^\s*\(/ || \$_ =~ /i\s*\$/);
218         (\$rho, \$theta) = _emake(\$_);
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');
236     return \$self;
237 }
239 sub new { &make }               # For backward compatibility only.
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 }
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 }
261 #
262 # pi
263 #
264 # The number defined as pi = 180 degrees
265 #
266 sub pi () { 4 * CORE::atan2(1, 1) }
268 #
269 # pi2
270 #
271 # The full circle
272 #
273 sub pi2 () { 2 * pi }
275 #
276 # pi4
277 #
278 # The full circle twice.
279 #
280 sub pi4 () { 4 * pi }
282 #
283 # pip2
284 #
285 # The quarter circle
286 #
287 sub pip2 () { pi / 2 }
289 #
290 # pip4
291 #
292 # The eighth circle.
293 #
294 sub pip4 () { pi / 4 }
296 #
297 # _uplog10
298 #
299 # Used in log10().
300 #
301 sub _uplog10 () { 1 / CORE::log(10) }
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 }
318 #
319 # _ip2
320 #
321 # Half of i.
322 #
323 sub _ip2 () { i / 2 }
325 #
326 # Attribute access/set routines
327 #
329 sub _cartesian {\$_->{c_dirty} ?
330                    \$_->_update_cartesian : \$_->{'cartesian'}}
331 sub _polar     {\$_->{p_dirty} ?
332                    \$_->_update_polar : \$_->{'polar'}}
334 sub _set_cartesian { \$_->{p_dirty}++; \$_->{c_dirty} = 0;
335                      \$_->{'cartesian'} = \$_ }
336 sub _set_polar     { \$_->{c_dirty}++; \$_->{p_dirty} = 0;
337                      \$_->{'polar'} = \$_ }
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 }
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 }
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 }
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);
401 }
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 }
433 #
434 # _divbyzero
435 #
436 # Die on division by zero.
437 #
438 sub _divbyzero {
439     my \$mess = "\$_: Division by zero.\n";
441     if (defined \$_) {
442         \$mess .= "(Because in the definition of \$_, the divisor ";
443         \$mess .= "\$_ " unless ("\$_" eq '0');
444         \$mess .= "is 0)\n";
445     }
447     my @up = caller(1);
449     \$mess .= "Died at \$up line \$up.\n";
451     die \$mess;
452 }
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 }
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 }
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 }
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 }
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 }
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 }
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                 \$_ = \$_;
594             } else {
595                 return CORE::abs(\$z);
596             }
597         }
598         if (defined \$rho) {
599             \$z->{'polar'} = [ \$rho, \${\$z->_polar} ];
600             \$z->{p_dirty} = 0;
601             \$z->{c_dirty} = 1;
602             return \$rho;
603         } else {
604             return \${\$z->_polar};
605         }
606 }
608 sub _theta {
609     my \$theta = \$_;
611     if    (\$\$theta >   pi()) { \$\$theta -= pi2 }
612     elsif (\$\$theta <= -pi()) { \$\$theta += pi2 }
613 }
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}, \$theta ];
626             \$z->{p_dirty} = 0;
627             \$z->{c_dirty} = 1;
628         } else {
629             \$theta = \${\$z->_polar};
630             _theta(\\$theta);
631         }
632         return \$theta;
633 }
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 }
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 }
677 #
679 #
680 # Die on bad root.
681 #
683     my \$mess = "Root '\$_' illegal, root rank must be positive integer.\n";
685     my @up = caller(1);
687     \$mess .= "Died at \$up line \$up.\n";
689     die \$mess;
690 }
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 }
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} ];
736             \$z->{c_dirty} = 0;
737             \$z->{p_dirty} = 1;
738         } else {
739             return \${\$z->_cartesian};
740         }
741 }
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}, \$Im ];
753             \$z->{c_dirty} = 0;
754             \$z->{p_dirty} = 1;
755         } else {
756             return \${\$z->_cartesian};
757         }
758 }
760 #
761 # rho
762 #
763 # Return or set rho(w).
764 #
765 sub rho {
766     Math::Complex::abs(@_);
767 }
769 #
770 # theta
771 #
772 # Return or set theta(w).
773 #
774 sub theta {
775     Math::Complex::arg(@_);
776 }
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 }
789 #
790 # _logofzero
791 #
792 # Die on logarithm of zero.
793 #
794 sub _logofzero {
795     my \$mess = "\$_: Logarithm of zero.\n";
797     if (defined \$_) {
798         \$mess .= "(Because in the definition of \$_, the argument ";
799         \$mess .= "\$_ " unless (\$_ eq '0');
800         \$mess .= "is 0)\n";
801     }
803     my @up = caller(1);
805     \$mess .= "Died at \$up line \$up.\n";
807     die \$mess;
808 }
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 }
828 #
829 # ln
830 #
831 # Alias for log().
832 #
833 sub ln { Math::Complex::log(@_) }
835 #
836 # log10
837 #
838 # Compute log10(z).
839 #
841 sub log10 {
842         return Math::Complex::log(\$_) * _uplog10;
843 }
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 }
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 }
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 }
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 }
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 }
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 }
928 #
929 # cosec
930 #
931 # Alias for csc().
932 #
933 sub cosec { Math::Complex::csc(@_) }
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 }
947 #
948 # cotan
949 #
950 # Alias for cot().
951 #
952 sub cotan { Math::Complex::cot(@_) }
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 = \$_;
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 }
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 = \$_;
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 }
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 }
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 }
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 }
1042 #
1043 # acosec
1044 #
1045 # Alias for acsc().
1046 #
1047 sub acosec { Math::Complex::acsc(@_) }
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 }
1064 #
1065 # acotan
1066 #
1067 # Alias for acot().
1068 #
1069 sub acotan { Math::Complex::acot(@_) }
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 }
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 }
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 }
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 }
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 }
1151 #
1152 # cosech
1153 #
1154 # Alias for csch().
1155 #
1156 sub cosech { Math::Complex::csch(@_) }
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 }
1173 #
1174 # cotanh
1175 #
1176 # Alias for coth().
1177 #
1178 sub cotanh { Math::Complex::coth(@_) }
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 }
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 }
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 }
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 }
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 }
1266 #
1267 # acosech
1268 #
1269 # Alias for acosh().
1270 #
1271 sub acosech { Math::Complex::acsch(@_) }
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 }
1290 #
1291 # acotanh
1292 #
1293 # Alias for acot().
1294 #
1295 sub acotanh { Math::Complex::acoth(@_) }
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 }
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;
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         }
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         }
1364         # Called as a class method
1365         %DISPLAY_FORMAT = %display_format;
1366         return
1367             wantarray ?
1368                 %DISPLAY_FORMAT :
1369                     \$DISPLAY_FORMAT{style};
1370 }
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;
1385         my \$style = \$z->display_format;
1387         \$style = \$DISPLAY_FORMAT{style} unless defined \$style;
1389         return \$z->_stringify_polar if \$style =~ /^p/i;
1390         return \$z->_stringify_cartesian;
1391 }
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);
1403         my %format = \$z->display_format;
1404         my \$format = \$format{format};
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         }
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         }
1438         my \$str = \$re;
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         }
1451         return \$str;
1452 }
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;
1465         my %format = \$z->display_format;
1466         my \$format = \$format{format};
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         }
1476         return "[\$r,\$theta]" if defined \$theta;
1478         #
1479         # Try to identify pi/n and friends.
1480         #
1482         \$t -= int(CORE::abs(\$t) / pi2) * pi2;
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         }
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         }
1503         return "[\$r,\$theta]";
1504 }
1506 sub Inf {
1507     return \$Inf;
1508 }
1510 1;
1511 __END__
1513 =pod
1517 Math::Complex - complex numbers and associated mathematical functions
1521         use Math::Complex;
1523         \$z = Math::Complex->make(5, 6);
1524         \$t = 4 - 3*i + \$z;
1525         \$j = cplxe(1, 2*pi/3);
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.
1534 If you wonder what complex numbers are, they were invented to be able to solve
1535 the following equation:
1537         x*x = -1
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.
1543 The arithmetics with pure imaginary numbers works just like you would expect
1544 it with real numbers... you just have to remember that
1546         i*i = -1
1548 so you have:
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
1556 Complex numbers are numbers that have both a real part and an imaginary
1557 part, and are usually noted:
1559         a + bi
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:
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
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
1573         z = a + bi
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.
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:
1583         [rho, theta]
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:
1589         rho * exp(i * theta)
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:
1594         a = rho * cos(theta)
1595         b = rho * sin(theta)
1597 which is also expressed by this formula:
1599         z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
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)>.
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>.
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.
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):
1625         sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
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:
1631         sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
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
1637         sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
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.
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?).
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.
1652          z = a + bi
1653         ~z = a - bi
1655 Simple... Now look:
1657         z * ~z = (a + bi) * (a - bi) = a*a + b*b
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:
1662         rho = abs(z) = sqrt(a*a + b*b)
1664 so
1666         z * ~z = abs(z) ** 2
1668 If z is a pure real number (i.e. C<b == 0>), then the above yields:
1670         a * a = abs(a) ** 2
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.
1681 Given the following notations:
1683         z1 = a + bi = r1 * exp(i * t1)
1684         z2 = c + di = r2 * exp(i * t2)
1685         z = <any complex or real number>
1687 the following (overloaded) operations are supported on complex numbers:
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.
1703 The definition used for complex arguments of atan2() is
1705        -i log((x + iy)/sqrt(x*x+y*y))
1707 Note that atan2(0, 0) is not well-defined.
1709 The following extra operations are supported on both real and complex
1710 numbers:
1712         Re(z) = a
1713         Im(z) = b
1714         arg(z) = t
1715         abs(z) = r
1717         cbrt(z) = z ** (1/3)
1718         log10(z) = log(z) / log(10)
1719         logn(z, n) = log(z) / log(n)
1721         tan(z) = sin(z) / cos(z)
1723         csc(z) = 1 / sin(z)
1724         sec(z) = 1 / cos(z)
1725         cot(z) = 1 / tan(z)
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))
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))
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))
1739         csch(z) = 1 / sinh(z)
1740         sech(z) = 1 / cosh(z)
1741         coth(z) = 1 / tanh(z)
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))
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))
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.
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:
1764         1 + j + j*j = 0;
1766 is a simple matter of writing:
1768         \$j = ((root(1, 3));
1770 The I<k>th root for C<z = [r,t]> is given by:
1772         (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
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>.
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.
1785 To create a complex number, use either:
1787         \$z = Math::Complex->make(3, 4);
1788         \$z = cplx(3, 4);
1790 if you know the cartesian form of the number, or
1792         \$z = 3 + 4*i;
1794 if you like. To create a number using the polar form, use either:
1796         \$z = Math::Complex->emake(5, pi/3);
1797         \$x = cplxe(5, pi/3);
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).
1803 It is possible to write:
1805         \$x = cplxe(-3, pi/4);
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).
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.
1815         \$z1 = cplx(-2,  1);
1816         \$z2 = cplx(\$z1, 4);
1818 The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also
1819 understand a single (string) argument of the forms
1821         2-3i
1822         -3i
1823         [2,3]
1824         [2,-3pi/4]
1825         
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.
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).
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>.
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.
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.
1851 For instance:
1853         use Math::Complex;
1855         Math::Complex::display_format('polar');
1856         \$j = (root(1, 3));
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"
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>.
1865 For the reverse of stringifying, see the C<make> and C<emake>.
1867 =head2 CHANGED IN PERL 5.6
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.
1873 The old display format style, which can have values C<"cartesian"> or
1874 C<"polar">, can be changed using the C<"style"> parameter.
1876         \$j->display_format(style => "polar");
1878 The one parameter calling convention also still works.
1880         \$j->display_format("polar");
1882 There are two new display parameters.
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>.
1889         # the \$j from the above example
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"
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.
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.
1909 is simple and almost transparent.
1911 Here are some examples:
1913         use Math::Complex;
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";
1919         \$z = -16 + 0*i;                 # Force it to be a complex
1920         print "sqrt(\$z) = ", sqrt(\$z), "\n";
1922         \$k = exp(i * 2*pi/3);
1923         print "\$j - \$k = ", \$j - \$k, "\n";
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.
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:
1935     use Math::Complex ':pi';
1936     \$third_of_circle = pi2 / 3;
1940 The floating point infinity can be exported as a subroutine Inf():
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;
1947 Note that the stringified form of infinity varies between platforms:
1948 it can be for example any of
1950    inf
1951    infinity
1952    INF
1953    1.#INF
1955 or it can be something else.
1957 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
1959 The division (/) and the following functions
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
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
1971         cot(0): Division by zero.
1972         (Because in the definition of cot(0), the divisor sin(0) is 0)
1973         Died at ...
1975 or
1977         atanh(-1): Logarithm of zero.
1978         Died at...
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.
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.
1995 =head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
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
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 ...
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... ;-)
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.
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.