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