From: Jarkko Hietaniemi Date: Thu, 2 Jan 1997 23:17:59 +0000 (+1200) Subject: Math::Complex update X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=commitdiff_plain;h=0c721ce2c118567e416384760f5ad175d956d33b;p=p5sagit%2Fp5-mst-13.2.git Math::Complex update --- diff --git a/lib/Math/Complex.pm b/lib/Math/Complex.pm index fce53f7..30194eb 100644 --- a/lib/Math/Complex.pm +++ b/lib/Math/Complex.pm @@ -1,34 +1,48 @@ # $RCSFile$ # # Complex numbers and associated mathematical functions -# -- Raphael Manfredi, Sept 1996 +# -- Raphael Manfredi, September 1996 +# -- Jarkko Hietaniemi, March 1997 require Exporter; package Math::Complex; @ISA = qw(Exporter); +use strict; + +use vars qw(@EXPORT $package $display + $pi $i $ilog10 $logn %logn); + @EXPORT = qw( pi i Re Im arg + sqrt exp log ln log10 logn cbrt root - tan cotan asin acos atan acotan - sinh cosh tanh cotanh asinh acosh atanh acotanh + tan + cosec csc sec cotan cot + asin acos atan + acosec acsc asec acotan acot + sinh cosh tanh + cosech csch sech cotanh coth + asinh acosh atanh + acosech acsch asech acotanh acoth cplx cplxe ); use overload - '+' => \&plus, - '-' => \&minus, - '*' => \&multiply, - '/' => \÷, + '+' => \&plus, + '-' => \&minus, + '*' => \&multiply, + '/' => \÷, '**' => \&power, '<=>' => \&spaceship, 'neg' => \&negate, - '~' => \&conjugate, + '~' => \&conjugate, 'abs' => \&abs, 'sqrt' => \&sqrt, 'exp' => \&exp, 'log' => \&log, 'sin' => \&sin, 'cos' => \&cos, + 'tan' => \&tan, 'atan2' => \&atan2, qw("" stringify); @@ -87,7 +101,7 @@ sub new { &make } # For backward compatibility only. # sub cplx { my ($re, $im) = @_; - return $package->make($re, $im); + return $package->make($re, defined $im ? $im : 0); } # @@ -98,7 +112,7 @@ sub cplx { # sub cplxe { my ($rho, $theta) = @_; - return $package->emake($rho, $theta); + return $package->emake($rho, defined $theta ? $theta : 0); } # @@ -129,8 +143,10 @@ sub i () { # 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] } @@ -168,8 +184,9 @@ sub update_polar { # sub plus { my ($z1, $z2, $regular) = @_; + $z2 = cplx($z2, 0) unless ref $z2; my ($re1, $im1) = @{$z1->cartesian}; - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2); + my ($re2, $im2) = @{$z2->cartesian}; unless (defined $regular) { $z1->set_cartesian([$re1 + $re2, $im1 + $im2]); return $z1; @@ -184,8 +201,9 @@ sub plus { # sub minus { my ($z1, $z2, $inverted) = @_; + $z2 = cplx($z2, 0) unless ref $z2; my ($re1, $im1) = @{$z1->cartesian}; - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2); + my ($re2, $im2) = @{$z2->cartesian}; unless (defined $inverted) { $z1->set_cartesian([$re1 - $re2, $im1 - $im2]); return $z1; @@ -203,7 +221,8 @@ sub minus { sub multiply { my ($z1, $z2, $regular) = @_; my ($r1, $t1) = @{$z1->polar}; - my ($r2, $t2) = ref $z2 ? @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi); + my ($r2, $t2) = ref $z2 ? + @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi); unless (defined $regular) { $z1->set_polar([$r1 * $r2, $t1 + $t2]); return $z1; @@ -212,6 +231,20 @@ sub multiply { } # +# divbyzero +# +# Die on division by zero. +# +sub divbyzero { + warn $package . '::' . "$_[0]: Division by zero.\n"; + warn "(Because in the definition of $_[0], $_[1] is 0)\n" + if (defined $_[1]); + my @up = caller(1); + my $dmess = "Died at $up[1] line $up[2].\n"; + die $dmess; +} + +# # (divide) # # Computes z1/z2. @@ -219,14 +252,20 @@ sub multiply { sub divide { my ($z1, $z2, $inverted) = @_; my ($r1, $t1) = @{$z1->polar}; - my ($r2, $t2) = ref $z2 ? @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi); + 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; } - return $inverted ? - (ref $z1)->emake($r2 / $r1, $t2 - $t1) : - (ref $z1)->emake($r1 / $r2, $t1 - $t2); + if ($inverted) { + divbyzero "$z2/0" if ($r1 == 0); + return (ref $z1)->emake($r2 / $r1, $t2 - $t1); + } else { + divbyzero "$z1/0" if ($r2 == 0); + return (ref $z1)->emake($r1 / $r2, $t1 - $t2); + } } # @@ -248,8 +287,9 @@ sub power { # sub spaceship { my ($z1, $z2, $inverted) = @_; + $z2 = cplx($z2, 0) unless ref $z2; my ($re1, $im1) = @{$z1->cartesian}; - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2); + my ($re2, $im2) = @{$z2->cartesian}; my $sgn = $inverted ? -1 : 1; return $sgn * ($re1 <=> $re2) if $re1 != $re2; return $sgn * ($im1 <=> $im2); @@ -292,6 +332,7 @@ sub conjugate { # sub abs { my ($z) = @_; + return abs($z) unless ref $z; my ($r, $t) = @{$z->polar}; return abs($r); } @@ -303,7 +344,7 @@ sub abs { # sub arg { my ($z) = @_; - return 0 unless ref $z; + return ($z < 0 ? pi : 0) unless ref $z; my ($r, $t) = @{$z->polar}; return $t; } @@ -311,10 +352,11 @@ sub arg { # # (sqrt) # -# Compute sqrt(z) (positive only). +# Compute sqrt(z). # sub sqrt { my ($z) = @_; + $z = cplx($z, 0) unless ref $z; my ($r, $t) = @{$z->polar}; return (ref $z)->emake(sqrt($r), $t/2); } @@ -322,11 +364,11 @@ sub sqrt { # # cbrt # -# Compute cbrt(z) (cubic root, primary only). +# Compute cbrt(z) (cubic root). # sub cbrt { my ($z) = @_; - return $z ** (1/3) unless ref $z; + return cplx($z, 0) ** (1/3) unless ref $z; my ($r, $t) = @{$z->polar}; return (ref $z)->emake($r**(1/3), $t/3); } @@ -389,6 +431,7 @@ 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); } @@ -400,21 +443,32 @@ sub exp { # sub log { my ($z) = @_; + $z = cplx($z, 0) unless ref $z; my ($r, $t) = @{$z->polar}; + my ($x, $y) = @{$z->cartesian}; + $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); } # +# ln +# +# Alias for log(). +# +sub ln { Math::Complex::log(@_) } + +# # log10 # # Compute log10(z). # sub log10 { my ($z) = @_; - $log10 = log(10) unless defined $log10; - return log($z) / $log10 unless ref $z; + my $ilog10 = 1 / log(10) unless defined $ilog10; + return log(cplx($z, 0)) * $ilog10 unless ref $z; my ($r, $t) = @{$z->polar}; - return (ref $z)->make(log($r) / $log10, $t / $log10); + return (ref $z)->make(log($r) * $ilog10, $t * $ilog10); } # @@ -424,9 +478,10 @@ 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) / log($n); + return log($z) / $logn; } # @@ -439,7 +494,8 @@ sub cos { 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); + return (ref $z)->make(cos($x) * ($ey + $ey_1)/2, + sin($x) * ($ey_1 - $ey)/2); } # @@ -452,7 +508,8 @@ sub sin { 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); + return (ref $z)->make(sin($x) * ($ey + $ey_1)/2, + cos($x) * ($ey - $ey_1)/2); } # @@ -462,29 +519,70 @@ sub sin { # sub tan { my ($z) = @_; - return sin($z) / cos($z); + my $cz = cos($z); + divbyzero "tan($z)", "cos($z)" if ($cz == 0); + return sin($z) / $cz; } # -# cotan +# 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; +} + # -# Computes cotan(z) = 1 / tan(z). +# cosec # -sub cotan { +# Alias for csc(). +# +sub cosec { Math::Complex::csc(@_) } + +# +# cot +# +# Computes cot(z) = 1 / tan(z). +# +sub cot { my ($z) = @_; - return cos($z) / sin($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) = @_; - my $cz = $z*$z - 1; - $cz = cplx($cz, 0) if !ref $cz && $cz < 0; # Force complex if <0 - return ~i * log($z + sqrt $cz); # ~i is -i + $z = cplx($z, 0) unless ref $z; + return ~i * log($z + (Re($z) * Im($z) > 0 ? 1 : -1) * sqrt($z*$z - 1)); } # @@ -494,43 +592,80 @@ sub acos { # sub asin { my ($z) = @_; - my $cz = 1 - $z*$z; - $cz = cplx($cz, 0) if !ref $cz && $cz < 0; # Force complex if <0 - return ~i * log(i * $z + sqrt $cz); # ~i is -i + $z = cplx($z, 0) unless ref $z; + return ~i * log(i * $z + sqrt(1 - $z*$z)); } # # atan # -# Computes the arc tagent atan(z) = i/2 log((i+z) / (i-z)). +# Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)). # sub atan { my ($z) = @_; - return i/2 * log((i + $z) / (i - $z)); + divbyzero "atan($z)", "i - $z" if ($z == i); + return i/2*log((i + $z) / (i - $z)); } # -# acotan +# asec +# +# Computes the arc secant asec(z) = acos(1 / z). +# +sub asec { + my ($z) = @_; + return acos(1 / $z); +} + +# +# acosec +# +# Computes the arc cosecant sec(z) = asin(1 / z). +# +sub acosec { + my ($z) = @_; + return asin(1 / $z); +} + +# +# acsc # -# Computes the arc cotangent acotan(z) = -i/2 log((i+z) / (z-i)) +# Alias for acosec(). +# +sub acsc { Math::Complex::acosec(@_) } + # -sub acotan { +# acot +# +# Computes the arc cotangent acot(z) = -i/2 log((i+z) / (z-i)) +# +sub acot { my ($z) = @_; + divbyzero "acot($z)", "$z - i" if ($z == i); return i/-2 * log((i + $z) / ($z - i)); } # +# 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 ($x, $y) = ref $z ? @{$z->cartesian} : ($z); + $z = cplx($z, 0) unless ref $z; + my ($x, $y) = @{$z->cartesian}; 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); + return (ref $z)->make(cos($y) * ($ex + $ex_1)/2, + sin($y) * ($ex - $ex_1)/2); } # @@ -540,11 +675,13 @@ sub cosh { # sub sinh { my ($z) = @_; - my ($x, $y) = ref $z ? @{$z->cartesian} : ($z); + $z = cplx($z, 0) unless ref $z; + my ($x, $y) = @{$z->cartesian}; 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); + return (ref $z)->make(cos($y) * ($ex - $ex_1)/2, + sin($y) * ($ex + $ex_1)/2); } # @@ -554,29 +691,70 @@ sub sinh { # sub tanh { my ($z) = @_; - return sinh($z) / cosh($z); + my $cz = cosh($z); + divbyzero "tanh($z)", "cosh($z)" if ($cz == 0); + return sinh($z) / $cz; } # -# cotanh +# 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). # -# Comptutes the hyperbolic cotangent cotanh(z) = cosh(z) / 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(@_) } + # -sub cotanh { +# coth +# +# Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z). +# +sub coth { my ($z) = @_; - return cosh($z) / sinh($z); + my $sz = sinh($z); + divbyzero "coth($z)", "sinh($z)" if ($sz == 0); + return cosh($z) / $sz; } # +# cotanh +# +# Alias for coth(). +# +sub cotanh { Math::Complex::coth(@_) } + +# # acosh # # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)). # sub acosh { my ($z) = @_; - my $cz = $z*$z - 1; - $cz = cplx($cz, 0) if !ref $cz && $cz < 0; # Force complex if <0 - return log($z + sqrt $cz); + $z = cplx($z, 0) unless ref $z; # asinh(-2) + return log($z + sqrt($z*$z - 1)); } # @@ -586,8 +764,8 @@ sub acosh { # sub asinh { my ($z) = @_; - my $cz = $z*$z + 1; # Already complex if <0 - return log($z + sqrt $cz); + $z = cplx($z, 0) unless ref $z; # asinh(-2) + return log($z + sqrt($z*$z + 1)); } # @@ -597,24 +775,62 @@ sub asinh { # sub atanh { my ($z) = @_; + $z = cplx($z, 0) unless ref $z; # atanh(-2) + divbyzero 'atanh(1)', "1 - $z" if ($z == 1); my $cz = (1 + $z) / (1 - $z); - $cz = cplx($cz, 0) if !ref $cz && $cz < 0; # Force complex if <0 return log($cz) / 2; } # -# acotanh +# asech +# +# Computes the hyperbolic arc secant asech(z) = acosh(1 / z). +# +sub asech { + my ($z) = @_; + divbyzero 'asech(0)', $z if ($z == 0); + return acosh(1 / $z); +} + +# +# acsch # -# Computes the arc hyperbolic cotangent acotanh(z) = 1/2 log((1+z) / (z-1)). +# Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z). # -sub acotanh { +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 arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)). +# +sub acoth { + my ($z) = @_; + $z = cplx($z, 0) unless ref $z; # acoth(-2) + divbyzero 'acoth(1)', "$z - 1" if ($z == 1); my $cz = (1 + $z) / ($z - 1); - $cz = cplx($cz, 0) if !ref $cz && $cz < 0; # Force complex if <0 return log($cz) / 2; } # +# acotanh +# +# Alias for acot(). +# +sub acotanh { Math::Complex::acoth(@_) } + +# # (atan2) # # Compute atan(z1/z2). @@ -622,7 +838,7 @@ sub acotanh { sub atan2 { my ($z1, $z2, $inverted) = @_; my ($re1, $im1) = @{$z1->cartesian}; - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2); + my ($re2, $im2) = @{$z2->cartesian}; my $tan; if (defined $inverted && $inverted) { # atan(z2/z1) return pi * ($re2 > 0 ? 1 : -1) if $re1 == 0 && $im1 == 0; @@ -653,7 +869,7 @@ sub display_format { if (ref $self) { # Called as a method $format = shift; - } else { # Regular procedure call + } else { # Regular procedure call $format = $self; undef $self; } @@ -709,7 +925,7 @@ sub stringify_cartesian { elsif ($y == -1) { $im = '-i' } elsif (abs($y) >= 1e-14) { $im = $y . "i" } - my $str; + my $str = ''; $str = $re if defined $re; $str .= "+$im" if defined $im; $str =~ s/\+-/-/; @@ -728,22 +944,24 @@ sub stringify_polar { my $z = shift; my ($r, $t) = @{$z->polar}; my $theta; + my $eps = 1e-14; - return '[0,0]' if $r <= 1e-14; + return '[0,0]' if $r <= $eps; my $tpi = 2 * pi; my $nt = $t / $tpi; $nt = ($nt - int($nt)) * $tpi; $nt += $tpi if $nt < 0; # Range [0, 2pi] - if (abs($nt) <= 1e-14) { $theta = 0 } - elsif (abs(pi-$nt) <= 1e-14) { $theta = 'pi' } + if (abs($nt) <= $eps) { $theta = 0 } + elsif (abs(pi-$nt) <= $eps) { $theta = 'pi' } if (defined $theta) { - $r = int($r + ($r < 0 ? -1 : 1) * 1e-14) - if int(abs($r)) != int(abs($r) + 1e-14); - $theta = int($theta + ($theta < 0 ? -1 : 1) * 1e-14) - if int(abs($theta)) != int(abs($theta) + 1e-14); + $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\]"; } @@ -756,18 +974,20 @@ sub stringify_polar { 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) <= 1e-14) { - $theta = ($nt < 0 ? '-':'').($k == 1 ? 'pi':"${k}pi").'/'.abs($n); + if (abs($kpi/$n - $nt) <= $eps) { + $theta = ($nt < 0 ? '-':''). + ($k == 1 ? 'pi':"${k}pi").'/'.abs($n); last; } } $theta = $nt unless defined $theta; - $r = int($r + ($r < 0 ? -1 : 1) * 1e-14) - if int(abs($r)) != int(abs($r) + 1e-14); - $theta = int($theta + ($theta < 0 ? -1 : 1) * 1e-14) - if int(abs($theta)) != int(abs($theta) + 1e-14); + $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)); return "\[$r,$theta\]"; } @@ -974,24 +1194,41 @@ numbers: logn(z, n) = log(z) / log(n) tan(z) = sin(z) / cos(z) - cotan(z) = 1 / tan(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 + sqrt(z*z-1)) atan(z) = i/2 * log((i+z) / (i-z)) - acotan(z) = -i/2 * log((i+z) / (z-i)) + + acsc(z) = asin(1 / z) + asec(z) = acos(1 / z) + acot(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) - cotanh(z) = 1 / tanh(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)) - acotanh(z) = 1/2 * log((1+z) / (z-1)) -The I function is available to compute all the Ith + acsch(z) = asinh(1 / z) + 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. + +The I function is available to compute all the I roots of some complex, where I is a strictly positive integer. There are exactly I such roots, returned as a list. Getting the number mathematicians call C such that: @@ -1006,10 +1243,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 operation 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 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 @@ -1027,9 +1265,9 @@ if you like. To create a number using the trigonometric 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). -(Mnmemonic: C is used as a notation for complex numbers in the trigonometric -form). +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). It is possible to write: @@ -1104,6 +1342,7 @@ 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. -=head1 AUTHOR +=head1 AUTHORS -Raphael Manfredi > + Raphael Manfredi > + Jarkko Hietaniemi > diff --git a/t/lib/complex.t b/t/lib/complex.t index 1ffd7d5..5dcac6c 100755 --- a/t/lib/complex.t +++ b/t/lib/complex.t @@ -3,7 +3,8 @@ # $RCSfile$ # # Regression tests for the new Math::Complex pacakge -# -- Raphael Manfredi, Sept 1996 +# -- Raphael Manfredi, Septemeber 1996 +# -- Jarkko Hietaniemi Manfredi, March 1997 BEGIN { chdir 't' if -d 't'; @INC = '../lib'; @@ -13,12 +14,13 @@ use Math::Complex; $test = 0; $| = 1; $script = ''; -$epsilon = 1e-10; +my $eps = 1e-4; # for example root() is quite bad while () { - next if /^#/ || /^\s*$/; - chop; - $set_test = 0; # Assume not a test over a set of values + s/^\s+//; + next if $_ eq '' || /^\#/; + chomp; + $test_set = 0; # Assume not a test over a set of values if (/^&(.*)/) { $op = $1; next; @@ -28,14 +30,16 @@ while () { next; } elsif (s/^\|//) { - $set_test = 1; # Requests we loop over the set... + $test_set = 1; # Requests we loop over the set... } my @args = split(/:/); - if ($set_test) { + if ($test_set == 1) { my $i; for ($i = 0; $i < @set; $i++) { - $target = $set[$i]; # complex number - $zvalue = $val[$i]; # textual value as found in set definition + # complex number + $target = $set[$i]; + # textual value as found in set definition + $zvalue = $val[$i]; test($zvalue, $target, @args); } } else { @@ -57,7 +61,7 @@ sub test { } if (defined $z) { $args = "'$op'"; # Really the value - $try = "abs(\$z0 - \$z1) <= 1e-10 ? \$z1 : \$z0"; + $try = "abs(\$z0 - \$z1) <= $eps ? \$z1 : \$z0"; $script .= "\$res = $try; "; $script .= "check($test, $args[0], \$res, \$z$#args, $args);\n"; } else { @@ -114,7 +118,15 @@ sub value { sub check { my ($test, $try, $got, $expected, @z) = @_; - if ("$got" eq "$expected" || ($expected =~ /^-?\d/ && $got == $expected)) { + +# print "# @_\n"; + + if ("$got" eq "$expected" + || + ($expected =~ /^-?\d/ && $got == $expected) + || + (abs($got - $expected) < $eps) + ) { print "ok $test\n"; } else { print "not ok $test\n"; @@ -163,10 +175,24 @@ __END__ [4,pi]:[2,pi/2]:[2,pi/2] [2,pi/2]:[4,pi]:[0.5,-(pi)/2] +&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) @@ -185,6 +211,7 @@ __END__ (3,4):(3,4):1 &sqrt +-9:(0,3) (-100,0):(0,10) (16,-30):(5,-3) @@ -209,46 +236,214 @@ __END__ |'z - ~z':'2*i*Im(z)' |'z * ~z':'abs(z) * abs(z)' -{ (4,3); [3,2]; (-3,4); (0,2); 3; 1; (-5, 0); [2,1] } +{ (2,3); [3,2]; (-3,2); (0,2); 3; 1.2; -3; (-3, 0); (-2, -1); [2,1] } -|'exp(z)':'exp(a) * exp(i * b)' +|'(root(z, 4))[1] ** 4':'z' +|'(root(z, 5))[3] ** 5':'z' +|'(root(z, 8))[7] ** 8':'z' |'abs(z)':'r' -|'sqrt(z) * sqrt(z)':'z' -|'sqrt(z)':'sqrt(r) * exp(i * t/2)' +|'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)' -|'log(z)':'log(r) + i*t' -|'sin(asin(z))':'z' |'cos(acos(z))':'z' -|'tan(atan(z))':'z' -|'cotan(acotan(z))':'z' |'cos(z) ** 2 + sin(z) ** 2':1 -|'cosh(z) ** 2 - sinh(z) ** 2':1 |'cos(z)':'cosh(i*z)' -|'cotan(z)':'1 / tan(z)' -|'cotanh(z)':'1 / tanh(z)' -|'i*sin(z)':'sinh(i*z)' -|'z**z':'exp(z * log(z))' -|'log(exp(z))':'z' +|'cosh(z) ** 2 - sinh(z) ** 2':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, 3)':'log(z) / log(3)' |'logn(z, 2)':'log(z) / log(2)' -|'(root(z, 4))[1] ** 4':'z' -|'(root(z, 8))[7] ** 8':'z' +|'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; (-1,0.5); (0,0.5); 0.5; (2, 0) } +{ (1,1); [1,0.5]; (-2, -1); 2; -3; (-1,0.5); (0,0.5); 0.5; (2, 0); (-1, -2) } -|'sinh(asinh(z))':'z' |'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' -|'cotanh(acotanh(z))':'z' -{ (0.2,-0.4); [1,0.5]; -1.2; (-1,0.5); (0,-0.5); 0.5; (1.1, 0) } +{ (0.2,-0.4); [1,0.5]; -1.2; (-1,0.5); 0.5; (1.1, 0) } -|'asin(sin(z))':'z' |'acos(cos(z)) ** 2':'z * z' -|'atan(tan(z))':'z' -|'asinh(sinh(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' +&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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 3):( 0.14694666622553, -0.23182380450040) +(-2, 3):( -0.14694666622553, -0.23182380450040) +(-2,-3):( -0.14694666622553, 0.23182380450040) +( 2,-3):( 0.14694666622553, 0.23182380450040) + +# eof