Re: [Help-smalltalk] Re: [PATCH] Fraction and Integer asFloat: answer nearest floating point

[Paolo Bonzini]

Thanks, I applied a pretty big patch with all your changes and mine.

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.  Rewrite #asFloat:.
        * kernel/LargeInt.st: Rewrite #asFloat:.  Add #trailingZeros.
        * kernel/SmallInt.st: Add #trailingZeros.
        * kernel/Integer.st: Add #trailingZeros.
        * kernel/MappedColl.st: Use map size as collection size.

        * tests/floatmath.st: Add accuracy tests.
        * tests/floatmath.ok: Adjust.
 
2007-01-27  Paolo Bonzini  <address@hidden>

        * lib-src/truncl.c: New.

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

        * libgst/prims.def: Use truncl and lrint to implement
        conversion from float to integer.


--- orig/NEWS
+++ mod/NEWS
@@ -27,6 +27,9 @@ o   #copyFrom:to: is uniformly 0-based f
 o   Fixed a garbage collection bug that typically occurred when installing
     GNU Smalltalk, or when launching the installed image.
 
+o   Fixed many floating point rounding bugs in LargeIntegers and Fractions,
+    thanks to Nicolas Cellier.
+
 o   gst-package honors the INSTALL command found by configure.
 
 o   gst-config does not "forget" to prefix the library directories with -L.


--- orig/configure.ac
+++ mod/configure.ac
@@ -204,7 +204,7 @@ AC_CHECK_LIB(m, atan)
 GST_REPLACE_POLL
 AC_REPLACE_FUNCS(putenv strdup strerror strsignal mkstemp getpagesize \
        getdtablesize strstr ftruncate floorl ceill sqrtl frexpl ldexpl asinl \
-       acosl atanl logl expl tanl sinl cosl strsep strpbrk)
+       acosl atanl logl expl tanl sinl cosl truncl strsep strpbrk)
 AC_CHECK_FUNCS_ONCE(gethostname memcpy memmove sighold uname sbrk usleep lstat 
\
        grantpt popen getrusage gettimeofday getcwd fork strchr \
        sigsetmask alarm select mprotect madvise nl_langinfo waitpid \


--- 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 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,33 +385,59 @@ storeOn: aStream
 
 asFloat: characterization
     "Answer the receiver converted to a Float"
-    | n d shift sign |
-    n := numerator.
-    d := denominator.
-    sign := n sign * d sign.
-
-    "Avoid answering NaNs and infinite values.
-     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 ].
-
-    d = 0 ifTrue: [
-       ^sign > 0
-           ifTrue: [ characterization infinity ]
-           ifFalse: [ characterization negativeInfinity ] ].
-    n = 0 ifTrue: [
-       ^sign > 0
-           ifTrue: [ characterization coerce: 0 ]
-           ifFalse: [ characterization negativeInfinity reciprocal ] ].
+    "Answer the receiver converted to a Float"
 
-    ^(characterization coerce: n)
-        / (characterization coerce: d)
-!
+    | n d sign hn hd hq nBits q q1 r exponent floatExponent | 
+    sign := numerator sign * denominator sign.
+    n := numerator abs.
+    d := denominator abs.
+    hn := n highBit.
+    hd := d highBit.
+
+    "If both numerator and denominator are represented exactly in floating
+     point number, then fastest thing to do is to use hardwired float division"
+    nBits := characterization precision + 1.
+    (hn < nBits and: [hd < nBits]) 
+        ifTrue: [^(characterization coerce: numerator) / (characterization 
coerce: denominator)].
+
+    "Try and obtain a mantissa with characterization precision + 1 bits by 
integer division.
+     Additional bit is a helper for rounding mode.
+     First guess is rough, we might get one more bit or one less"
+    exponent := hn - hd - nBits.
+    exponent > 0 
+        ifTrue: [d := d bitShift: exponent]
+        ifFalse: [n := n bitShift: exponent negated].
+    q := n quo: d.
+    r := n - (q * d).
+    hq := q highBit.
+
+    "check for gradual underflow, in which case we should use less bits"
+    floatExponent := exponent + hq.
+    floatExponent >= (characterization emin - 1) ifFalse: [nBits := nBits + 
floatExponent - characterization emin+1].
+
+    "Use exactly nBits"
+    hq > nBits 
+        ifTrue: 
+            [exponent := exponent + hq - nBits.
+            r := (q bitAnd: (1 bitShift: hq - nBits) - 1) * d + r.
+            q := q bitShift: nBits - hq].
+    hq < nBits 
+        ifTrue: 
+            [exponent := exponent + hq - nBits.
+            q1 := (r bitShift: nBits - hq) quo: d.
+            q := (q bitShift: nBits - hq) bitAnd: q1.
+            r := (r bitShift: nBits - hq) - (q1 * d)].
+
+    "check if we should round upward.
+    The case of exact half (q bitAnd: 1) = 1 & (r = 0)
+    will be handled by Integer>>asFloat:"
+    ((q bitAnd: 1) = 0 or: [r = 0]) ifFalse: [q := q + 1].
+
+    "build the Float"
+    ^(sign > 0
+        ifTrue: [characterization coerce: q]
+        ifFalse: [(characterization coerce: q) negated]) 
+            timesTwoPower: exponent!
 
 reduce
     "Reduce the fraction."


