help-smalltalk
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Help-smalltalk] BUG gst2.2: Float asExactFraction asFloat inexact


From: Paolo Bonzini
Subject: Re: [Help-smalltalk] BUG gst2.2: Float asExactFraction asFloat inexact
Date: Sat, 27 Jan 2007 16:51:58 +0100
User-agent: Thunderbird 1.5.0.9 (Macintosh/20061207)


Unfortunately, current implementation of both Largeinteger>>asFloat and Fraction>>asFloat will cause successive round off errors...

With 2.3.1 and your fix to asExactFraction, this works. But there are more bugs in general, and I tried to rewrite the methods to use timesTwoPower: as much as possible, and to avoid problems with denormals.

More over, denormalized number (gradual underflow) will convert back to zero...

| x |
x := 1.1d-310.
^x asExactFraction asFloat = 0

This is exposing a bug in the lexer.

st> 1.1d-310 = 0!
true
st> 1.1d-306 / 10000 = 0!
false

However, there is another problem, in that (1 << 1072) is not a representable floating point number, but its reciprocal is a (denormal) floating point number. Therefore, the problem is that for denormals asExactFraction returns a denominator that is not representable, and the bug is actually in Fraction>>#asFloat.

The attached patch fixes all these bugs, except the one in the lexer. I would like a double-check in Fraction>>#asFloat, maybe the rounding could be improved if the numerator or denominator do not fit in a Float.

If you can provide some regression tests for rounding errors in #asExactFraction and #asFloat, please do. Just write the code for the borderline cases you can think of, I'll check if the results are right since they may still be wrong with your outdated image.

Paolo
2007-01-27  Paolo Bonzini  <address@hidden>
            Nicolas Cellier  <address@hidden>

        * kernel/Float.st: Rewrite #truncated and #asExactFraction.
        * kernel/Fraction.st: Use #asFloatD in #hash.  Be careful about
        denormals in #asFloat:
        * kernel/LargeInt.st: Rewrite #asFloat:.
        * kernel/MappedColl.st: Use map size as collection size.

2007-01-27  Paolo Bonzini  <address@hidden>

        * libgst/prims.def: Use trunc{,f,l} and lrint{,f,l} to truncate
        floats.

--- orig/kernel/Float.st
+++ mod/kernel/Float.st
@@ -194,27 +194,30 @@ truncated
         ifFalse: [ float := self negated ].
 
     exponent := float exponent.
-    bytes := ByteArray new: exponent // 8 + 2.
-    float := float timesTwoPower: (exponent bitClear: 7) negated.
+    bytes := LargePositiveInteger new: (self class precision + 7) // 8 + 1.
+    float := float timesTwoPower: float class precision - exponent - 8.
 
