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