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