# -- Daniel S. Lewart Since Sep 1997
#
-require Exporter;
package Math::Complex;
-use 5.005_64;
+our($VERSION, @ISA, @EXPORT, %EXPORT_TAGS, $Inf);
+
+$VERSION = 1.34;
+
+BEGIN {
+ unless ($^O eq 'unicosmk') {
+ my $e = $!;
+ # We do want an arithmetic overflow, Inf INF inf Infinity:.
+ undef $Inf unless eval <<'EOE' and $Inf =~ /^inf(?:inity)?$/i;
+ local $SIG{FPE} = sub {die};
+ my $t = CORE::exp 30;
+ $Inf = CORE::exp $t;
+EOE
+ if (!defined $Inf) { # Try a different method
+ undef $Inf unless eval <<'EOE' and $Inf =~ /^inf(?:inity)?$/i;
+ local $SIG{FPE} = sub {die};
+ my $t = 1;
+ $Inf = $t + "1e99999999999999999999999999999999";
+EOE
+ }
+ $! = $e; # Clear ERANGE.
+ }
+ $Inf = "Inf" if !defined $Inf || !($Inf > 0); # Desperation.
+}
+
use strict;
-our($VERSION, @ISA, @EXPORT, %EXPORT_TAGS);
+my $i;
+my %LOGN;
-my ( $i, %logn );
+# Regular expression for floating point numbers.
+my $gre = qr'\s*([\+\-]?(?:(?:(?:\d+(?:_\d+)*(?:\.\d*(?:_\d+)*)?|\.\d+(?:_\d+)*)(?:[eE][\+\-]?\d+(?:_\d+)*)?)))';
-$VERSION = sprintf("%s", q$Id: Complex.pm,v 1.26 1998/11/01 00:00:00 dsl Exp $ =~ /(\d+\.\d+)/);
+require Exporter;
@ISA = qw(Exporter);
# Package "privates"
#
-my $package = 'Math::Complex'; # Package name
my %DISPLAY_FORMAT = ('style' => 'cartesian',
'polar_pretty_print' => 1);
my $eps = 1e-14; # Epsilon
-my $Inf;
-unless ($^O eq 'unicos') { # Unicos gets a fatal runtime error
- $Inf = CORE::exp(CORE::exp(30));
-}
-$Inf = "Inf" if !defined $Inf || !$Inf > 0;
-
#
# Object attributes (internal):
# cartesian [real, imaginary] -- cartesian form
die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n";
}
+sub _remake {
+ my $arg = shift;
+ my ($made, $p, $q);
+
+ if ($arg =~ /^(?:$gre)?$gre\s*i\s*$/) {
+ ($p, $q) = ($1 || 0, $2);
+ $made = 'cart';
+ } elsif ($arg =~ /^\s*\[\s*$gre\s*(?:,\s*$gre\s*)?\]\s*$/) {
+ ($p, $q) = ($1, $2 || 0);
+ $made = 'exp';
+ }
+
+ if ($made) {
+ $p =~ s/^\+//;
+ $q =~ s/^\+//;
+ }
+
+ return ($made, $p, $q);
+}
+
#
# ->make
#
sub make {
my $self = bless {}, shift;
my ($re, $im) = @_;
+ if (@_ == 1) {
+ my ($remade, $p, $q) = _remake($re);
+ if ($remade) {
+ if ($remade eq 'cart') {
+ ($re, $im) = ($p, $q);
+ } else {
+ return (ref $self)->emake($p, $q);
+ }
+ }
+ }
my $rre = ref $re;
if ( $rre ) {
if ( $rre eq ref $self ) {
_cannot_make("imaginary part", $rim);
}
}
+ _cannot_make("real part", $re) unless $re =~ /^$gre$/;
+ $im ||= 0;
+ _cannot_make("imaginary part", $im) unless $im =~ /^$gre$/;
$self->{'cartesian'} = [ $re, $im ];
$self->{c_dirty} = 0;
$self->{p_dirty} = 1;
sub emake {
my $self = bless {}, shift;
my ($rho, $theta) = @_;
+ if (@_ == 1) {
+ my ($remade, $p, $q) = _remake($rho);
+ if ($remade) {
+ if ($remade eq 'exp') {
+ ($rho, $theta) = ($p, $q);
+ } else {
+ return (ref $self)->make($p, $q);
+ }
+ }
+ }
my $rrh = ref $rho;
if ( $rrh ) {
if ( $rrh eq ref $self ) {
$rho = -$rho;
$theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
}
+ _cannot_make("rho", $rho) unless $rho =~ /^$gre$/;
+ $theta ||= 0;
+ _cannot_make("theta", $theta) unless $theta =~ /^$gre$/;
$self->{'polar'} = [$rho, $theta];
$self->{p_dirty} = 0;
$self->{c_dirty} = 1;
# 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 Im {
my ($z, $Im) = @_;
- return $z unless ref $z;
+ return 0 unless ref $z;
if (defined $Im) {
$z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ];
$z->{c_dirty} = 0;
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)
+ my $logn = $LOGN{$n};
+ $logn = $LOGN{$n} = CORE::log($n) unless defined $logn; # Cache log(n)
return &log($z) / $logn;
}
my $z = $_[0];
return CORE::atan2(CORE::sqrt(1-$z*$z), $z)
if (! ref $z) && CORE::abs($z) <= 1;
- my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
+ $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 $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 __PACKAGE__->make($u, $v);
+ return (ref $z)->make($u, $v);
}
#
my $z = $_[0];
return CORE::atan2($z, CORE::sqrt(1-$z*$z))
if (! ref $z) && CORE::abs($z) <= 1;
- my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
+ $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 $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 __PACKAGE__->make($u, $v);
+ return (ref $z)->make($u, $v);
}
#
return $ex ? ($ex + 1/$ex)/2 : $Inf;
}
my ($x, $y) = @{$z->cartesian};
- my $cy = CORE::cos($y);
- my $sy = CORE::cos($y);
$ex = CORE::exp($x);
my $ex_1 = $ex ? 1 / $ex : $Inf;
return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
my $sy = CORE::sin($y);
$ex = CORE::exp($x);
my $ex_1 = $ex ? 1 / $ex : $Inf;
- return (ref $z)->make($cy * ($ex - $ex_1)/2,
- $sy * ($ex + $ex_1)/2);
+ return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
+ CORE::sin($y) * ($ex + $ex_1)/2);
}
#
return cplx(0, CORE::atan2(CORE::sqrt(1 - $re*$re), $re))
if CORE::abs($re) < 1;
}
- my $s = &sqrt($z*$z - 1);
- my $t = $z + $s;
- $t = 1/(2*$s) if $t == 0 || $t && &abs(cosh(&log($t)) - $z) > $eps;
- return &log($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;
+ my $u = &log($t);
+ $u->Im(-$u->Im) if $re < 0 && $im == 0;
+ return $re < 0 ? -$u : $u;
}
#
my $t = $z + CORE::sqrt($z*$z + 1);
return CORE::log($t) if $t;
}
- my $s = &sqrt($z*$z + 1);
- my $t = $z + $s;
- # Try Taylor series if looking bad.
- $t = 1/(2*$s) if $t == 0 || $t && &abs(sinh(&log($t)) - $z) > $eps;
+ 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);
}
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;
- }
- } else { # Called as a class method
- if (@_ = 1) {
- $display_format{style} = $self;
- } else {
- my %new = @_;
- @display_format{keys %new} = values %new;
- }
- undef $self;
+ }
+ if (@_ == 1) {
+ $display_format{style} = shift;
+ } else {
+ my %new = @_;
+ @display_format{keys %new} = values %new;
}
- if (defined $self) {
+ if (ref $self) { # Called as an object method
$self->{display_format} = { %display_format };
return
wantarray ?
$self->{display_format}->{style};
}
+ # Called as a class method
%DISPLAY_FORMAT = %display_format;
return
wantarray ?
}
if ($y) {
- if ($y == 1) { $im = "" }
- elsif ($y == -1) { $im = "-" }
- elsif ($y =~ /^(NaN[QS]?)$/i) {
+ if ($y =~ /^(NaN[QS]?)$/i) {
$im = $y;
} else {
if ($y =~ /^-?$Inf$/oi) {
$im = $y;
} else {
- $im = defined $format ? sprintf($format, $y) : $y;
+ $im =
+ defined $format ?
+ sprintf($format, $y) :
+ ($y == 1 ? "" : ($y == -1 ? "-" : $y));
}
}
$im .= "i";
$t -= int(CORE::abs($t) / pit2) * pit2;
- if ($format{polar_pretty_print}) {
+ if ($format{polar_pretty_print} && $t) {
my ($a, $b);
- for $a (2, 3, 4, 6, 8, 12, 16, 24, 30, 32, 36, 48, 60, 64, 72) {
+ for $a (2..9) {
$b = $t * $a / pi;
- if (int($b) == $b) {
+ if ($b =~ /^-?\d+$/) {
$b = $b < 0 ? "-" : "" if CORE::abs($b) == 1;
$theta = "${b}pi/$a";
last;
__END__
=pod
+
=head1 NAME
Math::Complex - complex numbers and associated mathematical functions
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
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.
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);
+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]
+
+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.
+
=head1 STRINGIFICATION
When printed, a complex number is usually shown under its 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]),
+(where I<n> is a positive integer and I<k> an integer within [-9, +9]),
this is called I<polar pretty-printing>.
=head2 CHANGED IN PERL 5.6
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. (The one
-parameter calling convention also still works.)
+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 parts of the complex number(s). The
-default is C<undef>, which corresponds usually (this is somewhat
-system-dependent) to C<"%.15g">. You can revert to the default by
-setting the format string to C<undef>.
+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' => '%.6f');
+ $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. If
-you want to know the whole truth for a complex number, you must call
-both the class method and the object method:
+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
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
Note that because we are operating on approximations of real numbers,
these errors can happen when merely `too close' to the singularities
-listed above. For example C<tan(2*atan2(1,1)+1e-15)> will die of
-division by zero.
+listed above.
=head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS
=head1 AUTHORS
-Raphael Manfredi <F<Raphael_Manfredi@pobox.com>> and
-Jarkko Hietaniemi <F<jhi@iki.fi>>.
+Daniel S. Lewart <F<d-lewart@uiuc.edu>>
-Extensive patches by Daniel S. Lewart <F<d-lewart@uiuc.edu>>.
+Original authors Raphael Manfredi <F<Raphael_Manfredi@pobox.com>> and
+Jarkko Hietaniemi <F<jhi@iki.fi>>
=cut