From: cilibrar@... (Rudi Cilibrasi)
Date: 2003-04-18T17:52:11+09:00
Subject: Re: Roundoff problem with Float and Marshal

I have implemented the idea I gave earlier.  I had to use 2 bytes to
pass all my test cases.  There is a 2-bit space for other floating formats
also.  I hope this technique preserves portability as well as before,
yet also gives numerical stability which is essential for some types of
scientific computing, e.g. distributed and long-term simulations.
What do you think of this solution?

Rudi

--- marshal.c	2003-04-18 10:25:29.000000000 +0200
+++ marshal.c	2003-04-18 10:25:02.000000000 +0200
@@ -182,4 +182,59 @@ w_long(x, arg)
 }
 
+#define FT_OTHER 0
+#define FT_IEEE754 1
+#define FT_UNDEFINED 5
+
+static int checkFloatingType()
+{
+  static int floatingType = FT_UNDEFINED;
+  if (floatingType == FT_UNDEFINED) {
+    volatile double t = 1.234;
+    volatile unsigned long *b = (unsigned long *) &t;
+    if (b[0] == 0xc8b43958 && b[1] == 0x3ff3be76)
+      floatingType = FT_IEEE754;
+    /* ... add other checks here ... */
+    if (floatingType == FT_UNDEFINED)
+      floatingType = FT_OTHER;
+  }
+  return floatingType;
+}
+
+static double fixMantissaBits(double d, unsigned char *buf)
+{
+  int floatingType = checkFloatingType();
+  int encodedType = buf[1] >> 6;
+  if (floatingType != encodedType)
+    return d;
+  switch(floatingType) {
+    case FT_IEEE754:
+      ((unsigned char *) &d)[0] = buf[0];
+      ((unsigned char *) &d)[1] &= 0xc0;
+      ((unsigned char *) &d)[1] |= buf[1] & ~0xc0;
+      break;
+    default:
+      break;
+  }
+  return d;
+}
+
+/* Saves 2 bytes total of mantissa and floatingType */
+void saveMantissaBits(double d, unsigned char *buf)
+{
+  unsigned char result;
+  int floatingType = checkFloatingType();
+  switch(floatingType) {
+    case FT_IEEE754:
+      buf[0] = ((unsigned char *) &d)[0];
+      buf[1] = ((unsigned char *) &d)[1];
+      break;
+    default:
+      buf[0] = 1;
+      buf[1] = 1;
+      break;
+  }
+  buf[1] = (buf[1] & ~0xc0) | (floatingType << 6);
+}
+
 static void
 w_float(d, arg)
@@ -202,7 +257,8 @@ w_float(d, arg)
     else {
 	/* xxx: should not use system's sprintf(3) */
-	sprintf(buf, "%.16g", d);
+  saveMantissaBits(d, buf);
+	sprintf(buf+2, "%.17g", d);
     }
-    w_bytes(buf, strlen(buf), arg);
+    w_bytes(buf, strlen(buf+2)+2, arg);
 }
 
@@ -930,5 +986,6 @@ r_object0(arg, proc)
 	    }
 	    else {
-		d = strtod(RSTRING(str)->ptr, 0);
+		d = strtod(RSTRING(str)->ptr+2, 0);
+    d = fixMantissaBits(d,RSTRING(str)->ptr);
 	    }
 	    v = rb_float_new(d);