complex.pm

MSYS在windows下模拟了一个类unix的终端

#
# z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
#
sub root {
	my ($z, $n) = @_;
	_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;
	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;
	}
	return @root;
}

#
# Re
#
# Return or set Re(z).
#
sub Re {
	my ($z, $Re) = @_;
	return $z unless ref $z;
	if (defined $Re) {
	    $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ];
	    $z->{c_dirty} = 0;
	    $z->{p_dirty} = 1;
	} else {
	    return ${$z->cartesian}[0];
	}
}

#
# Im
#
# Return or set Im(z).
#
sub Im {
	my ($z, $Im) = @_;
	return 0 unless ref $z;
	if (defined $Im) {
	    $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ];
	    $z->{c_dirty} = 0;
	    $z->{p_dirty} = 1;
	} else {
	    return ${$z->cartesian}[1];
	}
}

#
# rho
#
# Return or set rho(w).
#
sub rho {
    Math::Complex::abs(@_);
}

#
# theta
#
# Return or set theta(w).
#
sub theta {
    Math::Complex::arg(@_);
}

#
# (exp)
#
# Computes exp(z).
#
sub exp {
	my ($z) = @_;
	my ($x, $y) = @{$z->cartesian};
	return (ref $z)->emake(CORE::exp($x), $y);
}

#
# _logofzero
#
# Die on logarithm of zero.
#
sub _logofzero {
    my $mess = "$_[0]: Logarithm of zero.\n";

    if (defined $_[1]) {
	$mess .= "(Because in the definition of $_[0], the argument ";
	$mess .= "$_[1] " unless ($_[1] eq '0');
	$mess .= "is 0)\n";
    }

    my @up = caller(1);

    $mess .= "Died at $up[1] line $up[2].\n";

    die $mess;
}

#
# (log)
#
# Compute log(z).
#
sub log {
	my ($z) = @_;
	unless (ref $z) {
	    _logofzero("log") if $z == 0;
	    return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi);
	}
	my ($r, $t) = @{$z->polar};
	_logofzero("log") if $r == 0;
	if    ($t >   pi()) { $t -= pit2 }
	elsif ($t <= -pi()) { $t += pit2 }
	return (ref $z)->make(CORE::log($r), $t);
}

#
# ln
#
# Alias for log().
#
sub ln { Math::Complex::log(@_) }

#
# log10
#
# Compute log10(z).
#

sub log10 {
	return Math::Complex::log($_[0]) * uplog10;
}

#
# logn
#
# Compute logn(z,n) = log(z) / log(n)
#
sub logn {
	my ($z, $n) = @_;
	$z = cplx($z, 0) unless ref $z;
	my $logn = $LOGN{$n};
	$logn = $LOGN{$n} = CORE::log($n) unless defined $logn;	# Cache log(n)
	return &log($z) / $logn;
}

#
# (cos)
#
# Compute cos(z) = (exp(iz) + exp(-iz))/2.
#
sub cos {
	my ($z) = @_;
	return CORE::cos($z) unless ref $z;
	my ($x, $y) = @{$z->cartesian};
	my $ey = CORE::exp($y);
	my $sx = CORE::sin($x);
	my $cx = CORE::cos($x);
	my $ey_1 = $ey ? 1 / $ey : $Inf;
	return (ref $z)->make($cx * ($ey + $ey_1)/2,
			      $sx * ($ey_1 - $ey)/2);
}

#
# (sin)
#
# Compute sin(z) = (exp(iz) - exp(-iz))/2.
#
sub sin {
	my ($z) = @_;
	return CORE::sin($z) unless ref $z;
	my ($x, $y) = @{$z->cartesian};
	my $ey = CORE::exp($y);
	my $sx = CORE::sin($x);
	my $cx = CORE::cos($x);
	my $ey_1 = $ey ? 1 / $ey : $Inf;
	return (ref $z)->make($sx * ($ey + $ey_1)/2,
			      $cx * ($ey - $ey_1)/2);
}

#
# tan
#
# Compute tan(z) = sin(z) / cos(z).
#
sub tan {
	my ($z) = @_;
	my $cz = &cos($z);
	_divbyzero "tan($z)", "cos($z)" if $cz == 0;
	return &sin($z) / $cz;
}

#
# sec
#
# Computes the secant sec(z) = 1 / cos(z).
#
sub sec {
	my ($z) = @_;
	my $cz = &cos($z);
	_divbyzero "sec($z)", "cos($z)" if ($cz == 0);
	return 1 / $cz;
}

