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