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