#
# csc
#
# Computes the cosecant csc(z) = 1 / sin(z).
#
sub csc {
	my ($z) = @_;
	my $sz = &sin($z);
	_divbyzero "csc($z)", "sin($z)" if ($sz == 0);
	return 1 / $sz;
}

#
# cosec
#
# Alias for csc().
#
sub cosec { Math::Complex::csc(@_) }

#
# cot
#
# Computes cot(z) = cos(z) / sin(z).
#
sub cot {
	my ($z) = @_;
	my $sz = &sin($z);
	_divbyzero "cot($z)", "sin($z)" if ($sz == 0);
	return &cos($z) / $sz;
}

#
# cotan
#
# Alias for cot().
#
sub cotan { Math::Complex::cot(@_) }

#
# acos
#
# Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
#
sub acos {
	my $z = $_[0];
	return CORE::atan2(CORE::sqrt(1-$z*$z), $z)
	    if (! ref $z) && CORE::abs($z) <= 1;
	$z = cplx($z, 0) unless ref $z;
	my ($x, $y) = @{$z->cartesian};
	return 0 if $x == 1 && $y == 0;
	my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
	my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
	my $alpha = ($t1 + $t2)/2;
	my $beta  = ($t1 - $t2)/2;
	$alpha = 1 if $alpha < 1;
	if    ($beta >  1) { $beta =  1 }
	elsif ($beta < -1) { $beta = -1 }
	my $u = CORE::atan2(CORE::sqrt(1-$beta*$beta), $beta);
	my $v = CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
	$v = -$v if $y > 0 || ($y == 0 && $x < -1);
	return (ref $z)->make($u, $v);
}

#
# asin
#
# Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
#
sub asin {
	my $z = $_[0];
	return CORE::atan2($z, CORE::sqrt(1-$z*$z))
	    if (! ref $z) && CORE::abs($z) <= 1;
	$z = cplx($z, 0) unless ref $z;
	my ($x, $y) = @{$z->cartesian};
	return 0 if $x == 0 && $y == 0;
	my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y);
	my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y);
	my $alpha = ($t1 + $t2)/2;
	my $beta  = ($t1 - $t2)/2;
	$alpha = 1 if $alpha < 1;
	if    ($beta >  1) { $beta =  1 }
	elsif ($beta < -1) { $beta = -1 }
	my $u =  CORE::atan2($beta, CORE::sqrt(1-$beta*$beta));
	my $v = -CORE::log($alpha + CORE::sqrt($alpha*$alpha-1));
	$v = -$v if $y > 0 || ($y == 0 && $x < -1);
	return (ref $z)->make($u, $v);
}

#
# atan
#
# Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
#
sub atan {
	my ($z) = @_;
	return CORE::atan2($z, 1) unless ref $z;
	my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
	return 0 if $x == 0 && $y == 0;
	_divbyzero "atan(i)"  if ( $z == i);
	_logofzero "atan(-i)" if (-$z == i); # -i is a bad file test...
	my $log = &log((i + $z) / (i - $z));
	return ip2 * $log;
}

#
# asec
#
# Computes the arc secant asec(z) = acos(1 / z).
#
sub asec {
	my ($z) = @_;
	_divbyzero "asec($z)", $z if ($z == 0);
	return acos(1 / $z);
}

#
# acsc
#
# Computes the arc cosecant acsc(z) = asin(1 / z).
#
sub acsc {
	my ($z) = @_;
	_divbyzero "acsc($z)", $z if ($z == 0);
	return asin(1 / $z);
}

#
# acosec
#
# Alias for acsc().
#
sub acosec { Math::Complex::acsc(@_) }

#
# acot
#
# Computes the arc cotangent acot(z) = atan(1 / z)
#
sub acot {
	my ($z) = @_;
	_divbyzero "acot(0)"  if $z == 0;
	return ($z >= 0) ? CORE::atan2(1, $z) : CORE::atan2(-1, -$z)
	    unless ref $z;
	_divbyzero "acot(i)"  if ($z - i == 0);
	_logofzero "acot(-i)" if ($z + i == 0);
	return atan(1 / $z);
}

#
# acotan
#
# Alias for acot().
#
sub acotan { Math::Complex::acot(@_) }

#
# cosh
#
# Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
#
sub cosh {
	my ($z) = @_;
	my $ex;
	unless (ref $z) {
	    $ex = CORE::exp($z);
	    return $ex ? ($ex + 1/$ex)/2 : $Inf;
	}
	my ($x, $y) = @{$z->cartesian};
	$ex = CORE::exp($x);
	my $ex_1 = $ex ? 1 / $ex : $Inf;
	return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
			      CORE::sin($y) * ($ex - $ex_1)/2);
}

