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