package Math::Complex;
-our($VERSION, @ISA, @EXPORT, %EXPORT_TAGS, $Inf);
+use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $Inf);
-$VERSION = 1.31;
+$VERSION = 1.36;
BEGIN {
unless ($^O eq 'unicosmk') {
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);
sqrt log ln
log10 logn cbrt root
cplx cplxe
+ atan2
),
@trig);
+my @pi = qw(pi pi2 pi4 pip2 pip4);
+
+@EXPORT_OK = @pi;
+
%EXPORT_TAGS = (
'trig' => [@trig],
+ 'pi' => [@pi],
);
use overload
- '+' => \&plus,
- '-' => \&minus,
- '*' => \&multiply,
- '/' => \÷,
- '**' => \&power,
- '==' => \&numeq,
- '<=>' => \&spaceship,
- 'neg' => \&negate,
- '~' => \&conjugate,
+ '+' => \&_plus,
+ '-' => \&_minus,
+ '*' => \&_multiply,
+ '/' => \&_divide,
+ '**' => \&_power,
+ '==' => \&_numeq,
+ '<=>' => \&_spaceship,
+ 'neg' => \&_negate,
+ '~' => \&_conjugate,
'abs' => \&abs,
'sqrt' => \&sqrt,
'exp' => \&exp,
'cos' => \&cos,
'tan' => \&tan,
'atan2' => \&atan2,
- qw("" stringify);
+ '""' => \&_stringify;
#
# Package "privates"
# 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)
+# bn_cartesian
+# bnc_dirty
#
# Die on bad *make() arguments.
sub _cannot_make {
- die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n";
+ 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);
}
#
# Create a new complex number (cartesian form)
#
sub make {
- my $self = bless {}, shift;
- my ($re, $im) = @_;
- my $rre = ref $re;
- if ( $rre ) {
- if ( $rre eq ref $self ) {
- $re = Re($re);
- } else {
- _cannot_make("real part", $rre);
- }
- }
- my $rim = ref $im;
- if ( $rim ) {
- if ( $rim eq ref $self ) {
- $im = Im($im);
- } else {
- _cannot_make("imaginary part", $rim);
- }
- }
- $self->{'cartesian'} = [ $re, $im ];
- $self->{c_dirty} = 0;
- $self->{p_dirty} = 1;
- $self->display_format('cartesian');
- return $self;
+ 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;
}
#
# Create a new complex number (exponential form)
#
sub emake {
- my $self = bless {}, shift;
- my ($rho, $theta) = @_;
- my $rrh = ref $rho;
- if ( $rrh ) {
- if ( $rrh eq ref $self ) {
- $rho = rho($rho);
- } else {
- _cannot_make("rho", $rrh);
- }
- }
- my $rth = ref $theta;
- if ( $rth ) {
- if ( $rth eq ref $self ) {
- $theta = theta($theta);
- } else {
- _cannot_make("theta", $rth);
- }
- }
+ 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();
}
- $self->{'polar'} = [$rho, $theta];
- $self->{p_dirty} = 0;
- $self->{c_dirty} = 1;
- $self->display_format('polar');
- return $self;
+ }
+ 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.
# This avoids the burden of writing Math::Complex->make(re, im).
#
sub cplx {
- my ($re, $im) = @_;
- return __PACKAGE__->make($re, defined $im ? $im : 0);
+ return __PACKAGE__->make(@_);
}
#
# This avoids the burden of writing Math::Complex->emake(rho, theta).
#
sub cplxe {
- my ($rho, $theta) = @_;
- return __PACKAGE__->emake($rho, defined $theta ? $theta : 0);
+ return __PACKAGE__->emake(@_);
}
#
sub pi () { 4 * CORE::atan2(1, 1) }
#
-# pit2
+# pi2
#
# The full circle
#
-sub pit2 () { 2 * pi }
+sub pi2 () { 2 * pi }
+
+#
+# pi4
+#
+# The full circle twice.
+#
+sub pi4 () { 4 * pi }
#
# pip2
sub pip2 () { pi / 2 }
#
-# deg1
+# pip4
#
-# One degree in radians, used in stringify_polar.
+# The eighth circle.
#
-
-sub deg1 () { pi / 180 }
+sub pip4 () { pi / 4 }
#
-# uplog10
+# _uplog10
#
# Used in log10().
#
-sub uplog10 () { 1 / CORE::log(10) }
+sub _uplog10 () { 1 / CORE::log(10) }
#
# i
}
#
-# ip2
+# _ip2
#
# Half of i.
#
-sub ip2 () { i / 2 }
+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 _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]->{'cartesian'} = $_[1] }
-sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
+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
+# ->_update_cartesian
#
# Recompute and return the cartesian form, given accurate polar form.
#
-sub update_cartesian {
+sub _update_cartesian {
my $self = shift;
my ($r, $t) = @{$self->{'polar'}};
$self->{c_dirty} = 0;
#
#
-# ->update_polar
+# ->_update_polar
#
# Recompute and return the polar form, given accurate cartesian form.
#
-sub update_polar {
+sub _update_polar {
my $self = shift;
my ($x, $y) = @{$self->{'cartesian'}};
$self->{p_dirty} = 0;
}
#
-# (plus)
+# (_plus)
#
# Computes z1+z2.
#
-sub plus {
+sub _plus {
my ($z1, $z2, $regular) = @_;
- my ($re1, $im1) = @{$z1->cartesian};
+ my ($re1, $im1) = @{$z1->_cartesian};
$z2 = cplx($z2) unless ref $z2;
- my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
+ my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
unless (defined $regular) {
- $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
+ $z1->_set_cartesian([$re1 + $re2, $im1 + $im2]);
return $z1;
}
return (ref $z1)->make($re1 + $re2, $im1 + $im2);
}
#
-# (minus)
+# (_minus)
#
# Computes z1-z2.
#
-sub minus {
+sub _minus {
my ($z1, $z2, $inverted) = @_;
- my ($re1, $im1) = @{$z1->cartesian};
+ my ($re1, $im1) = @{$z1->_cartesian};
$z2 = cplx($z2) unless ref $z2;
- my ($re2, $im2) = @{$z2->cartesian};
+ my ($re2, $im2) = @{$z2->_cartesian};
unless (defined $inverted) {
- $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
+ $z1->_set_cartesian([$re1 - $re2, $im1 - $im2]);
return $z1;
}
return $inverted ?
}
#
-# (multiply)
+# (_multiply)
#
# Computes z1*z2.
#
-sub multiply {
+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 ($r1, $t1) = @{$z1->_polar};
+ my ($r2, $t2) = @{$z2->_polar};
my $t = $t1 + $t2;
- if ($t > pi()) { $t -= pit2 }
- elsif ($t <= -pi()) { $t += pit2 }
+ if ($t > pi()) { $t -= pi2 }
+ elsif ($t <= -pi()) { $t += pi2 }
unless (defined $regular) {
- $z1->set_polar([$r1 * $r2, $t]);
+ $z1->_set_polar([$r1 * $r2, $t]);
return $z1;
}
return (ref $z1)->emake($r1 * $r2, $t);
} else {
- my ($x1, $y1) = @{$z1->cartesian};
+ my ($x1, $y1) = @{$z1->_cartesian};
if (ref $z2) {
- my ($x2, $y2) = @{$z2->cartesian};
+ 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);
}
#
-# (divide)
+# (_divide)
#
# Computes z1/z2.
#
-sub divide {
+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 ($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 -= pit2 }
- elsif ($t <= -pi()) { $t += pit2 }
+ 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 -= pit2 }
- elsif ($t <= -pi()) { $t += pit2 }
+ 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};
+ ($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};
+ my ($x1, $y1) = @{$z1->_cartesian};
if (ref $z2) {
- ($x2, $y2) = @{$z2->cartesian};
+ ($x2, $y2) = @{$z2->_cartesian};
$d = $x2*$x2 + $y2*$y2;
_divbyzero "$z1/0" if $d == 0;
my $u = ($x1*$x2 + $y1*$y2)/$d;
}
#
-# (power)
+# (_power)
#
# Computes z1**z2 = exp(z2 * log z1)).
#
-sub power {
+sub _power {
my ($z1, $z2, $inverted) = @_;
if ($inverted) {
return 1 if $z1 == 0 || $z2 == 1;
# 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;
+ cplx(@{$w->_cartesian}) : $w;
}
#
-# (spaceship)
+# (_spaceship)
#
# Computes z1 <=> z2.
# Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i.
#
-sub spaceship {
+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 ($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)
+# (_numeq)
#
# Computes z1 == z2.
#
-# (Required in addition to spaceship() because of NaNs.)
-sub numeq {
+# (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);
+ 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)
+# (_negate)
#
# Computes -z.
#
-sub negate {
+sub _negate {
my ($z) = @_;
if ($z->{c_dirty}) {
- my ($r, $t) = @{$z->polar};
+ my ($r, $t) = @{$z->_polar};
$t = ($t <= 0) ? $t + pi : $t - pi;
return (ref $z)->emake($r, $t);
}
- my ($re, $im) = @{$z->cartesian};
+ my ($re, $im) = @{$z->_cartesian};
return (ref $z)->make(-$re, -$im);
}
#
-# (conjugate)
+# (_conjugate)
#
-# Compute complex's conjugate.
+# Compute complex's _conjugate.
#
-sub conjugate {
+sub _conjugate {
my ($z) = @_;
if ($z->{c_dirty}) {
- my ($r, $t) = @{$z->polar};
+ my ($r, $t) = @{$z->_polar};
return (ref $z)->emake($r, -$t);
}
- my ($re, $im) = @{$z->cartesian};
+ my ($re, $im) = @{$z->_cartesian};
return (ref $z)->make($re, -$im);
}
}
}
if (defined $rho) {
- $z->{'polar'} = [ $rho, ${$z->polar}[1] ];
+ $z->{'polar'} = [ $rho, ${$z->_polar}[1] ];
$z->{p_dirty} = 0;
$z->{c_dirty} = 1;
return $rho;
} else {
- return ${$z->polar}[0];
+ return ${$z->_polar}[0];
}
}
sub _theta {
my $theta = $_[0];
- if ($$theta > pi()) { $$theta -= pit2 }
- elsif ($$theta <= -pi()) { $$theta += pit2 }
+ if ($$theta > pi()) { $$theta -= pi2 }
+ elsif ($$theta <= -pi()) { $$theta += pi2 }
}
#
return $z unless ref $z;
if (defined $theta) {
_theta(\$theta);
- $z->{'polar'} = [ ${$z->polar}[0], $theta ];
+ $z->{'polar'} = [ ${$z->_polar}[0], $theta ];
$z->{p_dirty} = 0;
$z->{c_dirty} = 1;
} else {
- $theta = ${$z->polar}[1];
+ $theta = ${$z->_polar}[1];
_theta(\$theta);
}
return $theta;
#
sub sqrt {
my ($z) = @_;
- my ($re, $im) = ref $z ? @{$z->cartesian} : ($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};
+ my ($r, $t) = @{$z->_polar};
return (ref $z)->emake(CORE::sqrt($r), $t/2);
}
-CORE::exp(CORE::log(-$z)/3) :
($z > 0 ? CORE::exp(CORE::log($z)/3): 0)
unless ref $z;
- my ($r, $t) = @{$z->polar};
+ my ($r, $t) = @{$z->_polar};
return 0 if $r == 0;
return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3);
}
# Die on bad root.
#
sub _rootbad {
- my $mess = "Root $_[0] illegal, root rank must be positive integer.\n";
+ my $mess = "Root '$_[0]' illegal, root rank must be positive integer.\n";
my @up = caller(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) = @_;
+ 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 @root;
- my $k;
- my $theta_inc = pit2 / $n;
+ @{$z->_polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
+ my $theta_inc = pi2 / $n;
my $rho = $r ** (1/$n);
- my $theta;
my $cartesian = ref $z && $z->{c_dirty} == 0;
- for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
- my $w = cplxe($rho, $theta);
- # Yes, $cartesian is loop invariant.
- push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
+ 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;
}
- return @root;
}
#
my ($z, $Re) = @_;
return $z unless ref $z;
if (defined $Re) {
- $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ];
+ $z->{'cartesian'} = [ $Re, ${$z->_cartesian}[1] ];
$z->{c_dirty} = 0;
$z->{p_dirty} = 1;
} else {
- return ${$z->cartesian}[0];
+ return ${$z->_cartesian}[0];
}
}
my ($z, $Im) = @_;
return 0 unless ref $z;
if (defined $Im) {
- $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ];
+ $z->{'cartesian'} = [ ${$z->_cartesian}[0], $Im ];
$z->{c_dirty} = 0;
$z->{p_dirty} = 1;
} else {
- return ${$z->cartesian}[1];
+ return ${$z->_cartesian}[1];
}
}
#
sub exp {
my ($z) = @_;
- my ($x, $y) = @{$z->cartesian};
+ my ($x, $y) = @{$z->_cartesian};
return (ref $z)->emake(CORE::exp($x), $y);
}
_logofzero("log") if $z == 0;
return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
}
- my ($r, $t) = @{$z->polar};
+ my ($r, $t) = @{$z->_polar};
_logofzero("log") if $r == 0;
- if ($t > pi()) { $t -= pit2 }
- elsif ($t <= -pi()) { $t += pit2 }
+ if ($t > pi()) { $t -= pi2 }
+ elsif ($t <= -pi()) { $t += pi2 }
return (ref $z)->make(CORE::log($r), $t);
}
#
sub log10 {
- return Math::Complex::log($_[0]) * uplog10;
+ return Math::Complex::log($_[0]) * _uplog10;
}
#
sub cos {
my ($z) = @_;
return CORE::cos($z) unless ref $z;
- my ($x, $y) = @{$z->cartesian};
+ my ($x, $y) = @{$z->_cartesian};
my $ey = CORE::exp($y);
my $sx = CORE::sin($x);
my $cx = CORE::cos($x);
sub sin {
my ($z) = @_;
return CORE::sin($z) unless ref $z;
- my ($x, $y) = @{$z->cartesian};
+ my ($x, $y) = @{$z->_cartesian};
my $ey = CORE::exp($y);
my $sx = CORE::sin($x);
my $cx = CORE::cos($x);
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};
+ 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);
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};
+ 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);
sub atan {
my ($z) = @_;
return CORE::atan2($z, 1) unless ref $z;
- my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
+ 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;
+ return _ip2 * $log;
}
#
$ex = CORE::exp($z);
return $ex ? ($ex + 1/$ex)/2 : $Inf;
}
- my ($x, $y) = @{$z->cartesian};
+ 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,
$ex = CORE::exp($z);
return $ex ? ($ex - 1/$ex)/2 : "-$Inf";
}
- my ($x, $y) = @{$z->cartesian};
+ my ($x, $y) = @{$z->_cartesian};
my $cy = CORE::cos($y);
my $sy = CORE::sin($y);
$ex = CORE::exp($x);
unless (ref $z) {
$z = cplx($z, 0);
}
- my ($re, $im) = @{$z->cartesian};
+ my ($re, $im) = @{$z->_cartesian};
if ($im == 0) {
return CORE::log($re + CORE::sqrt($re*$re - 1))
if $re >= 1;
#
# (atan2)
#
-# Compute atan(z1/z2).
+# 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) = @{$z1->cartesian};
+ ($re1, $im1) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
+ ($re2, $im2) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
} else {
- ($re1, $im1) = @{$z1->cartesian};
- ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
+ ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0);
+ ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0);
}
- if ($im2 == 0) {
- return CORE::atan2($re1, $re2) if $im1 == 0;
- return ($im1<=>0) * pip2 if $re2 == 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 ));
}
- my $w = atan($z1/$z2);
- my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
- $u += pi if $re2 < 0;
- $u -= pit2 if $u > pi;
- return cplx($u, $v);
+ return CORE::atan2($re1, $re2);
}
#
}
#
-# (stringify)
+# (_stringify)
#
# Show nicely formatted complex number under its cartesian or polar form,
# depending on the current display format:
# . Otherwise, use the generic current default for all complex numbers,
# which is a package global variable.
#
-sub stringify {
+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;
+ return $z->_stringify_polar if $style =~ /^p/i;
+ return $z->_stringify_cartesian;
}
#
-# ->stringify_cartesian
+# ->_stringify_cartesian
#
# Stringify as a cartesian representation 'a+bi'.
#
-sub stringify_cartesian {
+sub _stringify_cartesian {
my $z = shift;
- my ($x, $y) = @{$z->cartesian};
+ my ($x, $y) = @{$z->_cartesian};
my ($re, $im);
my %format = $z->display_format;
#
-# ->stringify_polar
+# ->_stringify_polar
#
# Stringify as a polar representation '[r,t]'.
#
-sub stringify_polar {
+sub _stringify_polar {
my $z = shift;
- my ($r, $t) = @{$z->polar};
+ my ($r, $t) = @{$z->_polar};
my $theta;
my %format = $z->display_format;
# Try to identify pi/n and friends.
#
- $t -= int(CORE::abs($t) / pit2) * pit2;
+ $t -= int(CORE::abs($t) / pi2) * pi2;
if ($format{polar_pretty_print} && $t) {
my ($a, $b);
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> will be
-noted C<abs(z)>.
+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<theta> is zero or I<pi>.
+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
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 an horizontal bar above the number, or C<~z> here.
+with a horizontal bar above the number, or C<~z> here.
z = a + bi
~z = a - bi
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(z1, z2) = atan(z1/z2)
+ 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:
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 also mutators. The C<cbrt>
+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.
(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> 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
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
-either the C<make> or C<emake>: the appropriate component of
+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);
-=head1 STRINGIFICATION
+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.
+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
(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
$j->arg(2); # (the last two aka rho, theta)
# can be used also as mutators.
+=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;
+
=head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
The division (/) and the following functions
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 the
+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.
+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
=head1 BUGS
Saying C<use Math::Complex;> exports many mathematical routines in the
-caller environment and even overrides some (C<sqrt>, C<log>).
+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
=head1 AUTHORS
-Raphael Manfredi <F<Raphael_Manfredi@pobox.com>> and
-Jarkko Hietaniemi <F<jhi@iki.fi>>.
-
-Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.
+Daniel S. Lewart <F<lewart!at!uiuc.edu>>
+Jarkko Hietaniemi <F<jhi!at!iki.fi>>
+Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>
=cut
projects / p5sagit/p5-mst-13.2.git / blobdiff