This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
perl 5.000
[perl5.git] / lib / Math / Complex.pm
CommitLineData
a0d0e21e
LW
1#
2# Perl5 Package for complex numbers
3#
4# 1994 by David Nadler
5# Coding know-how provided by Tom Christiansen, Tim Bunce, and Larry Wall
6# sqrt() added by Tom Christiansen; beware should have two roots,
7# but only returns one. (use wantarray?)
8#
9#
10# The functions "Re", "Im", and "arg" are provided.
11# "~" is used as the conjugation operator and "abs" is overloaded.
12#
13# Transcendental functions overloaded: so far only sin, cos, and exp.
14#
15
16package Math::Complex;
17
18require Exporter;
19
20@ISA = ('Exporter');
21
22# just to make use happy
23
24%OVERLOAD= (
25 '+' => sub { my($x1,$y1,$x2,$y2) = (@{$_[0]},@{$_[1]});
26 bless [ $x1+$x2, $y1+$y2];
27 },
28
29 '-' => sub { my($x1,$y1,$x2,$y2) = (@{$_[0]},@{$_[1]});
30 bless [ $x1-$x2, $y1-$y2];
31 },
32
33 '*' => sub { my($x1,$y1,$x2,$y2) = (@{$_[0]},@{$_[1]});
34 bless [ $x1*$x2-$y1*$y2,$x1*$y2+$x2*$y1];
35 },
36
37 '/' => sub { my($x1,$y1,$x2,$y2) = (@{$_[0]},@{$_[1]});
38 my $q = $x2*$x2+$y2*$y2;
39 bless [($x1*$x2+$y1*$y2)/$q, ($y1*$x2-$y2*$x1)/$q];
40 },
41
42 'neg' => sub { my($x,$y) = @{$_[0]}; bless [ -$x, -$y];
43 },
44
45 '~' => sub { my($x,$y) = @{$_[0]}; bless [ $x, -$y];
46 },
47
48 'abs' => sub { my($x,$y) = @{$_[0]}; sqrt $x*$x+$y*$y;
49 },
50
51 'cos' => sub { my($x,$y) = @{$_[0]};
52 my ($ab,$c,$s) = (exp $y, cos $x, sin $x);
53 my $abr = 1/(2*$ab); $ab /= 2;
54 bless [ ($abr+$ab)*$c, ($abr-$ab)*$s];
55 },
56
57 'sin' => sub { my($x,$y) = @{$_[0]};
58 my ($ab,$c,$s) = (exp $y, cos $x, sin $x);
59 my $abr = 1/(2*$ab); $ab /= 2;
60 bless [ (-$abr-$ab)*$s, ($abr-$ab)*$c];
61 },
62
63 'exp' => sub { my($x,$y) = @{$_[0]};
64 my ($ab,$c,$s) = (exp $x, cos $y, sin $y);
65 bless [ $ab*$c, $ab*$s ];
66 },
67
68 'sqrt' => sub {
69 my($zr,$zi) = @{$_[0]};
70 my ($x, $y, $r, $w);
71 my $c = new Math::Complex (0,0);
72 if (($zr == 0) && ($zi == 0)) {
73 # nothing, $c already set
74 }
75 else {
76 $x = abs($zr);
77 $y = abs($zi);
78 if ($x >= $y) {
79 $r = $y/$x;
80 $w = sqrt($x) * sqrt(0.5*(1.0+sqrt(1.0+$r*$r)));
81 }
82 else {
83 $r = $x/$y;
84 $w = sqrt($y) * sqrt($y) * sqrt(0.5*($r+sqrt(1.0+$r*$r)));
85 }
86 if ( $zr >= 0) {
87 @$c = ($w, $zi/(2 * $w) );
88 }
89 else {
90 $c->[1] = ($zi >= 0) ? $w : -$w;
91 $c->[0] = $zi/(2.0* $c->[1]);
92 }
93 }
94 return $c;
95 },
96
97 qw("" stringify)
98);
99
100sub new {
101 shift;
102 my @C = @_;
103 bless \@C;
104}
105
106sub Re {
107 my($x,$y) = @{$_[0]};
108 $x;
109}
110
111sub Im {
112 my($x,$y) = @{$_[0]};
113 $y;
114}
115
116sub arg {
117 my($x,$y) = @{$_[0]};
118 atan2($y,$x);
119}
120
121sub stringify {
122 my($x,$y) = @{$_[0]};
123 my($re,$im);
124
125 $re = $x if ($x);
126 if ($y == 1) {$im = 'i';}
127 elsif ($y == -1){$im = '-i';}
128 elsif ($y) {$im = "${y}i"; }
129
130 local $_ = $re.'+'.$im;
131 s/\+-/-/;
132 s/^\+//;
133 s/[\+-]$//;
134 $_ = 0 if ($_ eq '');
135 return $_;
136}