diff options
Diffstat (limited to 'ext')
-rw-r--r-- | ext/Math-Complex/lib/Math/Complex.pm | 2093 | ||||
-rw-r--r-- | ext/Math-Complex/lib/Math/Trig.pm | 768 | ||||
-rw-r--r-- | ext/Math-Complex/t/Complex.t | 1136 | ||||
-rw-r--r-- | ext/Math-Complex/t/Trig.t | 393 | ||||
-rw-r--r-- | ext/Math-Complex/t/underbar.t | 27 |
5 files changed, 0 insertions, 4417 deletions
diff --git a/ext/Math-Complex/lib/Math/Complex.pm b/ext/Math-Complex/lib/Math/Complex.pm deleted file mode 100644 index 8475a2b5d0..0000000000 --- a/ext/Math-Complex/lib/Math/Complex.pm +++ /dev/null @@ -1,2093 +0,0 @@ -# -# 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 strict; - -use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $Inf $ExpInf); - -$VERSION = 1.56; - -use Config; - -BEGIN { - my %DBL_MAX = - ( - 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} = { }; - 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 = 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; - -@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 - ); - -@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); - -@EXPORT_OK = @pi; - -%EXPORT_TAGS = ( - 'trig' => [@trig], - 'pi' => [@pi], -); - -use overload - '+' => \&_plus, - '-' => \&_minus, - '*' => \&_multiply, - '/' => \&_divide, - '**' => \&_power, - '==' => \&_numeq, - '<=>' => \&_spaceship, - 'neg' => \&_negate, - '~' => \&_conjugate, - 'abs' => \&abs, - 'sqrt' => \&sqrt, - 'exp' => \&exp, - 'log' => \&log, - 'sin' => \&sin, - 'cos' => \&cos, - 'tan' => \&tan, - '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; - $q =~ s/^\+//; - $q =~ s/^(-?)inf$/"${1}9**9**9"/e; - } - - 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; - $q =~ s/^(-?)inf$/"${1}9**9**9"/e; - } - - return ($p, $q); -} - -# -# ->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, $theta) 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>> - -=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/ext/Math-Complex/lib/Math/Trig.pm b/ext/Math-Complex/lib/Math/Trig.pm deleted file mode 100644 index b7767bebcc..0000000000 --- a/ext/Math-Complex/lib/Math/Trig.pm +++ /dev/null @@ -1,768 +0,0 @@ -# -# Trigonometric functions, mostly inherited from Math::Complex. -# -- Jarkko Hietaniemi, since April 1997 -# -- Raphael Manfredi, September 1996 (indirectly: because of Math::Complex) -# - -require Exporter; -package Math::Trig; - -use 5.005; -use strict; - -use Math::Complex 1.56; -use Math::Complex qw(:trig :pi); - -use vars qw($VERSION $PACKAGE @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS); - -@ISA = qw(Exporter); - -$VERSION = 1.20; - -my @angcnv = qw(rad2deg rad2grad - deg2rad deg2grad - grad2rad grad2deg); - -my @areal = qw(asin_real acos_real); - -@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); - -@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 - -%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 $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1); - - my $lat0 = pip2 - $phi0; - my $lat1 = pip2 - $phi1; - - my $direction = - acos_real((sin($lat1) - sin($lat0) * cos($distance)) / - (cos($lat0) * sin($distance))); - - $direction = pi2 - $direction - if sin($theta1 - $theta0) < 0; - - return rad2rad($direction); -} - -*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>> and -Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>. - -=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/ext/Math-Complex/t/Complex.t b/ext/Math-Complex/t/Complex.t deleted file mode 100644 index 687d6220df..0000000000 --- a/ext/Math-Complex/t/Complex.t +++ /dev/null @@ -1,1136 +0,0 @@ -#!./perl - -# -# Regression tests for the Math::Complex pacakge -# -- Raphael Manfredi since Sep 1996 -# -- Jarkko Hietaniemi since Mar 1997 -# -- Daniel S. Lewart since Sep 1997 - -BEGIN { - if ($ENV{PERL_CORE}) { - chdir 't' if -d 't'; - #@INC = '../lib'; - } -} - -use Math::Complex 1.54; - -use vars qw($VERSION); - -$VERSION = 1.92; - -my ($args, $op, $target, $test, $test_set, $try, $val, $zvalue, @set, @val); - -$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. - -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); -my $inf = 9**9**9; -'; - -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' => '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]" - -{ (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/ext/Math-Complex/t/Trig.t b/ext/Math-Complex/t/Trig.t deleted file mode 100644 index ae0b0fd1d6..0000000000 --- a/ext/Math-Complex/t/Trig.t +++ /dev/null @@ -1,393 +0,0 @@ -#!./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 - -BEGIN { - if ($ENV{PERL_CORE}) { - chdir 't' if -d 't'; - #@INC = '../lib'; - } -} - -BEGIN { - eval { require Test::More }; - if ($@) { - # We are willing to lose testing in e.g. 5.00504. - print "1..0 # No Test::More, skipping\n"; - exit(0); - } else { - import Test::More; - } -} - -plan(tests => 153); - -use Math::Trig 1.18; -use Math::Trig 1.18 qw(:pi Inf); - -my $pip2 = pi / 2; - -use strict; - -use vars qw($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)); -} - -print "# Infinity\n"; - -my $BigDouble = 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), 1.3441e+43, 1e-3)); -ok(near(sech(100), 7.4402e-44, 1e-3)); -ok(near(cosh(100), 1.3441e+43, 1e-3)); -ok(near(csch(100), 7.4402e-44, 1e-3)); -ok(near(tanh(100), 1)); -ok(near(coth(100), 1)); - -ok(near(sinh(-100), -1.3441e+43, 1e-3)); -ok(near(sech(-100), 7.4402e-44, 1e-3)); -ok(near(cosh(-100), 1.3441e+43, 1e-3)); -ok(near(csch(-100), -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/ext/Math-Complex/t/underbar.t b/ext/Math-Complex/t/underbar.t deleted file mode 100644 index 643e86655c..0000000000 --- a/ext/Math-Complex/t/underbar.t +++ /dev/null @@ -1,27 +0,0 @@ -# -# Tests that the standard Perl 5 functions that we override -# that operate on the $_ will work correctly [perl #62412] -# - -use Test::More; - -my @f = qw(abs cos exp log sin sqrt); - -plan tests => scalar @f; - -use strict; - -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); -} - |