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