complex.pm

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

#
# Complex numbers and associated mathematical functions
# -- Raphael Manfredi	Since Sep 1996
# -- Jarkko Hietaniemi	Since Mar 1997
# -- Daniel S. Lewart	Since Sep 1997
#

package Math::Complex;

our($VERSION, @ISA, @EXPORT, %EXPORT_TAGS, $Inf);

$VERSION = 1.31;

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;

my $i;
my %LOGN;

require Exporter;

@ISA = qw(Exporter);

my @trig = qw(
	      pi
	      tan
	      csc cosec sec cot cotan
	      asin acos atan
	      acsc acosec asec acot acotan
	      sinh cosh tanh
	      csch cosech sech coth cotanh
	      asinh acosh atanh
	      acsch acosech asech acoth acotanh
	     );

@EXPORT = (qw(
	     i Re Im rho theta arg
	     sqrt log ln
	     log10 logn cbrt root
	     cplx cplxe
	     ),
	   @trig);

%EXPORT_TAGS = (
    'trig' => [@trig],
);

use overload
	'+'	=> \&plus,
	'-'	=> \&minus,
	'*'	=> \&multiply,
	'/'	=> \&divide,
	'**'	=> \&power,
	'=='	=> \&numeq,
	'<=>'	=> \&spaceship,
	'neg'	=> \&negate,
	'~'	=> \&conjugate,
	'abs'	=> \&abs,
	'sqrt'	=> \&sqrt,
	'exp'	=> \&exp,
	'log'	=> \&log,
	'sin'	=> \&sin,
	'cos'	=> \&cos,
	'tan'	=> \&tan,
	'atan2'	=> \&atan2,
	qw("" stringify);

#
# Package "privates"
#

my %DISPLAY_FORMAT = ('style' => 'cartesian',
		      'polar_pretty_print' => 1);
my $eps            = 1e-14;		# Epsilon

#
# Object attributes (internal):
#	cartesian	[real, imaginary] -- cartesian form
#	polar		[rho, theta] -- polar form
#	c_dirty		cartesian form not up-to-date
#	p_dirty		polar form not up-to-date
#	display		display format (package's global when not set)
#

# Die on bad *make() arguments.

sub _cannot_make {
    die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n";
}

#
# ->make
#
# Create a new complex number (cartesian form)
#
sub make {
	my $self = bless {}, shift;
	my ($re, $im) = @_;
	my $rre = ref $re;
	if ( $rre ) {
	    if ( $rre eq ref $self ) {
		$re = Re($re);
	    } else {
		_cannot_make("real part", $rre);
	    }
	}
	my $rim = ref $im;
	if ( $rim ) {
	    if ( $rim eq ref $self ) {
		$im = Im($im);
	    } else {
		_cannot_make("imaginary part", $rim);
	    }
	}
	$self->{'cartesian'} = [ $re, $im ];
	$self->{c_dirty} = 0;
	$self->{p_dirty} = 1;
	$self->display_format('cartesian');
	return $self;
}

#
# ->emake
#
# Create a new complex number (exponential form)
#
sub emake {
	my $self = bless {}, shift;
	my ($rho, $theta) = @_;
	my $rrh = ref $rho;
	if ( $rrh ) {
	    if ( $rrh eq ref $self ) {
		$rho = rho($rho);
	    } else {
		_cannot_make("rho", $rrh);
	    }
	}
	my $rth = ref $theta;
	if ( $rth ) {
	    if ( $rth eq ref $self ) {
		$theta = theta($theta);
	    } else {
		_cannot_make("theta", $rth);
	    }
	}
	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;
}

sub new { &make }		# For backward compatibility only.

#
# cplx
#
# Creates a complex number from a (re, im) tuple.
# This avoids the burden of writing Math::Complex->make(re, im).
#
sub cplx {
	my ($re, $im) = @_;
	return __PACKAGE__->make($re, defined $im ? $im : 0);
}

#
# cplxe
#
# Creates a complex number from a (rho, theta) tuple.
# This avoids the burden of writing Math::Complex->emake(rho, theta).
#
sub cplxe {
	my ($rho, $theta) = @_;
	return __PACKAGE__->emake($rho, defined $theta ? $theta : 0);
}

#
# pi
#
# The number defined as pi = 180 degrees
#
sub pi () { 4 * CORE::atan2(1, 1) }

#
# pit2
#
# The full circle
#
sub pit2 () { 2 * pi }

#
# pip2
#
# The quarter circle
#
sub pip2 () { pi / 2 }

#
# deg1
#
# One degree in radians, used in stringify_polar.
#

