← Index
Performance Profile   « block view • line view • sub view »
For /wise/base/deliv/dev/bin/getfix
  Run on Thu May 20 15:30:03 2010
Reported on Thu May 20 16:25:29 2010

File/opt/wise/lib/perl5/5.10.0/Math/Complex.pm
Statements Executed40
Total Time0.007931 seconds

Subroutines — ordered by exclusive time
Calls P F Exclusive
Time		 Inclusive
Time		 Subroutine
1114.6e-54.6e-5Math::Complex::atan2
1111.3e-51.3e-5Math::Complex::sqrt
00000Math::Complex::BEGIN
00000Math::Complex::Im
00000Math::Complex::Re
00000Math::Complex::_cannot_make
00000Math::Complex::_cartesian
00000Math::Complex::_conjugate
00000Math::Complex::_divbyzero
00000Math::Complex::_divide
00000Math::Complex::_emake
00000Math::Complex::_ip2
00000Math::Complex::_logofzero
00000Math::Complex::_make
00000Math::Complex::_minus
00000Math::Complex::_multiply
00000Math::Complex::_negate
00000Math::Complex::_numeq
00000Math::Complex::_plus
00000Math::Complex::_polar
00000Math::Complex::_power
00000Math::Complex::_rootbad
00000Math::Complex::_set_cartesian
00000Math::Complex::_set_polar
00000Math::Complex::_spaceship
00000Math::Complex::_stringify
00000Math::Complex::_stringify_cartesian
00000Math::Complex::_stringify_polar
00000Math::Complex::_theta
00000Math::Complex::_update_cartesian
00000Math::Complex::_update_polar
00000Math::Complex::abs
00000Math::Complex::acos
00000Math::Complex::acosec
00000Math::Complex::acosech
00000Math::Complex::acosh
00000Math::Complex::acot
00000Math::Complex::acotan
00000Math::Complex::acotanh
00000Math::Complex::acoth
00000Math::Complex::acsc
00000Math::Complex::acsch
00000Math::Complex::arg
00000Math::Complex::asec
00000Math::Complex::asech
00000Math::Complex::asin
00000Math::Complex::asinh
00000Math::Complex::atan
00000Math::Complex::atanh
00000Math::Complex::cbrt
00000Math::Complex::cos
00000Math::Complex::cosec
00000Math::Complex::cosech
00000Math::Complex::cosh
00000Math::Complex::cot
00000Math::Complex::cotan
00000Math::Complex::cotanh
00000Math::Complex::coth
00000Math::Complex::cplx
00000Math::Complex::cplxe
00000Math::Complex::csc
00000Math::Complex::csch
00000Math::Complex::display_format
00000Math::Complex::emake
00000Math::Complex::exp
00000Math::Complex::i
00000Math::Complex::ln
00000Math::Complex::log
00000Math::Complex::log10
00000Math::Complex::logn
00000Math::Complex::make
00000Math::Complex::new
00000Math::Complex::rho
00000Math::Complex::root
00000Math::Complex::sec
00000Math::Complex::sech
00000Math::Complex::sin
00000Math::Complex::sinh
00000Math::Complex::tan
00000Math::Complex::tanh
00000Math::Complex::theta

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

Report produced by the NYTProf 2.07 Perl profiler, developed by Tim Bunce and Adam Kaplan.