#
# sinh
#
# Computes the hyperbolic sine sinh(z) = (exp(z) - exp(-z))/2.
#
sub sinh {
	my ($z) = @_;
	my $ex;
	unless (ref $z) {
	    return 0 if $z == 0;
	    $ex = CORE::exp($z);
	    return $ex ? ($ex - 1/$ex)/2 : "-$Inf";
	}
	my ($x, $y) = @{$z->cartesian};
	my $cy = CORE::cos($y);
	my $sy = CORE::sin($y);
	$ex = CORE::exp($x);
	my $ex_1 = $ex ? 1 / $ex : $Inf;
	return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
			      CORE::sin($y) * ($ex + $ex_1)/2);
}

#
# tanh
#
# Computes the hyperbolic tangent tanh(z) = sinh(z) / cosh(z).
#
sub tanh {
	my ($z) = @_;
	my $cz = cosh($z);
	_divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
	return sinh($z) / $cz;
}

#
# sech
#
# Computes the hyperbolic secant sech(z) = 1 / cosh(z).
#
sub sech {
	my ($z) = @_;
	my $cz = cosh($z);
	_divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
	return 1 / $cz;
}

#
# csch
#
# Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
#
sub csch {
	my ($z) = @_;
	my $sz = sinh($z);
	_divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
	return 1 / $sz;
}

#
# cosech
#
# Alias for csch().
#
sub cosech { Math::Complex::csch(@_) }

#
# coth
#
# Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
#
sub coth {
	my ($z) = @_;
	my $sz = sinh($z);
	_divbyzero "coth($z)", "sinh($z)" if $sz == 0;
	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) = @_;
	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))
#
sub asinh {
	my ($z) = @_;
	unless (ref $z) {
	    my $t = $z + CORE::sqrt($z*$z + 1);
	    return CORE::log($t) if $t;
	}
	my $t = &sqrt($z * $z + 1) + $z;
	# Try Taylor if looking bad (this usually means that
	# $z was large negative, therefore the sqrt is really
	# close to abs(z), summing that with z...)
	$t = 1/(2 * $z) - 1/(8 * $z**3) + 1/(16 * $z**5) - 5/(128 * $z**7)
	    if $t == 0;
	return &log($t);
}

#
# atanh
#
# Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
#
sub atanh {
	my ($z) = @_;
	unless (ref $z) {
	    return CORE::log((1 + $z)/(1 - $z))/2 if CORE::abs($z) < 1;
	    $z = cplx($z, 0);
	}
	_divbyzero 'atanh(1)',  "1 - $z" if (1 - $z == 0);
	_logofzero 'atanh(-1)'           if (1 + $z == 0);
	return 0.5 * &log((1 + $z) / (1 - $z));
}

#
# asech
#
# Computes the 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 hyperbolic arc cosecant acsch(z) = asinh(1 / z).
#
sub acsch {
	my ($z) = @_;
	_divbyzero 'acsch(0)', $z if ($z == 0);
	return asinh(1 / $z);
}

#
# acosech
#
# Alias for acosh().
#
sub acosech { Math::Complex::acsch(@_) }

#
# acoth
#
# Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
#
sub acoth {
	my ($z) = @_;
	_divbyzero 'acoth(0)'            if ($z == 0);
	unless (ref $z) {
	    return CORE::log(($z + 1)/($z - 1))/2 if CORE::abs($z) > 1;
	    $z = cplx($z, 0);
	}
	_divbyzero 'acoth(1)',  "$z - 1" if ($z - 1 == 0);
	_logofzero 'acoth(-1)', "1 + $z" if (1 + $z == 0);
	return &log((1 + $z) / ($z - 1)) / 2;
}

#
# acotanh
#
# Alias for acot().
#
sub acotanh { Math::Complex::acoth(@_) }

#
# (atan2)
#
# Compute atan(z1/z2).
#
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};
	} else {
	    ($re1, $im1) = @{$z1->cartesian};
	    ($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;
	}
	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 (get if no argument) the display format for all complex numbers that
# don't happen to have overridden it via ->display_format
#
# When called as an object method, this actually sets the display format for
# the current object.
#
# Valid object formats are 'c' and 'p' for cartesian and polar. The first
# letter is used actually, so the type can be fully spelled out for clarity.
#
sub display_format {
	my $self  = shift;
	my %display_format = %DISPLAY_FORMAT;