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