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