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