--- orig/kernel/Integer.st
+++ mod/kernel/Integer.st
@@ -139,12 +139,18 @@ clearBit: index
 !
 
 noMask: anInteger
-    "True if no 1 bits in anInteger are 1 in the receiver"
+    "Answer true if no 1 bits in anInteger are 1 in the receiver."
     ^(self bitAnd: anInteger) = 0
 !
 
+trailingZeros
+    "Return the index of the lowest order 1 bit of the receiver,
+     minus 1."
+    self subclassResponsibility
+!
+
 highBit
-    "Return the index of the highest order 1 bit of the receiver"
+    "Return the index of the highest order 1 bit of the receiver."
     self subclassResponsibility
 !
 


--- orig/kernel/LargeInt.st
+++ mod/kernel/LargeInt.st
@@ -493,6 +493,17 @@ negated
 
 !LargeInteger methodsFor: 'bit operations'!
 
+trailingZeros
+    "Answer the number of trailing zero bits in the receiver"
+
+    | each |
+    1 to: self size do: [ :index |
+        (each := self digitAt: index) = 0
+            ifFalse: [ ^index * 8 - 8 + (TrailingZeros at: each) ].
+    ].
+    ^self highBit - 1
+!
+
 bitAnd: aNumber
     "Answer the receiver ANDed with aNumber"
     | newBytes |
@@ -977,17 +988,32 @@ sign
 asFloat: characterization
     "Answer the receiver converted to a Float"
 
-    | adjust result byte |
+    | nTruncatedBits mantissa exponent mask trailingBits inexact carry | 
+
+    "Check for number bigger than maximum mantissa"
+    nTruncatedBits := self highBit - characterization precision.
+    nTruncatedBits <= 0 
+        ifTrue: [^self fastAsFloat: characterization ].
+
+    mantissa := self bitShift: nTruncatedBits negated.
+    exponent := nTruncatedBits.
+
+    "Apply IEEE 754 round to nearest even default rounding mode"
+    carry := self bitAt: nTruncatedBits.
+    (carry = 1 and: [mantissa odd or: [self trailingZeros >= nTruncatedBits]]) 
+        ifTrue: [mantissa := mantissa + 1].
+    ^(characterization coerce: mantissa) timesTwoPower: exponent!
+
+fastAsFloat: characterization
+    "Conversion can be exact, construct Float by successive mul add operations"
+    | result byte |
     byte := characterization coerce: 256.
-    adjust := byte raisedToInteger: self size - 1.
-    result := 0.
+    result := characterization coerce: 0.
 
     self size to: 1 by: -1 do: [ :index |
-       result := (self at: index) * adjust + result.
-       adjust := adjust / byte.
+       result := result * byte + (self at: index).
     ].
-    ^result
-!
+    ^result!
 
 mostSignificantByte
     "Private - Answer the value of the most significant byte"


--- 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/kernel/SmallInt.st
+++ mod/kernel/SmallInt.st
@@ -94,6 +94,24 @@ generality
 
 !SmallInteger methodsFor: 'bit arithmetic'!
 
