Re: Roundoff problem with Float and Marshal

Hi,

At Fri, 18 Apr 2003 12:20:01 +0900,
Yukihiro Matsumoto wrote:
> It's the limitation of base 10 representation of float values that
> Marshal uses.  We need to use binary representation like XDR to solve
> this.  I'm not sure if it is worth for the cost of its implementation
> and portability.  Float tend to be inexact, after all.

I'd tried C99 style hexadecimal float.  I expect this is
portable (except for depending on frexp() and ldexp()), but
could not be loaded by current version.

But still, strtod() problem remains.


Index: marshal.c
===================================================================
RCS file: /pub/cvs/ruby/src/ruby/marshal.c,v
retrieving revision 1.84
diff -u -2 -p -r1.84 marshal.c
--- marshal.c	9 Apr 2003 05:08:25 -0000	1.84
+++ marshal.c	17 Apr 2003 14:01:58 -0000
@@ -44,5 +44,5 @@ shortlen(len, ds)
 
 #define MARSHAL_MAJOR   4
-#define MARSHAL_MINOR   8
+#define MARSHAL_MINOR   9
 
 #define TYPE_NIL	'0'
@@ -201,6 +201,39 @@ w_float(d, arg)
     }
     else {
-	/* xxx: should not use system's sprintf(3) */
-	sprintf(buf, "%.16g", d);
+	char *p = buf;
+	double i;
+	int e;
+
+	if (d < 0) {
+	    d = -d;
+	    *p++ = '-';
+	}
+	strcpy(p, "0x");
+	p += 2;
+	if (d > 0x10000 || d < 1.0) {
+	    d = frexp(d, &e);
+	    *p++ = '0';
+	}
+	else {
+	    d = modf(d, &i);
+	    sprintf(p, "%x", (int)i);
+	    p += strlen(p);
+	    e = 0;
+	}
+	if (d > 0) {
+	    *p++ = '.';
+	    do {
+		d = modf(d * 0x10000, &i);
+		sprintf(p, "%.4x", (int)i);
+		p += 4;
+	    } while (d > 0);
+	}
+	while (p[-1] == '0') --p;
+	if (e) {
+	    sprintf(p, "p%+d", e);
+	    p += strlen(p);
+	}
+	w_bytes(buf, p - buf, arg);
+	return;
     }
     w_bytes(buf, strlen(buf), arg);
@@ -778,4 +811,52 @@ r_symbol(arg)
 }
 
+static double
+r_float(str)
+    VALUE str;
+{
+    int s = 0, n;
+    unsigned long fr;
+    double d, t = 0.0;
+    char tmp[SIZEOF_LONG*CHAR_BIT/4+1], *e;
+    const char *ptr = RSTRING(str)->ptr;
+    const char *p = ptr;
+
+    if (*ptr == '-') {
+	s = -1;
+	p++;
+    }
+    if (strcmp(ptr, "nan") == 0) {
+	return t / t;
+    }
+    else if (strcmp(p, "inf") == 0) {
+	return (s ? -1.0 : 1.0) / t;
+    }
+    else if (strncmp(p, "0x", 2) != 0) {
+	return strtod(ptr, 0);
+    }
+    p = strchr(ptr = p, 'p');
+    n = p ? atoi(p + 1) : 0;
+    e = strchr(ptr, '.');
+    n += 4 * (e ? e - ptr : p ? p - ptr : strlen(ptr));
+    tmp[sizeof tmp - 1] = 0;
+    d = 0;
+    while (*ptr) {
+	strncpy(tmp, ptr, sizeof tmp - 1);
+	fr = strtoul(tmp, &e, 16);
+	d += ldexp((double)fr, n -= 4 * (e - tmp));
+	ptr += e - tmp;
+	if (*ptr == '.') {
+	    ptr++;
+	}
+	else if (*ptr == 'p') {
+	    break;
+	}
+    }
+    if (s) {
+	d = -d;
+    }
+    return d;
+}
+
 static char*
 r_unique(arg)
@@ -916,23 +997,6 @@ r_object0(arg, proc)
 
       case TYPE_FLOAT:
-	{
-	    double d, t = 0.0;
-	    VALUE str = r_bytes(arg);
-
-	    if (strcmp(RSTRING(str)->ptr, "nan") == 0) {
-		d = t / t;
-	    }
-	    else if (strcmp(RSTRING(str)->ptr, "inf") == 0) {
-		d = 1.0 / t;
-	    }
-	    else if (strcmp(RSTRING(str)->ptr, "-inf") == 0) {
-		d = -1.0 / t;
-	    }
-	    else {
-		d = strtod(RSTRING(str)->ptr, 0);
-	    }
-	    v = rb_float_new(d);
-	    r_regist(v, arg);
-	}
+	v = rb_float_new(r_float(r_bytes(arg)));
+	r_regist(v, arg);
 	break;
 


-- 
Nobu Nakada