From: Nobuyoshi Nakada <nobu@...>
Date: 2009-05-16T06:23:53+09:00
Subject: [ruby-core:23465] Re: [Feature #1408] 0.1.to_r not equal to  (1/10)

Hi,

At Fri, 1 May 2009 21:12:52 +0900,
Roger Pack wrote in [ruby-core:23345]:
> True--however the (current) code for String#to_s attempts to
> determine whether the floating point number "is the
> equivalent default for the rounded value" (i.e. if it round
> trips).

What about this?


Index: rational.c
===================================================================
--- rational.c	(revision 23433)
+++ rational.c	(working copy)
@@ -1286,4 +1286,5 @@ integer_to_r(VALUE self)
 }
 
+#if 0
 static void
 float_decode_internal(VALUE self, VALUE *rf, VALUE *rn)
@@ -1299,5 +1300,4 @@ float_decode_internal(VALUE self, VALUE 
 }
 
-#if 0
 static VALUE
 float_decode(VALUE self)
@@ -1310,11 +1310,82 @@ float_decode(VALUE self)
 #endif
 
+#if FLT_RADIX == 2 && SIZEOF_BDIGITS * 2 * CHAR_BIT > DBL_MANT_DIG
+# ifdef HAVE_LONG_LONG
+#   define BDIGITDBL2NUM(x) ULL2NUM(x)
+# else
+#   define BDIGITDBL2NUM(x) ULONG2NUM(x) 
+# endif
+#else
+# define NEEDS_FDIV
+static ID id_fdiv;
+fun2(fdiv)
+#endif
+
+static VALUE
+float_r_round(double a, double f, int n)
+{
+    int i, r;
+#ifdef BDIGITDBL2NUM
+    BDIGIT_DBL fn = (BDIGIT_DBL)fabs(f);
+    BDIGIT_DBL d1 = (BDIGIT_DBL)1 << -n, d2 = d1;
+    BDIGIT_DBL rv = d1 % fn;
+    VALUE b, d;
+    if (rv < 10) {
+	for (i = 1, r = (int)rv; i <= r; ++i) {
+	    if ((double)fn / --d2 != a) break;
+	    if (fn % (d1 = d2) == 0) break;
+	}
+    }
+    else if ((rv = fn - rv) && rv < 10) {
+	for (i = 1, r = (int)rv; i <= r; ++i) {
+	    if ((double)fn / ++d2 != a) break;
+	    if (fn % (d1 = d2) == 0) break;
+	}
+    }
+    b = BDIGITDBL2NUM(fn);
+    d = BDIGITDBL2NUM(d1);
+    if (f < 0) b = f_negate(b);
+#else
+    VALUE d2, fn, rv;
+    VALUE b = rb_dbl2big(f);
+    VALUE d = rb_big_pow(rb_uint2big(FLT_RADIX), INT2FIX(-n));
+    if (FIXNUM_P(d)) {
+	d = rb_uint2big(FIX2LONG(d));
+    }
+    d2 = d;
+    fn = f_abs(b);
+    rv = rb_big_modulo(d, fn);
+    if (FIXNUM_P(rv) && (r = FIX2LONG(rv)) < 10) {
+	for (i = 1; i <= r; ++i) {
+	    d2 = f_sub(d2, INT2FIX(1));
+	    if (RFLOAT_VALUE(f_fdiv(fn, d2)) != a) break;
+	    if (f_mod(fn, d = d2) == INT2FIX(0)) break;
+	}
+    }
+    else if (FIXNUM_P(rv = f_sub(fn, rv)) && (r = FIX2LONG(rv)) < 10) {
+	for (i = 1; i <= r; ++i) {
+	    d2 = f_add(d2, INT2FIX(1));
+	    if (RFLOAT_VALUE(f_fdiv(fn, d2)) != a) break;
+	    if (f_mod(fn, d = d2) == INT2FIX(0)) break;
+	}
+    }
+#endif
+    return rb_rational_new(b, d);
+}
+
 static VALUE
 float_to_r(VALUE self)
 {
-    VALUE f, n;
+    double a, f;
+    int n;
 
-    float_decode_internal(self, &f, &n);
-    return f_mul(f, f_expt(INT2FIX(FLT_RADIX), n));
+    a = RFLOAT_VALUE(self);
+    f = frexp(a, &n);
+    f = ldexp(f, DBL_MANT_DIG);
+    n -= DBL_MANT_DIG;
+    if (n <= DBL_MANT_DIG && f != 0) {
+	return float_r_round(a, f, n);
+    }
+    return f_mul(rb_dbl2big(f), f_expt(INT2FIX(FLT_RADIX), INT2FIX(n)));
 }
 
@@ -1569,4 +1640,7 @@ Init_Rational(void)
     id_to_s = rb_intern("to_s");
     id_truncate = rb_intern("truncate");
+#ifdef NEEDS_FDIV
+    id_fdiv = rb_intern("fdiv");
+#endif
 
     ml = (long)(log(DBL_MAX) / log(2.0) - 1);


-- 
Nobu Nakada