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