Skip to content

Commit f6fc109

Browse files
committed
change the Round function to combat machine rounding error
1 parent d16b4d0 commit f6fc109

1 file changed

Lines changed: 20 additions & 3 deletions

File tree

macros/core/PGauxiliaryFunctions.pl

Lines changed: 20 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -175,12 +175,29 @@ =head2 Round
175175
176176
Round(1.789,2) returns 1.79
177177
178+
This macro assumes an input x having lots of consecutive 9s such as 1.4999999999999997 is already
179+
the result of machine rounding error, and was meant to be 1.5. And so will round up to 2, not down
180+
to 1. Do not use this macro if a problem might work with honest numbers like 1.4999999999999997
181+
and need to round it.
182+
178183
=cut
179184

180-
# Round contributed bt Mark Schmitt 3-6-03
181185
sub Round {
182-
if (@_ == 1) { $_[0] > 0 ? int $_[0] + 0.5 : int $_[0] - 0.5 }
183-
elsif (@_ == 2) { $_[0] > 0 ? Round($_[0] * 10**$_[1]) / 10**$_[1] : Round($_[0] * 10**$_[1]) / 10**$_[1] }
186+
my ($x, $n) = @_;
187+
$n = 0 unless $n;
188+
my $e = (split(/E/, sprintf("%E", $x)))[1] + 0; # exponent for $x
189+
my $s = ($x < 0 ? -1 : 1); # the sign of $x
190+
my $N = $e + $n; # number of digits to retain
191+
$x += $s * 10**($e - 15); # adjust for repeated 9s
192+
return sprintf("%.${N}E", $x) + 0 unless $N == -1; # round the adjusted value
193+
194+
# For zero original digits, we add a digit just above the first one
195+
# in $x, round that, then remove the added digit, getting 0 if $x
196+
# didn't round up, or 1 in the proper place if it did. This means
197+
# that 0.005 rounds to .01, for example, when 2 digits are requested.
198+
199+
my $d = $s * 10**($e + 1);
200+
return sprintf("%.0E", $x + $d) - $d;
184201
}
185202

186203
=head2 lcm

0 commit comments

Comments
 (0)