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