+trailingZeros
+    "Return the index of the trailing zero bits in the receiver"
+
+    | n bit |
+    self = 0 ifTrue: [ ^0 ].
+
+    n := self.
+    bit := 0.
+
+    (n bitAnd: 16r3FFFFFFF) = 0 ifTrue: [
+       bit := bit + 30. n := n bitShift: -30 ].
+
+    (n bitAnd: 16rFFFF) = 0 ifTrue: [ bit := bit + 16. n := n bitShift: -16 ].
+    (n bitAnd: 16rFF) = 0 ifTrue: [ bit := bit + 8. n := n bitShift: -8 ].
+    (n bitAnd: 16rF) = 0 ifTrue: [ bit := bit + 4. n := n bitShift: -4 ].
+    (n bitAnd: 16r3) = 0 ifTrue: [ bit := bit + 2. n := n bitShift: -2 ].
+    ^bit + (1 - (n bitAnd: 16r1))!
+
 highBit
     "Return the index of the highest order 1 bit of the receiver"
 


--- 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;
        }
     }


--- orig/tests/floatmath.ok
+++ mod/tests/floatmath.ok
@@ -1,248 +1,74 @@
-
-Execution begins...
-returned value is 3.10000
-
-Execution begins...
-returned value is 3.45000
-
-Execution begins...
-returned value is 30000.0
-
-Execution begins...
-returned value is 34500.0
-
-Execution begins...
-returned value is 7.70000
-
-Execution begins...
-returned value is -8.62000
-
-Execution begins...
-returned value is false
-
-Execution begins...
-returned value is true
-
-Execution begins...
-returned value is true
-
-Execution begins...
-returned value is false
-
-Execution begins...
-returned value is true
-
-Execution begins...
-returned value is true
-
-Execution begins...
-returned value is true
-
-Execution begins...
-returned value is true
-
-Execution begins...
-returned value is 0.00180000
-
-Execution begins...
-returned value is 11250.0
-
-Execution begins...
-returned value is 3
-
-Execution begins...
-returned value is 0.141593
-
-Execution begins...
-returned value is 12
-
-Execution begins...
-returned value is 720.000
-
-Execution begins...
-returned value is 2.81250
-
-Execution begins...
-returned value is inf
-
-Execution begins...
-returned value is 'Inf'
-
-Execution begins...
-returned value is -inf
-
-Execution begins...
-returned value is '-Inf'
-
-Execution begins...
-returned value is nan
-
-Execution begins...
-returned value is 'NaN'
-
-Execution begins...
-returned value is '0.0'
-
-Execution begins...
-returned value is '-0.0'
-
-Execution begins...
-returned value is true
-
-Execution begins...
-returned value is false
-
-Execution begins...
-returned value is nan
-
-Execution begins...
-returned value is nan
-
-Execution begins...
-returned value is nan
-
-Execution begins...
-returned value is 5.00000
-
-Execution begins...
-returned value is 5.00000
-
-Execution begins...
-returned value is 7.00000
-
-Execution begins...
-returned value is 0.00000
-
-Execution begins...
-returned value is -0.00000
-
-Execution begins...
-returned value is -0.00000
-
-Execution begins...
-returned value is -0.00000
-
-Execution begins...
-returned value is nan
-
-Execution begins...
-returned value is nan
-
-Execution begins...
-returned value is nan
-
-Execution begins...
-returned value is 7.00000
-
-Execution begins...
-returned value is 7.00000
-
-Execution begins...
-returned value is 7.00000
-
-Execution begins...
-returned value is 0.00000
-
-Execution begins...
-returned value is 0.00000
-
-Execution begins...
-returned value is 0.00000
-
-Execution begins...
-returned value is -0.00000
-
-Execution begins...
-returned value is true
-
-Execution begins...
-returned value is true
-
-Execution begins...
-returned value is true
-
-Execution begins...
-returned value is true
-
-Execution begins...
-returned value is true
-
-Execution begins...
-returned value is true
-
-Execution begins...
-returned value is true
-
-Execution begins...
-returned value is true
-
-Execution begins...
-(-0.0 -0.0 0.0 0.0 true 0.0 0.0 true )
-(-0.0 0.0 0.0 0.0 true -0.0 -0.0 true )
-(0.0 -0.0 -0.0 -0.0 true 0.0 0.0 true )
-(0.0 0.0 0.0 0.0 true 0.0 0.0 true )
-returned value is Array new: 4 "<0>"
-
-Execution begins...
-(-0.0 -0.0 0.0 0.0 true 0.0 0.0 true )
-(-0.0 0.0 0.0 0.0 true -0.0 -0.0 true )
-(0.0 -0.0 -0.0 -0.0 true 0.0 0.0 true )
-(0.0 0.0 0.0 0.0 true 0.0 0.0 true )
-returned value is Array new: 4 "<0>"
-
-Execution begins...
-(-0.0 -0.0 0.0 0.0 true 0.0 0.0 true )
-(0.0 0.0 0.0 0.0 true 0.0 0.0 true )
-returned value is Array new: 2 "<0>"
-
-Execution begins...
-(-0.0 0.0 true true true )
-(0.0 -0.0 true true true )
-returned value is Array new: 2 "<0>"
-
-Execution begins...
-returned value is Float
-
-Execution begins...
-true->1.0e16
-returned value is 1.00000e+16
-
-Execution begins...
-true->1.2345e16
-returned value is 1.23450e+16
-
-Execution begins...
-true->10.0
-returned value is 10.0000
-
-Execution begins...
-true->17.7674749
-returned value is 17.7675
-
-Execution begins...
-true->0.12345
-returned value is 0.123450
-
-Execution begins...
-true->0.0000000012345
-returned value is 1.23450d-09
-
-Execution begins...
-true->1.2345e-9
-returned value is 1.23450e-09
-
-Execution begins...
-true->0.83205029433784
-returned value is 0.832050
-
-Execution begins...
-true->0.832050294337844
-returned value is 0.832050
-
-Execution begins...
-true->0.55470019622523
-returned value is 0.554700
-
-Execution begins...
-true->0.554700196225229
-returned value is 0.554700
+*** floatmath.ok       Wed Nov  8 11:54:47 2006
+--- /Volumes/disk0s8/devel/gst/+build/tests/floatmath.log      Sun Jan 28 
22:31:28 2007
+***************
+*** 204,213 ****
+--- 204,221 ----
+  returned value is Float
+  
+  Execution begins...
++ true->10000000000000000.0
++ returned value is 1.00000d+16
++ 
++ Execution begins...
+  true->1.0e16
+  returned value is 1.00000e+16
+  
+  Execution begins...
++ true->12345000000000000.0
++ returned value is 1.23450d+16
++ 
++ Execution begins...
+  true->1.2345e16
+  returned value is 1.23450e+16
+  
+***************
+*** 246,248 ****
+--- 254,301 ----
+  Execution begins...
+  true->0.554700196225229
+  returned value is 0.554700
++ 
++ Execution begins...
++ returned value is Float class
++ 
++ Execution begins...
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ returned value is FloatD
++ 
++ Execution begins...
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ true
++ returned value is FloatE


