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