File:
[LON-CAPA] /
loncom /
xml /
LCMathComplex.pm
Revision
1.1:
download - view:
text,
annotated -
select for diffs
Tue Oct 29 21:01:21 2013 UTC (10 years, 6 months ago) by
raeburn
Branches:
MAIN
CVS tags:
version_2_11_2_uiuc,
version_2_11_2_msu,
version_2_11_2_educog,
version_2_11_2,
version_2_11_1,
version_2_11_0_RC3,
version_2_11_0_RC2,
version_2_11_0,
HEAD
- Bug 6629.
- Re-enable access to functionality in Math::Complex in Safe Space by
copying Math::Complex 1.55 to LONCAPA::LCMathComplex and eliminating use
of Config.pm
- Whenever ./UPDATE is run to install or update LON-CAPA, the code which
sets $nvsize in the standard Math::Complex script will be run in
lcmathcomplex.piml and the value of $nvsize will be set to the
appropriate value: 4, 8, 10, 12 or 16 in LONCAPA/LCMathComplex.pm.
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