This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
Missing 'require' in auto-generated .pm by h2xs
[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#
8c03c583
JH
514# _logofzero
515#
516# Die on division by zero.
517#
518sub _logofzero {
519 my $mess = "$_[0]: Logarithm of zero.\n";
520
521 if (defined $_[1]) {
522 $mess .= "(Because in the definition of $_[0], the argument ";
523 $mess .= "$_[1] " unless ($_[1] eq '0');
524 $mess .= "is 0)\n";
525 }
526
527 my @up = caller(1);
528
529 $mess .= "Died at $up[1] line $up[2].\n";
530
531 die $mess;
532}
533
534#
66730be0
RM
535# (log)
536#
537# Compute log(z).
538#
539sub log {
540 my ($z) = @_;
0c721ce2 541 $z = cplx($z, 0) unless ref $z;
0c721ce2 542 my ($x, $y) = @{$z->cartesian};
5cd24f17 543 my ($r, $t) = @{$z->polar};
0c721ce2
JH
544 $t -= 2 * pi if ($t > pi() and $x < 0);
545 $t += 2 * pi if ($t < -pi() and $x < 0);
66730be0
RM
546 return (ref $z)->make(log($r), $t);
547}
548
549#
0c721ce2
JH
550# ln
551#
552# Alias for log().
553#
554sub ln { Math::Complex::log(@_) }
555
556#
66730be0
RM
557# log10
558#
559# Compute log10(z).
560#
5cd24f17 561
66730be0
RM
562sub log10 {
563 my ($z) = @_;
5cd24f17
PP
564
565 return log(cplx($z, 0)) * log10inv unless ref $z;
66730be0 566 my ($r, $t) = @{$z->polar};
5cd24f17 567 return (ref $z)->make(log($r) * log10inv, $t * log10inv);
66730be0
RM
568}
569
570#
571# logn
572#
573# Compute logn(z,n) = log(z) / log(n)
574#
575sub logn {
576 my ($z, $n) = @_;
0c721ce2 577 $z = cplx($z, 0) unless ref $z;
66730be0
RM
578 my $logn = $logn{$n};
579 $logn = $logn{$n} = log($n) unless defined $logn; # Cache log(n)
0c721ce2 580 return log($z) / $logn;
66730be0
RM
581}
582
583#
584# (cos)
585#
586# Compute cos(z) = (exp(iz) + exp(-iz))/2.
587#
588sub cos {
589 my ($z) = @_;
5cd24f17 590 $z = cplx($z, 0) unless ref $z;
66730be0
RM
591 my ($x, $y) = @{$z->cartesian};
592 my $ey = exp($y);
593 my $ey_1 = 1 / $ey;
0c721ce2
JH
594 return (ref $z)->make(cos($x) * ($ey + $ey_1)/2,
595 sin($x) * ($ey_1 - $ey)/2);
66730be0
RM
596}
597
598#
599# (sin)
600#
601# Compute sin(z) = (exp(iz) - exp(-iz))/2.
602#
603sub sin {
604 my ($z) = @_;
5cd24f17 605 $z = cplx($z, 0) unless ref $z;
66730be0
RM
606 my ($x, $y) = @{$z->cartesian};
607 my $ey = exp($y);
608 my $ey_1 = 1 / $ey;
0c721ce2
JH
609 return (ref $z)->make(sin($x) * ($ey + $ey_1)/2,
610 cos($x) * ($ey - $ey_1)/2);
66730be0
RM
611}
612
613#
614# tan
615#
616# Compute tan(z) = sin(z) / cos(z).
617#
618sub tan {
619 my ($z) = @_;
0c721ce2 620 my $cz = cos($z);
0e505df1 621 _divbyzero "tan($z)", "cos($z)" if ($cz == 0);
0c721ce2 622 return sin($z) / $cz;
66730be0
RM
623}
624
625#
0c721ce2
JH
626# sec
627#
628# Computes the secant sec(z) = 1 / cos(z).
629#
630sub sec {
631 my ($z) = @_;
632 my $cz = cos($z);
0e505df1 633 _divbyzero "sec($z)", "cos($z)" if ($cz == 0);
0c721ce2
JH
634 return 1 / $cz;
635}
636
637#
638# csc
639#
640# Computes the cosecant csc(z) = 1 / sin(z).
641#
642sub csc {
643 my ($z) = @_;
644 my $sz = sin($z);
0e505df1 645 _divbyzero "csc($z)", "sin($z)" if ($sz == 0);
0c721ce2
JH
646 return 1 / $sz;
647}
648
66730be0 649#
0c721ce2 650# cosec
66730be0 651#
0c721ce2
JH
652# Alias for csc().
653#
654sub cosec { Math::Complex::csc(@_) }
655
656#
657# cot
658#
659# Computes cot(z) = 1 / tan(z).
660#
661sub cot {
66730be0 662 my ($z) = @_;
0c721ce2 663 my $sz = sin($z);
0e505df1 664 _divbyzero "cot($z)", "sin($z)" if ($sz == 0);
0c721ce2 665 return cos($z) / $sz;
66730be0
RM
666}
667
668#
0c721ce2
JH
669# cotan
670#
671# Alias for cot().
672#
673sub cotan { Math::Complex::cot(@_) }
674
675#
66730be0
RM
676# acos
677#
678# Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
679#
680sub acos {
681 my ($z) = @_;
0c721ce2 682 $z = cplx($z, 0) unless ref $z;
8c03c583
JH
683 my ($re, $im) = @{$z->cartesian};
684 return atan2(sqrt(1 - $re * $re), $re)
685 if ($im == 0 and abs($re) <= 1.0);
686 my $acos = ~i * log($z + sqrt($z*$z - 1));
687 if ($im == 0 ||
688 (abs($re) < 1 && abs($im) < 1) ||
689 (abs($re) > 1 && abs($im) > 1
690 && !($re > 1 && $im > 1)
691 && !($re < -1 && $im < -1))) {
692 # this rule really, REALLY, must be simpler
693 return -$acos;
694 }
695 return $acos;
66730be0
RM
696}
697
698#
699# asin
700#
701# Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
702#
703sub asin {
704 my ($z) = @_;
0c721ce2 705 $z = cplx($z, 0) unless ref $z;
8c03c583
JH
706 my ($re, $im) = @{$z->cartesian};
707 return atan2($re, sqrt(1 - $re * $re))
708 if ($im == 0 and abs($re) <= 1.0);
0c721ce2 709 return ~i * log(i * $z + sqrt(1 - $z*$z));
66730be0
RM
710}
711
712#
713# atan
714#
0c721ce2 715# Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
66730be0
RM
716#
717sub atan {
718 my ($z) = @_;
5cd24f17 719 $z = cplx($z, 0) unless ref $z;
8c03c583
JH
720 _divbyzero "atan(i)" if ( $z == i);
721 _divbyzero "atan(-i)" if (-$z == i);
0c721ce2 722 return i/2*log((i + $z) / (i - $z));
a0d0e21e
LW
723}
724
66730be0 725#
0c721ce2
JH
726# asec
727#
728# Computes the arc secant asec(z) = acos(1 / z).
729#
730sub asec {
731 my ($z) = @_;
0e505df1 732 _divbyzero "asec($z)", $z if ($z == 0);
8c03c583
JH
733 $z = cplx($z, 0) unless ref $z;
734 my ($re, $im) = @{$z->cartesian};
735 if ($im == 0 && abs($re) >= 1.0) {
736 my $ire = 1 / $re;
737 return atan2(sqrt(1 - $ire * $ire), $ire);
738 }
739 my $asec = acos(1 / $z);
740 return ~$asec if $re < 0 && $re > -1 && $im == 0;
741 return -$asec if $im && !($re > 0 && $im > 0) && !($re < 0 && $im < 0);
742 return $asec;
0c721ce2
JH
743}
744
745#
5cd24f17 746# acsc
0c721ce2 747#
8c03c583 748# Computes the arc cosecant acsc(z) = asin(1 / z).
0c721ce2 749#
5cd24f17 750sub acsc {
0c721ce2 751 my ($z) = @_;
0e505df1 752 _divbyzero "acsc($z)", $z if ($z == 0);
8c03c583
JH
753 $z = cplx($z, 0) unless ref $z;
754 my ($re, $im) = @{$z->cartesian};
755 if ($im == 0 && abs($re) >= 1.0) {
756 my $ire = 1 / $re;
757 return atan2($ire, sqrt(1 - $ire * $ire));
758 }
759 my $acsc = asin(1 / $z);
760 return ~$acsc if $re < 0 && $re > -1 && $im == 0;
761 return $acsc;
0c721ce2
JH
762}
763
764#
5cd24f17 765# acosec
66730be0 766#
5cd24f17 767# Alias for acsc().
0c721ce2 768#
5cd24f17 769sub acosec { Math::Complex::acsc(@_) }
0c721ce2 770
66730be0 771#
0c721ce2
JH
772# acot
773#
8c03c583 774# Computes the arc cotangent acot(z) = atan(1 / z)
0c721ce2
JH
775#
776sub acot {
66730be0 777 my ($z) = @_;
8c03c583 778 _divbyzero "acot($z)" if ($z == 0);
5cd24f17 779 $z = cplx($z, 0) unless ref $z;
8c03c583
JH
780 _divbyzero "acot(i)", if ( $z == i);
781 _divbyzero "acot(-i)" if (-$z == i);
782 return atan(1 / $z);
66730be0
RM
783}
784
785#
0c721ce2
JH
786# acotan
787#
788# Alias for acot().
789#
790sub acotan { Math::Complex::acot(@_) }
791
792#
66730be0
RM
793# cosh
794#
795# Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
796#
797sub cosh {
798 my ($z) = @_;
0e505df1
JH
799 my $real;
800 unless (ref $z) {
801 $z = cplx($z, 0);
802 $real = 1;
803 }
804 my ($x, $y) = @{$z->cartesian};
66730be0
RM
805 my $ex = exp($x);
806 my $ex_1 = 1 / $ex;
f4837644 807 return cplx(0.5 * ($ex + $ex_1), 0) if $real;
0c721ce2
JH
808 return (ref $z)->make(cos($y) * ($ex + $ex_1)/2,
809 sin($y) * ($ex - $ex_1)/2);
66730be0
RM
810}
811
812#
813# sinh
814#
815# Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
816#
817sub sinh {
818 my ($z) = @_;
0e505df1
JH
819 my $real;
820 unless (ref $z) {
821 $z = cplx($z, 0);
822 $real = 1;
823 }
824 my ($x, $y) = @{$z->cartesian};
66730be0
RM
825 my $ex = exp($x);
826 my $ex_1 = 1 / $ex;
f4837644 827 return cplx(0.5 * ($ex - $ex_1), 0) if $real;
0c721ce2
JH
828 return (ref $z)->make(cos($y) * ($ex - $ex_1)/2,
829 sin($y) * ($ex + $ex_1)/2);
66730be0
RM
830}
831
832#
833# tanh
834#
835# Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
836#
837sub tanh {
838 my ($z) = @_;
0c721ce2 839 my $cz = cosh($z);
0e505df1 840 _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
0c721ce2 841 return sinh($z) / $cz;
66730be0
RM
842}
843
844#
0c721ce2
JH
845# sech
846#
847# Computes the hyperbolic secant sech(z) = 1 / cosh(z).
848#
849sub sech {
850 my ($z) = @_;
851 my $cz = cosh($z);
0e505df1 852 _divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
0c721ce2
JH
853 return 1 / $cz;
854}
855
856#
857# csch
858#
859# Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
66730be0 860#
0c721ce2
JH
861sub csch {
862 my ($z) = @_;
863 my $sz = sinh($z);
0e505df1 864 _divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
0c721ce2
JH
865 return 1 / $sz;
866}
867
868#
869# cosech
870#
871# Alias for csch().
872#
873sub cosech { Math::Complex::csch(@_) }
874
66730be0 875#
0c721ce2
JH
876# coth
877#
878# Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
879#
880sub coth {
66730be0 881 my ($z) = @_;
0c721ce2 882 my $sz = sinh($z);
0e505df1 883 _divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
0c721ce2 884 return cosh($z) / $sz;
66730be0
RM
885}
886
887#
0c721ce2
JH
888# cotanh
889#
890# Alias for coth().
891#
892sub cotanh { Math::Complex::coth(@_) }
893
894#
66730be0
RM
895# acosh
896#
8c03c583 897# Computes the arc hyperbolic cosine acosh(z) = log(z +- sqrt(z*z-1)).
66730be0
RM
898#
899sub acosh {
900 my ($z) = @_;
5cd24f17 901 $z = cplx($z, 0) unless ref $z;
8c03c583
JH
902 my ($re, $im) = @{$z->cartesian};
903 return log($re + sqrt(cplx($re*$re - 1, 0)))
904 if ($im == 0 && $re < 0);
0c721ce2 905 return log($z + sqrt($z*$z - 1));
66730be0
RM
906}
907
908#
909# asinh
910#
911# Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1))
912#
913sub asinh {
914 my ($z) = @_;
5cd24f17 915 $z = cplx($z, 0) unless ref $z;
0c721ce2 916 return log($z + sqrt($z*$z + 1));
66730be0
RM
917}
918
919#
920# atanh
921#
922# Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
923#
924sub atanh {
925 my ($z) = @_;
8c03c583
JH
926 _divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
927 _logofzero 'atanh(-1)' if ($z == -1);
5cd24f17 928 $z = cplx($z, 0) unless ref $z;
8c03c583
JH
929 my ($re, $im) = @{$z->cartesian};
930 if ($im == 0 && $re > 1) {
931 return cplx(atanh(1 / $re), pi/2);
932 }
933 return log((1 + $z) / (1 - $z)) / 2;
66730be0
RM
934}
935
936#
0c721ce2
JH
937# asech
938#
939# Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
940#
941sub asech {
942 my ($z) = @_;
0e505df1 943 _divbyzero 'asech(0)', $z if ($z == 0);
8c03c583
JH
944 $z = cplx($z, 0) unless ref $z;
945 my ($re, $im) = @{$z->cartesian};
946 if ($im == 0 && $re < 0) {
947 my $ire = 1 / $re;
948 return log($ire + sqrt(cplx($ire*$ire - 1, 0)));
949 }
0c721ce2
JH
950 return acosh(1 / $z);
951}
952
953#
954# acsch
66730be0 955#
0c721ce2 956# Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
66730be0 957#
0c721ce2 958sub acsch {
66730be0 959 my ($z) = @_;
0e505df1 960 _divbyzero 'acsch(0)', $z if ($z == 0);
0c721ce2
JH
961 return asinh(1 / $z);
962}
963
964#
965# acosech
966#
967# Alias for acosh().
968#
969sub acosech { Math::Complex::acsch(@_) }
970
971#
972# acoth
973#
974# Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
975#
976sub acoth {
977 my ($z) = @_;
8c03c583
JH
978 _divbyzero 'acoth(1)', "$z - 1" if ($z == 1);
979 _logofzero 'acoth(-1)' if ($z == -1);
5cd24f17 980 $z = cplx($z, 0) unless ref $z;
8c03c583
JH
981 my ($re, $im) = @{$z->cartesian};
982 if ($im == 0 and abs($re) < 1) {
983 return cplx(acoth(1/$re) , pi/2);
984 }
985 return log((1 + $z) / ($z - 1)) / 2;
66730be0
RM
986}
987
988#
0c721ce2
JH
989# acotanh
990#
991# Alias for acot().
992#
993sub acotanh { Math::Complex::acoth(@_) }
994
995#
66730be0
RM
996# (atan2)
997#
998# Compute atan(z1/z2).
999#
1000sub atan2 {
1001 my ($z1, $z2, $inverted) = @_;
5cd24f17
PP
1002 my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
1003 my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
66730be0
RM
1004 my $tan;
1005 if (defined $inverted && $inverted) { # atan(z2/z1)
1006 return pi * ($re2 > 0 ? 1 : -1) if $re1 == 0 && $im1 == 0;
1007 $tan = $z2 / $z1;
1008 } else {
1009 return pi * ($re1 > 0 ? 1 : -1) if $re2 == 0 && $im2 == 0;
1010 $tan = $z1 / $z2;
1011 }
1012 return atan($tan);
1013}
1014
1015#
1016# display_format
1017# ->display_format
1018#
1019# Set (fetch if no argument) display format for all complex numbers that
1020# don't happen to have overrriden it via ->display_format
1021#
1022# When called as a method, this actually sets the display format for
1023# the current object.
1024#
1025# Valid object formats are 'c' and 'p' for cartesian and polar. The first
1026# letter is used actually, so the type can be fully spelled out for clarity.
1027#
1028sub display_format {
1029 my $self = shift;
1030 my $format = undef;
1031
1032 if (ref $self) { # Called as a method
1033 $format = shift;
0c721ce2 1034 } else { # Regular procedure call
66730be0
RM
1035 $format = $self;
1036 undef $self;
1037 }
1038
1039 if (defined $self) {
1040 return defined $self->{display} ? $self->{display} : $display
1041 unless defined $format;
1042 return $self->{display} = $format;
1043 }
1044
1045 return $display unless defined $format;
1046 return $display = $format;
1047}
1048
1049#
1050# (stringify)
1051#
1052# Show nicely formatted complex number under its cartesian or polar form,
1053# depending on the current display format:
1054#
1055# . If a specific display format has been recorded for this object, use it.
1056# . Otherwise, use the generic current default for all complex numbers,
1057# which is a package global variable.
1058#
a0d0e21e 1059sub stringify {
66730be0
RM
1060 my ($z) = shift;
1061 my $format;
1062
1063 $format = $display;
1064 $format = $z->{display} if defined $z->{display};
1065
1066 return $z->stringify_polar if $format =~ /^p/i;
1067 return $z->stringify_cartesian;
1068}
1069
1070#
1071# ->stringify_cartesian
1072#
1073# Stringify as a cartesian representation 'a+bi'.
1074#
1075sub stringify_cartesian {
1076 my $z = shift;
1077 my ($x, $y) = @{$z->cartesian};
1078 my ($re, $im);
1079
55497cff
PP
1080 $x = int($x + ($x < 0 ? -1 : 1) * 1e-14)
1081 if int(abs($x)) != int(abs($x) + 1e-14);
1082 $y = int($y + ($y < 0 ? -1 : 1) * 1e-14)
1083 if int(abs($y)) != int(abs($y) + 1e-14);
1084
66730be0
RM
1085 $re = "$x" if abs($x) >= 1e-14;
1086 if ($y == 1) { $im = 'i' }
1087 elsif ($y == -1) { $im = '-i' }
40da2db3 1088 elsif (abs($y) >= 1e-14) { $im = $y . "i" }
66730be0 1089
0c721ce2 1090 my $str = '';
66730be0
RM
1091 $str = $re if defined $re;
1092 $str .= "+$im" if defined $im;
1093 $str =~ s/\+-/-/;
1094 $str =~ s/^\+//;
1095 $str = '0' unless $str;
1096
1097 return $str;
1098}
1099
1100#
1101# ->stringify_polar
1102#
1103# Stringify as a polar representation '[r,t]'.
1104#
1105sub stringify_polar {
1106 my $z = shift;
1107 my ($r, $t) = @{$z->polar};
1108 my $theta;
0c721ce2 1109 my $eps = 1e-14;
66730be0 1110
0c721ce2 1111 return '[0,0]' if $r <= $eps;
a0d0e21e 1112
66730be0
RM
1113 my $tpi = 2 * pi;
1114 my $nt = $t / $tpi;
1115 $nt = ($nt - int($nt)) * $tpi;
1116 $nt += $tpi if $nt < 0; # Range [0, 2pi]
a0d0e21e 1117
0c721ce2
JH
1118 if (abs($nt) <= $eps) { $theta = 0 }
1119 elsif (abs(pi-$nt) <= $eps) { $theta = 'pi' }
66730be0 1120
55497cff 1121 if (defined $theta) {
0c721ce2
JH
1122 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1123 if int(abs($r)) != int(abs($r) + $eps);
1124 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1125 if ($theta ne 'pi' and
1126 int(abs($theta)) != int(abs($theta) + $eps));
55497cff
PP
1127 return "\[$r,$theta\]";
1128 }
66730be0
RM
1129
1130 #
1131 # Okay, number is not a real. Try to identify pi/n and friends...
1132 #
1133
1134 $nt -= $tpi if $nt > pi;
1135 my ($n, $k, $kpi);
1136
1137 for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
1138 $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
0c721ce2
JH
1139 if (abs($kpi/$n - $nt) <= $eps) {
1140 $theta = ($nt < 0 ? '-':'').
1141 ($k == 1 ? 'pi':"${k}pi").'/'.abs($n);
66730be0
RM
1142 last;
1143 }
1144 }
1145
1146 $theta = $nt unless defined $theta;
1147
0c721ce2
JH
1148 $r = int($r + ($r < 0 ? -1 : 1) * $eps)
1149 if int(abs($r)) != int(abs($r) + $eps);
1150 $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
1151 if ($theta !~ m(^-?\d*pi/\d+$) and
1152 int(abs($theta)) != int(abs($theta) + $eps));
55497cff 1153
66730be0 1154 return "\[$r,$theta\]";
a0d0e21e 1155}
a5f75d66
AD
1156
11571;
1158__END__
1159
1160=head1 NAME
1161
66730be0 1162Math::Complex - complex numbers and associated mathematical functions
a5f75d66
AD
1163
1164=head1 SYNOPSIS
1165
66730be0 1166 use Math::Complex;
5aabfad6 1167
66730be0
RM
1168 $z = Math::Complex->make(5, 6);
1169 $t = 4 - 3*i + $z;
1170 $j = cplxe(1, 2*pi/3);
a5f75d66
AD
1171
1172=head1 DESCRIPTION
1173
66730be0
RM
1174This package lets you create and manipulate complex numbers. By default,
1175I<Perl> limits itself to real numbers, but an extra C<use> statement brings
1176full complex support, along with a full set of mathematical functions
1177typically associated with and/or extended to complex numbers.
1178
1179If you wonder what complex numbers are, they were invented to be able to solve
1180the following equation:
1181
1182 x*x = -1
1183
1184and by definition, the solution is noted I<i> (engineers use I<j> instead since
1185I<i> usually denotes an intensity, but the name does not matter). The number
1186I<i> is a pure I<imaginary> number.
1187
1188The arithmetics with pure imaginary numbers works just like you would expect
1189it with real numbers... you just have to remember that
1190
1191 i*i = -1
1192
1193so you have:
1194
1195 5i + 7i = i * (5 + 7) = 12i
1196 4i - 3i = i * (4 - 3) = i
1197 4i * 2i = -8
1198 6i / 2i = 3
1199 1 / i = -i
1200
1201Complex numbers are numbers that have both a real part and an imaginary
1202part, and are usually noted:
1203
1204 a + bi
1205
1206where C<a> is the I<real> part and C<b> is the I<imaginary> part. The
1207arithmetic with complex numbers is straightforward. You have to
1208keep track of the real and the imaginary parts, but otherwise the
1209rules used for real numbers just apply:
1210
1211 (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i
1212 (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i
1213
1214A graphical representation of complex numbers is possible in a plane
1215(also called the I<complex plane>, but it's really a 2D plane).
1216The number
1217
1218 z = a + bi
1219
1220is the point whose coordinates are (a, b). Actually, it would
1221be the vector originating from (0, 0) to (a, b). It follows that the addition
1222of two complex numbers is a vectorial addition.
1223
1224Since there is a bijection between a point in the 2D plane and a complex
1225number (i.e. the mapping is unique and reciprocal), a complex number
1226can also be uniquely identified with polar coordinates:
1227
1228 [rho, theta]
1229
1230where C<rho> is the distance to the origin, and C<theta> the angle between
1231the vector and the I<x> axis. There is a notation for this using the
1232exponential form, which is:
1233
1234 rho * exp(i * theta)
1235
1236where I<i> is the famous imaginary number introduced above. Conversion
1237between this form and the cartesian form C<a + bi> is immediate:
1238
1239 a = rho * cos(theta)
1240 b = rho * sin(theta)
1241
1242which is also expressed by this formula:
1243
1244 z = rho * exp(i * theta) = rho * (cos theta + i * sin theta)
1245
1246In other words, it's the projection of the vector onto the I<x> and I<y>
1247axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta>
1248the I<argument> of the complex number. The I<norm> of C<z> will be
1249noted C<abs(z)>.
1250
1251The polar notation (also known as the trigonometric
1252representation) is much more handy for performing multiplications and
1253divisions of complex numbers, whilst the cartesian notation is better
1254suited for additions and substractions. Real numbers are on the I<x>
1255axis, and therefore I<theta> is zero.
1256
1257All the common operations that can be performed on a real number have
1258been defined to work on complex numbers as well, and are merely
1259I<extensions> of the operations defined on real numbers. This means
1260they keep their natural meaning when there is no imaginary part, provided
1261the number is within their definition set.
1262
1263For instance, the C<sqrt> routine which computes the square root of
1264its argument is only defined for positive real numbers and yields a
1265positive real number (it is an application from B<R+> to B<R+>).
1266If we allow it to return a complex number, then it can be extended to
1267negative real numbers to become an application from B<R> to B<C> (the
1268set of complex numbers):
1269
1270 sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i
1271
1272It can also be extended to be an application from B<C> to B<C>,
1273whilst its restriction to B<R> behaves as defined above by using
1274the following definition:
1275
1276 sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2)
1277
1278Indeed, a negative real number can be noted C<[x,pi]>
1279(the modulus I<x> is always positive, so C<[x,pi]> is really C<-x>, a
1280negative number)
1281and the above definition states that
1282
1283 sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i
1284
1285which is exactly what we had defined for negative real numbers above.
a5f75d66 1286
66730be0
RM
1287All the common mathematical functions defined on real numbers that
1288are extended to complex numbers share that same property of working
1289I<as usual> when the imaginary part is zero (otherwise, it would not
1290be called an extension, would it?).
a5f75d66 1291
66730be0
RM
1292A I<new> operation possible on a complex number that is
1293the identity for real numbers is called the I<conjugate>, and is noted
1294with an horizontal bar above the number, or C<~z> here.
a5f75d66 1295
66730be0
RM
1296 z = a + bi
1297 ~z = a - bi
a5f75d66 1298
66730be0 1299Simple... Now look:
a5f75d66 1300
66730be0 1301 z * ~z = (a + bi) * (a - bi) = a*a + b*b
a5f75d66 1302
66730be0
RM
1303We saw that the norm of C<z> was noted C<abs(z)> and was defined as the
1304distance to the origin, also known as:
a5f75d66 1305
66730be0 1306 rho = abs(z) = sqrt(a*a + b*b)
a5f75d66 1307
66730be0
RM
1308so
1309
1310 z * ~z = abs(z) ** 2
1311
1312If z is a pure real number (i.e. C<b == 0>), then the above yields:
1313
1314 a * a = abs(a) ** 2
1315
1316which is true (C<abs> has the regular meaning for real number, i.e. stands
1317for the absolute value). This example explains why the norm of C<z> is
1318noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet
1319is the regular C<abs> we know when the complex number actually has no
1320imaginary part... This justifies I<a posteriori> our use of the C<abs>
1321notation for the norm.
1322
1323=head1 OPERATIONS
1324
1325Given the following notations:
1326
1327 z1 = a + bi = r1 * exp(i * t1)
1328 z2 = c + di = r2 * exp(i * t2)
1329 z = <any complex or real number>
1330
1331the following (overloaded) operations are supported on complex numbers:
1332
1333 z1 + z2 = (a + c) + i(b + d)
1334 z1 - z2 = (a - c) + i(b - d)
1335 z1 * z2 = (r1 * r2) * exp(i * (t1 + t2))
1336 z1 / z2 = (r1 / r2) * exp(i * (t1 - t2))
1337 z1 ** z2 = exp(z2 * log z1)
1338 ~z1 = a - bi
1339 abs(z1) = r1 = sqrt(a*a + b*b)
1340 sqrt(z1) = sqrt(r1) * exp(i * t1/2)
1341 exp(z1) = exp(a) * exp(i * b)
1342 log(z1) = log(r1) + i*t1
1343 sin(z1) = 1/2i (exp(i * z1) - exp(-i * z1))
1344 cos(z1) = 1/2 (exp(i * z1) + exp(-i * z1))
1345 abs(z1) = r1
1346 atan2(z1, z2) = atan(z1/z2)
1347
1348The following extra operations are supported on both real and complex
1349numbers:
1350
1351 Re(z) = a
1352 Im(z) = b
1353 arg(z) = t
1354
1355 cbrt(z) = z ** (1/3)
1356 log10(z) = log(z) / log(10)
1357 logn(z, n) = log(z) / log(n)
1358
1359 tan(z) = sin(z) / cos(z)
0c721ce2 1360
5aabfad6
PP
1361 csc(z) = 1 / sin(z)
1362 sec(z) = 1 / cos(z)
0c721ce2 1363 cot(z) = 1 / tan(z)
66730be0
RM
1364
1365 asin(z) = -i * log(i*z + sqrt(1-z*z))
1366 acos(z) = -i * log(z + sqrt(z*z-1))
1367 atan(z) = i/2 * log((i+z) / (i-z))
0c721ce2 1368
5aabfad6
PP
1369 acsc(z) = asin(1 / z)
1370 asec(z) = acos(1 / z)
8c03c583 1371 acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i))
66730be0
RM
1372
1373 sinh(z) = 1/2 (exp(z) - exp(-z))
1374 cosh(z) = 1/2 (exp(z) + exp(-z))
0c721ce2
JH
1375 tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
1376
5aabfad6
PP
1377 csch(z) = 1 / sinh(z)
1378 sech(z) = 1 / cosh(z)
0c721ce2 1379 coth(z) = 1 / tanh(z)
66730be0
RM
1380
1381 asinh(z) = log(z + sqrt(z*z+1))
1382 acosh(z) = log(z + sqrt(z*z-1))
1383 atanh(z) = 1/2 * log((1+z) / (1-z))
66730be0 1384
5aabfad6
PP
1385 acsch(z) = asinh(1 / z)
1386 asech(z) = acosh(1 / z)
0c721ce2
JH
1387 acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
1388
1389I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>, I<coth>,
1390I<acosech>, I<acotanh>, have aliases I<ln>, I<cosec>, I<cotan>,
1391I<acosec>, I<acotan>, I<cosech>, I<cotanh>, I<acosech>, I<acotanh>,
1392respectively.
1393
1394The I<root> function is available to compute all the I<n>
66730be0
RM
1395roots of some complex, where I<n> is a strictly positive integer.
1396There are exactly I<n> such roots, returned as a list. Getting the
1397number mathematicians call C<j> such that:
1398
1399 1 + j + j*j = 0;
1400
1401is a simple matter of writing:
1402
1403 $j = ((root(1, 3))[1];
1404
1405The I<k>th root for C<z = [r,t]> is given by:
1406
1407 (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
1408
f4837644
JH
1409The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
1410order to ensure its restriction to real numbers is conform to what you
1411would expect, the comparison is run on the real part of the complex
1412number first, and imaginary parts are compared only when the real
1413parts match.
66730be0
RM
1414
1415=head1 CREATION
1416
1417To create a complex number, use either:
1418
1419 $z = Math::Complex->make(3, 4);
1420 $z = cplx(3, 4);
1421
1422if you know the cartesian form of the number, or
1423
1424 $z = 3 + 4*i;
1425
1426if you like. To create a number using the trigonometric form, use either:
1427
1428 $z = Math::Complex->emake(5, pi/3);
1429 $x = cplxe(5, pi/3);
1430
0c721ce2
JH
1431instead. The first argument is the modulus, the second is the angle
1432(in radians, the full circle is 2*pi). (Mnmemonic: C<e> is used as a
1433notation for complex numbers in the trigonometric form).
66730be0
RM
1434
1435It is possible to write:
1436
1437 $x = cplxe(-3, pi/4);
1438
1439but that will be silently converted into C<[3,-3pi/4]>, since the modulus
1440must be positive (it represents the distance to the origin in the complex
1441plane).
1442
1443=head1 STRINGIFICATION
1444
1445When printed, a complex number is usually shown under its cartesian
1446form I<a+bi>, but there are legitimate cases where the polar format
1447I<[r,t]> is more appropriate.
1448
1449By calling the routine C<Math::Complex::display_format> and supplying either
1450C<"polar"> or C<"cartesian">, you override the default display format,
1451which is C<"cartesian">. Not supplying any argument returns the current
1452setting.
1453
1454This default can be overridden on a per-number basis by calling the
1455C<display_format> method instead. As before, not supplying any argument
1456returns the current display format for this number. Otherwise whatever you
1457specify will be the new display format for I<this> particular number.
1458
1459For instance:
1460
1461 use Math::Complex;
1462
1463 Math::Complex::display_format('polar');
1464 $j = ((root(1, 3))[1];
1465 print "j = $j\n"; # Prints "j = [1,2pi/3]
1466 $j->display_format('cartesian');
1467 print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i"
1468
1469The polar format attempts to emphasize arguments like I<k*pi/n>
1470(where I<n> is a positive integer and I<k> an integer within [-9,+9]).
1471
1472=head1 USAGE
1473
1474Thanks to overloading, the handling of arithmetics with complex numbers
1475is simple and almost transparent.
1476
1477Here are some examples:
1478
1479 use Math::Complex;
1480
1481 $j = cplxe(1, 2*pi/3); # $j ** 3 == 1
1482 print "j = $j, j**3 = ", $j ** 3, "\n";
1483 print "1 + j + j**2 = ", 1 + $j + $j**2, "\n";
1484
1485 $z = -16 + 0*i; # Force it to be a complex
1486 print "sqrt($z) = ", sqrt($z), "\n";
1487
1488 $k = exp(i * 2*pi/3);
1489 print "$j - $k = ", $j - $k, "\n";
a5f75d66 1490
5cd24f17 1491=head1 ERRORS DUE TO DIVISION BY ZERO
5aabfad6
PP
1492
1493The division (/) and the following functions
1494
1495 tan
1496 sec
1497 csc
1498 cot
5cd24f17
PP
1499 asec
1500 acsc
5aabfad6
PP
1501 atan
1502 acot
1503 tanh
1504 sech
1505 csch
1506 coth
1507 atanh
1508 asech
1509 acsch
1510 acoth
1511
1512cannot be computed for all arguments because that would mean dividing
8c03c583
JH
1513by zero or taking logarithm of zero. These situations cause fatal
1514runtime errors looking like this
5aabfad6
PP
1515
1516 cot(0): Division by zero.
5cd24f17 1517 (Because in the definition of cot(0), the divisor sin(0) is 0)
5aabfad6
PP
1518 Died at ...
1519
8c03c583
JH
1520or
1521
1522 atanh(-1): Logarithm of zero.
1523 Died at...
1524
1525For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
1526C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the
1527C<atanh>, C<acoth>, the argument cannot be C<1> (one). For the
1528C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one). For the
1529C<atan>, C<acot>, the argument cannot be C<i> (the imaginary unit).
1530For the C<atan>, C<acoth>, the argument cannot be C<-i> (the negative
1531imaginary unit). For the C<tan>, C<sec>, C<tanh>, C<sech>, the
1532argument cannot be I<pi/2 + k * pi>, where I<k> is any integer.
5cd24f17 1533
a5f75d66
AD
1534=head1 BUGS
1535
5cd24f17
PP
1536Saying C<use Math::Complex;> exports many mathematical routines in the
1537caller environment and even overrides some (C<sin>, C<cos>, C<sqrt>,
1538C<log>, C<exp>). This is construed as a feature by the Authors,
1539actually... ;-)
66730be0
RM
1540
1541The code is not optimized for speed, although we try to use the cartesian
1542form for addition-like operators and the trigonometric form for all
1543multiplication-like operators.
1544
1545The arg() routine does not ensure the angle is within the range [-pi,+pi]
1546(a side effect caused by multiplication and division using the trigonometric
1547representation).
a5f75d66 1548
66730be0
RM
1549All routines expect to be given real or complex numbers. Don't attempt to
1550use BigFloat, since Perl has currently no rule to disambiguate a '+'
1551operation (for instance) between two overloaded entities.
a5f75d66 1552
0c721ce2 1553=head1 AUTHORS
a5f75d66 1554
ace5de91
GS
1555Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>> and
1556Jarkko Hietaniemi <F<jhi@iki.fi>>.
5cd24f17
PP
1557
1558=cut
1559
1560# eof