summaryrefslogtreecommitdiff
path: root/lib/Math/Complex.pm
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Math/Complex.pm')
-rw-r--r--lib/Math/Complex.pm136
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 $_;
+}