| CC: | jkeenan [...] cpan.org |
| Subject: | [RT 116767] Reduce risk of overflow/underflow in _update_polar(). |
The following issue was reported in Perl 5's bug queue at rt.perl.org.
Since Math::Complex is listed (in Porting/Maintainers.pl) as being
maintained primarily on CPAN, I am forwarding this issue to this queue.
--jkeenan
From Peter John Acklam:
Compute $rho using a numerically better formula that does not
overflow or underflow during intermediate computations.
---
cpan/Math-Complex/lib/Math/Complex.pm | 24 ++++++++++++++++++++++--
1 files changed, 22 insertions(+), 2 deletions(-)
diff --git a/cpan/Math-Complex/lib/Math/Complex.pm
b/cpan/Math-Complex/lib/Math/Complex.pm
index 19fb164..1c498f3 100644
--- a/cpan/Math-Complex/lib/Math/Complex.pm
+++ b/cpan/Math-Complex/lib/Math/Complex.pm
@@ -420,8 +420,28 @@ sub _update_polar {
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)];
+
+ # To reduce risk of overflow and underflow, use
+ #
+ # / abs(x) * sqrt(1 + (y/x)**2) if x > y
+ # sqrt(x**2 + y**2) = | abs(y) * sqrt(1 + (x/y)**2) if x < y
+ # \ abs(x) * sqrt(2) if x = y
+
+ my $ax = CORE::abs($x);
+ my $ay = CORE::abs($y);
+ my $rho;
+ if ($ax > $ay) {
+ my $c = $ay/$ax;
+ $rho = $ax * CORE::sqrt(1 + $c * $c);
+ } elsif ($ax < $ay) {
+ my $c = $ax/$ay;
+ $rho = $ay * CORE::sqrt(1 + $c * $c);
+ } else {
+ $rho = $ax * CORE::sqrt(2);
+ }
+ my $theta = CORE::atan2($y, $x);
+
+ return $self->{'polar'} = [$rho, $theta];
}
#
--
1.7.9