-    bytes size - 1 to: 1 by: -1 do: [ :i |
-         bytes at: i put: float truncated.
-         float := float fractionPart timesTwoPower: 8
+    1 to: bytes size - 1 do: [ :i |
+         bytes digitAt: i put: (float fractionPart timesTwoPower: 8) truncated.
+         float := float integerPart timesTwoPower: -8.
     ].
 
-    ^positive
-        ifTrue: [ (LargeInteger from: bytes) ]
-        ifFalse: [ (LargeInteger from: bytes) negated ]
+    bytes := bytes bitShift: (exponent - float class precision).
+    positive ifFalse: [ bytes := bytes negated ].
+    ^bytes
 !
 
 asExactFraction
     "Convert the receiver into a fraction with optimal approximation,
      but with usually huge terms."
 
+    | shift mantissa |
     self checkCoercion.
-
-    ^(self timesTwoPower: self exponent + self class precision) truncated /
-     (1    bitShift:      self exponent + self class precision)
+    shift := self exponent negated + self class precision.
+    mantissa := (self timesTwoPower: shift) truncated.
+    ^shift negative
+        ifTrue: [mantissa * (1 bitShift: shift negated)]
+        ifFalse: [mantissa / (1 bitShift: shift)]
 !
 
 asFraction


--- orig/kernel/Fraction.st
+++ mod/kernel/Fraction.st
@@ -288,7 +288,7 @@ truncated
 hash
     "Answer an hash value for the receiver"
     denominator = 1 ifTrue: [ ^numerator hash ].
-    ^(numerator asFloatD / denominator asFloatD) hash
+    ^self asFloatD hash
 ! !
 
 
@@ -385,7 +385,7 @@ storeOn: aStream
 
 asFloat: characterization
     "Answer the receiver converted to a Float"
-    | n d shift sign |
+    | n d nShift dShift sign x |
     n := numerator.
     d := denominator.
     sign := n sign * d sign.
@@ -394,11 +394,15 @@ asFloat: characterization
      1e1800 asFloat / (1e1799 + 1) asFloat = NaN, but
      (1e1800 / (1e1799 + 1)) asFloat must be 10."
 
-    shift := (characterization emax - numerator highBit).
-    shift := shift min: (characterization emax - denominator highBit).
-    shift < 0 ifTrue: [
-       "Lose some more precision, but we MUST avoid infinites and NaNs!"
-       shift := shift - 10. n := n bitShift: shift. d := d bitShift: shift ].
+    nShift := (characterization emax - numerator highBit).
+    nShift < 0
+       ifTrue: [ n := n bitShift: nShift ]
+       ifFalse: [ nShift := 0 ].
+
+    dShift := (characterization emax - denominator highBit).
+    dShift < 0
+       ifTrue: [ d := d bitShift: dShift ]
+       ifFalse: [ dShift := 0 ].
 
     d = 0 ifTrue: [
        ^sign > 0
@@ -409,8 +413,8 @@ asFloat: characterization
            ifTrue: [ characterization coerce: 0 ]
            ifFalse: [ characterization negativeInfinity reciprocal ] ].
 
-    ^(characterization coerce: n)
-        / (characterization coerce: d)
+    x := (characterization coerce: n) / (characterization coerce: d).
+    ^x timesTwoPower: dShift - nShift
 !
 
 reduce


--- orig/kernel/LargeInt.st
+++ mod/kernel/LargeInt.st
@@ -977,14 +988,16 @@ sign
 asFloat: characterization
     "Answer the receiver converted to a Float"
 
-    | adjust result byte |
-    byte := characterization coerce: 256.
-    adjust := byte raisedToInteger: self size - 1.
+    | result byte shift |
+    shift := 8 * self size.
     result := 0.
 
     self size to: 1 by: -1 do: [ :index |
-       result := (self at: index) * adjust + result.
-       adjust := adjust / byte.
+       byte := self at: index.
+       byte = 0 ifFalse: [
+            byte := characterization coerce: byte.
+           result := (byte timesTwoPower: shift) + result ].
+       shift := shift - 8
     ].
     ^result
 !


--- orig/kernel/MappedColl.st
+++ mod/kernel/MappedColl.st
@@ -89,7 +89,7 @@ at: key put: value
 
 size
     "Answer the receiver's size"
-    ^domain size
+    ^map size
 !
 
 add: anObject


--- orig/libgst/prims.def
+++ mod/libgst/prims.def
@@ -1531,14 +1531,9 @@ primitive VMpr_FloatD_truncated [succeed
   if COMMON (RECEIVER_IS_CLASS (oop1, _gst_floatd_class))
     {
       double oopValue = FLOATD_OOP_VALUE (oop1);
-      if COMMON (oopValue >= 0.0 && oopValue <= MAX_ST_INT)
+      if COMMON (oopValue >= MIN_ST_INT && oopValue <= MAX_ST_INT)
        {
-         PUSH_INT ((intptr_t) (oopValue + 0.000000000000005));
-         PRIM_SUCCEEDED;
-       }
-      else if COMMON (oopValue < 0.0 && oopValue >= MIN_ST_INT)
-       {
-         PUSH_INT ((intptr_t) (oopValue - 0.000000000000005));
+         PUSH_INT (lrint (trunc (oopValue)));
          PRIM_SUCCEEDED;
        }
     }
@@ -1732,14 +1727,9 @@ primitive VMpr_FloatE_truncated [succeed
   if COMMON (RECEIVER_IS_CLASS (oop1, _gst_floate_class))
     {
       double oopValue = FLOATE_OOP_VALUE (oop1);
-      if COMMON (oopValue >= 0.0 && oopValue <= MAX_ST_INT)
+      if COMMON (oopValue >= MIN_ST_INT && oopValue <= MAX_ST_INT)
        {
-         PUSH_INT ((intptr_t) (oopValue + 0.000000000000005));
-         PRIM_SUCCEEDED;
-       }
-      else if COMMON (oopValue < 0.0 && oopValue >= MIN_ST_INT)
-       {
-         PUSH_INT ((intptr_t) (oopValue - 0.000000000000005));
+         PUSH_INT (lrintf (truncf (oopValue)));
          PRIM_SUCCEEDED;
        }
     }
@@ -1933,14 +1923,9 @@ primitive VMpr_FloatQ_truncated [succeed
   if COMMON (RECEIVER_IS_CLASS (oop1, _gst_floatq_class))
     {
       long double oopValue = FLOATQ_OOP_VALUE (oop1);
-      if COMMON (oopValue >= 0.0 && oopValue <= MAX_ST_INT)
-       {
-         PUSH_INT ((intptr_t) (oopValue + 0.000000000000005));
-         PRIM_SUCCEEDED;
-       }
-      else if COMMON (oopValue < 0.0 && oopValue >= MIN_ST_INT)
+      if COMMON (oopValue >= MIN_ST_INT && oopValue <= MAX_ST_INT)
        {
-         PUSH_INT ((intptr_t) (oopValue - 0.000000000000005));
+         PUSH_INT (lrintl (truncl (oopValue)));
          PRIM_SUCCEEDED;
        }
     }




reply via email to

[Prev in Thread] Current Thread [Next in Thread]