--- orig/tests/floatmath.st
+++ mod/tests/floatmath.st
@@ -8,7 +8,7 @@
 
 "======================================================================
 |
-| Copyright (C) 1988, 1989, 1999, 2006  Free Software Foundation.
+| Copyright (C) 1988, 1989, 1999, 2006, 2007  Free Software Foundation.
 | Written by Steve Byrne
 |
 | This file is part of GNU Smalltalk.
@@ -168,7 +168,9 @@ test
     (((Behavior evaluate: self printString) = self) -> self) printNl. 
 ! !
 
+1d16 test!
 1e16 test!
+1.2345d16 test!
 1.2345e16 test!
 10.0 test!
 (20 - 2.2325251) test!
@@ -179,3 +181,55 @@ test
 0.832050294337844 test!
 0.55470019622523 test!
 0.554700196225229 test!
+
+
+"Fun with rounding"
+
+!Float class methodsFor: 'testing'!
+
+assert: aBoolean
+    aBoolean ifFalse: [ self halt ] ifTrue: [ aBoolean printNl ]!
+
+test
+    | p |
+    p := 1 bitShift: self precision - 1.
+    self assert: (self coerce: p+0+(1/4)) asExactFraction = (p+0).
+    self assert: (self coerce: p+0+(1/2)) asExactFraction = (p+0).
+    self assert: (self coerce: p+0+(3/4)) asExactFraction = (p+1).
+    self assert: (self coerce: p+1+(1/4)) asExactFraction = (p+1).
+    self assert: (self coerce: p+1+(1/2)) asExactFraction = (p+2).
+    self assert: (self coerce: p+1+(3/4)) asExactFraction = (p+2).
+    
+    self assert: ((self emin - self precision - 1 to: self emax - 1) 
allSatisfy: [:i |
+        p := (self coerce: 1) timesTwoPower: i.
+       (self coerce: p asExactFraction) = p]).
+    
+    self assert: ((1 to: 1 + self precision - self emin) allSatisfy: [:i |
+        p := (self coerce: 1) timesTwoPower: i negated.
+       (self coerce: (1 bitShift: i) reciprocal negated) = p negated]).
+    
+    "check for negative zero"
+    p := 1 bitShift: 1 + self precision - self emin.
+    self assert: (self coerce: p reciprocal) positive.
+    self assert: (self coerce: p reciprocal negated) negative.
+    
+    "check for infinity"
+    p := 1 bitShift: self emax + 1.
+    self assert: (self coerce: p) = self infinity.
+    self assert: (self coerce: p negated) = self negativeInfinity.
+    
+    p := 1 bitShift: 1 + self precision - self emin.
+    self assert: (self coerce: p / 3) = self infinity.
+    self assert: (self coerce: p / -3) = self negativeInfinity.
+    
+    "check for non infinity/nan"
+    p := 1 bitShift: self emax + 1.
+    self assert: (self coerce: p / 3) isFinite.
+    self assert: (self coerce: p / -3) isFinite.
+    
+    p := 1 bitShift: 1 + self precision - self emin.
+    self assert: (self coerce: 3 / p) isFinite.
+    self assert: (self coerce: -3 / p) isFinite! !
+
+FloatD test!
+FloatE test!