sub deg1 () { pi / 180 }

#
# uplog10
#
# Used in log10().
#
sub uplog10 () { 1 / CORE::log(10) }

#
# i
#
# The number defined as i*i = -1;
#
sub i () {
        return $i if ($i);
	$i = bless {};
	$i->{'cartesian'} = [0, 1];
	$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
#

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] }

#
# ->update_cartesian
#
# Recompute and return the cartesian form, given accurate polar form.
#
sub update_cartesian {
	my $self = shift;
	my ($r, $t) = @{$self->{'polar'}};
	$self->{c_dirty} = 0;
	return $self->{'cartesian'} = [$r * CORE::cos($t), $r * CORE::sin($t)];
}

#
#
# ->update_polar
#
# Recompute and return the polar form, given accurate cartesian form.
#
sub update_polar {
	my $self = shift;
	my ($x, $y) = @{$self->{'cartesian'}};
	$self->{p_dirty} = 0;
	return $self->{'polar'} = [0, 0] if $x == 0 && $y == 0;
	return $self->{'polar'} = [CORE::sqrt($x*$x + $y*$y),
				   CORE::atan2($y, $x)];
}

#
# (plus)
#
# Computes z1+z2.
#
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]);
		return $z1;
	}
	return (ref $z1)->make($re1 + $re2, $im1 + $im2);
}

#
# (minus)
#
# Computes z1-z2.
#
sub minus {
	my ($z1, $z2, $inverted) = @_;
	my ($re1, $im1) = @{$z1->cartesian};
	$z2 = cplx($z2) unless ref $z2;
	my ($re2, $im2) = @{$z2->cartesian};
	unless (defined $inverted) {
		$z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
		return $z1;
	}
	return $inverted ?
		(ref $z1)->make($re2 - $re1, $im2 - $im1) :
		(ref $z1)->make($re1 - $re2, $im1 - $im2);

}

#
# (multiply)
#
# Computes z1*z2.
#
sub multiply {
        my ($z1, $z2, $regular) = @_;
	if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
	    # if both polar better use polar to avoid rounding errors
	    my ($r1, $t1) = @{$z1->polar};
	    my ($r2, $t2) = @{$z2->polar};
	    my $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);
	    }
	}
}

#
# _divbyzero
#
# Die on division by zero.
#
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 .= "is 0)\n";
    }

    my @up = caller(1);

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

    die $mess;
}

#
# (divide)
#
# Computes z1/z2.
#
sub divide {
	my ($z1, $z2, $inverted) = @_;
	if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) {
	    # if both polar better use polar to avoid rounding errors
	    my ($r1, $t1) = @{$z1->polar};
	    my ($r2, $t2) = @{$z2->polar};
	    my $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 {
	    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);
		}
	    }
	}
}

#
# (power)
#
# Computes z1**z2 = exp(z2 * log z1)).
#
sub power {
	my ($z1, $z2, $inverted) = @_;
	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.
#
sub spaceship {
	my ($z1, $z2, $inverted) = @_;
	my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
	my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
	my $sgn = $inverted ? -1 : 1;
	return $sgn * ($re1 <=> $re2) if $re1 != $re2;
	return $sgn * ($im1 <=> $im2);
}

#
# (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.
#
sub negate {
	my ($z) = @_;
	if ($z->{c_dirty}) {
		my ($r, $t) = @{$z->polar};
		$t = ($t <= 0) ? $t + pi : $t - pi;
		return (ref $z)->emake($r, $t);
	}
	my ($re, $im) = @{$z->cartesian};
	return (ref $z)->make(-$re, -$im);
}

#
# (conjugate)
#
# Compute complex's conjugate.
#
sub conjugate {
	my ($z) = @_;
	if ($z->{c_dirty}) {
		my ($r, $t) = @{$z->polar};
		return (ref $z)->emake($r, -$t);
	}
	my ($re, $im) = @{$z->cartesian};
	return (ref $z)->make($re, -$im);
}

#
# (abs)
#
# Compute or set complex's norm (rho).
#
sub abs {
	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 or set complex's argument (theta).
#
sub arg {
	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;
}

#
# (sqrt)
#
# 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) = @_;
	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(CORE::sqrt($r), $t/2);
}

#
# cbrt
#
# Compute cbrt(z) (cubic root).
#
# Why are we not returning three values?  The same answer as for sqrt().
#
sub cbrt {
	my ($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 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;
}

#
# root
#
# Computes all nth root for z, returning an array whose size is n.
# `n' must be a positive integer.
#
# The roots are given by (for k = 0..n-1):