diff options
Diffstat (limited to 'lib/Math/Complex.pm')
-rw-r--r-- | lib/Math/Complex.pm | 136 |
1 files changed, 136 insertions, 0 deletions
diff --git a/lib/Math/Complex.pm b/lib/Math/Complex.pm new file mode 100644 index 0000000000..a5a40b2486 --- /dev/null +++ b/lib/Math/Complex.pm @@ -0,0 +1,136 @@ +# +# Perl5 Package for complex numbers +# +# 1994 by David Nadler +# Coding know-how provided by Tom Christiansen, Tim Bunce, and Larry Wall +# sqrt() added by Tom Christiansen; beware should have two roots, +# but only returns one. (use wantarray?) +# +# +# The functions "Re", "Im", and "arg" are provided. +# "~" is used as the conjugation operator and "abs" is overloaded. +# +# Transcendental functions overloaded: so far only sin, cos, and exp. +# + +package Math::Complex; + +require Exporter; + +@ISA = ('Exporter'); + +# just to make use happy + +%OVERLOAD= ( + '+' => sub { my($x1,$y1,$x2,$y2) = (@{$_[0]},@{$_[1]}); + bless [ $x1+$x2, $y1+$y2]; + }, + + '-' => sub { my($x1,$y1,$x2,$y2) = (@{$_[0]},@{$_[1]}); + bless [ $x1-$x2, $y1-$y2]; + }, + + '*' => sub { my($x1,$y1,$x2,$y2) = (@{$_[0]},@{$_[1]}); + bless [ $x1*$x2-$y1*$y2,$x1*$y2+$x2*$y1]; + }, + + '/' => sub { my($x1,$y1,$x2,$y2) = (@{$_[0]},@{$_[1]}); + my $q = $x2*$x2+$y2*$y2; + bless [($x1*$x2+$y1*$y2)/$q, ($y1*$x2-$y2*$x1)/$q]; + }, + + 'neg' => sub { my($x,$y) = @{$_[0]}; bless [ -$x, -$y]; + }, + + '~' => sub { my($x,$y) = @{$_[0]}; bless [ $x, -$y]; + }, + + 'abs' => sub { my($x,$y) = @{$_[0]}; sqrt $x*$x+$y*$y; + }, + + 'cos' => sub { my($x,$y) = @{$_[0]}; + my ($ab,$c,$s) = (exp $y, cos $x, sin $x); + my $abr = 1/(2*$ab); $ab /= 2; + bless [ ($abr+$ab)*$c, ($abr-$ab)*$s]; + }, + + 'sin' => sub { my($x,$y) = @{$_[0]}; + my ($ab,$c,$s) = (exp $y, cos $x, sin $x); + my $abr = 1/(2*$ab); $ab /= 2; + bless [ (-$abr-$ab)*$s, ($abr-$ab)*$c]; + }, + + 'exp' => sub { my($x,$y) = @{$_[0]}; + my ($ab,$c,$s) = (exp $x, cos $y, sin $y); + bless [ $ab*$c, $ab*$s ]; + }, + + 'sqrt' => sub { + my($zr,$zi) = @{$_[0]}; + my ($x, $y, $r, $w); + my $c = new Math::Complex (0,0); + if (($zr == 0) && ($zi == 0)) { + # nothing, $c already set + } + else { + $x = abs($zr); + $y = abs($zi); + if ($x >= $y) { + $r = $y/$x; + $w = sqrt($x) * sqrt(0.5*(1.0+sqrt(1.0+$r*$r))); + } + else { + $r = $x/$y; + $w = sqrt($y) * sqrt($y) * sqrt(0.5*($r+sqrt(1.0+$r*$r))); + } + if ( $zr >= 0) { + @$c = ($w, $zi/(2 * $w) ); + } + else { + $c->[1] = ($zi >= 0) ? $w : -$w; + $c->[0] = $zi/(2.0* $c->[1]); + } + } + return $c; + }, + + qw("" stringify) +); + +sub new { + shift; + my @C = @_; + bless \@C; +} + +sub Re { + my($x,$y) = @{$_[0]}; + $x; +} + +sub Im { + my($x,$y) = @{$_[0]}; + $y; +} + +sub arg { + my($x,$y) = @{$_[0]}; + atan2($y,$x); +} + +sub stringify { + my($x,$y) = @{$_[0]}; + my($re,$im); + + $re = $x if ($x); + if ($y == 1) {$im = 'i';} + elsif ($y == -1){$im = '-i';} + elsif ($y) {$im = "${y}i"; } + + local $_ = $re.'+'.$im; + s/\+-/-/; + s/^\+//; + s/[\+-]$//; + $_ = 0 if ($_ eq ''); + return $_; +} |