diff options
author | Karl Williamson <khw@cpan.org> | 2022-09-30 21:28:19 -0600 |
---|---|---|
committer | James E Keenan <jkeenan@cpan.org> | 2022-10-01 09:36:13 -0400 |
commit | 74faf852d82b5401b5159457e3c856b5c314a785 (patch) | |
tree | f56a581502c27cfa0251d6162bfd70685b622c98 /dist/Math-Complex | |
parent | 8a234390ad95e31a35481824eece8519e7407c81 (diff) | |
download | perl-74faf852d82b5401b5159457e3c856b5c314a785.tar.gz |
Move Math-Complex from cpan/ to dist/
This module is now being maintained by p5p now.
Diffstat (limited to 'dist/Math-Complex')
-rw-r--r-- | dist/Math-Complex/lib/Math/Complex.pm | 2132 | ||||
-rw-r--r-- | dist/Math-Complex/lib/Math/Trig.pm | 761 | ||||
-rw-r--r-- | dist/Math-Complex/t/Complex.t | 1160 | ||||
-rw-r--r-- | dist/Math-Complex/t/Trig.t | 387 | ||||
-rw-r--r-- | dist/Math-Complex/t/underbar.t | 28 |
5 files changed, 4468 insertions, 0 deletions
diff --git a/dist/Math-Complex/lib/Math/Complex.pm b/dist/Math-Complex/lib/Math/Complex.pm new file mode 100644 index 0000000000..6cab2689bd --- /dev/null +++ b/dist/Math-Complex/lib/Math/Complex.pm @@ -0,0 +1,2132 @@ +# +# Complex numbers and associated mathematical functions +# -- Raphael Manfredi Since Sep 1996 +# -- Jarkko Hietaniemi Since Mar 1997 +# -- Daniel S. Lewart Since Sep 1997 +# + +package Math::Complex; + +{ use 5.006; } +use strict; + +our $VERSION = 1.59_02; + +use Config; + +our ($Inf, $ExpInf); +our ($vax_float, $has_inf, $has_nan); + +BEGIN { + $vax_float = (pack("d",1) =~ /^[\x80\x10]\x40/); + $has_inf = !$vax_float; + $has_nan = !$vax_float; + + unless ($has_inf) { + # For example in vax, there is no Inf, + # and just mentioning the DBL_MAX (1.70141183460469229e+38) + # causes SIGFPE. + + # These are pretty useless without a real infinity, + # but setting them makes for less warnings about their + # undefined values. + $Inf = "Inf"; + $ExpInf = "Inf"; + return; + } + + my %DBL_MAX = # These are IEEE 754 maxima. + ( + 4 => '1.70141183460469229e+38', + 8 => '1.7976931348623157e+308', + # AFAICT the 10, 12, and 16-byte long doubles + # all have the same maximum. + 10 => '1.1897314953572317650857593266280070162E+4932', + 12 => '1.1897314953572317650857593266280070162E+4932', + 16 => '1.1897314953572317650857593266280070162E+4932', + ); + + my $nvsize = $Config{nvsize} || + ($Config{uselongdouble} && $Config{longdblsize}) || + $Config{doublesize}; + die "Math::Complex: Could not figure out nvsize\n" + unless defined $nvsize; + die "Math::Complex: Cannot not figure out max nv (nvsize = $nvsize)\n" + unless defined $DBL_MAX{$nvsize}; + my $DBL_MAX = eval $DBL_MAX{$nvsize}; + die "Math::Complex: Could not figure out max nv (nvsize = $nvsize)\n" + unless defined $DBL_MAX; + my $BIGGER_THAN_THIS = 1e30; # Must find something bigger than this. + if ($^O eq 'unicosmk') { + $Inf = $DBL_MAX; + } else { + local $SIG{FPE} = sub { }; + local $!; + # We do want an arithmetic overflow, Inf INF inf Infinity. + for my $t ( + 'exp(99999)', # Enough even with 128-bit long doubles. + 'inf', + 'Inf', + 'INF', + 'infinity', + 'Infinity', + 'INFINITY', + '1e99999', + ) { + local $^W = 0; + my $i = eval "$t+1.0"; + if (defined $i && $i > $BIGGER_THAN_THIS) { + $Inf = $i; + last; + } + } + $Inf = $DBL_MAX unless defined $Inf; # Oh well, close enough. + die "Math::Complex: Could not get Infinity" + unless $Inf > $BIGGER_THAN_THIS; + $ExpInf = eval 'exp(99999)'; + } + # print "# On this machine, Inf = '$Inf'\n"; +} + +use Scalar::Util qw(set_prototype); + +use warnings; +no warnings 'syntax'; # To avoid the (_) warnings. + +BEGIN { + # For certain functions that we override, in 5.10 or better + # we can set a smarter prototype that will handle the lexical $_ + # (also a 5.10+ feature). + if ($] >= 5.010000) { + set_prototype \&abs, '_'; + set_prototype \&cos, '_'; + set_prototype \&exp, '_'; + set_prototype \&log, '_'; + set_prototype \&sin, '_'; + set_prototype \&sqrt, '_'; + } +} + +my $i; +my %LOGN; + +# Regular expression for floating point numbers. +# These days we could use Scalar::Util::lln(), I guess. +my $gre = qr'\s*([\+\-]?(?:(?:(?:\d+(?:_\d+)*(?:\.\d*(?:_\d+)*)?|\.\d+(?:_\d+)*)(?:[eE][\+\-]?\d+(?:_\d+)*)?))|inf)'i; + +require Exporter; + +our @ISA = qw(Exporter); + +my @trig = qw( + pi + tan + csc cosec sec cot cotan + asin acos atan + acsc acosec asec acot acotan + sinh cosh tanh + csch cosech sech coth cotanh + asinh acosh atanh + acsch acosech asech acoth acotanh + ); + +our @EXPORT = (qw( + i Re Im rho theta arg + sqrt log ln + log10 logn cbrt root + cplx cplxe + atan2 + ), + @trig); + +my @pi = qw(pi pi2 pi4 pip2 pip4 Inf); + +our @EXPORT_OK = @pi; + +our %EXPORT_TAGS = ( + 'trig' => [@trig], + 'pi' => [@pi], +); + +use overload + '=' => \&_copy, + '+=' => \&_plus, + '+' => \&_plus, + '-=' => \&_minus, + '-' => \&_minus, + '*=' => \&_multiply, + '*' => \&_multiply, + '/=' => \&_divide, + '/' => \&_divide, + '**=' => \&_power, + '**' => \&_power, + '==' => \&_numeq, + '<=>' => \&_spaceship, + 'neg' => \&_negate, + '~' => \&_conjugate, + 'abs' => \&abs, + 'sqrt' => \&sqrt, + 'exp' => \&exp, + 'log' => \&log, + 'sin' => \&sin, + 'cos' => \&cos, + 'atan2' => \&atan2, + '""' => \&_stringify; + +# +# Package "privates" +# + +my %DISPLAY_FORMAT = ('style' => 'cartesian', + 'polar_pretty_print' => 1); +my $eps = 1e-14; # Epsilon + +# +# Object attributes (internal): +# cartesian [real, imaginary] -- cartesian form +# polar [rho, theta] -- polar form +# c_dirty cartesian form not up-to-date +# p_dirty polar form not up-to-date +# display display format (package's global when not set) +# + +# Die on bad *make() arguments. + +sub _cannot_make { + die "@{[(caller(1))[3]]}: Cannot take $_[0] of '$_[1]'.\n"; +} + +sub _make { + my $arg = shift; + my ($p, $q); + + if ($arg =~ /^$gre$/) { + ($p, $q) = ($1, 0); + } elsif ($arg =~ /^(?:$gre)?$gre\s*i\s*$/) { + ($p, $q) = ($1 || 0, $2); + } elsif ($arg =~ /^\s*\(\s*$gre\s*(?:,\s*$gre\s*)?\)\s*$/) { + ($p, $q) = ($1, $2 || 0); + } + + if (defined $p) { + $p =~ s/^\+//; + $p =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf; + $q =~ s/^\+//; + $q =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf; + } + + return ($p, $q); +} + +sub _emake { + my $arg = shift; + my ($p, $q); + + if ($arg =~ /^\s*\[\s*$gre\s*(?:,\s*$gre\s*)?\]\s*$/) { + ($p, $q) = ($1, $2 || 0); + } elsif ($arg =~ m!^\s*\[\s*$gre\s*(?:,\s*([-+]?\d*\s*)?pi(?:/\s*(\d+))?\s*)?\]\s*$!) { + ($p, $q) = ($1, ($2 eq '-' ? -1 : ($2 || 1)) * pi() / ($3 || 1)); + } elsif ($arg =~ /^\s*\[\s*$gre\s*\]\s*$/) { + ($p, $q) = ($1, 0); + } elsif ($arg =~ /^\s*$gre\s*$/) { + ($p, $q) = ($1, 0); + } + + if (defined $p) { + $p =~ s/^\+//; + $q =~ s/^\+//; + $p =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf; + $q =~ s/^(-?)inf$/"${1}9**9**9"/e if $has_inf; + } + + return ($p, $q); +} + +sub _copy { + my $self = shift; + my $clone = {%$self}; + if ($self->{'cartesian'}) { + $clone->{'cartesian'} = [@{$self->{'cartesian'}}]; + } + if ($self->{'polar'}) { + $clone->{'polar'} = [@{$self->{'polar'}}]; + } + bless $clone,__PACKAGE__; + return $clone; +} + +# +# ->make +# +# Create a new complex number (cartesian form) +# +sub make { + my $self = bless {}, shift; + my ($re, $im); + if (@_ == 0) { + ($re, $im) = (0, 0); + } elsif (@_ == 1) { + return (ref $self)->emake($_[0]) + if ($_[0] =~ /^\s*\[/); + ($re, $im) = _make($_[0]); + } elsif (@_ == 2) { + ($re, $im) = @_; + } + if (defined $re) { + _cannot_make("real part", $re) unless $re =~ /^$gre$/; + } + $im ||= 0; + _cannot_make("imaginary part", $im) unless $im =~ /^$gre$/; + $self->_set_cartesian([$re, $im ]); + $self->display_format('cartesian'); + + return $self; +} + +# +# ->emake +# +# Create a new complex number (exponential form) +# +sub emake { + my $self = bless {}, shift; + my ($rho, $theta); + if (@_ == 0) { + ($rho, $theta) = (0, 0); + } elsif (@_ == 1) { + return (ref $self)->make($_[0]) + if ($_[0] =~ /^\s*\(/ || $_[0] =~ /i\s*$/); + ($rho, $theta) = _emake($_[0]); + } elsif (@_ == 2) { + ($rho, $theta) = @_; + } + if (defined $rho && defined $theta) { + if ($rho < 0) { + $rho = -$rho; + $theta = ($theta <= 0) ? $theta + pi() : $theta - pi(); + } + } + if (defined $rho) { + _cannot_make("rho", $rho) unless $rho =~ /^$gre$/; + } + $theta ||= 0; + _cannot_make("theta", $theta) unless $theta =~ /^$gre$/; + $self->_set_polar([$rho, $theta]); + $self->display_format('polar'); + + return $self; +} + +sub new { &make } # For backward compatibility only. + +# +# cplx +# +# Creates a complex number from a (re, im) tuple. +# This avoids the burden of writing Math::Complex->make(re, im). +# +sub cplx { + return __PACKAGE__->make(@_); +} + +# +# cplxe +# +# Creates a complex number from a (rho, theta) tuple. +# This avoids the burden of writing Math::Complex->emake(rho, theta). +# +sub cplxe { + return __PACKAGE__->emake(@_); +} + +# +# pi +# +# The number defined as pi = 180 degrees +# +sub pi () { 4 * CORE::atan2(1, 1) } + +# +# pi2 +# +# The full circle +# +sub pi2 () { 2 * pi } + +# +# pi4 +# +# The full circle twice. +# +sub pi4 () { 4 * pi } + +# +# pip2 +# +# The quarter circle +# +sub pip2 () { pi / 2 } + +# +# pip4 +# +# The eighth circle. +# +sub pip4 () { pi / 4 } + +# +# _uplog10 +# +# Used in log10(). +# +sub _uplog10 () { 1 / CORE::log(10) } + +# +# i +# +# The number defined as i*i = -1; +# +sub i () { + return $i if ($i); + $i = bless {}; + $i->{'cartesian'} = [0, 1]; + $i->{'polar'} = [1, pip2]; + $i->{c_dirty} = 0; + $i->{p_dirty} = 0; + return $i; +} + +# +# _ip2 +# +# Half of i. +# +sub _ip2 () { i / 2 } + +# +# Attribute access/set routines +# + +sub _cartesian {$_[0]->{c_dirty} ? + $_[0]->_update_cartesian : $_[0]->{'cartesian'}} +sub _polar {$_[0]->{p_dirty} ? + $_[0]->_update_polar : $_[0]->{'polar'}} + +sub _set_cartesian { $_[0]->{p_dirty}++; $_[0]->{c_dirty} = 0; + $_[0]->{'cartesian'} = $_[1] } +sub _set_polar { $_[0]->{c_dirty}++; $_[0]->{p_dirty} = 0; + $_[0]->{'polar'} = $_[1] } + +# +# ->_update_cartesian +# +# Recompute and return the cartesian form, given accurate polar form. +# +sub _update_cartesian { + my $self = shift; + my ($r, $t) = @{$self->{'polar'}}; + $self->{c_dirty} = 0; + return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)]; +} + +# +# +# ->_update_polar +# +# Recompute and return the polar form, given accurate cartesian form. +# +sub _update_polar { + my $self = shift; + my ($x, $y) = @{$self->{'cartesian'}}; + $self->{p_dirty} = 0; + return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0; + return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y), + CORE::atan2($y, $x)]; +} + +# +# (_plus) +# +# Computes z1+z2. +# +sub _plus { + my ($z1, $z2, $regular) = @_; + my ($re1, $im1) = @{$z1->_cartesian}; + $z2 = cplx($z2) unless ref $z2; + my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0); + unless (defined $regular) { + $z1->_set_cartesian([$re1 + $re2, $im1 + $im2]); + return $z1; + } + return (ref $z1)->make($re1 + $re2, $im1 + $im2); +} + +# +# (_minus) +# +# Computes z1-z2. +# +sub _minus { + my ($z1, $z2, $inverted) = @_; + my ($re1, $im1) = @{$z1->_cartesian}; + $z2 = cplx($z2) unless ref $z2; + my ($re2, $im2) = @{$z2->_cartesian}; + unless (defined $inverted) { + $z1->_set_cartesian([$re1 - $re2, $im1 - $im2]); + return $z1; + } + return $inverted ? + (ref $z1)->make($re2 - $re1, $im2 - $im1) : + (ref $z1)->make($re1 - $re2, $im1 - $im2); + +} + +# +# (_multiply) +# +# Computes z1*z2. +# +sub _multiply { + my ($z1, $z2, $regular) = @_; + if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) { + # if both polar better use polar to avoid rounding errors + my ($r1, $t1) = @{$z1->_polar}; + my ($r2, $t2) = @{$z2->_polar}; + my $t = $t1 + $t2; + if ($t > pi()) { $t -= pi2 } + elsif ($t <= -pi()) { $t += pi2 } + unless (defined $regular) { + $z1->_set_polar([$r1 * $r2, $t]); + return $z1; + } + return (ref $z1)->emake($r1 * $r2, $t); + } else { + my ($x1, $y1) = @{$z1->_cartesian}; + if (ref $z2) { + my ($x2, $y2) = @{$z2->_cartesian}; + return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2); + } else { + return (ref $z1)->make($x1*$z2, $y1*$z2); + } + } +} + +# +# _divbyzero +# +# Die on division by zero. +# +sub _divbyzero { + my $mess = "$_[0]: Division by zero.\n"; + + if (defined $_[1]) { + $mess .= "(Because in the definition of $_[0], the divisor "; + $mess .= "$_[1] " unless ("$_[1]" eq '0'); + $mess .= "is 0)\n"; + } + + my @up = caller(1); + + $mess .= "Died at $up[1] line $up[2].\n"; + + die $mess; +} + +# +# (_divide) +# +# Computes z1/z2. +# +sub _divide { + my ($z1, $z2, $inverted) = @_; + if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) { + # if both polar better use polar to avoid rounding errors + my ($r1, $t1) = @{$z1->_polar}; + my ($r2, $t2) = @{$z2->_polar}; + my $t; + if ($inverted) { + _divbyzero "$z2/0" if ($r1 == 0); + $t = $t2 - $t1; + if ($t > pi()) { $t -= pi2 } + elsif ($t <= -pi()) { $t += pi2 } + return (ref $z1)->emake($r2 / $r1, $t); + } else { + _divbyzero "$z1/0" if ($r2 == 0); + $t = $t1 - $t2; + if ($t > pi()) { $t -= pi2 } + elsif ($t <= -pi()) { $t += pi2 } + return (ref $z1)->emake($r1 / $r2, $t); + } + } else { + my ($d, $x2, $y2); + if ($inverted) { + ($x2, $y2) = @{$z1->_cartesian}; + $d = $x2*$x2 + $y2*$y2; + _divbyzero "$z2/0" if $d == 0; + return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d); + } else { + my ($x1, $y1) = @{$z1->_cartesian}; + if (ref $z2) { + ($x2, $y2) = @{$z2->_cartesian}; + $d = $x2*$x2 + $y2*$y2; + _divbyzero "$z1/0" if $d == 0; + my $u = ($x1*$x2 + $y1*$y2)/$d; + my $v = ($y1*$x2 - $x1*$y2)/$d; + return (ref $z1)->make($u, $v); + } else { + _divbyzero "$z1/0" if $z2 == 0; + return (ref $z1)->make($x1/$z2, $y1/$z2); + } + } + } +} + +# +# (_power) +# +# Computes z1**z2 = exp(z2 * log z1)). +# +sub _power { + my ($z1, $z2, $inverted) = @_; + if ($inverted) { + return 1 if $z1 == 0 || $z2 == 1; + return 0 if $z2 == 0 && Re($z1) > 0; + } else { + return 1 if $z2 == 0 || $z1 == 1; + return 0 if $z1 == 0 && Re($z2) > 0; + } + my $w = $inverted ? &exp($z1 * &log($z2)) + : &exp($z2 * &log($z1)); + # If both arguments cartesian, return cartesian, else polar. + return $z1->{c_dirty} == 0 && + (not ref $z2 or $z2->{c_dirty} == 0) ? + cplx(@{$w->_cartesian}) : $w; +} + +# +# (_spaceship) +# +# Computes z1 <=> z2. +# Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i. +# +sub _spaceship { + my ($z1, $z2, $inverted) = @_; + my ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0); + my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0); + my $sgn = $inverted ? -1 : 1; + return $sgn * ($re1 <=> $re2) if $re1 != $re2; + return $sgn * ($im1 <=> $im2); +} + +# +# (_numeq) +# +# Computes z1 == z2. +# +# (Required in addition to _spaceship() because of NaNs.) +sub _numeq { + my ($z1, $z2, $inverted) = @_; + my ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0); + my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0); + return $re1 == $re2 && $im1 == $im2 ? 1 : 0; +} + +# +# (_negate) +# +# Computes -z. +# +sub _negate { + my ($z) = @_; + if ($z->{c_dirty}) { + my ($r, $t) = @{$z->_polar}; + $t = ($t <= 0) ? $t + pi : $t - pi; + return (ref $z)->emake($r, $t); + } + my ($re, $im) = @{$z->_cartesian}; + return (ref $z)->make(-$re, -$im); +} + +# +# (_conjugate) +# +# Compute complex's _conjugate. +# +sub _conjugate { + my ($z) = @_; + if ($z->{c_dirty}) { + my ($r, $t) = @{$z->_polar}; + return (ref $z)->emake($r, -$t); + } + my ($re, $im) = @{$z->_cartesian}; + return (ref $z)->make($re, -$im); +} + +# +# (abs) +# +# Compute or set complex's norm (rho). +# +sub abs { + my ($z, $rho) = @_ ? @_ : $_; + unless (ref $z) { + if (@_ == 2) { + $_[0] = $_[1]; + } else { + return CORE::abs($z); + } + } + if (defined $rho) { + $z->{'polar'} = [ $rho, ${$z->_polar}[1] ]; + $z->{p_dirty} = 0; + $z->{c_dirty} = 1; + return $rho; + } else { + return ${$z->_polar}[0]; + } +} + +sub _theta { + my $theta = $_[0]; + + if ($$theta > pi()) { $$theta -= pi2 } + elsif ($$theta <= -pi()) { $$theta += pi2 } +} + +# +# arg +# +# Compute or set complex's argument (theta). +# +sub arg { + my ($z, $theta) = @_; + return $z unless ref $z; + if (defined $theta) { + _theta(\$theta); + $z->{'polar'} = [ ${$z->_polar}[0], $theta ]; + $z->{p_dirty} = 0; + $z->{c_dirty} = 1; + } else { + $theta = ${$z->_polar}[1]; + _theta(\$theta); + } + return $theta; +} + +# +# (sqrt) +# +# Compute sqrt(z). +# +# It is quite tempting to use wantarray here so that in list context +# sqrt() would return the two solutions. This, however, would +# break things like +# +# print "sqrt(z) = ", sqrt($z), "\n"; +# +# The two values would be printed side by side without no intervening +# whitespace, quite confusing. +# Therefore if you want the two solutions use the root(). +# +sub sqrt { + my ($z) = @_ ? $_[0] : $_; + my ($re, $im) = ref $z ? @{$z->_cartesian} : ($z, 0); + return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re) + if $im == 0; + my ($r, $t) = @{$z->_polar}; + return (ref $z)->emake(CORE::sqrt($r), $t/2); +} + +# +# cbrt +# +# Compute cbrt(z) (cubic root). +# +# Why are we not returning three values? The same answer as for sqrt(). +# +sub cbrt { + my ($z) = @_; + return $z < 0 ? + -CORE::exp(CORE::log(-$z)/3) : + ($z > 0 ? CORE::exp(CORE::log($z)/3): 0) + unless ref $z; + my ($r, $t) = @{$z->_polar}; + return 0 if $r == 0; + return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3); +} + +# +# _rootbad +# +# Die on bad root. +# +sub _rootbad { + my $mess = "Root '$_[0]' illegal, root rank must be positive integer.\n"; + + my @up = caller(1); + + $mess .= "Died at $up[1] line $up[2].\n"; + + die $mess; +} + +# +# root +# +# Computes all nth root for z, returning an array whose size is n. +# `n' must be a positive integer. +# +# The roots are given by (for k = 0..n-1): +# +# z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n)) +# +sub root { + my ($z, $n, $k) = @_; + _rootbad($n) if ($n < 1 or int($n) != $n); + my ($r, $t) = ref $z ? + @{$z->_polar} : (CORE::abs($z), $z >= 0 ? 0 : pi); + my $theta_inc = pi2 / $n; + my $rho = $r ** (1/$n); + my $cartesian = ref $z && $z->{c_dirty} == 0; + if (@_ == 2) { + my @root; + for (my $i = 0, my $theta = $t / $n; + $i < $n; + $i++, $theta += $theta_inc) { + my $w = cplxe($rho, $theta); + # Yes, $cartesian is loop invariant. + push @root, $cartesian ? cplx(@{$w->_cartesian}) : $w; + } + return @root; + } elsif (@_ == 3) { + my $w = cplxe($rho, $t / $n + $k * $theta_inc); + return $cartesian ? cplx(@{$w->_cartesian}) : $w; + } +} + +# +# Re +# +# Return or set Re(z). +# +sub Re { + my ($z, $Re) = @_; + return $z unless ref $z; + if (defined $Re) { + $z->{'cartesian'} = [ $Re, ${$z->_cartesian}[1] ]; + $z->{c_dirty} = 0; + $z->{p_dirty} = 1; + } else { + return ${$z->_cartesian}[0]; + } +} + +# +# Im +# +# Return or set Im(z). +# +sub Im { + my ($z, $Im) = @_; + return 0 unless ref $z; + if (defined $Im) { + $z->{'cartesian'} = [ ${$z->_cartesian}[0], $Im ]; + $z->{c_dirty} = 0; + $z->{p_dirty} = 1; + } else { + return ${$z->_cartesian}[1]; + } +} + +# +# rho +# +# Return or set rho(w). +# +sub rho { + Math::Complex::abs(@_); +} + +# +# theta +# +# Return or set theta(w). +# +sub theta { + Math::Complex::arg(@_); +} + +# +# (exp) +# +# Computes exp(z). +# +sub exp { + my ($z) = @_ ? @_ : $_; + return CORE::exp($z) unless ref $z; + my ($x, $y) = @{$z->_cartesian}; + return (ref $z)->emake(CORE::exp($x), $y); +} + +# +# _logofzero +# +# Die on logarithm of zero. +# +sub _logofzero { + my $mess = "$_[0]: Logarithm of zero.\n"; + + if (defined $_[1]) { + $mess .= "(Because in the definition of $_[0], the argument "; + $mess .= "$_[1] " unless ($_[1] eq '0'); + $mess .= "is 0)\n"; + } + + my @up = caller(1); + + $mess .= "Died at $up[1] line $up[2].\n"; + + die $mess; +} + +# +# (log) +# +# Compute log(z). +# +sub log { + my ($z) = @_ ? @_ : $_; + unless (ref $z) { + _logofzero("log") if $z == 0; + return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi); + } + my ($r, $t) = @{$z->_polar}; + _logofzero("log") if $r == 0; + if ($t > pi()) { $t -= pi2 } + elsif ($t <= -pi()) { $t += pi2 } + return (ref $z)->make(CORE::log($r), $t); +} + +# +# ln +# +# Alias for log(). +# +sub ln { Math::Complex::log(@_) } + +# +# log10 +# +# Compute log10(z). +# + +sub log10 { + return Math::Complex::log($_[0]) * _uplog10; +} + +# +# logn +# +# Compute logn(z,n) = log(z) / log(n) +# +sub logn { + my ($z, $n) = @_; + $z = cplx($z, 0) unless ref $z; + my $logn = $LOGN{$n}; + $logn = $LOGN{$n} = CORE::log($n) unless defined $logn; # Cache log(n) + return &log($z) / $logn; +} + +# +# (cos) +# +# Compute cos(z) = (exp(iz) + exp(-iz))/2. +# +sub cos { + my ($z) = @_ ? @_ : $_; + return CORE::cos($z) unless ref $z; + my ($x, $y) = @{$z->_cartesian}; + my $ey = CORE::exp($y); + my $sx = CORE::sin($x); + my $cx = CORE::cos($x); + my $ey_1 = $ey ? 1 / $ey : Inf(); + return (ref $z)->make($cx * ($ey + $ey_1)/2, + $sx * ($ey_1 - $ey)/2); +} + +# +# (sin) +# +# Compute sin(z) = (exp(iz) - exp(-iz))/2. +# +sub sin { + my ($z) = @_ ? @_ : $_; + return CORE::sin($z) unless ref $z; + my ($x, $y) = @{$z->_cartesian}; + my $ey = CORE::exp($y); + my $sx = CORE::sin($x); + my $cx = CORE::cos($x); + my $ey_1 = $ey ? 1 / $ey : Inf(); + return (ref $z)->make($sx * ($ey + $ey_1)/2, + $cx * ($ey - $ey_1)/2); +} + +# +# tan +# +# Compute tan(z) = sin(z) / cos(z). +# +sub tan { + my ($z) = @_; + my $cz = &cos($z); + _divbyzero "tan($z)", "cos($z)" if $cz == 0; + return &sin($z) / $cz; +} + +# +# sec +# +# Computes the secant sec(z) = 1 / cos(z). +# +sub sec { + my ($z) = @_; + my $cz = &cos($z); + _divbyzero "sec($z)", "cos($z)" if ($cz == 0); + return 1 / $cz; +} + +# +# csc +# +# Computes the cosecant csc(z) = 1 / sin(z). +# +sub csc { + my ($z) = @_; + my $sz = &sin($z); + _divbyzero "csc($z)", "sin($z)" if ($sz == 0); + return 1 / $sz; +} + +# +# cosec +# +# Alias for csc(). +# +sub cosec { Math::Complex::csc(@_) } + +# +# cot +# +# Computes cot(z) = cos(z) / sin(z). +# +sub cot { + my ($z) = @_; + my $sz = &sin($z); + _divbyzero "cot($z)", "sin($z)" if ($sz == 0); + return &cos($z) / $sz; +} + +# +# cotan +# +# Alias for cot(). +# +sub cotan { Math::Complex::cot(@_) } + +# +# acos +# +# Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)). +# +sub acos { + my $z = $_[0]; + return CORE::atan2(CORE::sqrt(1-$z*$z), $z) + if (! ref $z) && CORE::abs($z) <= 1; + $z = cplx($z, 0) unless ref $z; + my ($x, $y) = @{$z->_cartesian}; + return 0 if $x == 1 && $y == 0; + my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y); + my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y); + my $alpha = ($t1 + $t2)/2; + my $beta = ($t1 - $t2)/2; + $alpha = 1 if $alpha < 1; + if ($beta > 1) { $beta = 1 } + elsif ($beta < -1) { $beta = -1 } + my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta); + my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1)); + $v = -$v if $y > 0 || ($y == 0 && $x < -1); + return (ref $z)->make($u, $v); +} + +# +# asin +# +# Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)). +# +sub asin { + my $z = $_[0]; + return CORE::atan2($z, CORE::sqrt(1-$z*$z)) + if (! ref $z) && CORE::abs($z) <= 1; + $z = cplx($z, 0) unless ref $z; + my ($x, $y) = @{$z->_cartesian}; + return 0 if $x == 0 && $y == 0; + my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y); + my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y); + my $alpha = ($t1 + $t2)/2; + my $beta = ($t1 - $t2)/2; + $alpha = 1 if $alpha < 1; + if ($beta > 1) { $beta = 1 } + elsif ($beta < -1) { $beta = -1 } + my $u = CORE::atan2($beta, CORE::sqrt(1-$beta*$beta)); + my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1)); + $v = -$v if $y > 0 || ($y == 0 && $x < -1); + return (ref $z)->make($u, $v); +} + +# +# atan +# +# Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)). +# +sub atan { + my ($z) = @_; + return CORE::atan2($z, 1) unless ref $z; + my ($x, $y) = ref $z ? @{$z->_cartesian} : ($z, 0); + return 0 if $x == 0 && $y == 0; + _divbyzero "atan(i)" if ( $z == i); + _logofzero "atan(-i)" if (-$z == i); # -i is a bad file test... + my $log = &log((i + $z) / (i - $z)); + return _ip2 * $log; +} + +# +# asec +# +# Computes the arc secant asec(z) = acos(1 / z). +# +sub asec { + my ($z) = @_; + _divbyzero "asec($z)", $z if ($z == 0); + return acos(1 / $z); +} + +# +# acsc +# +# Computes the arc cosecant acsc(z) = asin(1 / z). +# +sub acsc { + my ($z) = @_; + _divbyzero "acsc($z)", $z if ($z == 0); + return asin(1 / $z); +} + +# +# acosec +# +# Alias for acsc(). +# +sub acosec { Math::Complex::acsc(@_) } + +# +# acot +# +# Computes the arc cotangent acot(z) = atan(1 / z) +# +sub acot { + my ($z) = @_; + _divbyzero "acot(0)" if $z == 0; + return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z) + unless ref $z; + _divbyzero "acot(i)" if ($z - i == 0); + _logofzero "acot(-i)" if ($z + i == 0); + return atan(1 / $z); +} + +# +# acotan +# +# Alias for acot(). +# +sub acotan { Math::Complex::acot(@_) } + +# +# cosh +# +# Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2. +# +sub cosh { + my ($z) = @_; + my $ex; + unless (ref $z) { + $ex = CORE::exp($z); + return $ex ? ($ex == $ExpInf ? Inf() : ($ex + 1/$ex)/2) : Inf(); + } + my ($x, $y) = @{$z->_cartesian}; + $ex = CORE::exp($x); + my $ex_1 = $ex ? 1 / $ex : Inf(); + return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2, + CORE::sin($y) * ($ex - $ex_1)/2); +} + +# +# sinh +# +# Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2. +# +sub sinh { + my ($z) = @_; + my $ex; + unless (ref $z) { + return 0 if $z == 0; + $ex = CORE::exp($z); + return $ex ? ($ex == $ExpInf ? Inf() : ($ex - 1/$ex)/2) : -Inf(); + } + my ($x, $y) = @{$z->_cartesian}; + my $cy = CORE::cos($y); + my $sy = CORE::sin($y); + $ex = CORE::exp($x); + my $ex_1 = $ex ? 1 / $ex : Inf(); + return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2, + CORE::sin($y) * ($ex + $ex_1)/2); +} + +# +# tanh +# +# Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z). +# +sub tanh { + my ($z) = @_; + my $cz = cosh($z); + _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0); + my $sz = sinh($z); + return 1 if $cz == $sz; + return -1 if $cz == -$sz; + return $sz / $cz; +} + +# +# sech +# +# Computes the hyperbolic secant sech(z) = 1 / cosh(z). +# +sub sech { + my ($z) = @_; + my $cz = cosh($z); + _divbyzero "sech($z)", "cosh($z)" if ($cz == 0); + return 1 / $cz; +} + +# +# csch +# +# Computes the hyperbolic cosecant csch(z) = 1 / sinh(z). +# +sub csch { + my ($z) = @_; + my $sz = sinh($z); + _divbyzero "csch($z)", "sinh($z)" if ($sz == 0); + return 1 / $sz; +} + +# +# cosech +# +# Alias for csch(). +# +sub cosech { Math::Complex::csch(@_) } + +# +# coth +# +# Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z). +# +sub coth { + my ($z) = @_; + my $sz = sinh($z); + _divbyzero "coth($z)", "sinh($z)" if $sz == 0; + my $cz = cosh($z); + return 1 if $cz == $sz; + return -1 if $cz == -$sz; + return $cz / $sz; +} + +# +# cotanh +# +# Alias for coth(). +# +sub cotanh { Math::Complex::coth(@_) } + +# +# acosh +# +# Computes the area/inverse hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)). +# +sub acosh { + my ($z) = @_; + unless (ref $z) { + $z = cplx($z, 0); + } + my ($re, $im) = @{$z->_cartesian}; + if ($im == 0) { + return CORE::log($re + CORE::sqrt($re*$re - 1)) + if $re >= 1; + return cplx(0, CORE::atan2(CORE::sqrt(1 - $re*$re), $re)) + if CORE::abs($re) < 1; + } + my $t = &sqrt($z * $z - 1) + $z; + # Try Taylor if looking bad (this usually means that + # $z was large negative, therefore the sqrt is really + # close to abs(z), summing that with z...) + $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7) + if $t == 0; + my $u = &log($t); + $u->Im(-$u->Im) if $re < 0 && $im == 0; + return $re < 0 ? -$u : $u; +} + +# +# asinh +# +# Computes the area/inverse hyperbolic sine asinh(z) = log(z + sqrt(z*z+1)) +# +sub asinh { + my ($z) = @_; + unless (ref $z) { + my $t = $z + CORE::sqrt($z*$z + 1); + return CORE::log($t) if $t; + } + my $t = &sqrt($z * $z + 1) + $z; + # Try Taylor if looking bad (this usually means that + # $z was large negative, therefore the sqrt is really + # close to abs(z), summing that with z...) + $t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7) + if $t == 0; + return &log($t); +} + +# +# atanh +# +# Computes the area/inverse hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)). +# +sub atanh { + my ($z) = @_; + unless (ref $z) { + return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1; + $z = cplx($z, 0); + } + _divbyzero 'atanh(1)', "1 - $z" if (1 - $z == 0); + _logofzero 'atanh(-1)' if (1 + $z == 0); + return 0.5 * &log((1 + $z) / (1 - $z)); +} + +# +# asech +# +# Computes the area/inverse hyperbolic secant asech(z) = acosh(1 / z). +# +sub asech { + my ($z) = @_; + _divbyzero 'asech(0)', "$z" if ($z == 0); + return acosh(1 / $z); +} + +# +# acsch +# +# Computes the area/inverse hyperbolic cosecant acsch(z) = asinh(1 / z). +# +sub acsch { + my ($z) = @_; + _divbyzero 'acsch(0)', $z if ($z == 0); + return asinh(1 / $z); +} + +# +# acosech +# +# Alias for acosh(). +# +sub acosech { Math::Complex::acsch(@_) } + +# +# acoth +# +# Computes the area/inverse hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)). +# +sub acoth { + my ($z) = @_; + _divbyzero 'acoth(0)' if ($z == 0); + unless (ref $z) { + return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1; + $z = cplx($z, 0); + } + _divbyzero 'acoth(1)', "$z - 1" if ($z - 1 == 0); + _logofzero 'acoth(-1)', "1 + $z" if (1 + $z == 0); + return &log((1 + $z) / ($z - 1)) / 2; +} + +# +# acotanh +# +# Alias for acot(). +# +sub acotanh { Math::Complex::acoth(@_) } + +# +# (atan2) +# +# Compute atan(z1/z2), minding the right quadrant. +# +sub atan2 { + my ($z1, $z2, $inverted) = @_; + my ($re1, $im1, $re2, $im2); + if ($inverted) { + ($re1, $im1) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0); + ($re2, $im2) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0); + } else { + ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0); + ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0); + } + if ($im1 || $im2) { + # In MATLAB the imaginary parts are ignored. + # warn "atan2: Imaginary parts ignored"; + # http://documents.wolfram.com/mathematica/functions/ArcTan + # NOTE: Mathematica ArcTan[x,y] while atan2(y,x) + my $s = $z1 * $z1 + $z2 * $z2; + _divbyzero("atan2") if $s == 0; + my $i = &i; + my $r = $z2 + $z1 * $i; + return -$i * &log($r / &sqrt( $s )); + } + return CORE::atan2($re1, $re2); +} + +# +# display_format +# ->display_format +# +# Set (get if no argument) the display format for all complex numbers that +# don't happen to have overridden it via ->display_format +# +# When called as an object method, this actually sets the display format for +# the current object. +# +# Valid object formats are 'c' and 'p' for cartesian and polar. The first +# letter is used actually, so the type can be fully spelled out for clarity. +# +sub display_format { + my $self = shift; + my %display_format = %DISPLAY_FORMAT; + + if (ref $self) { # Called as an object method + if (exists $self->{display_format}) { + my %obj = %{$self->{display_format}}; + @display_format{keys %obj} = values %obj; + } + } + if (@_ == 1) { + $display_format{style} = shift; + } else { + my %new = @_; + @display_format{keys %new} = values %new; + } + + if (ref $self) { # Called as an object method + $self->{display_format} = { %display_format }; + return + wantarray ? + %{$self->{display_format}} : + $self->{display_format}->{style}; + } + + # Called as a class method + %DISPLAY_FORMAT = %display_format; + return + wantarray ? + %DISPLAY_FORMAT : + $DISPLAY_FORMAT{style}; +} + +# +# (_stringify) +# +# Show nicely formatted complex number under its cartesian or polar form, +# depending on the current display format: +# +# . If a specific display format has been recorded for this object, use it. +# . Otherwise, use the generic current default for all complex numbers, +# which is a package global variable. +# +sub _stringify { + my ($z) = shift; + + my $style = $z->display_format; + + $style = $DISPLAY_FORMAT{style} unless defined $style; + + return $z->_stringify_polar if $style =~ /^p/i; + return $z->_stringify_cartesian; +} + +# +# ->_stringify_cartesian +# +# Stringify as a cartesian representation 'a+bi'. +# +sub _stringify_cartesian { + my $z = shift; + my ($x, $y) = @{$z->_cartesian}; + my ($re, $im); + + my %format = $z->display_format; + my $format = $format{format}; + + if ($x) { + if ($x =~ /^NaN[QS]?$/i) { + $re = $x; + } else { + if ($x =~ /^-?\Q$Inf\E$/oi) { + $re = $x; + } else { + $re = defined $format ? sprintf($format, $x) : $x; + } + } + } else { + undef $re; + } + + if ($y) { + if ($y =~ /^(NaN[QS]?)$/i) { + $im = $y; + } else { + if ($y =~ /^-?\Q$Inf\E$/oi) { + $im = $y; + } else { + $im = + defined $format ? + sprintf($format, $y) : + ($y == 1 ? "" : ($y == -1 ? "-" : $y)); + } + } + $im .= "i"; + } else { + undef $im; + } + + my $str = $re; + + if (defined $im) { + if ($y < 0) { + $str .= $im; + } elsif ($y > 0 || $im =~ /^NaN[QS]?i$/i) { + $str .= "+" if defined $re; + $str .= $im; + } + } elsif (!defined $re) { + $str = "0"; + } + + return $str; +} + + +# +# ->_stringify_polar +# +# Stringify as a polar representation '[r,t]'. +# +sub _stringify_polar { + my $z = shift; + my ($r, $t) = @{$z->_polar}; + my $theta; + + my %format = $z->display_format; + my $format = $format{format}; + + if ($t =~ /^NaN[QS]?$/i || $t =~ /^-?\Q$Inf\E$/oi) { + $theta = $t; + } elsif ($t == pi) { + $theta = "pi"; + } elsif ($r == 0 || $t == 0) { + $theta = defined $format ? sprintf($format, $t) : $t; + } + + return "[$r,$theta]" if defined $theta; + + # + # Try to identify pi/n and friends. + # + + $t -= int(CORE::abs($t) / pi2) * pi2; + + if ($format{polar_pretty_print} && $t) { + my ($a, $b); + for $a (2..9) { + $b = $t * $a / pi; + if ($b =~ /^-?\d+$/) { + $b = $b < 0 ? "-" : "" if CORE::abs($b) == 1; + $theta = "${b}pi/$a"; + last; + } + } + } + + if (defined $format) { + $r = sprintf($format, $r); + $theta = sprintf($format, $t) unless defined $theta; + } else { + $theta = $t unless defined $theta; + } + + return "[$r,$theta]"; +} + +sub Inf { + return $Inf; +} + +1; +__END__ + +=pod + +=head1 NAME + +Math::Complex - complex numbers and associated mathematical functions + +=head1 SYNOPSIS + + use Math::Complex; + + $z = Math::Complex->make(5, 6); + $t = 4 - 3*i + $z; + $j = cplxe(1, 2*pi/3); + +=head1 DESCRIPTION + +This package lets you create and manipulate complex numbers. By default, +I<Perl> limits itself to real numbers, but an extra C<use> statement brings +full complex support, along with a full set of mathematical functions +typically associated with and/or extended to complex numbers. + +If you wonder what complex numbers are, they were invented to be able to solve +the following equation: + + x*x = -1 + +and by definition, the solution is noted I<i> (engineers use I<j> instead since +I<i> usually denotes an intensity, but the name does not matter). The number +I<i> is a pure I<imaginary> number. + +The arithmetics with pure imaginary numbers works just like you would expect +it with real numbers... you just have to remember that + + i*i = -1 + +so you have: + + 5i + 7i = i * (5 + 7) = 12i + 4i - 3i = i * (4 - 3) = i + 4i * 2i = -8 + 6i / 2i = 3 + 1 / i = -i + +Complex numbers are numbers that have both a real part and an imaginary +part, and are usually noted: + + a + bi + +where C<a> is the I<real> part and C<b> is the I<imaginary> part. The +arithmetic with complex numbers is straightforward. You have to +keep track of the real and the imaginary parts, but otherwise the +rules used for real numbers just apply: + + (4 + 3i) + (5 - 2i) = (4 + 5) + i(3 - 2) = 9 + i + (2 + i) * (4 - i) = 2*4 + 4i -2i -i*i = 8 + 2i + 1 = 9 + 2i + +A graphical representation of complex numbers is possible in a plane +(also called the I<complex plane>, but it's really a 2D plane). +The number + + z = a + bi + +is the point whose coordinates are (a, b). Actually, it would +be the vector originating from (0, 0) to (a, b). It follows that the addition +of two complex numbers is a vectorial addition. + +Since there is a bijection between a point in the 2D plane and a complex +number (i.e. the mapping is unique and reciprocal), a complex number +can also be uniquely identified with polar coordinates: + + [rho, theta] + +where C<rho> is the distance to the origin, and C<theta> the angle between +the vector and the I<x> axis. There is a notation for this using the +exponential form, which is: + + rho * exp(i * theta) + +where I<i> is the famous imaginary number introduced above. Conversion +between this form and the cartesian form C<a + bi> is immediate: + + a = rho * cos(theta) + b = rho * sin(theta) + +which is also expressed by this formula: + + z = rho * exp(i * theta) = rho * (cos theta + i * sin theta) + +In other words, it's the projection of the vector onto the I<x> and I<y> +axes. Mathematicians call I<rho> the I<norm> or I<modulus> and I<theta> +the I<argument> of the complex number. The I<norm> of C<z> is +marked here as C<abs(z)>. + +The polar notation (also known as the trigonometric representation) is +much more handy for performing multiplications and divisions of +complex numbers, whilst the cartesian notation is better suited for +additions and subtractions. Real numbers are on the I<x> axis, and +therefore I<y> or I<theta> is zero or I<pi>. + +All the common operations that can be performed on a real number have +been defined to work on complex numbers as well, and are merely +I<extensions> of the operations defined on real numbers. This means +they keep their natural meaning when there is no imaginary part, provided +the number is within their definition set. + +For instance, the C<sqrt> routine which computes the square root of +its argument is only defined for non-negative real numbers and yields a +non-negative real number (it is an application from B<R+> to B<R+>). +If we allow it to return a complex number, then it can be extended to +negative real numbers to become an application from B<R> to B<C> (the +set of complex numbers): + + sqrt(x) = x >= 0 ? sqrt(x) : sqrt(-x)*i + +It can also be extended to be an application from B<C> to B<C>, +whilst its restriction to B<R> behaves as defined above by using +the following definition: + + sqrt(z = [r,t]) = sqrt(r) * exp(i * t/2) + +Indeed, a negative real number can be noted C<[x,pi]> (the modulus +I<x> is always non-negative, so C<[x,pi]> is really C<-x>, a negative +number) and the above definition states that + + sqrt([x,pi]) = sqrt(x) * exp(i*pi/2) = [sqrt(x),pi/2] = sqrt(x)*i + +which is exactly what we had defined for negative real numbers above. +The C<sqrt> returns only one of the solutions: if you want the both, +use the C<root> function. + +All the common mathematical functions defined on real numbers that +are extended to complex numbers share that same property of working +I<as usual> when the imaginary part is zero (otherwise, it would not +be called an extension, would it?). + +A I<new> operation possible on a complex number that is +the identity for real numbers is called the I<conjugate>, and is noted +with a horizontal bar above the number, or C<~z> here. + + z = a + bi + ~z = a - bi + +Simple... Now look: + + z * ~z = (a + bi) * (a - bi) = a*a + b*b + +We saw that the norm of C<z> was noted C<abs(z)> and was defined as the +distance to the origin, also known as: + + rho = abs(z) = sqrt(a*a + b*b) + +so + + z * ~z = abs(z) ** 2 + +If z is a pure real number (i.e. C<b == 0>), then the above yields: + + a * a = abs(a) ** 2 + +which is true (C<abs> has the regular meaning for real number, i.e. stands +for the absolute value). This example explains why the norm of C<z> is +noted C<abs(z)>: it extends the C<abs> function to complex numbers, yet +is the regular C<abs> we know when the complex number actually has no +imaginary part... This justifies I<a posteriori> our use of the C<abs> +notation for the norm. + +=head1 OPERATIONS + +Given the following notations: + + z1 = a + bi = r1 * exp(i * t1) + z2 = c + di = r2 * exp(i * t2) + z = <any complex or real number> + +the following (overloaded) operations are supported on complex numbers: + + z1 + z2 = (a + c) + i(b + d) + z1 - z2 = (a - c) + i(b - d) + z1 * z2 = (r1 * r2) * exp(i * (t1 + t2)) + z1 / z2 = (r1 / r2) * exp(i * (t1 - t2)) + z1 ** z2 = exp(z2 * log z1) + ~z = a - bi + abs(z) = r1 = sqrt(a*a + b*b) + sqrt(z) = sqrt(r1) * exp(i * t/2) + exp(z) = exp(a) * exp(i * b) + log(z) = log(r1) + i*t + sin(z) = 1/2i (exp(i * z1) - exp(-i * z)) + cos(z) = 1/2 (exp(i * z1) + exp(-i * z)) + atan2(y, x) = atan(y / x) # Minding the right quadrant, note the order. + +The definition used for complex arguments of atan2() is + + -i log((x + iy)/sqrt(x*x+y*y)) + +Note that atan2(0, 0) is not well-defined. + +The following extra operations are supported on both real and complex +numbers: + + Re(z) = a + Im(z) = b + arg(z) = t + abs(z) = r + + cbrt(z) = z ** (1/3) + log10(z) = log(z) / log(10) + logn(z, n) = log(z) / log(n) + + tan(z) = sin(z) / cos(z) + + csc(z) = 1 / sin(z) + sec(z) = 1 / cos(z) + cot(z) = 1 / tan(z) + + asin(z) = -i * log(i*z + sqrt(1-z*z)) + acos(z) = -i * log(z + i*sqrt(1-z*z)) + atan(z) = i/2 * log((i+z) / (i-z)) + + acsc(z) = asin(1 / z) + asec(z) = acos(1 / z) + acot(z) = atan(1 / z) = -i/2 * log((i+z) / (z-i)) + + sinh(z) = 1/2 (exp(z) - exp(-z)) + cosh(z) = 1/2 (exp(z) + exp(-z)) + tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z)) + + csch(z) = 1 / sinh(z) + sech(z) = 1 / cosh(z) + coth(z) = 1 / tanh(z) + + asinh(z) = log(z + sqrt(z*z+1)) + acosh(z) = log(z + sqrt(z*z-1)) + atanh(z) = 1/2 * log((1+z) / (1-z)) + + acsch(z) = asinh(1 / z) + asech(z) = acosh(1 / z) + acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1)) + +I<arg>, I<abs>, I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>, +I<coth>, I<acosech>, I<acotanh>, have aliases I<rho>, I<theta>, I<ln>, +I<cosec>, I<cotan>, I<acosec>, I<acotan>, I<cosech>, I<cotanh>, +I<acosech>, I<acotanh>, respectively. C<Re>, C<Im>, C<arg>, C<abs>, +C<rho>, and C<theta> can be used also as mutators. The C<cbrt> +returns only one of the solutions: if you want all three, use the +C<root> function. + +The I<root> function is available to compute all the I<n> +roots of some complex, where I<n> is a strictly positive integer. +There are exactly I<n> such roots, returned as a list. Getting the +number mathematicians call C<j> such that: + + 1 + j + j*j = 0; + +is a simple matter of writing: + + $j = ((root(1, 3))[1]; + +The I<k>th root for C<z = [r,t]> is given by: + + (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n) + +You can return the I<k>th root directly by C<root(z, n, k)>, +indexing starting from I<zero> and ending at I<n - 1>. + +The I<spaceship> numeric comparison operator, E<lt>=E<gt>, is also +defined. In order to ensure its restriction to real numbers is conform +to what you would expect, the comparison is run on the real part of +the complex number first, and imaginary parts are compared only when +the real parts match. + +=head1 CREATION + +To create a complex number, use either: + + $z = Math::Complex->make(3, 4); + $z = cplx(3, 4); + +if you know the cartesian form of the number, or + + $z = 3 + 4*i; + +if you like. To create a number using the polar form, use either: + + $z = Math::Complex->emake(5, pi/3); + $x = cplxe(5, pi/3); + +instead. The first argument is the modulus, the second is the angle +(in radians, the full circle is 2*pi). (Mnemonic: C<e> is used as a +notation for complex numbers in the polar form). + +It is possible to write: + + $x = cplxe(-3, pi/4); + +but that will be silently converted into C<[3,-3pi/4]>, since the +modulus must be non-negative (it represents the distance to the origin +in the complex plane). + +It is also possible to have a complex number as either argument of the +C<make>, C<emake>, C<cplx>, and C<cplxe>: the appropriate component of +the argument will be used. + + $z1 = cplx(-2, 1); + $z2 = cplx($z1, 4); + +The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also +understand a single (string) argument of the forms + + 2-3i + -3i + [2,3] + [2,-3pi/4] + [2] + +in which case the appropriate cartesian and exponential components +will be parsed from the string and used to create new complex numbers. +The imaginary component and the theta, respectively, will default to zero. + +The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also +understand the case of no arguments: this means plain zero or (0, 0). + +=head1 DISPLAYING + +When printed, a complex number is usually shown under its cartesian +style I<a+bi>, but there are legitimate cases where the polar style +I<[r,t]> is more appropriate. The process of converting the complex +number into a string that can be displayed is known as I<stringification>. + +By calling the class method C<Math::Complex::display_format> and +supplying either C<"polar"> or C<"cartesian"> as an argument, you +override the default display style, which is C<"cartesian">. Not +supplying any argument returns the current settings. + +This default can be overridden on a per-number basis by calling the +C<display_format> method instead. As before, not supplying any argument +returns the current display style for this number. Otherwise whatever you +specify will be the new display style for I<this> particular number. + +For instance: + + use Math::Complex; + + Math::Complex::display_format('polar'); + $j = (root(1, 3))[1]; + print "j = $j\n"; # Prints "j = [1,2pi/3]" + $j->display_format('cartesian'); + print "j = $j\n"; # Prints "j = -0.5+0.866025403784439i" + +The polar style attempts to emphasize arguments like I<k*pi/n> +(where I<n> is a positive integer and I<k> an integer within [-9, +9]), +this is called I<polar pretty-printing>. + +For the reverse of stringifying, see the C<make> and C<emake>. + +=head2 CHANGED IN PERL 5.6 + +The C<display_format> class method and the corresponding +C<display_format> object method can now be called using +a parameter hash instead of just a one parameter. + +The old display format style, which can have values C<"cartesian"> or +C<"polar">, can be changed using the C<"style"> parameter. + + $j->display_format(style => "polar"); + +The one parameter calling convention also still works. + + $j->display_format("polar"); + +There are two new display parameters. + +The first one is C<"format">, which is a sprintf()-style format string +to be used for both numeric parts of the complex number(s). The is +somewhat system-dependent but most often it corresponds to C<"%.15g">. +You can revert to the default by setting the C<format> to C<undef>. + + # the $j from the above example + + $j->display_format('format' => '%.5f'); + print "j = $j\n"; # Prints "j = -0.50000+0.86603i" + $j->display_format('format' => undef); + print "j = $j\n"; # Prints "j = -0.5+0.86603i" + +Notice that this affects also the return values of the +C<display_format> methods: in list context the whole parameter hash +will be returned, as opposed to only the style parameter value. +This is a potential incompatibility with earlier versions if you +have been calling the C<display_format> method in list context. + +The second new display parameter is C<"polar_pretty_print">, which can +be set to true or false, the default being true. See the previous +section for what this means. + +=head1 USAGE + +Thanks to overloading, the handling of arithmetics with complex numbers +is simple and almost transparent. + +Here are some examples: + + use Math::Complex; + + $j = cplxe(1, 2*pi/3); # $j ** 3 == 1 + print "j = $j, j**3 = ", $j ** 3, "\n"; + print "1 + j + j**2 = ", 1 + $j + $j**2, "\n"; + + $z = -16 + 0*i; # Force it to be a complex + print "sqrt($z) = ", sqrt($z), "\n"; + + $k = exp(i * 2*pi/3); + print "$j - $k = ", $j - $k, "\n"; + + $z->Re(3); # Re, Im, arg, abs, + $j->arg(2); # (the last two aka rho, theta) + # can be used also as mutators. + +=head1 CONSTANTS + +=head2 PI + +The constant C<pi> and some handy multiples of it (pi2, pi4, +and pip2 (pi/2) and pip4 (pi/4)) are also available if separately +exported: + + use Math::Complex ':pi'; + $third_of_circle = pi2 / 3; + +=head2 Inf + +The floating point infinity can be exported as a subroutine Inf(): + + use Math::Complex qw(Inf sinh); + my $AlsoInf = Inf() + 42; + my $AnotherInf = sinh(1e42); + print "$AlsoInf is $AnotherInf\n" if $AlsoInf == $AnotherInf; + +Note that the stringified form of infinity varies between platforms: +it can be for example any of + + inf + infinity + INF + 1.#INF + +or it can be something else. + +Also note that in some platforms trying to use the infinity in +arithmetic operations may result in Perl crashing because using +an infinity causes SIGFPE or its moral equivalent to be sent. +The way to ignore this is + + local $SIG{FPE} = sub { }; + +=head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO + +The division (/) and the following functions + + log ln log10 logn + tan sec csc cot + atan asec acsc acot + tanh sech csch coth + atanh asech acsch acoth + +cannot be computed for all arguments because that would mean dividing +by zero or taking logarithm of zero. These situations cause fatal +runtime errors looking like this + + cot(0): Division by zero. + (Because in the definition of cot(0), the divisor sin(0) is 0) + Died at ... + +or + + atanh(-1): Logarithm of zero. + Died at... + +For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>, +C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the +logarithmic functions and the C<atanh>, C<acoth>, the argument cannot +be C<1> (one). For the C<atanh>, C<acoth>, the argument cannot be +C<-1> (minus one). For the C<atan>, C<acot>, the argument cannot be +C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument +cannot be C<-i> (the negative imaginary unit). For the C<tan>, +C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k> +is any integer. atan2(0, 0) is undefined, and if the complex arguments +are used for atan2(), a division by zero will happen if z1**2+z2**2 == 0. + +Note that because we are operating on approximations of real numbers, +these errors can happen when merely `too close' to the singularities +listed above. + +=head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS + +The C<make> and C<emake> accept both real and complex arguments. +When they cannot recognize the arguments they will die with error +messages like the following + + Math::Complex::make: Cannot take real part of ... + Math::Complex::make: Cannot take real part of ... + Math::Complex::emake: Cannot take rho of ... + Math::Complex::emake: Cannot take theta of ... + +=head1 BUGS + +Saying C<use Math::Complex;> exports many mathematical routines in the +caller environment and even overrides some (C<sqrt>, C<log>, C<atan2>). +This is construed as a feature by the Authors, actually... ;-) + +All routines expect to be given real or complex numbers. Don't attempt to +use BigFloat, since Perl has currently no rule to disambiguate a '+' +operation (for instance) between two overloaded entities. + +In Cray UNICOS there is some strange numerical instability that results +in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast. Beware. +The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex. +Whatever it is, it does not manifest itself anywhere else where Perl runs. + +=head1 SEE ALSO + +L<Math::Trig> + +=head1 AUTHORS + +Daniel S. Lewart <F<lewart!at!uiuc.edu>>, +Jarkko Hietaniemi <F<jhi!at!iki.fi>>, +Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>, +Zefram <zefram@fysh.org> + +=head1 LICENSE + +This library is free software; you can redistribute it and/or modify +it under the same terms as Perl itself. + +=cut + +1; + +# eof diff --git a/dist/Math-Complex/lib/Math/Trig.pm b/dist/Math-Complex/lib/Math/Trig.pm new file mode 100644 index 0000000000..1d9612a41c --- /dev/null +++ b/dist/Math-Complex/lib/Math/Trig.pm @@ -0,0 +1,761 @@ +# +# Trigonometric functions, mostly inherited from Math::Complex. +# -- Jarkko Hietaniemi, since April 1997 +# -- Raphael Manfredi, September 1996 (indirectly: because of Math::Complex) +# + +package Math::Trig; + +{ use 5.006; } +use strict; + +use Math::Complex 1.59; +use Math::Complex qw(:trig :pi); +require Exporter; + +our @ISA = qw(Exporter); + +our $VERSION = 1.23; + +my @angcnv = qw(rad2deg rad2grad + deg2rad deg2grad + grad2rad grad2deg); + +my @areal = qw(asin_real acos_real); + +our @EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}}, + @angcnv, @areal); + +my @rdlcnv = qw(cartesian_to_cylindrical + cartesian_to_spherical + cylindrical_to_cartesian + cylindrical_to_spherical + spherical_to_cartesian + spherical_to_cylindrical); + +my @greatcircle = qw( + great_circle_distance + great_circle_direction + great_circle_bearing + great_circle_waypoint + great_circle_midpoint + great_circle_destination + ); + +my @pi = qw(pi pi2 pi4 pip2 pip4); + +our @EXPORT_OK = (@rdlcnv, @greatcircle, @pi, 'Inf'); + +# See e.g. the following pages: +# http://www.movable-type.co.uk/scripts/LatLong.html +# http://williams.best.vwh.net/avform.htm + +our %EXPORT_TAGS = ('radial' => [ @rdlcnv ], + 'great_circle' => [ @greatcircle ], + 'pi' => [ @pi ]); + +sub _DR () { pi2/360 } +sub _RD () { 360/pi2 } +sub _DG () { 400/360 } +sub _GD () { 360/400 } +sub _RG () { 400/pi2 } +sub _GR () { pi2/400 } + +# +# Truncating remainder. +# + +sub _remt ($$) { + # Oh yes, POSIX::fmod() would be faster. Possibly. If it is available. + $_[0] - $_[1] * int($_[0] / $_[1]); +} + +# +# Angle conversions. +# + +sub rad2rad($) { _remt($_[0], pi2) } + +sub deg2deg($) { _remt($_[0], 360) } + +sub grad2grad($) { _remt($_[0], 400) } + +sub rad2deg ($;$) { my $d = _RD * $_[0]; $_[1] ? $d : deg2deg($d) } + +sub deg2rad ($;$) { my $d = _DR * $_[0]; $_[1] ? $d : rad2rad($d) } + +sub grad2deg ($;$) { my $d = _GD * $_[0]; $_[1] ? $d : deg2deg($d) } + +sub deg2grad ($;$) { my $d = _DG * $_[0]; $_[1] ? $d : grad2grad($d) } + +sub rad2grad ($;$) { my $d = _RG * $_[0]; $_[1] ? $d : grad2grad($d) } + +sub grad2rad ($;$) { my $d = _GR * $_[0]; $_[1] ? $d : rad2rad($d) } + +# +# acos and asin functions which always return a real number +# + +sub acos_real { + return 0 if $_[0] >= 1; + return pi if $_[0] <= -1; + return acos($_[0]); +} + +sub asin_real { + return &pip2 if $_[0] >= 1; + return -&pip2 if $_[0] <= -1; + return asin($_[0]); +} + +sub cartesian_to_spherical { + my ( $x, $y, $z ) = @_; + + my $rho = sqrt( $x * $x + $y * $y + $z * $z ); + + return ( $rho, + atan2( $y, $x ), + $rho ? acos_real( $z / $rho ) : 0 ); +} + +sub spherical_to_cartesian { + my ( $rho, $theta, $phi ) = @_; + + return ( $rho * cos( $theta ) * sin( $phi ), + $rho * sin( $theta ) * sin( $phi ), + $rho * cos( $phi ) ); +} + +sub spherical_to_cylindrical { + my ( $x, $y, $z ) = spherical_to_cartesian( @_ ); + + return ( sqrt( $x * $x + $y * $y ), $_[1], $z ); +} + +sub cartesian_to_cylindrical { + my ( $x, $y, $z ) = @_; + + return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z ); +} + +sub cylindrical_to_cartesian { + my ( $rho, $theta, $z ) = @_; + + return ( $rho * cos( $theta ), $rho * sin( $theta ), $z ); +} + +sub cylindrical_to_spherical { + return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) ); +} + +sub great_circle_distance { + my ( $theta0, $phi0, $theta1, $phi1, $rho ) = @_; + + $rho = 1 unless defined $rho; # Default to the unit sphere. + + my $lat0 = pip2 - $phi0; + my $lat1 = pip2 - $phi1; + + return $rho * + acos_real( cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) + + sin( $lat0 ) * sin( $lat1 ) ); +} + +sub great_circle_direction { + my ( $theta0, $phi0, $theta1, $phi1 ) = @_; + + my $lat0 = pip2 - $phi0; + my $lat1 = pip2 - $phi1; + + return rad2rad(pi2 - + atan2(sin($theta0-$theta1) * cos($lat1), + cos($lat0) * sin($lat1) - + sin($lat0) * cos($lat1) * cos($theta0-$theta1))); +} + +*great_circle_bearing = \&great_circle_direction; + +sub great_circle_waypoint { + my ( $theta0, $phi0, $theta1, $phi1, $point ) = @_; + + $point = 0.5 unless defined $point; + + my $d = great_circle_distance( $theta0, $phi0, $theta1, $phi1 ); + + return undef if $d == pi; + + my $sd = sin($d); + + return ($theta0, $phi0) if $sd == 0; + + my $A = sin((1 - $point) * $d) / $sd; + my $B = sin( $point * $d) / $sd; + + my $lat0 = pip2 - $phi0; + my $lat1 = pip2 - $phi1; + + my $x = $A * cos($lat0) * cos($theta0) + $B * cos($lat1) * cos($theta1); + my $y = $A * cos($lat0) * sin($theta0) + $B * cos($lat1) * sin($theta1); + my $z = $A * sin($lat0) + $B * sin($lat1); + + my $theta = atan2($y, $x); + my $phi = acos_real($z); + + return ($theta, $phi); +} + +sub great_circle_midpoint { + great_circle_waypoint(@_[0..3], 0.5); +} + +sub great_circle_destination { + my ( $theta0, $phi0, $dir0, $dst ) = @_; + + my $lat0 = pip2 - $phi0; + + my $phi1 = asin_real(sin($lat0)*cos($dst) + + cos($lat0)*sin($dst)*cos($dir0)); + + my $theta1 = $theta0 + atan2(sin($dir0)*sin($dst)*cos($lat0), + cos($dst)-sin($lat0)*sin($phi1)); + + my $dir1 = great_circle_bearing($theta1, $phi1, $theta0, $phi0) + pi; + + $dir1 -= pi2 if $dir1 > pi2; + + return ($theta1, $phi1, $dir1); +} + +1; + +__END__ +=pod + +=head1 NAME + +Math::Trig - trigonometric functions + +=head1 SYNOPSIS + + use Math::Trig; + + $x = tan(0.9); + $y = acos(3.7); + $z = asin(2.4); + + $halfpi = pi/2; + + $rad = deg2rad(120); + + # Import constants pi2, pip2, pip4 (2*pi, pi/2, pi/4). + use Math::Trig ':pi'; + + # Import the conversions between cartesian/spherical/cylindrical. + use Math::Trig ':radial'; + + # Import the great circle formulas. + use Math::Trig ':great_circle'; + +=head1 DESCRIPTION + +C<Math::Trig> defines many trigonometric functions not defined by the +core Perl which defines only the C<sin()> and C<cos()>. The constant +B<pi> is also defined as are a few convenience functions for angle +conversions, and I<great circle formulas> for spherical movement. + +=head1 TRIGONOMETRIC FUNCTIONS + +The tangent + +=over 4 + +=item B<tan> + +=back + +The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot +are aliases) + +B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan> + +The arcus (also known as the inverse) functions of the sine, cosine, +and tangent + +B<asin>, B<acos>, B<atan> + +The principal value of the arc tangent of y/x + +B<atan2>(y, x) + +The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc +and acotan/acot are aliases). Note that atan2(0, 0) is not well-defined. + +B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan> + +The hyperbolic sine, cosine, and tangent + +B<sinh>, B<cosh>, B<tanh> + +The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch +and cotanh/coth are aliases) + +B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh> + +The area (also known as the inverse) functions of the hyperbolic +sine, cosine, and tangent + +B<asinh>, B<acosh>, B<atanh> + +The area cofunctions of the hyperbolic sine, cosine, and tangent +(acsch/acosech and acoth/acotanh are aliases) + +B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh> + +The trigonometric constant B<pi> and some of handy multiples +of it are also defined. + +B<pi, pi2, pi4, pip2, pip4> + +=head2 ERRORS DUE TO DIVISION BY ZERO + +The following functions + + acoth + acsc + acsch + asec + asech + atanh + cot + coth + csc + csch + sec + sech + tan + tanh + +cannot be computed for all arguments because that would mean dividing +by zero or taking logarithm of zero. These situations cause fatal +runtime errors looking like this + + cot(0): Division by zero. + (Because in the definition of cot(0), the divisor sin(0) is 0) + Died at ... + +or + + atanh(-1): Logarithm of zero. + Died at... + +For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>, +C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the +C<atanh>, C<acoth>, the argument cannot be C<1> (one). For the +C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one). For the +C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k * +pi>, where I<k> is any integer. + +Note that atan2(0, 0) is not well-defined. + +=head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS + +Please note that some of the trigonometric functions can break out +from the B<real axis> into the B<complex plane>. For example +C<asin(2)> has no definition for plain real numbers but it has +definition for complex numbers. + +In Perl terms this means that supplying the usual Perl numbers (also +known as scalars, please see L<perldata>) as input for the +trigonometric functions might produce as output results that no more +are simple real numbers: instead they are complex numbers. + +The C<Math::Trig> handles this by using the C<Math::Complex> package +which knows how to handle complex numbers, please see L<Math::Complex> +for more information. In practice you need not to worry about getting +complex numbers as results because the C<Math::Complex> takes care of +details like for example how to display complex numbers. For example: + + print asin(2), "\n"; + +should produce something like this (take or leave few last decimals): + + 1.5707963267949-1.31695789692482i + +That is, a complex number with the real part of approximately C<1.571> +and the imaginary part of approximately C<-1.317>. + +=head1 PLANE ANGLE CONVERSIONS + +(Plane, 2-dimensional) angles may be converted with the following functions. + +=over + +=item deg2rad + + $radians = deg2rad($degrees); + +=item grad2rad + + $radians = grad2rad($gradians); + +=item rad2deg + + $degrees = rad2deg($radians); + +=item grad2deg + + $degrees = grad2deg($gradians); + +=item deg2grad + + $gradians = deg2grad($degrees); + +=item rad2grad + + $gradians = rad2grad($radians); + +=back + +The full circle is 2 I<pi> radians or I<360> degrees or I<400> gradians. +The result is by default wrapped to be inside the [0, {2pi,360,400}[ circle. +If you don't want this, supply a true second argument: + + $zillions_of_radians = deg2rad($zillions_of_degrees, 1); + $negative_degrees = rad2deg($negative_radians, 1); + +You can also do the wrapping explicitly by rad2rad(), deg2deg(), and +grad2grad(). + +=over 4 + +=item rad2rad + + $radians_wrapped_by_2pi = rad2rad($radians); + +=item deg2deg + + $degrees_wrapped_by_360 = deg2deg($degrees); + +=item grad2grad + + $gradians_wrapped_by_400 = grad2grad($gradians); + +=back + +=head1 RADIAL COORDINATE CONVERSIONS + +B<Radial coordinate systems> are the B<spherical> and the B<cylindrical> +systems, explained shortly in more detail. + +You can import radial coordinate conversion functions by using the +C<:radial> tag: + + use Math::Trig ':radial'; + + ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z); + ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z); + ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z); + ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z); + ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi); + ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi); + +B<All angles are in radians>. + +=head2 COORDINATE SYSTEMS + +B<Cartesian> coordinates are the usual rectangular I<(x, y, z)>-coordinates. + +Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional +coordinates which define a point in three-dimensional space. They are +based on a sphere surface. The radius of the sphere is B<rho>, also +known as the I<radial> coordinate. The angle in the I<xy>-plane +(around the I<z>-axis) is B<theta>, also known as the I<azimuthal> +coordinate. The angle from the I<z>-axis is B<phi>, also known as the +I<polar> coordinate. The North Pole is therefore I<0, 0, rho>, and +the Gulf of Guinea (think of the missing big chunk of Africa) I<0, +pi/2, rho>. In geographical terms I<phi> is latitude (northward +positive, southward negative) and I<theta> is longitude (eastward +positive, westward negative). + +B<BEWARE>: some texts define I<theta> and I<phi> the other way round, +some texts define the I<phi> to start from the horizontal plane, some +texts use I<r> in place of I<rho>. + +Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional +coordinates which define a point in three-dimensional space. They are +based on a cylinder surface. The radius of the cylinder is B<rho>, +also known as the I<radial> coordinate. The angle in the I<xy>-plane +(around the I<z>-axis) is B<theta>, also known as the I<azimuthal> +coordinate. The third coordinate is the I<z>, pointing up from the +B<theta>-plane. + +=head2 3-D ANGLE CONVERSIONS + +Conversions to and from spherical and cylindrical coordinates are +available. Please notice that the conversions are not necessarily +reversible because of the equalities like I<pi> angles being equal to +I<-pi> angles. + +=over 4 + +=item cartesian_to_cylindrical + + ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z); + +=item cartesian_to_spherical + + ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z); + +=item cylindrical_to_cartesian + + ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z); + +=item cylindrical_to_spherical + + ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z); + +Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>. + +=item spherical_to_cartesian + + ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi); + +=item spherical_to_cylindrical + + ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi); + +Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>. + +=back + +=head1 GREAT CIRCLE DISTANCES AND DIRECTIONS + +A great circle is section of a circle that contains the circle +diameter: the shortest distance between two (non-antipodal) points on +the spherical surface goes along the great circle connecting those two +points. + +=head2 great_circle_distance + +You can compute spherical distances, called B<great circle distances>, +by importing the great_circle_distance() function: + + use Math::Trig 'great_circle_distance'; + + $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]); + +The I<great circle distance> is the shortest distance between two +points on a sphere. The distance is in C<$rho> units. The C<$rho> is +optional, it defaults to 1 (the unit sphere), therefore the distance +defaults to radians. + +If you think geographically the I<theta> are longitudes: zero at the +Greenwhich meridian, eastward positive, westward negative -- and the +I<phi> are latitudes: zero at the North Pole, northward positive, +southward negative. B<NOTE>: this formula thinks in mathematics, not +geographically: the I<phi> zero is at the North Pole, not at the +Equator on the west coast of Africa (Bay of Guinea). You need to +subtract your geographical coordinates from I<pi/2> (also known as 90 +degrees). + + $distance = great_circle_distance($lon0, pi/2 - $lat0, + $lon1, pi/2 - $lat1, $rho); + +=head2 great_circle_direction + +The direction you must follow the great circle (also known as I<bearing>) +can be computed by the great_circle_direction() function: + + use Math::Trig 'great_circle_direction'; + + $direction = great_circle_direction($theta0, $phi0, $theta1, $phi1); + +=head2 great_circle_bearing + +Alias 'great_circle_bearing' for 'great_circle_direction' is also available. + + use Math::Trig 'great_circle_bearing'; + + $direction = great_circle_bearing($theta0, $phi0, $theta1, $phi1); + +The result of great_circle_direction is in radians, zero indicating +straight north, pi or -pi straight south, pi/2 straight west, and +-pi/2 straight east. + +=head2 great_circle_destination + +You can inversely compute the destination if you know the +starting point, direction, and distance: + + use Math::Trig 'great_circle_destination'; + + # $diro is the original direction, + # for example from great_circle_bearing(). + # $distance is the angular distance in radians, + # for example from great_circle_distance(). + # $thetad and $phid are the destination coordinates, + # $dird is the final direction at the destination. + + ($thetad, $phid, $dird) = + great_circle_destination($theta, $phi, $diro, $distance); + +or the midpoint if you know the end points: + +=head2 great_circle_midpoint + + use Math::Trig 'great_circle_midpoint'; + + ($thetam, $phim) = + great_circle_midpoint($theta0, $phi0, $theta1, $phi1); + +The great_circle_midpoint() is just a special case of + +=head2 great_circle_waypoint + + use Math::Trig 'great_circle_waypoint'; + + ($thetai, $phii) = + great_circle_waypoint($theta0, $phi0, $theta1, $phi1, $way); + +Where the $way is a value from zero ($theta0, $phi0) to one ($theta1, +$phi1). Note that antipodal points (where their distance is I<pi> +radians) do not have waypoints between them (they would have an an +"equator" between them), and therefore C<undef> is returned for +antipodal points. If the points are the same and the distance +therefore zero and all waypoints therefore identical, the first point +(either point) is returned. + +The thetas, phis, direction, and distance in the above are all in radians. + +You can import all the great circle formulas by + + use Math::Trig ':great_circle'; + +Notice that the resulting directions might be somewhat surprising if +you are looking at a flat worldmap: in such map projections the great +circles quite often do not look like the shortest routes -- but for +example the shortest possible routes from Europe or North America to +Asia do often cross the polar regions. (The common Mercator projection +does B<not> show great circles as straight lines: straight lines in the +Mercator projection are lines of constant bearing.) + +=head1 EXAMPLES + +To calculate the distance between London (51.3N 0.5W) and Tokyo +(35.7N 139.8E) in kilometers: + + use Math::Trig qw(great_circle_distance deg2rad); + + # Notice the 90 - latitude: phi zero is at the North Pole. + sub NESW { deg2rad($_[0]), deg2rad(90 - $_[1]) } + my @L = NESW( -0.5, 51.3); + my @T = NESW(139.8, 35.7); + my $km = great_circle_distance(@L, @T, 6378); # About 9600 km. + +The direction you would have to go from London to Tokyo (in radians, +straight north being zero, straight east being pi/2). + + use Math::Trig qw(great_circle_direction); + + my $rad = great_circle_direction(@L, @T); # About 0.547 or 0.174 pi. + +The midpoint between London and Tokyo being + + use Math::Trig qw(great_circle_midpoint); + + my @M = great_circle_midpoint(@L, @T); + +or about 69 N 89 E, in the frozen wastes of Siberia. + +B<NOTE>: you B<cannot> get from A to B like this: + + Dist = great_circle_distance(A, B) + Dir = great_circle_direction(A, B) + C = great_circle_destination(A, Dist, Dir) + +and expect C to be B, because the bearing constantly changes when +going from A to B (except in some special case like the meridians or +the circles of latitudes) and in great_circle_destination() one gives +a B<constant> bearing to follow. + +=head2 CAVEAT FOR GREAT CIRCLE FORMULAS + +The answers may be off by few percentages because of the irregular +(slightly aspherical) form of the Earth. The errors are at worst +about 0.55%, but generally below 0.3%. + +=head2 Real-valued asin and acos + +For small inputs asin() and acos() may return complex numbers even +when real numbers would be enough and correct, this happens because of +floating-point inaccuracies. You can see these inaccuracies for +example by trying theses: + + print cos(1e-6)**2+sin(1e-6)**2 - 1,"\n"; + printf "%.20f", cos(1e-6)**2+sin(1e-6)**2,"\n"; + +which will print something like this + + -1.11022302462516e-16 + 0.99999999999999988898 + +even though the expected results are of course exactly zero and one. +The formulas used to compute asin() and acos() are quite sensitive to +this, and therefore they might accidentally slip into the complex +plane even when they should not. To counter this there are two +interfaces that are guaranteed to return a real-valued output. + +=over 4 + +=item asin_real + + use Math::Trig qw(asin_real); + + $real_angle = asin_real($input_sin); + +Return a real-valued arcus sine if the input is between [-1, 1], +B<inclusive> the endpoints. For inputs greater than one, pi/2 +is returned. For inputs less than minus one, -pi/2 is returned. + +=item acos_real + + use Math::Trig qw(acos_real); + + $real_angle = acos_real($input_cos); + +Return a real-valued arcus cosine if the input is between [-1, 1], +B<inclusive> the endpoints. For inputs greater than one, zero +is returned. For inputs less than minus one, pi is returned. + +=back + +=head1 BUGS + +Saying C<use Math::Trig;> exports many mathematical routines in the +caller environment and even overrides some (C<sin>, C<cos>). This is +construed as a feature by the Authors, actually... ;-) + +The code is not optimized for speed, especially because we use +C<Math::Complex> and thus go quite near complex numbers while doing +the computations even when the arguments are not. This, however, +cannot be completely avoided if we want things like C<asin(2)> to give +an answer instead of giving a fatal runtime error. + +Do not attempt navigation using these formulas. + +L<Math::Complex> + +=head1 AUTHORS + +Jarkko Hietaniemi <F<jhi!at!iki.fi>>, +Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>, +Zefram <zefram@fysh.org> + +=head1 LICENSE + +This library is free software; you can redistribute it and/or modify +it under the same terms as Perl itself. + +=cut + +# eof diff --git a/dist/Math-Complex/t/Complex.t b/dist/Math-Complex/t/Complex.t new file mode 100644 index 0000000000..c4fd96f8bd --- /dev/null +++ b/dist/Math-Complex/t/Complex.t @@ -0,0 +1,1160 @@ +#!./perl + +# +# Regression tests for the Math::Complex pacakge +# -- Raphael Manfredi since Sep 1996 +# -- Jarkko Hietaniemi since Mar 1997 +# -- Daniel S. Lewart since Sep 1997 + +use strict; +use warnings; + +use Math::Complex 1.54; + +# they are used later in the test and not exported by Math::Complex +*_stringify_cartesian = \&Math::Complex::_stringify_cartesian; +*_stringify_polar = \&Math::Complex::_stringify_polar; + +our $vax_float = (pack("d",1) =~ /^[\x80\x10]\x40/); +our $has_inf = !$vax_float; + +my ($args, $op, $target, $test, $test_set, $try, $val, $zvalue, @set, @val); +my ($bad, $z); + +$test = 0; +$| = 1; +my @script = ( + 'my ($res, $s0,$s1,$s2,$s3,$s4,$s5,$s6,$s7,$s8,$s9,$s10,$z0,$z1,$z2);' . + "\n\n" +); +my $eps = 1e-13; + +if ($^O eq 'unicos') { # For some reason root() produces very inaccurate + $eps = 1e-10; # results in Cray UNICOS, and occasionally also +} # cos(), sin(), cosh(), sinh(). The division + # of doubles is the current suspect. + +$test++; +push @script, "{ my \$t=$test; ".q{ + my $a = Math::Complex->new(1); + my $b = $a; + $a += 2; + print "not " unless "$a" eq "3" && "$b" eq "1"; + print "ok $t\n"; +}."}"; + +while (<DATA>) { + s/^\s+//; + next if $_ eq '' || /^\#/; + chomp; + $test_set = 0; # Assume not a test over a set of values + if (/^&(.+)/) { + $op = $1; + next; + } + elsif (/^\{(.+)\}/) { + set($1, \@set, \@val); + next; + } + elsif (s/^\|//) { + $test_set = 1; # Requests we loop over the set... + } + my @args = split(/:/); + if ($test_set == 1) { + my $i; + for ($i = 0; $i < @set; $i++) { + # complex number + $target = $set[$i]; + # textual value as found in set definition + $zvalue = $val[$i]; + test($zvalue, $target, @args); + } + } else { + test($op, undef, @args); + } +} + +# + +sub test_mutators { + my $op; + + $test++; +push(@script, <<'EOT'); +{ + my $z = cplx( 1, 1); + $z->Re(2); + $z->Im(3); + print "# $test Re(z) = ",$z->Re(), " Im(z) = ", $z->Im(), " z = $z\n"; + print 'not ' unless Re($z) == 2 and Im($z) == 3; +EOT + push(@script, qq(print "ok $test\\n"}\n)); + + $test++; +push(@script, <<'EOT'); +{ + my $z = cplx( 1, 1); + $z->abs(3 * sqrt(2)); + print "# $test Re(z) = ",$z->Re(), " Im(z) = ", $z->Im(), " z = $z\n"; + print 'not ' unless (abs($z) - 3 * sqrt(2)) < $eps and + (arg($z) - pi / 4 ) < $eps and + (Re($z) - 3 ) < $eps and + (Im($z) - 3 ) < $eps; +EOT + push(@script, qq(print "ok $test\\n"}\n)); + + $test++; +push(@script, <<'EOT'); +{ + my $z = cplx( 1, 1); + $z->arg(-3 / 4 * pi); + print "# $test Re(z) = ",$z->Re(), " Im(z) = ", $z->Im(), " z = $z\n"; + print 'not ' unless (arg($z) + 3 / 4 * pi) < $eps and + (abs($z) - sqrt(2) ) < $eps and + (Re($z) + 1 ) < $eps and + (Im($z) + 1 ) < $eps; +EOT + push(@script, qq(print "ok $test\\n"}\n)); +} + +test_mutators(); + +my $constants = ' +my $i = cplx(0, 1); +my $pi = cplx(pi, 0); +my $pii = cplx(0, pi); +my $pip2 = cplx(pi/2, 0); +my $pip4 = cplx(pi/4, 0); +my $zero = cplx(0, 0); +'; + +if ($has_inf) { + $constants .= <<'EOF'; +my $inf = 9**9**9; +EOF +} + +push(@script, $constants); + + +# test the divbyzeros + +sub test_dbz { + for my $op (@_) { + $test++; + push(@script, <<EOT); + eval '$op'; + (\$bad) = (\$@ =~ /(.+)/); + print "# $test op = $op divbyzero? \$bad...\n"; + print 'not ' unless (\$@ =~ /Division by zero/); +EOT + push(@script, qq(print "ok $test\\n";\n)); + } +} + +# test the logofzeros + +sub test_loz { + for my $op (@_) { + $test++; + push(@script, <<EOT); + eval '$op'; + (\$bad) = (\$@ =~ /(.+)/); + print "# $test op = $op logofzero? \$bad...\n"; + print 'not ' unless (\$@ =~ /Logarithm of zero/); +EOT + push(@script, qq(print "ok $test\\n";\n)); + } +} + +test_dbz( + 'i/0', + 'acot(0)', + 'acot(+$i)', +# 'acoth(-1)', # Log of zero. + 'acoth(0)', + 'acoth(+1)', + 'acsc(0)', + 'acsch(0)', + 'asec(0)', + 'asech(0)', + 'atan($i)', +# 'atanh(-1)', # Log of zero. + 'atanh(+1)', + 'cot(0)', + 'coth(0)', + 'csc(0)', + 'csch(0)', + 'atan(cplx(0, 1), cplx(1, 0))', + ); + +test_loz( + 'log($zero)', + 'atan(-$i)', + 'acot(-$i)', + 'atanh(-1)', + 'acoth(-1)', + ); + +# test the bad roots + +sub test_broot { + for my $op (@_) { + $test++; + push(@script, <<EOT); + eval 'root(2, $op)'; + (\$bad) = (\$@ =~ /(.+)/); + print "# $test op = $op badroot? \$bad...\n"; + print 'not ' unless (\$@ =~ /root rank must be/); +EOT + push(@script, qq(print "ok $test\\n";\n)); + } +} + +test_broot(qw(-3 -2.1 0 0.99)); + +sub test_display_format { + $test++; + push @script, <<EOS; + print "# package display_format cartesian?\n"; + print "not " unless Math::Complex->display_format eq 'cartesian'; + print "ok $test\n"; +EOS + + push @script, <<EOS; + my \$j = (root(1,3))[1]; + + \$j->display_format('polar'); +EOS + + $test++; + push @script, <<EOS; + print "# j display_format polar?\n"; + print "not " unless \$j->display_format eq 'polar'; + print "ok $test\n"; +EOS + + $test++; + push @script, <<EOS; + print "# j = \$j\n"; + print "not " unless "\$j" eq "[1,2pi/3]"; + print "ok $test\n"; + + my %display_format; + + %display_format = \$j->display_format; +EOS + + $test++; + push @script, <<EOS; + print "# display_format{style} polar?\n"; + print "not " unless \$display_format{style} eq 'polar'; + print "ok $test\n"; +EOS + + $test++; + push @script, <<EOS; + print "# keys %display_format == 2?\n"; + print "not " unless keys %display_format == 2; + print "ok $test\n"; + + \$j->display_format('style' => 'cartesian', 'format' => '%.5f'); +EOS + + $test++; + push @script, <<EOS; + print "# j = \$j\n"; + print "not " unless "\$j" eq "-0.50000+0.86603i"; + print "ok $test\n"; + + %display_format = \$j->display_format; +EOS + + $test++; + push @script, <<EOS; + print "# display_format{format} %.5f?\n"; + print "not " unless \$display_format{format} eq '%.5f'; + print "ok $test\n"; +EOS + + $test++; + push @script, <<EOS; + print "# keys %display_format == 3?\n"; + print "not " unless keys %display_format == 3; + print "ok $test\n"; + + \$j->display_format('format' => undef); +EOS + + $test++; + push @script, <<EOS; + print "# j = \$j\n"; + print "not " unless "\$j" =~ /^-0(?:\\.5(?:0000\\d+)?|\\.49999\\d+)\\+0.86602540\\d+i\$/; + print "ok $test\n"; + + \$j->display_format('style' => 'polar', 'polar_pretty_print' => 0); +EOS + + $test++; + push @script, <<EOS; + print "# j = \$j\n"; + print "not " unless "\$j" =~ /^\\[1,2\\.09439510\\d+\\]\$/; + print "ok $test\n"; + + \$j->display_format('style' => 'polar', 'format' => "%.4g"); +EOS + + $test++; + push @script, <<EOS; + print "# j = \$j\n"; + print "not " unless "\$j" =~ /^\\[1,2\\.094\\]\$/; + print "ok $test\n"; + + \$j->display_format('style' => 'cartesian', 'format' => '(%.5g)'); +EOS + + $test++; + push @script, <<EOS; + print "# j = \$j\n"; + print "not " unless "\$j" eq "(-0.5)+(0.86603)i"; + print "ok $test\n"; +EOS + + $test++; + push @script, <<EOS; + print "# j display_format cartesian?\n"; + print "not " unless \$j->display_format eq 'cartesian'; + print "ok $test\n"; +EOS +} + +test_display_format(); + +sub test_remake { + $test++; + push @script, <<EOS; + print "# remake 2+3i\n"; + \$z = cplx('2+3i'); + print "not " unless \$z == Math::Complex->make(2,3); + print "ok $test\n"; +EOS + + $test++; + push @script, <<EOS; + print "# make 3i\n"; + \$z = Math::Complex->make('3i'); + print "not " unless \$z == cplx(0,3); + print "ok $test\n"; +EOS + + $test++; + push @script, <<EOS; + print "# emake [2,3]\n"; + \$z = Math::Complex->emake('[2,3]'); + print "not " unless \$z == cplxe(2,3); + print "ok $test\n"; +EOS + + $test++; + push @script, <<EOS; + print "# make (2,3)\n"; + \$z = Math::Complex->make('(2,3)'); + print "not " unless \$z == cplx(2,3); + print "ok $test\n"; +EOS + + $test++; + push @script, <<EOS; + print "# emake [2,3pi/8]\n"; + \$z = Math::Complex->emake('[2,3pi/8]'); + print "not " unless \$z == cplxe(2,3*\$pi/8); + print "ok $test\n"; +EOS + + $test++; + push @script, <<EOS; + print "# emake [2]\n"; + \$z = Math::Complex->emake('[2]'); + print "not " unless \$z == cplxe(2); + print "ok $test\n"; +EOS +} + +sub test_no_args { + push @script, <<'EOS'; +{ + print "# cplx, cplxe, make, emake without arguments\n"; +EOS + + $test++; + push @script, <<EOS; + my \$z0 = cplx(); + print ((\$z0->Re() == 0) ? "ok $test\n" : "not ok $test\n"); +EOS + + $test++; + push @script, <<EOS; + print ((\$z0->Im() == 0) ? "ok $test\n" : "not ok $test\n"); +EOS + + $test++; + push @script, <<EOS; + my \$z1 = cplxe(); + print ((\$z1->rho() == 0) ? "ok $test\n" : "not ok $test\n"); +EOS + + $test++; + push @script, <<EOS; + print ((\$z1->theta() == 0) ? "ok $test\n" : "not ok $test\n"); +EOS + + $test++; + push @script, <<EOS; + my \$z2 = Math::Complex->make(); + print ((\$z2->Re() == 0) ? "ok $test\n" : "not ok $test\n"); +EOS + + $test++; + push @script, <<EOS; + print ((\$z2->Im() == 0) ? "ok $test\n" : "not ok $test\n"); +EOS + + $test++; + push @script, <<EOS; + my \$z3 = Math::Complex->emake(); + print ((\$z3->rho() == 0) ? "ok $test\n" : "not ok $test\n"); +EOS + + $test++; + push @script, <<EOS; + print ((\$z3->theta() == 0) ? "ok $test\n" : "not ok $test\n"); +} +EOS +} + +sub test_atan2 { + push @script, <<'EOS'; +print "# atan2() with some real arguments\n"; +EOS + my @real = (-1, 0, 1); + for my $x (@real) { + for my $y (@real) { + next if $x == 0 && $y == 0; + $test++; + push @script, <<EOS; +print ((Math::Complex::atan2($y, $x) == CORE::atan2($y, $x)) ? "ok $test\n" : "not ok $test\n"); +EOS + } + } + push @script, <<'EOS'; + print "# atan2() with some complex arguments\n"; +EOS + $test++; + push @script, <<EOS; + print (abs(atan2(0, cplx(0, 1))) < $eps ? "ok $test\n" : "not ok $test\n"); +EOS + $test++; + push @script, <<EOS; + print (abs(atan2(cplx(0, 1), 0) - \$pip2) < $eps ? "ok $test\n" : "not ok $test\n"); +EOS + $test++; + push @script, <<EOS; + print (abs(atan2(cplx(0, 1), cplx(0, 1)) - \$pip4) < $eps ? "ok $test\n" : "not ok $test\n"); +EOS + $test++; + push @script, <<EOS; + print (abs(atan2(cplx(0, 1), cplx(1, 1)) - cplx(0.553574358897045, 0.402359478108525)) < $eps ? "ok $test\n" : "not ok $test\n"); +EOS +} + +sub test_decplx { +} + +test_remake(); + +test_no_args(); + +test_atan2(); + +test_decplx(); + +print "1..$test\n"; +#print @script, "\n"; +eval join '', @script; +die $@ if $@; + +sub abop { + my ($op) = @_; + + push(@script, qq(print "# $op=\n";)); +} + +sub test { + my ($op, $z, @args) = @_; + my ($baop) = 0; + $test++; + my $i; + $baop = 1 if ($op =~ s/;=$//); + for ($i = 0; $i < @args; $i++) { + $val = value($args[$i]); + push @script, "\$z$i = $val;\n"; + } + if (defined $z) { + $args = "'$op'"; # Really the value + $try = "abs(\$z0 - \$z1) <= $eps ? \$z1 : \$z0"; + push @script, "\$res = $try; "; + push @script, "check($test, $args[0], \$res, \$z$#args, $args);\n"; + } else { + my ($try, $args); + if (@args == 2) { + $try = "$op \$z0"; + $args = "'$args[0]'"; + } else { + $try = ($op =~ /^\w/) ? "$op(\$z0, \$z1)" : "\$z0 $op \$z1"; + $args = "'$args[0]', '$args[1]'"; + } + push @script, "\$res = $try; "; + push @script, "check($test, '$try', \$res, \$z$#args, $args);\n"; + if (@args > 2 and $baop) { # binary assignment ops + $test++; + # check the op= works + push @script, <<EOB; +{ + my \$za = cplx(ref \$z0 ? \@{\$z0->_cartesian} : (\$z0, 0)); + + my (\$z1r, \$z1i) = ref \$z1 ? \@{\$z1->_cartesian} : (\$z1, 0); + + my \$zb = cplx(\$z1r, \$z1i); + + \$za $op= \$zb; + my (\$zbr, \$zbi) = \@{\$zb->_cartesian}; + + check($test, '\$z0 $op= \$z1', \$za, \$z$#args, $args); +EOB + $test++; + # check that the rhs has not changed + push @script, qq(print "not " unless (\$zbr == \$z1r and \$zbi == \$z1i);); + push @script, qq(print "ok $test\\n";\n); + push @script, "}\n"; + } + } +} + +sub set { + my ($set, $setref, $valref) = @_; + @{$setref} = (); + @{$valref} = (); + my @set = split(/;\s*/, $set); + my @res; + my $i; + for ($i = 0; $i < @set; $i++) { + push(@{$valref}, $set[$i]); + my $val = value($set[$i]); + push @script, "\$s$i = $val;\n"; + push @{$setref}, "\$s$i"; + } +} + +sub value { + local ($_) = @_; + if (/^\s*\((.*),(.*)\)/) { + return "cplx($1,$2)"; + } + elsif (/^\s*([\-\+]?(?:\d+(\.\d+)?|\.\d+)(?:[e[\-\+]\d+])?)/) { + return "cplx($1,0)"; + } + elsif (/^\s*\[(.*),(.*)\]/) { + return "cplxe($1,$2)"; + } + elsif (/^\s*'(.*)'/) { + my $ex = $1; + $ex =~ s/\bz\b/$target/g; + $ex =~ s/\br\b/abs($target)/g; + $ex =~ s/\bt\b/arg($target)/g; + $ex =~ s/\ba\b/Re($target)/g; + $ex =~ s/\bb\b/Im($target)/g; + return $ex; + } + elsif (/^\s*"(.*)"/) { + return "\"$1\""; + } + return $_; +} + +sub check { + my ($test, $try, $got, $expected, @z) = @_; + + print "# @_\n"; + + if ("$got" eq "$expected" + || + ($expected =~ /^-?\d/ && $got == $expected) + || + (abs(Math::Complex->make($got) - Math::Complex->make($expected)) < $eps) + || + (abs($got - $expected) < $eps) + ) { + print "ok $test\n"; + } else { + print "not ok $test\n"; + my $args = (@z == 1) ? "z = $z[0]" : "z0 = $z[0], z1 = $z[1]"; + print "# '$try' expected: '$expected' got: '$got' for $args\n"; + } +} + +sub addsq { + my ($z1, $z2) = @_; + return ($z1 + i*$z2) * ($z1 - i*$z2); +} + +sub subsq { + my ($z1, $z2) = @_; + return ($z1 + $z2) * ($z1 - $z2); +} + +__END__ +&+;= +(3,4):(3,4):(6,8) +(-3,4):(3,-4):(0,0) +(3,4):-3:(0,4) +1:(4,2):(5,2) +[2,0]:[2,pi]:(0,0) + +&++ +(2,1):(3,1) + +&-;= +(2,3):(-2,-3) +[2,pi/2]:[2,-(pi)/2] +2:[2,0]:(0,0) +[3,0]:2:(1,0) +3:(4,5):(-1,-5) +(4,5):3:(1,5) +(2,1):(3,5):(-1,-4) + +&-- +(1,2):(0,2) +[2,pi]:[3,pi] + +&*;= +(0,1):(0,1):(-1,0) +(4,5):(1,0):(4,5) +[2,2*pi/3]:(1,0):[2,2*pi/3] +2:(0,1):(0,2) +(0,1):3:(0,3) +(0,1):(4,1):(-1,4) +(2,1):(4,-1):(9,2) + +&/;= +(3,4):(3,4):(1,0) +(4,-5):1:(4,-5) +1:(0,1):(0,-1) +(0,6):(0,2):(3,0) +(9,2):(4,-1):(2,1) +[4,pi]:[2,pi/2]:[2,pi/2] +[2,pi/2]:[4,pi]:[0.5,-(pi)/2] + +&**;= +(2,0):(3,0):(8,0) +(3,0):(2,0):(9,0) +(2,3):(4,0):(-119,-120) +(0,0):(1,0):(0,0) +(0,0):(2,3):(0,0) +(1,0):(0,0):(1,0) +(1,0):(1,0):(1,0) +(1,0):(2,3):(1,0) +(2,3):(0,0):(1,0) +(2,3):(1,0):(2,3) +(0,0):(0,0):(1,0) + +&Re +(3,4):3 +(-3,4):-3 +[1,pi/2]:0 + +&Im +(3,4):4 +(3,-4):-4 +[1,pi/2]:1 + +&abs +(3,4):5 +(-3,4):5 + +&arg +[2,0]:0 +[-2,0]:pi + +&~ +(4,5):(4,-5) +(-3,4):(-3,-4) +[2,pi/2]:[2,-(pi)/2] + +&< +(3,4):(1,2):0 +(3,4):(3,2):0 +(3,4):(3,8):1 +(4,4):(5,129):1 + +&== +(3,4):(4,5):0 +(3,4):(3,5):0 +(3,4):(2,4):0 +(3,4):(3,4):1 + +&sqrt +-9:(0,3) +(-100,0):(0,10) +(16,-30):(5,-3) + +&_stringify_cartesian +(-100,0):"-100" +(0,1):"i" +(4,-3):"4-3i" +(4,0):"4" +(-4,0):"-4" +(-2,4):"-2+4i" +(-2,-1):"-2-i" + +&_stringify_polar +[-1, 0]:"[1,pi]" +[1, pi/3]:"[1,pi/3]" +[6, -2*pi/3]:"[6,-2pi/3]" +[0.5, -9*pi/11]:"[0.5,-9pi/11]" +[1, 0.5]:"[1, 0.5]" + +{ (4,3); [3,2]; (-3,4); (0,2); [2,1] } + +|'z + ~z':'2*Re(z)' +|'z - ~z':'2*i*Im(z)' +|'z * ~z':'abs(z) * abs(z)' + +{ (0.5, 0); (-0.5, 0); (2,3); [3,2]; (-3,2); (0,2); 3; 1.2; (-3, 0); (-2, -1); [2,1] } + +|'(root(z, 4))[1] ** 4':'z' +|'(root(z, 5))[3] ** 5':'z' +|'(root(z, 8))[7] ** 8':'z' +|'(root(z, 8, 0)) ** 8':'z' +|'(root(z, 8, 7)) ** 8':'z' +|'abs(z)':'r' +|'acot(z)':'acotan(z)' +|'acsc(z)':'acosec(z)' +|'acsc(z)':'asin(1 / z)' +|'asec(z)':'acos(1 / z)' +|'cbrt(z)':'cbrt(r) * exp(i * t/3)' +|'cos(acos(z))':'z' +|'addsq(cos(z), sin(z))':1 +|'cos(z)':'cosh(i*z)' +|'subsq(cosh(z), sinh(z))':1 +|'cot(acot(z))':'z' +|'cot(z)':'1 / tan(z)' +|'cot(z)':'cotan(z)' +|'csc(acsc(z))':'z' +|'csc(z)':'1 / sin(z)' +|'csc(z)':'cosec(z)' +|'exp(log(z))':'z' +|'exp(z)':'exp(a) * exp(i * b)' +|'ln(z)':'log(z)' +|'log(exp(z))':'z' +|'log(z)':'log(r) + i*t' +|'log10(z)':'log(z) / log(10)' +|'logn(z, 2)':'log(z) / log(2)' +|'logn(z, 3)':'log(z) / log(3)' +|'sec(asec(z))':'z' +|'sec(z)':'1 / cos(z)' +|'sin(asin(z))':'z' +|'sin(i * z)':'i * sinh(z)' +|'sqrt(z) * sqrt(z)':'z' +|'sqrt(z)':'sqrt(r) * exp(i * t/2)' +|'tan(atan(z))':'z' +|'z**z':'exp(z * log(z))' + +{ (1,1); [1,0.5]; (-2, -1); 2; -3; (-1,0.5); (0,0.5); 0.5; (2, 0); (-1, -2) } + +|'cosh(acosh(z))':'z' +|'coth(acoth(z))':'z' +|'coth(z)':'1 / tanh(z)' +|'coth(z)':'cotanh(z)' +|'csch(acsch(z))':'z' +|'csch(z)':'1 / sinh(z)' +|'csch(z)':'cosech(z)' +|'sech(asech(z))':'z' +|'sech(z)':'1 / cosh(z)' +|'sinh(asinh(z))':'z' +|'tanh(atanh(z))':'z' + +{ (0.2,-0.4); [1,0.5]; -1.2; (-1,0.5); 0.5; (1.1, 0) } + +|'acos(cos(z)) ** 2':'z * z' +|'acosh(cosh(z)) ** 2':'z * z' +|'acoth(z)':'acotanh(z)' +|'acoth(z)':'atanh(1 / z)' +|'acsch(z)':'acosech(z)' +|'acsch(z)':'asinh(1 / z)' +|'asech(z)':'acosh(1 / z)' +|'asin(sin(z))':'z' +|'asinh(sinh(z))':'z' +|'atan(tan(z))':'z' +|'atanh(tanh(z))':'z' + +&log +(-2.0,0):( 0.69314718055995, 3.14159265358979) +(-1.0,0):( 0 , 3.14159265358979) +(-0.5,0):( -0.69314718055995, 3.14159265358979) +( 0.5,0):( -0.69314718055995, 0 ) +( 1.0,0):( 0 , 0 ) +( 2.0,0):( 0.69314718055995, 0 ) + +&log +( 2, 3):( 1.28247467873077, 0.98279372324733) +(-2, 3):( 1.28247467873077, 2.15879893034246) +(-2,-3):( 1.28247467873077, -2.15879893034246) +( 2,-3):( 1.28247467873077, -0.98279372324733) + +&sin +(-2.0,0):( -0.90929742682568, 0 ) +(-1.0,0):( -0.84147098480790, 0 ) +(-0.5,0):( -0.47942553860420, 0 ) +( 0.0,0):( 0 , 0 ) +( 0.5,0):( 0.47942553860420, 0 ) +( 1.0,0):( 0.84147098480790, 0 ) +( 2.0,0):( 0.90929742682568, 0 ) + +&sin +( 2, 3):( 9.15449914691143, -4.16890695996656) +(-2, 3):( -9.15449914691143, -4.16890695996656) +(-2,-3):( -9.15449914691143, 4.16890695996656) +( 2,-3):( 9.15449914691143, 4.16890695996656) + +&cos +(-2.0,0):( -0.41614683654714, 0 ) +(-1.0,0):( 0.54030230586814, 0 ) +(-0.5,0):( 0.87758256189037, 0 ) +( 0.0,0):( 1 , 0 ) +( 0.5,0):( 0.87758256189037, 0 ) +( 1.0,0):( 0.54030230586814, 0 ) +( 2.0,0):( -0.41614683654714, 0 ) + +&cos +( 2, 3):( -4.18962569096881, -9.10922789375534) +(-2, 3):( -4.18962569096881, 9.10922789375534) +(-2,-3):( -4.18962569096881, -9.10922789375534) +( 2,-3):( -4.18962569096881, 9.10922789375534) + +&tan +(-2.0,0):( 2.18503986326152, 0 ) +(-1.0,0):( -1.55740772465490, 0 ) +(-0.5,0):( -0.54630248984379, 0 ) +( 0.0,0):( 0 , 0 ) +( 0.5,0):( 0.54630248984379, 0 ) +( 1.0,0):( 1.55740772465490, 0 ) +( 2.0,0):( -2.18503986326152, 0 ) + +&tan +( 2, 3):( -0.00376402564150, 1.00323862735361) +(-2, 3):( 0.00376402564150, 1.00323862735361) +(-2,-3):( 0.00376402564150, -1.00323862735361) +( 2,-3):( -0.00376402564150, -1.00323862735361) + +&sec +(-2.0,0):( -2.40299796172238, 0 ) +(-1.0,0):( 1.85081571768093, 0 ) +(-0.5,0):( 1.13949392732455, 0 ) +( 0.0,0):( 1 , 0 ) +( 0.5,0):( 1.13949392732455, 0 ) +( 1.0,0):( 1.85081571768093, 0 ) +( 2.0,0):( -2.40299796172238, 0 ) + +&sec +( 2, 3):( -0.04167496441114, 0.09061113719624) +(-2, 3):( -0.04167496441114, -0.09061113719624) +(-2,-3):( -0.04167496441114, 0.09061113719624) +( 2,-3):( -0.04167496441114, -0.09061113719624) + +&csc +(-2.0,0):( -1.09975017029462, 0 ) +(-1.0,0):( -1.18839510577812, 0 ) +(-0.5,0):( -2.08582964293349, 0 ) +( 0.5,0):( 2.08582964293349, 0 ) +( 1.0,0):( 1.18839510577812, 0 ) +( 2.0,0):( 1.09975017029462, 0 ) + +&csc +( 2, 3):( 0.09047320975321, 0.04120098628857) +(-2, 3):( -0.09047320975321, 0.04120098628857) +(-2,-3):( -0.09047320975321, -0.04120098628857) +( 2,-3):( 0.09047320975321, -0.04120098628857) + +&cot +(-2.0,0):( 0.45765755436029, 0 ) +(-1.0,0):( -0.64209261593433, 0 ) +(-0.5,0):( -1.83048772171245, 0 ) +( 0.5,0):( 1.83048772171245, 0 ) +( 1.0,0):( 0.64209261593433, 0 ) +( 2.0,0):( -0.45765755436029, 0 ) + +&cot +( 2, 3):( -0.00373971037634, -0.99675779656936) +(-2, 3):( 0.00373971037634, -0.99675779656936) +(-2,-3):( 0.00373971037634, 0.99675779656936) +( 2,-3):( -0.00373971037634, 0.99675779656936) + +&asin +(-2.0,0):( -1.57079632679490, 1.31695789692482) +(-1.0,0):( -1.57079632679490, 0 ) +(-0.5,0):( -0.52359877559830, 0 ) +( 0.0,0):( 0 , 0 ) +( 0.5,0):( 0.52359877559830, 0 ) +( 1.0,0):( 1.57079632679490, 0 ) +( 2.0,0):( 1.57079632679490, -1.31695789692482) + +&asin +( 2, 3):( 0.57065278432110, 1.98338702991654) +(-2, 3):( -0.57065278432110, 1.98338702991654) +(-2,-3):( -0.57065278432110, -1.98338702991654) +( 2,-3):( 0.57065278432110, -1.98338702991654) + +&acos +(-2.0,0):( 3.14159265358979, -1.31695789692482) +(-1.0,0):( 3.14159265358979, 0 ) +(-0.5,0):( 2.09439510239320, 0 ) +( 0.0,0):( 1.57079632679490, 0 ) +( 0.5,0):( 1.04719755119660, 0 ) +( 1.0,0):( 0 , 0 ) +( 2.0,0):( 0 , 1.31695789692482) + +&acos +( 2, 3):( 1.00014354247380, -1.98338702991654) +(-2, 3):( 2.14144911111600, -1.98338702991654) +(-2,-3):( 2.14144911111600, 1.98338702991654) +( 2,-3):( 1.00014354247380, 1.98338702991654) + +&atan +(-2.0,0):( -1.10714871779409, 0 ) +(-1.0,0):( -0.78539816339745, 0 ) +(-0.5,0):( -0.46364760900081, 0 ) +( 0.0,0):( 0 , 0 ) +( 0.5,0):( 0.46364760900081, 0 ) +( 1.0,0):( 0.78539816339745, 0 ) +( 2.0,0):( 1.10714871779409, 0 ) + +&atan +( 2, 3):( 1.40992104959658, 0.22907268296854) +(-2, 3):( -1.40992104959658, 0.22907268296854) +(-2,-3):( -1.40992104959658, -0.22907268296854) +( 2,-3):( 1.40992104959658, -0.22907268296854) + +&asec +(-2.0,0):( 2.09439510239320, 0 ) +(-1.0,0):( 3.14159265358979, 0 ) +(-0.5,0):( 3.14159265358979, -1.31695789692482) +( 0.5,0):( 0 , 1.31695789692482) +( 1.0,0):( 0 , 0 ) +( 2.0,0):( 1.04719755119660, 0 ) + +&asec +( 2, 3):( 1.42041072246703, 0.23133469857397) +(-2, 3):( 1.72118193112276, 0.23133469857397) +(-2,-3):( 1.72118193112276, -0.23133469857397) +( 2,-3):( 1.42041072246703, -0.23133469857397) + +&acsc +(-2.0,0):( -0.52359877559830, 0 ) +(-1.0,0):( -1.57079632679490, 0 ) +(-0.5,0):( -1.57079632679490, 1.31695789692482) +( 0.5,0):( 1.57079632679490, -1.31695789692482) +( 1.0,0):( 1.57079632679490, 0 ) +( 2.0,0):( 0.52359877559830, 0 ) + +&acsc +( 2, 3):( 0.15038560432786, -0.23133469857397) +(-2, 3):( -0.15038560432786, -0.23133469857397) +(-2,-3):( -0.15038560432786, 0.23133469857397) +( 2,-3):( 0.15038560432786, 0.23133469857397) + +&acot +(-2.0,0):( -0.46364760900081, 0 ) +(-1.0,0):( -0.78539816339745, 0 ) +(-0.5,0):( -1.10714871779409, 0 ) +( 0.5,0):( 1.10714871779409, 0 ) +( 1.0,0):( 0.78539816339745, 0 ) +( 2.0,0):( 0.46364760900081, 0 ) + +&acot +( 2, 3):( 0.16087527719832, -0.22907268296854) +(-2, 3):( -0.16087527719832, -0.22907268296854) +(-2,-3):( -0.16087527719832, 0.22907268296854) +( 2,-3):( 0.16087527719832, 0.22907268296854) + +&sinh +(-2.0,0):( -3.62686040784702, 0 ) +(-1.0,0):( -1.17520119364380, 0 ) +(-0.5,0):( -0.52109530549375, 0 ) +( 0.0,0):( 0 , 0 ) +( 0.5,0):( 0.52109530549375, 0 ) +( 1.0,0):( 1.17520119364380, 0 ) +( 2.0,0):( 3.62686040784702, 0 ) + +&sinh +( 2, 3):( -3.59056458998578, 0.53092108624852) +(-2, 3):( 3.59056458998578, 0.53092108624852) +(-2,-3):( 3.59056458998578, -0.53092108624852) +( 2,-3):( -3.59056458998578, -0.53092108624852) + +&cosh +(-2.0,0):( 3.76219569108363, 0 ) +(-1.0,0):( 1.54308063481524, 0 ) +(-0.5,0):( 1.12762596520638, 0 ) +( 0.0,0):( 1 , 0 ) +( 0.5,0):( 1.12762596520638, 0 ) +( 1.0,0):( 1.54308063481524, 0 ) +( 2.0,0):( 3.76219569108363, 0 ) + +&cosh +( 2, 3):( -3.72454550491532, 0.51182256998738) +(-2, 3):( -3.72454550491532, -0.51182256998738) +(-2,-3):( -3.72454550491532, 0.51182256998738) +( 2,-3):( -3.72454550491532, -0.51182256998738) + +&tanh +(-2.0,0):( -0.96402758007582, 0 ) +(-1.0,0):( -0.76159415595576, 0 ) +(-0.5,0):( -0.46211715726001, 0 ) +( 0.0,0):( 0 , 0 ) +( 0.5,0):( 0.46211715726001, 0 ) +( 1.0,0):( 0.76159415595576, 0 ) +( 2.0,0):( 0.96402758007582, 0 ) + +&tanh +( 2, 3):( 0.96538587902213, -0.00988437503832) +(-2, 3):( -0.96538587902213, -0.00988437503832) +(-2,-3):( -0.96538587902213, 0.00988437503832) +( 2,-3):( 0.96538587902213, 0.00988437503832) + +&sech +(-2.0,0):( 0.26580222883408, 0 ) +(-1.0,0):( 0.64805427366389, 0 ) +(-0.5,0):( 0.88681888397007, 0 ) +( 0.0,0):( 1 , 0 ) +( 0.5,0):( 0.88681888397007, 0 ) +( 1.0,0):( 0.64805427366389, 0 ) +( 2.0,0):( 0.26580222883408, 0 ) + +&sech +( 2, 3):( -0.26351297515839, -0.03621163655877) +(-2, 3):( -0.26351297515839, 0.03621163655877) +(-2,-3):( -0.26351297515839, -0.03621163655877) +( 2,-3):( -0.26351297515839, 0.03621163655877) + +&csch +(-2.0,0):( -0.27572056477178, 0 ) +(-1.0,0):( -0.85091812823932, 0 ) +(-0.5,0):( -1.91903475133494, 0 ) +( 0.5,0):( 1.91903475133494, 0 ) +( 1.0,0):( 0.85091812823932, 0 ) +( 2.0,0):( 0.27572056477178, 0 ) + +&csch +( 2, 3):( -0.27254866146294, -0.04030057885689) +(-2, 3):( 0.27254866146294, -0.04030057885689) +(-2,-3):( 0.27254866146294, 0.04030057885689) +( 2,-3):( -0.27254866146294, 0.04030057885689) + +&coth +(-2.0,0):( -1.03731472072755, 0 ) +(-1.0,0):( -1.31303528549933, 0 ) +(-0.5,0):( -2.16395341373865, 0 ) +( 0.5,0):( 2.16395341373865, 0 ) +( 1.0,0):( 1.31303528549933, 0 ) +( 2.0,0):( 1.03731472072755, 0 ) + +&coth +( 2, 3):( 1.03574663776500, 0.01060478347034) +(-2, 3):( -1.03574663776500, 0.01060478347034) +(-2,-3):( -1.03574663776500, -0.01060478347034) +( 2,-3):( 1.03574663776500, -0.01060478347034) + +&asinh +(-2.0,0):( -1.44363547517881, 0 ) +(-1.0,0):( -0.88137358701954, 0 ) +(-0.5,0):( -0.48121182505960, 0 ) +( 0.0,0):( 0 , 0 ) +( 0.5,0):( 0.48121182505960, 0 ) +( 1.0,0):( 0.88137358701954, 0 ) +( 2.0,0):( 1.44363547517881, 0 ) + +&asinh +( 2, 3):( 1.96863792579310, 0.96465850440760) +(-2, 3):( -1.96863792579310, 0.96465850440761) +(-2,-3):( -1.96863792579310, -0.96465850440761) +( 2,-3):( 1.96863792579310, -0.96465850440760) + +&acosh +(-2.0,0):( 1.31695789692482, 3.14159265358979) +(-1.0,0):( 0, 3.14159265358979) +(-0.5,0):( 0, 2.09439510239320) +( 0.0,0):( 0, 1.57079632679490) +( 0.5,0):( 0, 1.04719755119660) +( 1.0,0):( 0 , 0 ) +( 2.0,0):( 1.31695789692482, 0 ) + +&acosh +( 2, 3):( 1.98338702991654, 1.00014354247380) +(-2, 3):( 1.98338702991653, 2.14144911111600) +(-2,-3):( 1.98338702991653, -2.14144911111600) +( 2,-3):( 1.98338702991654, -1.00014354247380) + +&atanh +(-2.0,0):( -0.54930614433405, 1.57079632679490) +(-0.5,0):( -0.54930614433405, 0 ) +( 0.0,0):( 0 , 0 ) +( 0.5,0):( 0.54930614433405, 0 ) +( 2.0,0):( 0.54930614433405, 1.57079632679490) + +&atanh +( 2, 3):( 0.14694666622553, 1.33897252229449) +(-2, 3):( -0.14694666622553, 1.33897252229449) +(-2,-3):( -0.14694666622553, -1.33897252229449) +( 2,-3):( 0.14694666622553, -1.33897252229449) + +&asech +(-2.0,0):( 0 , 2.09439510239320) +(-1.0,0):( 0 , 3.14159265358979) +(-0.5,0):( 1.31695789692482, 3.14159265358979) +( 0.5,0):( 1.31695789692482, 0 ) +( 1.0,0):( 0 , 0 ) +( 2.0,0):( 0 , 1.04719755119660) + +&asech +( 2, 3):( 0.23133469857397, -1.42041072246703) +(-2, 3):( 0.23133469857397, -1.72118193112276) +(-2,-3):( 0.23133469857397, 1.72118193112276) +( 2,-3):( 0.23133469857397, 1.42041072246703) + +&acsch +(-2.0,0):( -0.48121182505960, 0 ) +(-1.0,0):( -0.88137358701954, 0 ) +(-0.5,0):( -1.44363547517881, 0 ) +( 0.5,0):( 1.44363547517881, 0 ) +( 1.0,0):( 0.88137358701954, 0 ) +( 2.0,0):( 0.48121182505960, 0 ) + +&acsch +( 2, 3):( 0.15735549884499, -0.22996290237721) +(-2, 3):( -0.15735549884499, -0.22996290237721) +(-2,-3):( -0.15735549884499, 0.22996290237721) +( 2,-3):( 0.15735549884499, 0.22996290237721) + +&acoth +(-2.0,0):( -0.54930614433405, 0 ) +(-0.5,0):( -0.54930614433405, 1.57079632679490) +( 0.5,0):( 0.54930614433405, 1.57079632679490) +( 2.0,0):( 0.54930614433405, 0 ) + +&acoth +( 2, 3):( 0.14694666622553, -0.23182380450040) +(-2, 3):( -0.14694666622553, -0.23182380450040) +(-2,-3):( -0.14694666622553, 0.23182380450040) +( 2,-3):( 0.14694666622553, 0.23182380450040) + +# eof diff --git a/dist/Math-Complex/t/Trig.t b/dist/Math-Complex/t/Trig.t new file mode 100644 index 0000000000..a9a12556b6 --- /dev/null +++ b/dist/Math-Complex/t/Trig.t @@ -0,0 +1,387 @@ +#!./perl + +# +# Regression tests for the Math::Trig package +# +# The tests here are quite modest as the Math::Complex tests exercise +# these interfaces quite vigorously. +# +# -- Jarkko Hietaniemi, April 1997 + +use strict; +use warnings; +use Test::More tests => 153; + +use Math::Trig 1.18; +use Math::Trig 1.18 qw(:pi Inf); + +our $vax_float = (pack("d",1) =~ /^[\x80\x10]\x40/); +our $has_inf = !$vax_float; + +my $pip2 = pi / 2; + +use strict; + +our($x, $y, $z); + +my $eps = 1e-11; + +if ($^O eq 'unicos') { # See lib/Math/Complex.pm and t/lib/complex.t. + $eps = 1e-10; +} + +sub near { + my $e = defined $_[2] ? $_[2] : $eps; + my $d = $_[1] ? abs($_[0]/$_[1] - 1) : abs($_[0]); + print "# near? $_[0] $_[1] : $d : $e\n"; + $_[1] ? ($d < $e) : abs($_[0]) < $e; +} + +print "# Sanity checks\n"; + +ok(near(sin(1), 0.841470984807897)); +ok(near(cos(1), 0.54030230586814)); +ok(near(tan(1), 1.5574077246549)); + +ok(near(sec(1), 1.85081571768093)); +ok(near(csc(1), 1.18839510577812)); +ok(near(cot(1), 0.642092615934331)); + +ok(near(asin(1), 1.5707963267949)); +ok(near(acos(1), 0)); +ok(near(atan(1), 0.785398163397448)); + +ok(near(asec(1), 0)); +ok(near(acsc(1), 1.5707963267949)); +ok(near(acot(1), 0.785398163397448)); + +ok(near(sinh(1), 1.1752011936438)); +ok(near(cosh(1), 1.54308063481524)); +ok(near(tanh(1), 0.761594155955765)); + +ok(near(sech(1), 0.648054273663885)); +ok(near(csch(1), 0.850918128239322)); +ok(near(coth(1), 1.31303528549933)); + +ok(near(asinh(1), 0.881373587019543)); +ok(near(acosh(1), 0)); +ok(near(atanh(0.9), 1.47221948958322)); # atanh(1.0) would be an error. + +ok(near(asech(0.9), 0.467145308103262)); +ok(near(acsch(2), 0.481211825059603)); +ok(near(acoth(2), 0.549306144334055)); + +print "# Basics\n"; + +$x = 0.9; +ok(near(tan($x), sin($x) / cos($x))); + +ok(near(sinh(2), 3.62686040784702)); + +ok(near(acsch(0.1), 2.99822295029797)); + +$x = asin(2); +is(ref $x, 'Math::Complex'); + +# avoid using Math::Complex here +$x =~ /^([^-]+)(-[^i]+)i$/; +($y, $z) = ($1, $2); +ok(near($y, 1.5707963267949)); +ok(near($z, -1.31695789692482)); + +ok(near(deg2rad(90), pi/2)); + +ok(near(rad2deg(pi), 180)); + +use Math::Trig ':radial'; + +{ + my ($r,$t,$z) = cartesian_to_cylindrical(1,1,1); + + ok(near($r, sqrt(2))); + ok(near($t, deg2rad(45))); + ok(near($z, 1)); + + ($x,$y,$z) = cylindrical_to_cartesian($r, $t, $z); + + ok(near($x, 1)); + ok(near($y, 1)); + ok(near($z, 1)); + + ($r,$t,$z) = cartesian_to_cylindrical(1,1,0); + + ok(near($r, sqrt(2))); + ok(near($t, deg2rad(45))); + ok(near($z, 0)); + + ($x,$y,$z) = cylindrical_to_cartesian($r, $t, $z); + + ok(near($x, 1)); + ok(near($y, 1)); + ok(near($z, 0)); +} + +{ + my ($r,$t,$f) = cartesian_to_spherical(1,1,1); + + ok(near($r, sqrt(3))); + ok(near($t, deg2rad(45))); + ok(near($f, atan2(sqrt(2), 1))); + + ($x,$y,$z) = spherical_to_cartesian($r, $t, $f); + + ok(near($x, 1)); + ok(near($y, 1)); + ok(near($z, 1)); + + ($r,$t,$f) = cartesian_to_spherical(1,1,0); + + ok(near($r, sqrt(2))); + ok(near($t, deg2rad(45))); + ok(near($f, deg2rad(90))); + + ($x,$y,$z) = spherical_to_cartesian($r, $t, $f); + + ok(near($x, 1)); + ok(near($y, 1)); + ok(near($z, 0)); +} + +{ + my ($r,$t,$z) = cylindrical_to_spherical(spherical_to_cylindrical(1,1,1)); + + ok(near($r, 1)); + ok(near($t, 1)); + ok(near($z, 1)); + + ($r,$t,$z) = spherical_to_cylindrical(cylindrical_to_spherical(1,1,1)); + + ok(near($r, 1)); + ok(near($t, 1)); + ok(near($z, 1)); +} + +{ + use Math::Trig 'great_circle_distance'; + + ok(near(great_circle_distance(0, 0, 0, pi/2), pi/2)); + + ok(near(great_circle_distance(0, 0, pi, pi), pi)); + + # London to Tokyo. + my @L = (deg2rad(-0.5), deg2rad(90 - 51.3)); + my @T = (deg2rad(139.8), deg2rad(90 - 35.7)); + + my $km = great_circle_distance(@L, @T, 6378); + + ok(near($km, 9605.26637021388)); +} + +{ + my $R2D = 57.295779513082320876798154814169; + + sub frac { $_[0] - int($_[0]) } + + my $lotta_radians = deg2rad(1E+20, 1); + ok(near($lotta_radians, 1E+20/$R2D)); + + my $negat_degrees = rad2deg(-1E20, 1); + ok(near($negat_degrees, -1E+20*$R2D)); + + my $posit_degrees = rad2deg(-10000, 1); + ok(near($posit_degrees, -10000*$R2D)); +} + +{ + use Math::Trig 'great_circle_direction'; + + ok(near(great_circle_direction(0, 0, 0, pi/2), pi)); + +# Retired test: Relies on atan2(0, 0), which is not portable. +# ok(near(great_circle_direction(0, 0, pi, pi), -pi()/2)); + + my @London = (deg2rad( -0.167), deg2rad(90 - 51.3)); + my @Tokyo = (deg2rad( 139.5), deg2rad(90 - 35.7)); + my @Berlin = (deg2rad ( 13.417), deg2rad(90 - 52.533)); + my @Paris = (deg2rad ( 2.333), deg2rad(90 - 48.867)); + + ok(near(rad2deg(great_circle_direction(@London, @Tokyo)), + 31.791945393073)); + + ok(near(rad2deg(great_circle_direction(@Tokyo, @London)), + 336.069766430326)); + + ok(near(rad2deg(great_circle_direction(@Berlin, @Paris)), + 246.800348034667)); + + ok(near(rad2deg(great_circle_direction(@Paris, @Berlin)), + 58.2079877553156)); + + use Math::Trig 'great_circle_bearing'; + + ok(near(rad2deg(great_circle_bearing(@Paris, @Berlin)), + 58.2079877553156)); + + use Math::Trig 'great_circle_waypoint'; + use Math::Trig 'great_circle_midpoint'; + + my ($lon, $lat); + + ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.0); + + ok(near($lon, $London[0])); + + ok(near($lat, $London[1])); + + ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 1.0); + + ok(near($lon, $Tokyo[0])); + + ok(near($lat, $Tokyo[1])); + + ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.5); + + ok(near($lon, 1.55609593577679)); # 89.16 E + + ok(near($lat, 0.36783532946162)); # 68.93 N + + ($lon, $lat) = great_circle_midpoint(@London, @Tokyo); + + ok(near($lon, 1.55609593577679)); # 89.16 E + + ok(near($lat, 0.367835329461615)); # 68.93 N + + ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.25); + + ok(near($lon, 0.516073562850837)); # 29.57 E + + ok(near($lat, 0.400231313403387)); # 67.07 N + + ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.75); + + ok(near($lon, 2.17494903805952)); # 124.62 E + + ok(near($lat, 0.617809294053591)); # 54.60 N + + use Math::Trig 'great_circle_destination'; + + my $dir1 = great_circle_direction(@London, @Tokyo); + my $dst1 = great_circle_distance(@London, @Tokyo); + + ($lon, $lat) = great_circle_destination(@London, $dir1, $dst1); + + ok(near($lon, $Tokyo[0])); + + ok(near($lat, $pip2 - $Tokyo[1])); + + my $dir2 = great_circle_direction(@Tokyo, @London); + my $dst2 = great_circle_distance(@Tokyo, @London); + + ($lon, $lat) = great_circle_destination(@Tokyo, $dir2, $dst2); + + ok(near($lon, $London[0])); + + ok(near($lat, $pip2 - $London[1])); + + my $dir3 = (great_circle_destination(@London, $dir1, $dst1))[2]; + + ok(near($dir3, 2.69379263839118)); # about 154.343 deg + + my $dir4 = (great_circle_destination(@Tokyo, $dir2, $dst2))[2]; + + ok(near($dir4, 3.6993902625701)); # about 211.959 deg + + ok(near($dst1, $dst2)); +} + +SKIP: { +# With netbsd-vax (or any vax) there is neither Inf, nor 1e40. +skip("different float range", 42) if $vax_float; +skip("no inf", 42) unless $has_inf; + +print "# Infinity\n"; + +my $BigDouble = eval '1e40'; + +# E.g. netbsd-alpha core dumps on Inf arith without this. +local $SIG{FPE} = sub { }; + +ok(Inf() > $BigDouble); # This passes in netbsd-alpha. +ok(Inf() + $BigDouble > $BigDouble); # This coredumps in netbsd-alpha. +ok(Inf() + $BigDouble == Inf()); +ok(Inf() - $BigDouble > $BigDouble); +ok(Inf() - $BigDouble == Inf()); +ok(Inf() * $BigDouble > $BigDouble); +ok(Inf() * $BigDouble == Inf()); +ok(Inf() / $BigDouble > $BigDouble); +ok(Inf() / $BigDouble == Inf()); + +ok(-Inf() < -$BigDouble); +ok(-Inf() + $BigDouble < $BigDouble); +ok(-Inf() + $BigDouble == -Inf()); +ok(-Inf() - $BigDouble < -$BigDouble); +ok(-Inf() - $BigDouble == -Inf()); +ok(-Inf() * $BigDouble < -$BigDouble); +ok(-Inf() * $BigDouble == -Inf()); +ok(-Inf() / $BigDouble < -$BigDouble); +ok(-Inf() / $BigDouble == -Inf()); + +print "# sinh/sech/cosh/csch/tanh/coth unto infinity\n"; + +ok(near(sinh(100), eval '1.3441e+43', 1e-3)); +ok(near(sech(100), eval '7.4402e-44', 1e-3)); +ok(near(cosh(100), eval '1.3441e+43', 1e-3)); +ok(near(csch(100), eval '7.4402e-44', 1e-3)); +ok(near(tanh(100), 1)); +ok(near(coth(100), 1)); + +ok(near(sinh(-100), eval '-1.3441e+43', 1e-3)); +ok(near(sech(-100), eval ' 7.4402e-44', 1e-3)); +ok(near(cosh(-100), eval ' 1.3441e+43', 1e-3)); +ok(near(csch(-100), eval '-7.4402e-44', 1e-3)); +ok(near(tanh(-100), -1)); +ok(near(coth(-100), -1)); + +cmp_ok(sinh(1e5), '==', Inf()); +cmp_ok(sech(1e5), '==', 0); +cmp_ok(cosh(1e5), '==', Inf()); +cmp_ok(csch(1e5), '==', 0); +cmp_ok(tanh(1e5), '==', 1); +cmp_ok(coth(1e5), '==', 1); + +cmp_ok(sinh(-1e5), '==', -Inf()); +cmp_ok(sech(-1e5), '==', 0); +cmp_ok(cosh(-1e5), '==', Inf()); +cmp_ok(csch(-1e5), '==', 0); +cmp_ok(tanh(-1e5), '==', -1); +cmp_ok(coth(-1e5), '==', -1); + +} + +print "# great_circle_distance with small angles\n"; + +for my $e (qw(1e-2 1e-3 1e-4 1e-5)) { + # Can't assume == 0 because of floating point fuzz, + # but let's hope for at least < $e. + cmp_ok(great_circle_distance(0, $e, 0, $e), '<', $e); +} + +print "# asin_real, acos_real\n"; + +is(acos_real(-2.0), pi); +is(acos_real(-1.0), pi); +is(acos_real(-0.5), acos(-0.5)); +is(acos_real( 0.0), acos( 0.0)); +is(acos_real( 0.5), acos( 0.5)); +is(acos_real( 1.0), 0); +is(acos_real( 2.0), 0); + +is(asin_real(-2.0), -&pip2); +is(asin_real(-1.0), -&pip2); +is(asin_real(-0.5), asin(-0.5)); +is(asin_real( 0.0), asin( 0.0)); +is(asin_real( 0.5), asin( 0.5)); +is(asin_real( 1.0), pip2); +is(asin_real( 2.0), pip2); + +# eof diff --git a/dist/Math-Complex/t/underbar.t b/dist/Math-Complex/t/underbar.t new file mode 100644 index 0000000000..809e8805a0 --- /dev/null +++ b/dist/Math-Complex/t/underbar.t @@ -0,0 +1,28 @@ +# +# Tests that the standard Perl 5 functions that we override +# that operate on the $_ will work correctly [perl #62412] +# + +use Test::More; + +use strict; +use warnings; + +my @f = qw(abs cos exp log sin sqrt); + +plan tests => scalar @f; + +use Math::Complex; + +my %CORE; + +for my $f (@f) { + local $_ = 0.5; + $CORE{$f} = eval "CORE::$f"; +} + +for my $f (@f) { + local $_ = 0.5; + is(eval "Math::Complex::$f", $CORE{$f}, $f); +} + |