X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=blobdiff_plain;f=lib%2FMath%2FComplex.pm;h=19d30b062e3ece2d6e193dbc0598c4b95de2f59e;hb=d1be9408a3c14848d30728674452e191ba5fffaa;hp=7d5a014203266a321ce414d2f57742395ecc99f7;hpb=b0399d3c2d4e16b936b7db5a1aa3f810e4a39864;p=p5sagit%2Fp5-mst-13.2.git diff --git a/lib/Math/Complex.pm b/lib/Math/Complex.pm index 7d5a014..19d30b0 100644 --- a/lib/Math/Complex.pm +++ b/lib/Math/Complex.pm @@ -1,26 +1,49 @@ -# $RCSFile$ # # Complex numbers and associated mathematical functions -# -- Raphael Manfredi, September 1996 -# -- Jarkko Hietaniemi, March-April 1997 +# -- Raphael Manfredi Since Sep 1996 +# -- Jarkko Hietaniemi Since Mar 1997 +# -- Daniel S. Lewart Since Sep 1997 +# -require Exporter; package Math::Complex; +our($VERSION, @ISA, @EXPORT, %EXPORT_TAGS, $Inf); + +$VERSION = 1.32; + +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; -use vars qw($VERSION @ISA - @EXPORT %EXPORT_TAGS - $package $display - $i $logn %logn); +my $i; +my %LOGN; -@ISA = qw(Exporter); +require Exporter; -$VERSION = 1.01; +@ISA = qw(Exporter); my @trig = qw( pi - sin cos tan + tan csc cosec sec cot cotan asin acos atan acsc acosec asec acot acotan @@ -31,8 +54,8 @@ my @trig = qw( ); @EXPORT = (qw( - i Re Im arg - sqrt exp log ln + i Re Im rho theta arg + sqrt log ln log10 logn cbrt root cplx cplxe ), @@ -48,6 +71,7 @@ use overload '*' => \&multiply, '/' => \÷, '**' => \&power, + '==' => \&numeq, '<=>' => \&spaceship, 'neg' => \&negate, '~' => \&conjugate, @@ -62,11 +86,12 @@ use overload qw("" stringify); # -# Package globals +# Package "privates" # -$package = 'Math::Complex'; # Package name -$display = 'cartesian'; # Default display format +my %DISPLAY_FORMAT = ('style' => 'cartesian', + 'polar_pretty_print' => 1); +my $eps = 1e-14; # Epsilon # # Object attributes (internal): @@ -77,6 +102,12 @@ $display = 'cartesian'; # Default display format # 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"; +} + # # ->make # @@ -85,9 +116,26 @@ $display = 'cartesian'; # Default display format sub make { my $self = bless {}, shift; my ($re, $im) = @_; - $self->{'cartesian'} = [$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; } @@ -99,10 +147,30 @@ sub make { sub emake { my $self = bless {}, shift; my ($rho, $theta) = @_; - $theta += pi() if $rho < 0; - $self->{'polar'} = [abs($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); + } + } + 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; } @@ -116,7 +184,7 @@ sub new { &make } # For backward compatibility only. # sub cplx { my ($re, $im) = @_; - return $package->make($re, defined $im ? $im : 0); + return __PACKAGE__->make($re, defined $im ? $im : 0); } # @@ -127,24 +195,44 @@ sub cplx { # sub cplxe { my ($rho, $theta) = @_; - return $package->emake($rho, defined $theta ? $theta : 0); + return __PACKAGE__->emake($rho, defined $theta ? $theta : 0); } # # pi # -# The number defined as 2 * pi = 360 degrees +# The number defined as pi = 180 degrees # +sub pi () { 4 * CORE::atan2(1, 1) } -use constant pi => 4 * atan2(1, 1); +# +# pit2 +# +# The full circle +# +sub pit2 () { 2 * pi } # -# log2inv +# pip2 # -# Used in log10(). +# The quarter circle # +sub pip2 () { pi / 2 } -use constant log10inv => 1 / log(10); +# +# deg1 +# +# One degree in radians, used in stringify_polar. +# + +sub deg1 () { pi / 180 } + +# +# uplog10 +# +# Used in log10(). +# +sub uplog10 () { 1 / CORE::log(10) } # # i @@ -155,13 +243,20 @@ sub i () { return $i if ($i); $i = bless {}; $i->{'cartesian'} = [0, 1]; - $i->{'polar'} = [1, pi/2]; + $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 # @@ -182,7 +277,7 @@ sub update_cartesian { my $self = shift; my ($r, $t) = @{$self->{'polar'}}; $self->{c_dirty} = 0; - return $self->{'cartesian'} = [$r * cos $t, $r * sin $t]; + return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)]; } # @@ -196,7 +291,8 @@ sub update_polar { my ($x, $y) = @{$self->{'cartesian'}}; $self->{p_dirty} = 0; return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0; - return $self->{'polar'} = [sqrt($x*$x + $y*$y), atan2($y, $x)]; + return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y), + CORE::atan2($y, $x)]; } # @@ -207,6 +303,7 @@ sub update_polar { 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]); @@ -223,7 +320,8 @@ sub plus { sub minus { my ($z1, $z2, $inverted) = @_; my ($re1, $im1) = @{$z1->cartesian}; - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); + $z2 = cplx($z2) unless ref $z2; + my ($re2, $im2) = @{$z2->cartesian}; unless (defined $inverted) { $z1->set_cartesian([$re1 - $re2, $im1 - $im2]); return $z1; @@ -231,6 +329,7 @@ sub minus { return $inverted ? (ref $z1)->make($re2 - $re1, $im2 - $im1) : (ref $z1)->make($re1 - $re2, $im1 - $im2); + } # @@ -239,48 +338,46 @@ sub minus { # Computes z1*z2. # sub multiply { - my ($z1, $z2, $regular) = @_; - my ($r1, $t1) = @{$z1->polar}; - my ($r2, $t2) = ref $z2 ? - @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi); - unless (defined $regular) { - $z1->set_polar([$r1 * $r2, $t1 + $t2]); + 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 -= pit2 } + elsif ($t <= -pi()) { $t += pit2 } + 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); + } } - return (ref $z1)->emake($r1 * $r2, $t1 + $t2); } # -# divbyzero +# _divbyzero # # Die on division by zero. # -sub divbyzero { +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 .= "$_[1] " unless ("$_[1]" eq '0'); $mess .= "is 0)\n"; } my @up = caller(1); - - $mess .= "Died at $up[1] line $up[2].\n"; - die $mess; -} - -# -# zerotozero -# -# Die on zero raised to the zeroth. -# -sub zerotozero { - my $mess = "The zero raised to the zeroth power is not defined.\n"; - - my @up = caller(1); - $mess .= "Died at $up[1] line $up[2].\n"; die $mess; @@ -293,20 +390,45 @@ sub zerotozero { # sub divide { my ($z1, $z2, $inverted) = @_; - my ($r1, $t1) = @{$z1->polar}; - my ($r2, $t2) = ref $z2 ? - @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi); - unless (defined $inverted) { - divbyzero "$z1/0" if ($r2 == 0); - $z1->set_polar([$r1 / $r2, $t1 - $t2]); - return $z1; - } - if ($inverted) { - divbyzero "$z2/0" if ($r1 == 0); - return (ref $z1)->emake($r2 / $r1, $t2 - $t1); + 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 -= pit2 } + elsif ($t <= -pi()) { $t += pit2 } + 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 } + return (ref $z1)->emake($r1 / $r2, $t); + } } else { - divbyzero "$z1/0" if ($r2 == 0); - return (ref $z1)->emake($r1 / $r2, $t1 - $t2); + 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); + } + } } } @@ -317,16 +439,26 @@ sub divide { # sub power { my ($z1, $z2, $inverted) = @_; - zerotozero if ($z1 == 0 and $z2 == 0); - return exp($z1 * log $z2) if defined $inverted && $inverted; - return exp($z2 * log $z1); + 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. +# Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i. # sub spaceship { my ($z1, $z2, $inverted) = @_; @@ -338,6 +470,19 @@ sub spaceship { } # +# (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. @@ -346,7 +491,8 @@ sub negate { my ($z) = @_; if ($z->{c_dirty}) { my ($r, $t) = @{$z->polar}; - return (ref $z)->emake($r, pi + $t); + $t = ($t <= 0) ? $t + pi : $t - pi; + return (ref $z)->emake($r, $t); } my ($re, $im) = @{$z->cartesian}; return (ref $z)->make(-$re, -$im); @@ -370,25 +516,52 @@ sub conjugate { # # (abs) # -# Compute complex's norm (rho). +# Compute or set complex's norm (rho). # sub abs { - my ($z) = @_; - return abs($z) unless ref $z; - my ($r, $t) = @{$z->polar}; - return abs($r); + 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 -= pit2 } + elsif ($$theta <= -pi()) { $$theta += pit2 } } # # arg # -# Compute complex's argument (theta). +# Compute or set complex's argument (theta). # sub arg { - my ($z) = @_; - return ($z < 0 ? pi : 0) unless ref $z; - my ($r, $t) = @{$z->polar}; - return $t; + 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; } # @@ -396,11 +569,23 @@ sub arg { # # 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) = @_; - $z = cplx($z, 0) unless ref $z; + 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(sqrt($r), $t/2); + return (ref $z)->emake(CORE::sqrt($r), $t/2); } # @@ -408,11 +593,32 @@ sub sqrt { # # Compute cbrt(z) (cubic root). # +# Why are we not returning three values? The same answer as for sqrt(). +# sub cbrt { my ($z) = @_; - return cplx($z, 0) ** (1/3) unless ref $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 (ref $z)->emake($r**(1/3), $t/3); + 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; } # @@ -427,17 +633,19 @@ sub cbrt { # sub root { my ($z, $n) = @_; - $n = int($n + 0.5); - return undef unless $n > 0; - my ($r, $t) = ref $z ? @{$z->polar} : (abs($z), $z >= 0 ? 0 : pi); + _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 = 2 * pi / $n; + my $theta_inc = pit2 / $n; my $rho = $r ** (1/$n); my $theta; - my $complex = ref($z) || $package; + my $cartesian = ref $z && $z->{c_dirty} == 0; for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) { - push(@root, $complex->emake($rho, $theta)); + my $w = cplxe($rho, $theta); + # Yes, $cartesian is loop invariant. + push @root, $cartesian ? cplx(@{$w->cartesian}) : $w; } return @root; } @@ -445,25 +653,53 @@ sub root { # # Re # -# Return Re(z). +# Return or set Re(z). # sub Re { - my ($z) = @_; + my ($z, $Re) = @_; return $z unless ref $z; - my ($re, $im) = @{$z->cartesian}; - return $re; + if (defined $Re) { + $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ]; + $z->{c_dirty} = 0; + $z->{p_dirty} = 1; + } else { + return ${$z->cartesian}[0]; + } } # # Im # -# Return Im(z). +# Return or set Im(z). # sub Im { - my ($z) = @_; + my ($z, $Im) = @_; return 0 unless ref $z; - my ($re, $im) = @{$z->cartesian}; - return $im; + 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(@_); } # @@ -473,9 +709,29 @@ sub Im { # sub exp { my ($z) = @_; - $z = cplx($z, 0) unless ref $z; my ($x, $y) = @{$z->cartesian}; - return (ref $z)->emake(exp($x), $y); + 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; } # @@ -485,12 +741,15 @@ sub exp { # sub log { my ($z) = @_; - $z = cplx($z, 0) unless ref $z; - my ($x, $y) = @{$z->cartesian}; + unless (ref $z) { + _logofzero("log") if $z == 0; + return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi); + } my ($r, $t) = @{$z->polar}; - $t -= 2 * pi if ($t > pi() and $x < 0); - $t += 2 * pi if ($t < -pi() and $x < 0); - return (ref $z)->make(log($r), $t); + _logofzero("log") if $r == 0; + if ($t > pi()) { $t -= pit2 } + elsif ($t <= -pi()) { $t += pit2 } + return (ref $z)->make(CORE::log($r), $t); } # @@ -507,11 +766,7 @@ sub ln { Math::Complex::log(@_) } # sub log10 { - my ($z) = @_; - - return log(cplx($z, 0)) * log10inv unless ref $z; - my ($r, $t) = @{$z->polar}; - return (ref $z)->make(log($r) * log10inv, $t * log10inv); + return Math::Complex::log($_[0]) * uplog10; } # @@ -522,9 +777,9 @@ sub log10 { sub logn { my ($z, $n) = @_; $z = cplx($z, 0) unless ref $z; - my $logn = $logn{$n}; - $logn = $logn{$n} = log($n) unless defined $logn; # Cache log(n) - return log($z) / $logn; + my $logn = $LOGN{$n}; + $logn = $LOGN{$n} = CORE::log($n) unless defined $logn; # Cache log(n) + return &log($z) / $logn; } # @@ -534,12 +789,14 @@ sub logn { # sub cos { my ($z) = @_; - $z = cplx($z, 0) unless ref $z; + return CORE::cos($z) unless ref $z; my ($x, $y) = @{$z->cartesian}; - my $ey = exp($y); - my $ey_1 = 1 / $ey; - return (ref $z)->make(cos($x) * ($ey + $ey_1)/2, - sin($x) * ($ey_1 - $ey)/2); + 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); } # @@ -549,12 +806,14 @@ sub cos { # sub sin { my ($z) = @_; - $z = cplx($z, 0) unless ref $z; + return CORE::sin($z) unless ref $z; my ($x, $y) = @{$z->cartesian}; - my $ey = exp($y); - my $ey_1 = 1 / $ey; - return (ref $z)->make(sin($x) * ($ey + $ey_1)/2, - cos($x) * ($ey - $ey_1)/2); + 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); } # @@ -564,9 +823,9 @@ sub sin { # sub tan { my ($z) = @_; - my $cz = cos($z); - divbyzero "tan($z)", "cos($z)" if ($cz == 0); - return sin($z) / $cz; + my $cz = &cos($z); + _divbyzero "tan($z)", "cos($z)" if $cz == 0; + return &sin($z) / $cz; } # @@ -576,8 +835,8 @@ sub tan { # sub sec { my ($z) = @_; - my $cz = cos($z); - divbyzero "sec($z)", "cos($z)" if ($cz == 0); + my $cz = &cos($z); + _divbyzero "sec($z)", "cos($z)" if ($cz == 0); return 1 / $cz; } @@ -588,8 +847,8 @@ sub sec { # sub csc { my ($z) = @_; - my $sz = sin($z); - divbyzero "csc($z)", "sin($z)" if ($sz == 0); + my $sz = &sin($z); + _divbyzero "csc($z)", "sin($z)" if ($sz == 0); return 1 / $sz; } @@ -603,13 +862,13 @@ sub cosec { Math::Complex::csc(@_) } # # cot # -# Computes cot(z) = 1 / tan(z). +# 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; + my $sz = &sin($z); + _divbyzero "cot($z)", "sin($z)" if ($sz == 0); + return &cos($z) / $sz; } # @@ -625,9 +884,23 @@ sub cotan { Math::Complex::cot(@_) } # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)). # sub acos { - my ($z) = @_; + 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; - return ~i * log($z + (Re($z) * Im($z) > 0 ? 1 : -1) * sqrt($z*$z - 1)); + 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); } # @@ -636,9 +909,23 @@ sub acos { # Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)). # sub asin { - my ($z) = @_; + 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; - return ~i * log(i * $z + sqrt(1 - $z*$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); } # @@ -648,9 +935,13 @@ sub asin { # sub atan { my ($z) = @_; - $z = cplx($z, 0) unless ref $z; - divbyzero "atan($z)", "i - $z" if ($z == i); - return i/2*log((i + $z) / (i - $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; } # @@ -660,18 +951,18 @@ sub atan { # sub asec { my ($z) = @_; - divbyzero "asec($z)", $z if ($z == 0); + _divbyzero "asec($z)", $z if ($z == 0); return acos(1 / $z); } # # acsc # -# Computes the arc cosecant sec(z) = asin(1 / z). +# Computes the arc cosecant acsc(z) = asin(1 / z). # sub acsc { my ($z) = @_; - divbyzero "acsc($z)", $z if ($z == 0); + _divbyzero "acsc($z)", $z if ($z == 0); return asin(1 / $z); } @@ -685,13 +976,16 @@ sub acosec { Math::Complex::acsc(@_) } # # acot # -# Computes the arc cotangent acot(z) = -i/2 log((i+z) / (z-i)) +# Computes the arc cotangent acot(z) = atan(1 / z) # sub acot { my ($z) = @_; - $z = cplx($z, 0) unless ref $z; - divbyzero "acot($z)", "$z - i" if ($z == i); - return i/-2 * log((i + $z) / ($z - i)); + _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); } # @@ -708,12 +1002,16 @@ sub acotan { Math::Complex::acot(@_) } # sub cosh { my ($z) = @_; - my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0); - my $ex = exp($x); - my $ex_1 = 1 / $ex; - return ($ex + $ex_1)/2 unless ref $z; - return (ref $z)->make(cos($y) * ($ex + $ex_1)/2, - sin($y) * ($ex - $ex_1)/2); + my $ex; + unless (ref $z) { + $ex = CORE::exp($z); + return $ex ? ($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); } # @@ -723,12 +1021,19 @@ sub cosh { # sub sinh { my ($z) = @_; - my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0); - my $ex = exp($x); - my $ex_1 = 1 / $ex; - return ($ex - $ex_1)/2 unless ref $z; - return (ref $z)->make(cos($y) * ($ex - $ex_1)/2, - sin($y) * ($ex + $ex_1)/2); + my $ex; + unless (ref $z) { + return 0 if $z == 0; + $ex = CORE::exp($z); + return $ex ? ($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); } # @@ -739,7 +1044,7 @@ sub sinh { sub tanh { my ($z) = @_; my $cz = cosh($z); - divbyzero "tanh($z)", "cosh($z)" if ($cz == 0); + _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0); return sinh($z) / $cz; } @@ -751,7 +1056,7 @@ sub tanh { sub sech { my ($z) = @_; my $cz = cosh($z); - divbyzero "sech($z)", "cosh($z)" if ($cz == 0); + _divbyzero "sech($z)", "cosh($z)" if ($cz == 0); return 1 / $cz; } @@ -763,7 +1068,7 @@ sub sech { sub csch { my ($z) = @_; my $sz = sinh($z); - divbyzero "csch($z)", "sinh($z)" if ($sz == 0); + _divbyzero "csch($z)", "sinh($z)" if ($sz == 0); return 1 / $sz; } @@ -782,7 +1087,7 @@ sub cosech { Math::Complex::csch(@_) } sub coth { my ($z) = @_; my $sz = sinh($z); - divbyzero "coth($z)", "sinh($z)" if ($sz == 0); + _divbyzero "coth($z)", "sinh($z)" if $sz == 0; return cosh($z) / $sz; } @@ -800,19 +1105,45 @@ sub cotanh { Math::Complex::coth(@_) } # sub acosh { my ($z) = @_; - $z = cplx($z, 0) unless ref $z; - return log($z + sqrt($z*$z - 1)); + 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 arc hyperbolic sine asinh(z) = log(z + sqrt(z*z-1)) +# Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z+1)) # sub asinh { my ($z) = @_; - $z = cplx($z, 0) unless ref $z; - return log($z + sqrt($z*$z + 1)); + 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); } # @@ -822,10 +1153,13 @@ sub asinh { # sub atanh { my ($z) = @_; - divbyzero 'atanh(1)', "1 - $z" if ($z == 1); - $z = cplx($z, 0) unless ref $z; - my $cz = (1 + $z) / (1 - $z); - return log($cz) / 2; + 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)); } # @@ -835,7 +1169,7 @@ sub atanh { # sub asech { my ($z) = @_; - divbyzero 'asech(0)', $z if ($z == 0); + _divbyzero 'asech(0)', "$z" if ($z == 0); return acosh(1 / $z); } @@ -846,7 +1180,7 @@ sub asech { # sub acsch { my ($z) = @_; - divbyzero 'acsch(0)', $z if ($z == 0); + _divbyzero 'acsch(0)', $z if ($z == 0); return asinh(1 / $z); } @@ -864,10 +1198,14 @@ sub acosech { Math::Complex::acsch(@_) } # sub acoth { my ($z) = @_; - divbyzero 'acoth(1)', "$z - 1" if ($z == 1); - $z = cplx($z, 0) unless ref $z; - my $cz = (1 + $z) / ($z - 1); - return log($cz) / 2; + _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; } # @@ -884,51 +1222,69 @@ sub acotanh { Math::Complex::acoth(@_) } # sub atan2 { my ($z1, $z2, $inverted) = @_; - my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0); - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); - my $tan; - if (defined $inverted && $inverted) { # atan(z2/z1) - return pi * ($re2 > 0 ? 1 : -1) if $re1 == 0 && $im1 == 0; - $tan = $z2 / $z1; + my ($re1, $im1, $re2, $im2); + if ($inverted) { + ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); + ($re2, $im2) = @{$z1->cartesian}; } else { - return pi * ($re1 > 0 ? 1 : -1) if $re2 == 0 && $im2 == 0; - $tan = $z1 / $z2; + ($re1, $im1) = @{$z1->cartesian}; + ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); } - return atan($tan); + if ($im2 == 0) { + return CORE::atan2($re1, $re2) if $im1 == 0; + return ($im1<=>0) * pip2 if $re2 == 0; + } + 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); } # # display_format # ->display_format # -# Set (fetch if no argument) display format for all complex numbers that -# don't happen to have overrriden it via ->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 a method, this actually sets the display format for +# 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 $format = undef; - - if (ref $self) { # Called as a method - $format = shift; - } else { # Regular procedure call - $format = $self; - undef $self; + 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 (defined $self) { - return defined $self->{display} ? $self->{display} : $display - unless defined $format; - return $self->{display} = $format; + if (ref $self) { # Called as an object method + $self->{display_format} = { %display_format }; + return + wantarray ? + %{$self->{display_format}} : + $self->{display_format}->{style}; } - return $display unless defined $format; - return $display = $format; + # Called as a class method + %DISPLAY_FORMAT = %display_format; + return + wantarray ? + %DISPLAY_FORMAT : + $DISPLAY_FORMAT{style}; } # @@ -943,12 +1299,12 @@ sub display_format { # sub stringify { my ($z) = shift; - my $format; - $format = $display; - $format = $z->{display} if defined $z->{display}; + my $style = $z->display_format; + + $style = $DISPLAY_FORMAT{style} unless defined $style; - return $z->stringify_polar if $format =~ /^p/i; + return $z->stringify_polar if $style =~ /^p/i; return $z->stringify_cartesian; } @@ -962,26 +1318,58 @@ sub stringify_cartesian { my ($x, $y) = @{$z->cartesian}; my ($re, $im); - $x = int($x + ($x < 0 ? -1 : 1) * 1e-14) - if int(abs($x)) != int(abs($x) + 1e-14); - $y = int($y + ($y < 0 ? -1 : 1) * 1e-14) - if int(abs($y)) != int(abs($y) + 1e-14); + my %format = $z->display_format; + my $format = $format{format}; + + if ($x) { + if ($x =~ /^NaN[QS]?$/i) { + $re = $x; + } else { + if ($x =~ /^-?$Inf$/oi) { + $re = $x; + } else { + $re = defined $format ? sprintf($format, $x) : $x; + } + } + } else { + undef $re; + } - $re = "$x" if abs($x) >= 1e-14; - if ($y == 1) { $im = 'i' } - elsif ($y == -1) { $im = '-i' } - elsif (abs($y) >= 1e-14) { $im = $y . "i" } + if ($y) { + if ($y =~ /^(NaN[QS]?)$/i) { + $im = $y; + } else { + if ($y =~ /^-?$Inf$/oi) { + $im = $y; + } else { + $im = + defined $format ? + sprintf($format, $y) : + ($y == 1 ? "" : ($y == -1 ? "-" : $y)); + } + } + $im .= "i"; + } else { + undef $im; + } - my $str = ''; - $str = $re if defined $re; - $str .= "+$im" if defined $im; - $str =~ s/\+-/-/; - $str =~ s/^\+//; - $str = '0' unless $str; + 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 # @@ -991,57 +1379,53 @@ sub stringify_polar { my $z = shift; my ($r, $t) = @{$z->polar}; my $theta; - my $eps = 1e-14; - return '[0,0]' if $r <= $eps; + my %format = $z->display_format; + my $format = $format{format}; - my $tpi = 2 * pi; - my $nt = $t / $tpi; - $nt = ($nt - int($nt)) * $tpi; - $nt += $tpi if $nt < 0; # Range [0, 2pi] - - if (abs($nt) <= $eps) { $theta = 0 } - elsif (abs(pi-$nt) <= $eps) { $theta = 'pi' } - - if (defined $theta) { - $r = int($r + ($r < 0 ? -1 : 1) * $eps) - if int(abs($r)) != int(abs($r) + $eps); - $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps) - if ($theta ne 'pi' and - int(abs($theta)) != int(abs($theta) + $eps)); - return "\[$r,$theta\]"; + if ($t =~ /^NaN[QS]?$/i || $t =~ /^-?$Inf$/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; + # - # Okay, number is not a real. Try to identify pi/n and friends... + # Try to identify pi/n and friends. # - $nt -= $tpi if $nt > pi; - my ($n, $k, $kpi); - - for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) { - $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5); - if (abs($kpi/$n - $nt) <= $eps) { - $theta = ($nt < 0 ? '-':''). - ($k == 1 ? 'pi':"${k}pi").'/'.abs($n); - last; + $t -= int(CORE::abs($t) / pit2) * pit2; + + 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; } + } } - $theta = $nt unless defined $theta; - - $r = int($r + ($r < 0 ? -1 : 1) * $eps) - if int(abs($r)) != int(abs($r) + $eps); - $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps) - if ($theta !~ m(^-?\d*pi/\d+$) and - int(abs($theta)) != int(abs($theta) + $eps)); + if (defined $format) { + $r = sprintf($format, $r); + $theta = sprintf($format, $theta) unless defined $theta; + } else { + $theta = $t unless defined $theta; + } - return "\[$r,$theta\]"; + return "[$r,$theta]"; } 1; __END__ +=pod + =head1 NAME Math::Complex - complex numbers and associated mathematical functions @@ -1049,7 +1433,7 @@ 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); @@ -1126,7 +1510,7 @@ between this form and the cartesian form C is immediate: which is also expressed by this formula: - z = rho * exp(i * theta) = rho * (cos theta + i * sin theta) + z = rho * exp(i * theta) = rho * (cos theta + i * sin theta) In other words, it's the projection of the vector onto the I and I axes. Mathematicians call I the I or I and I @@ -1136,8 +1520,8 @@ noted C. 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 substractions. Real numbers are on the I -axis, and therefore I is zero. +suited for additions and subtractions. Real numbers are on the I +axis, and therefore I is zero or I. 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 @@ -1146,8 +1530,8 @@ they keep their natural meaning when there is no imaginary part, provided the number is within their definition set. For instance, the C routine which computes the square root of -its argument is only defined for positive real numbers and yields a -positive real number (it is an application from B to B). +its argument is only defined for non-negative real numbers and yields a +non-negative real number (it is an application from B to B). If we allow it to return a complex number, then it can be extended to negative real numbers to become an application from B to B (the set of complex numbers): @@ -1160,14 +1544,15 @@ 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 is always positive, so C<[x,pi]> is really C<-x>, a -negative number) -and the above definition states that +Indeed, a negative real number can be noted C<[x,pi]> (the modulus +I 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 returns only one of the solutions: if you want the both, +use the C function. All the common mathematical functions defined on real numbers that are extended to complex numbers share that same property of working @@ -1176,7 +1561,7 @@ be called an extension, would it?). A I operation possible on a complex number that is the identity for real numbers is called the I, 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 @@ -1220,14 +1605,13 @@ the following (overloaded) operations are supported on complex numbers: z1 * z2 = (r1 * r2) * exp(i * (t1 + t2)) z1 / z2 = (r1 / r2) * exp(i * (t1 - t2)) z1 ** z2 = exp(z2 * log z1) - ~z1 = a - bi - abs(z1) = r1 = sqrt(a*a + b*b) - sqrt(z1) = sqrt(r1) * exp(i * t1/2) - exp(z1) = exp(a) * exp(i * b) - log(z1) = log(r1) + i*t1 - sin(z1) = 1/2i (exp(i * z1) - exp(-i * z1)) - cos(z1) = 1/2 (exp(i * z1) + exp(-i * z1)) - abs(z1) = r1 + ~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(z1, z2) = atan(z1/z2) The following extra operations are supported on both real and complex @@ -1236,6 +1620,7 @@ numbers: Re(z) = a Im(z) = b arg(z) = t + abs(z) = r cbrt(z) = z ** (1/3) log10(z) = log(z) / log(10) @@ -1248,12 +1633,12 @@ numbers: cot(z) = 1 / tan(z) asin(z) = -i * log(i*z + sqrt(1-z*z)) - acos(z) = -i * log(z + sqrt(z*z-1)) + 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) = -i/2 * log((i+z) / (z-i)) + 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)) @@ -1262,7 +1647,7 @@ numbers: 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)) @@ -1271,10 +1656,13 @@ numbers: asech(z) = acosh(1 / z) acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1)) -I, I, I, I, I, I, I, -I, I, have aliases I, I, I, -I, I, I, I, I, I, -respectively. +I, I, I, I, I, I, I, I, +I, I, I, have aliases I, I, I, +I, I, I, I, I, I, +I, I, respectively. C, C, C, C, +C, and C can be used also as mutators. The C +returns only one of the solutions: if you want all three, use the +C function. The I function is available to compute all the I roots of some complex, where I is a strictly positive integer. @@ -1291,11 +1679,11 @@ The Ith root for C is given by: (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n) -The I comparison operator 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. +The I comparison operator, E=E, 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 @@ -1308,51 +1696,98 @@ if you know the cartesian form of the number, or $z = 3 + 4*i; -if you like. To create a number using the trigonometric form, use either: +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). (Mnmemonic: C is used as a -notation for complex numbers in the trigonometric form). +(in radians, the full circle is 2*pi). (Mnemonic: C 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 positive (it represents the distance to the origin in the complex -plane). +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 +either the C or C: the appropriate component of +the argument will be used. + + $z1 = cplx(-2, 1); + $z2 = cplx($z1, 4); =head1 STRINGIFICATION When printed, a complex number is usually shown under its cartesian -form I, but there are legitimate cases where the polar format +style I, but there are legitimate cases where the polar style I<[r,t]> is more appropriate. -By calling the routine C and supplying either -C<"polar"> or C<"cartesian">, you override the default display format, -which is C<"cartesian">. Not supplying any argument returns the current -setting. +By calling the class method C 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 method instead. As before, not supplying any argument -returns the current display format for this number. Otherwise whatever you -specify will be the new display format for I particular number. +returns the current display style for this number. Otherwise whatever you +specify will be the new display style for I 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 = (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 format attempts to emphasize arguments like I -(where I is a positive integer and I an integer within [-9,+9]). +The polar style attempts to emphasize arguments like I +(where I is a positive integer and I an integer within [-9, +9]), +this is called I. + +=head2 CHANGED IN PERL 5.6 + +The C class method and the corresponding +C 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 to C. + + # 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 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 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 @@ -1373,65 +1808,82 @@ Here are some examples: $k = exp(i * 2*pi/3); print "$j - $k = ", $j - $k, "\n"; -=head1 ERRORS DUE TO DIVISION BY ZERO + $z->Re(3); # Re, Im, arg, abs, + $j->arg(2); # (the last two aka rho, theta) + # can be used also as mutators. + +=head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO The division (/) and the following functions - tan - sec - csc - cot - asec - acsc - atan - acot - tanh - sech - csch - coth - atanh - asech - acsch - acoth + 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. These situations cause fatal runtime errors looking like this +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 ... -For the C, C, C, C, C, C, C, -C, the argument cannot be C<0> (zero). For the C, -C, the argument cannot be C<1> (one). For the C, C, -the argument cannot be C (the imaginary unit). For the C, -C, C, C, the argument cannot be I, where -I is any integer. +or -=head1 BUGS + atanh(-1): Logarithm of zero. + Died at... -Saying C exports many mathematical routines in the -caller environment and even overrides some (C, C, C, -C, C). This is construed as a feature by the Authors, -actually... ;-) +For the C, C, C, C, C, C, C, +C, C, the argument cannot be C<0> (zero). For the +logarithmic functions and the C, C, the argument cannot +be C<1> (one). For the C, C, the argument cannot be +C<-1> (minus one). For the C, C, the argument cannot be +C (the imaginary unit). For the C, C, the argument +cannot be C<-i> (the negative imaginary unit). For the C, +C, C, the argument cannot be I, where I +is any integer. + +Note that because we are operating on approximations of real numbers, +these errors can happen when merely `too close' to the singularities +listed above. -The code is not optimized for speed, although we try to use the cartesian -form for addition-like operators and the trigonometric form for all -multiplication-like operators. +=head1 ERRORS DUE TO INDIGESTIBLE ARGUMENTS -The arg() routine does not ensure the angle is within the range [-pi,+pi] -(a side effect caused by multiplication and division using the trigonometric -representation). +The C and C 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 exports many mathematical routines in the +caller environment and even overrides some (C, C). +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 AUTHORS - Raphael Manfredi > - Jarkko Hietaniemi > +Raphael Manfredi > and +Jarkko Hietaniemi >. + +Extensive patches by Daniel S. Lewart >. =cut +1; + # eof