* added files

--- /dev/null
+++ 
/Volumes/disk0s8/devel/gst/,,address@hidden/new-files-archive/./lib-src/.arch-ids/truncl.c.id
@@ -0,0 +1 @@
+Paolo Bonzini <address@hidden> Sun Jan 28 22:17:40 2007 18819.0
--- /dev/null
+++ 
/Volumes/disk0s8/devel/gst/,,address@hidden/new-files-archive/./lib-src/truncl.c
@@ -0,0 +1,79 @@
+/******************************** -*- C -*- ****************************
+ *
+ *      Emulation for truncl
+ *
+ *
+ ***********************************************************************/
+
+/***********************************************************************
+ *
+ * Copyright 2007 Free Software Foundation, Inc.
+ * Written by Paolo Bonzini.
+ *
+ * This file is part of GNU Smalltalk.
+ *
+ * GNU Smalltalk is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the Free
+ * Software Foundation; either version 2, or (at your option) any later
+ * version.
+ *
+ * Linking GNU Smalltalk statically or dynamically with other modules is
+ * making a combined work based on GNU Smalltalk.  Thus, the terms and
+ * conditions of the GNU General Public License cover the whole
+ * combination.
+ *
+ * In addition, as a special exception, the Free Software Foundation
+ * give you permission to combine GNU Smalltalk with free software
+ * programs or libraries that are released under the GNU LGPL and with
+ * independent programs running under the GNU Smalltalk virtual machine.
+ *
+ * You may copy and distribute such a system following the terms of the
+ * GNU GPL for GNU Smalltalk and the licenses of the other code
+ * concerned, provided that you include the source code of that other
+ * code when and as the GNU GPL requires distribution of source code.
+ *
+ * Note that people who make modified versions of GNU Smalltalk are not
+ * obligated to grant this special exception for their modified
+ * versions; it is their choice whether to do so.  The GNU General
+ * Public License gives permission to release a modified version without
+ * this exception; this exception also makes it possible to release a
+ * modified version which carries forward this exception.
+ *
+ * GNU Smalltalk is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
+ * more details.
+ *
+ * You should have received a copy of the GNU General Public License along with
+ * GNU Smalltalk; see the file COPYING.  If not, write to the Free Software
+ * Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ ***********************************************************************/
+
+#include <float.h>
+
+#include "mathl.h"
+
+/* To compute the integer part of X, sum a big enough
+   integer so that the precision of the floating point
+   number is exactly 1.  */
+
+long double
+truncl(long double x)
+{
+  long double y;
+  if (x < 0.0L)
+    {
+      y = -(1.0L / LDBL_EPSILON - x - 1.0 / LDBL_EPSILON);
+      if (y < x)
+        y = y + 1.0L;
+    }
+  else
+    {
+      y = 1.0L / LDBL_EPSILON + x - 1.0 / LDBL_EPSILON;
+      if (y > x)
+        y = y - 1.0L;
+    }
+
+  return y;
+}