Perl rounding error again

Question

The example to show the problem:

• having a number 105;
• divide with 1000 (result 0.105)
• rouded to 2 decimal places should be: 0.11

Now, several scripts - based on answers to another questions:

This is mostly recommented and mostly upvoted solution is using printf.

use 5.014;
use warnings;
my $i = 105;
printf "%.2f\n", $i/1000;   #prints 0.10

but prints a wrong result. In the comment to https://stackoverflow.com/a/1838885 @Sinan Unur says (6 times upvoted comment):

Use sprintf("%.3f", $value) for mathematical purposes too.

but, it didn't works "sometimes"... like above.

The another recommented solution Math::BigFloat:

use 5.014;
use warnings;
use Math::BigFloat;
my $i = 105;

Math::BigFloat->precision(-2);
my $r = Math::BigFloat->new($i/1000);

say "$r";   #0.10 ;(

Wrong result too. Another recommened one bignum:

use 5.014;
use warnings;
use bignum ( p => -2 );
my $i = 105;
my $r = $i/1000;
say "$r";   #0.10 ;(

wrong again. ;(

Now the working ones:

use 5.014;
use warnings;
use Math::Round;
my $i = 105;
say nearest(0.01, $i/1000); #GREAT prints 0.11 :)

good result 0.11, however a comment here https://stackoverflow.com/a/571740 complains about it.

and finally another recommendation "by my own" function:

use 5.014;
use warnings;
my $i = 105;
my $f = $i/1000;

say myround($f,2);  # 0.11

sub myround {
    my($float, $prec) = @_;
    my $f = $float * (10**$prec);
    my $r = int($f + $f/abs($f*2));
    return $r/(10**$prec);
}

prints 0.11 too, but can't prove it's correctness.

For the reference I was read:

• How do you round a floating point number in Perl?
• In Perl, how can I limit the number of places after the decimal point but have no trailing zeroes?

• How do I set the floating point precision in Perl?

• and many-many others...

• and finally this too: http://www.exploringbinary.com/inconsistent-rounding-of-printed-floating-point-numbers/ what gives a an really good overall view to the problem.

I understand than it is common problem to all languages, but please, after all above reading - I still have this question:

What is the error-proof way in perl to round a floating point number to N decimal places - with mathematically correct way, e.g. what will round results like 105/1000 correctly to N decimal places without "surprises"...

Answer 1

You're expecting a specific behaviour when the number is exactly 0.105, but floating point errors mean you can't expect a number to be exactly what you think it is.

105/1000 is a periodic number in binary just like 1/3 is periodic in decimal.

105/1000
       ____________________
= 0.00011010111000010100011 (bin)

~ 0.00011010111000010100011110101110000101000111101011100001 (bin)

= 0.10499999999999999611421941381195210851728916168212890625

0.1049999... is less than 0.105, so it rounds to 0.10.

But even if you had 0.105 exactly, that would still round to 0.10 since sprintf rounds half to even. A better test is 155/1000

155/1000
       ____________________
= 0.00100111101011100001010 (bin)

~ 0.0010011110101110000101000111101011100001010001111010111 (bin)

= 0.1549999999999999988897769753748434595763683319091796875

0.155 should round to 0.16, but it rounds to 0.15 due to floating point error.

$ perl -E'$_ = 155; say sprintf("%.2f", $_/1000);'
0.15

$ perl -E'$_ = 155; say sprintf("%.0f", $_/10)/100;'
0.16

The second one works because 5/10 isn't periodic, and therein lies the solution. As Sinan Unur said, you can correct the error by using sprintf. But you have to round to an integer if you don't want to lose your work.

$ perl -E'
   $_ = 155/1000;

   $_ *= 1000;                # Move decimal point past significant.
   $_ = sprintf("%.0f", $_);  # Fix floating-point error.
   $_ /= 10;                  # 5/10 is not periodic
   $_ = sprintf("%.0f", $_);  # Do our rounding.
   $_ /= 100;                 # Restore decimal point.

   say;
'
0.16

That will fix the rounding error, allowing sprintf to properly round half to even.

0.105  =>  0.10
0.115  =>  0.12
0.125  =>  0.12
0.135  =>  0.14
0.145  =>  0.14
0.155  =>  0.16
0.165  =>  0.16

If you want to round half up instead, you'll need to using something other than sprintf to do the final rounding. Or you could add s/5\z/6/; before the division by 10.

---

But that's complicated.

The first sentence of the answer is key. You're expecting a specific behaviour when the number is exactly 0.105, but floating point errors mean you can't expect a number to be exactly what you think it is. The solution is to introduce a tolerance. That's what rounding using sprintf does, but it's a blunt tool.

use strict;
use warnings;
use feature qw( say );

use POSIX qw( ceil floor );

sub round_half_up {
   my ($num, $places, $tol) = @_;

   my $mul = 1; $mul *= 10 for 1..$places;
   my $sign = $num >= 0 ? +1 : -1;

   my $scaled = $num * $sign * $mul;
   my $frac = $scaled - int($scaled);

   if ($sign >= 0) {
      if ($frac < 0.5-$tol) {
         return floor($scaled) / $mul;
      } else {
         return ceil($scaled) / $mul;
      }
   } else {
      if ($frac < 0.5+$tol) {
         return -floor($scaled) / $mul;
      } else {
         return -ceil($scaled) / $mul;
      }
   }
}

say sprintf '%5.2f', round_half_up( 0.10510000, 2, 0.00001);  #  0.11
say sprintf '%5.2f', round_half_up( 0.10500001, 2, 0.00001);  #  0.11  Within tol
say sprintf '%5.2f', round_half_up( 0.10500000, 2, 0.00001);  #  0.11  Within tol
say sprintf '%5.2f', round_half_up( 0.10499999, 2, 0.00001);  #  0.11  Within tol
say sprintf '%5.2f', round_half_up( 0.10410000, 2, 0.00001);  #  0.10
say sprintf '%5.2f', round_half_up(-0.10410000, 2, 0.00001);  # -0.10
say sprintf '%5.2f', round_half_up(-0.10499999, 2, 0.00001);  # -0.10  Within tol
say sprintf '%5.2f', round_half_up(-0.10500000, 2, 0.00001);  # -0.10  Within tol
say sprintf '%5.2f', round_half_up(-0.10500001, 2, 0.00001);  # -0.10  Within tol
say sprintf '%5.2f', round_half_up(-0.10510000, 2, 0.00001);  # -0.11

There's probably existing solutions that work along the same lines.

Answer 2

In the old Integer math days of programming, we use to pretend to use decimal places:

N = 345
DISPLAY N        # Displays 345
DISPLAY (1.2) N  # Displays 3.45

We learned a valuable trick when attempting to round sales taxes correctly:

my $amount = 1.344;
my $amount_rounded = sprintf "%.2f", $amount + .005;
my $amount2 = 1.345;
my $amount_rounded2 = sprintf "%.2f", $amount2 + .005;
say "$amount_rounted   $amount_rounded2";  # prints 1.34 and 1.35

By adding in 1/2 of the precision, I display the rounding correctly. When the number is 1.344, adding .005 made it 1.349, and chopping off the last digit displays dip lays 1.344. When I do the same thing with 1.345, adding in .005 makes it 1.350 and removing the last digit displays it as 1.35.

You could do this with a subroutine that will return the rounded amount.

---

Interesting...

There is a PerlFAQ on this subject. It recommends simply using printf to get the correct results:

use strict;
use warnings;
use feature qw(say);

my $number = .105;
say "$number";
printf "%.2f\n", $number;   # Prints .10 which is incorrect
printf "%.2f\n", 3.1459;    # Prins 3.15 which is correct

For Pi, this works, but not for .105. However:

use strict;
use warnings;
use feature qw(say);

my $number = .1051;
say "$number";
printf "%.2f\n", $number;   # Prints .11 which is correct
printf "%.2f\n", 3.1459;    # Prints 3.15 which is correct

This looks like an issue with the way Perl stores .105 internally. Probably something like .10499999999 which would be correctly rounded downwards. I also noticed that Perl warns me about using round and rounding as possible future reserved words.

Answer 3

Your custom function should mostly work as expected. Here's how it works and how you can verify it's correct:

sub myround {
    my($float, $prec) = @_;

    # Prevent div-by-zero later on
    if ($float == 0) { return 0; }

    # Moves the decimal $prec places to the right
    # Example: $float = 1.234, $prec = 2
    #   $f = $float * 10^2;
    #   $f = $float * 100;
    #   $f = 123.4;
    my $f = $float * (10**$prec);

    # Round 0.5 away from zero using $f/abs($f*2)
    #   if $f is positive, "$f/abs($f*2)" becomes  0.5
    #   if $f is negative, "$f/abs($f*2)" becomes -0.5
    #   if $f is zero, we have a problem (hence the earlier if statement)
    # In our example:
    #   $f = 123.4 + (123.4 / (123.4 * 2));
    #   $f = 123.4 + (0.5);
    #   $f = 123.9;
    # Then we truncate to integer:
    #   $r = int(123.9);
    #   $f = 123;
    my $r = int($f + $f/abs($f*2));

    # Lastly, we shift the deciaml back to where it should be:
    #   $r / 10^2
    #   $r / 100
    #   123 / 100
    #   return 1.23;
    return $r/(10**$prec);
}

However, the following it will throw an error for $float = 0, so there's an additional if statement at the beginning.

The nice thing about the above function is that it's possible to round to negative decimal places, allowing you round to the left of the decimal. For example, myround(123, -2) will give 100.

Answer 4

I'd add use bignum to your original code example.

use 5.014;
use warnings;
use bignum;

my $i = 105;                # Parsed as Math::BigInt
my $r = $i / 1000;          # Overloaded division produces Math::BigFloat
say $r->ffround(-2, +inf);  # Avoid using printf and the resulting downgrade to common float.

This solves the error you made in your use Math::BigFloat example by parsing your numbers into objects imediately and not waiting for you to pass the results of a round off error into Math::BigFloat->new

Answer 5

According to my experience, Perl printf/sprintf uses wrong algorithm. I made this conclusion considering at least the following simple example:

# The same floating part for both numbers (*.8830 or *.8829) is expected in the rounded value, but it is different for some reason:
printf("%.4f\n", "7.88295"); # gives 7.8830
printf("%.4f\n", "8.88295"); # gives 8.8829

The integer part should not have any influence in this example, but it has. I got this result with Perl 5.8.8.