Project

General

Profile

Bug #14635

Float#round(n) returns a wrong result when n is big

Added by mame (Yusuke Endoh) 9 months ago. Updated 9 months ago.

Status:
Feedback
Priority:
Normal
Assignee:
-
Target version:
-
[ruby-core:86323]

Description

First of all, don't confuse that this is a usual floating-point error issue.

The following looks inconsistent:

3.0e-31           #=> 3.0e-31
3.0e-31.round(31) #=> 3.0000000000000003e-31

What it should be

A Float value is actually a range.

3.0e-31 represents a range of 0.299999999999999959315060e-30 .. 0.300000000000000003105637e-30 (the bounds are approximate). I call this range "A".
3.0000000000000003e-31 represents a range of 0.300000000000000003105637e-30 .. 0.300000000000000046896214e-30. I call this range "B".

x.round(31) should (1) multiple x with 10**31, (2) round it as an integer, and (3) divide it with 10**31.

In this case:

(1) 3.0e-31 * 10**31 is a range of 2.99999999999999959315060 .. 3.00000000000000003105637.
(2) The rounded result is 3, whichever value is chosen from the range above.
(3) 3.0 / 10**31 is within the range "A", not within the range "B", so the result should be 3.0e-31, not 3.0000000000000003e-31.

How the bug occurs

The reason why 3.0e-31.round(31) returns 3.0000000000000003e-31, is the implementation issue of Float#round. It does the following:

(1) f = pow(10, b)
(2) x = round(x * f) as an integer
(3) return x / f

However, a double variable f cannot represent pow(10, 31) precisely. In other words, the 10**31 must be handled as an integer, but the implementation handles it as an inexact floating-point value. This is the issue.

How to fix

The issue is simple, but it might be very difficult to fix. strtod handles a string "3.0e-31" correctly. So, by doing the same as strtod, this issue would be fixed. However, the strtod implementation looks very difficult, at least to me. Contribution from mathematician is welcome.
(Honestly, I don't want to see such a complication in the source code. Another simpler approach would be more preferable.)

References

This issue has been already reported in #5273 by marcandre. But the status of the ticket looks unclear; I cannot see how many issues remains. So, I created this ticket for just one bug that I could confirm.


Related issues

Related to Ruby trunk - Bug #5273: Float#round returns the wrong floats for higher precisionClosed

History

#1 Updated by mame (Yusuke Endoh) 9 months ago

  • Related to Bug #5273: Float#round returns the wrong floats for higher precision added

#2 Updated by mame (Yusuke Endoh) 9 months ago

  • Subject changed from Float#round sometimes returns a wrong result to Float#round(n) returns a wrong result when n is big

#3 [ruby-core:86331] Updated by mame (Yusuke Endoh) 9 months ago

I've found a much simpler solution: when n is big, it should first translate the float to a rational, then call Rational#round, and finally translate the resulting rational to a float. It is slow, only when n >= 23 for Float#round(n).

Currently, Rational#to_f has the same inaccuracy issue, which can be fixed by #14637. The following patch includes the hunk for #14637.

diff --git a/bignum.c b/bignum.c
index b4c7560034..fd5f385cac 100644
--- a/bignum.c
+++ b/bignum.c
@@ -6178,9 +6178,7 @@ rb_big_fdiv_double(VALUE x, VALUE y)
        return big_fdiv_int(x, rb_int2big(FIX2LONG(y)));
     }
     else if (RB_BIGNUM_TYPE_P(y)) {
-   dy = rb_big2dbl(y);
-   if (isinf(dx) || isinf(dy))
-       return big_fdiv_int(x, y);
+   return big_fdiv_int(x, y);
     }
     else if (RB_FLOAT_TYPE_P(y)) {
    dy = RFLOAT_VALUE(y);
diff --git a/internal.h b/internal.h
index 9b6a213151..0bf20b19b0 100644
--- a/internal.h
+++ b/internal.h
@@ -1689,6 +1689,7 @@ VALUE rb_cstr_to_rat(const char *, int);
 VALUE rb_rational_abs(VALUE self);
 VALUE rb_rational_cmp(VALUE self, VALUE other);
 VALUE rb_numeric_quo(VALUE x, VALUE y);
+VALUE rb_flo_round_by_rational(int argc, VALUE *argv, VALUE num);

 /* re.c */
 VALUE rb_reg_compile(VALUE str, int options, const char *sourcefile, int sourceline);
diff --git a/numeric.c b/numeric.c
index 01856c7f20..15b27e9132 100644
--- a/numeric.c
+++ b/numeric.c
@@ -2239,6 +2239,10 @@ flo_round(int argc, VALUE *argv, VALUE num)
    frexp(number, &binexp);
    if (float_round_overflow(ndigits, binexp)) return num;
    if (float_round_underflow(ndigits, binexp)) return DBL2NUM(0);
+   if (ndigits > DBL_MANT_DIG * log(2.0) / log(5.0)) {
+       /* In this case, pow(10, ndigits) cannot be accurate. */
+       return rb_flo_round_by_rational(argc, argv, num);
+   }
    f = pow(10, ndigits);
    x = ROUND_CALL(mode, round, (number, f));
    return DBL2NUM(x / f);
diff --git a/rational.c b/rational.c
index d88f50f886..01bb88d1ae 100644
--- a/rational.c
+++ b/rational.c
@@ -1533,6 +1533,13 @@ nurat_round_n(int argc, VALUE *argv, VALUE self)
     return f_round_common(argc, argv, self, round_func);
 }

+static VALUE float_to_r(VALUE self);
+VALUE
+rb_flo_round_by_rational(int argc, VALUE *argv, VALUE num)
+{
+    return nurat_to_f(nurat_round_n(argc, argv, float_to_r(num)));
+}
+
 static double
 nurat_to_double(VALUE self)
 {
@@ -2016,7 +2023,6 @@ integer_denominator(VALUE self)
     return INT2FIX(1);
 }

-static VALUE float_to_r(VALUE self);
 /*
  * call-seq:
  *    flo.numerator  ->  integer

Also available in: Atom PDF