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