--- a/LimitedPrecisionReal.st	Fri Jun 30 10:28:54 2006 +0200
+++ b/LimitedPrecisionReal.st	Fri Jun 30 11:52:52 2006 +0200
@@ -77,14 +77,36 @@
      Since floats have a limited precision, you usually loose bits when doing this
      (see Float decimalPrecision, LongFloat decimalPrecision."
 
-    |newFloat sign absVal|
+    |newFloat sign absVal completeAbsVal delta mask|
 
     newFloat := self zero.
     sign := anInteger sign.
-    absVal := anInteger abs.
     (sign ~~ 0) ifTrue:[
+        "be careful to do round to nearest even floating point value -
+         thanks to Nicolas Cellier for this code"
+
+        absVal := anInteger abs.
+        delta := absVal highBit - self precision.
+        delta > 0 ifTrue: [
+            completeAbsVal := absVal.
+            "eliminate insignificant/trailing bits"
+            absVal := absVal bitShift: delta negated.
+            "inexact := trailingBits ~= 0.
+             Round to nearest even"
+            (completeAbsVal bitAt:delta) ~= 0 ifTrue:[
+                mask := (1 bitShift:delta-1) - 1.
+                ((completeAbsVal bitAnd:mask) ~= 0 or:[absVal odd]) ifTrue:[
+                    absVal := absVal + 1
+                ].
+            ].
+        ] ifFalse: [
+            delta := 0
+        ].
         absVal digitLength to:1 by:-1 do:[:i |
-            newFloat := (newFloat * 256.0) + (absVal digitByteAt:i) asFloat
+            newFloat := (newFloat * 256) + (absVal digitByteAt:i).
+        ].
+        delta ~~ 0 ifTrue: [
+            newFloat := newFloat timesTwoPower:delta.
         ].
         (sign < 0) ifTrue:[
             newFloat := newFloat negated
@@ -92,12 +114,16 @@
     ].
     ^ newFloat
 
+
     "
+     ShortFloat fromInteger:2
+     12345678901234567890 asShortFloat            
+
      1234567890 asFloat                     
      1234567890 asFloat asInteger                    
      -1234567890 asFloat asInteger                    
 
-     12345678901234567890 asFloat             
+     12345678901234567890 asFloat storeString            
      12345678901234567890 asFloat asInteger   
      -12345678901234567890 asFloat asInteger   
 
@@ -111,7 +137,16 @@
 
      1234567890123456789012345678901234567890 asLongFloat           
      1234567890123456789012345678901234567890 asLongFloat asInteger  
-     -1234567890123456789012345678901234567890 asLongFloat asInteger  
+     -1234567890123456789012345678901234567890 asLongFloat asInteger
+
+     'this test is on 65 bits'.
+     self assert: 16r1FFFFFFFFFFFF0801 asDouble ~= 16r1FFFFFFFFFFFF0800 asDouble.
+     'this test is on 64 bits'.
+     self assert: 16r1FFFFFFFFFFFF0802 asDouble ~= 16r1FFFFFFFFFFFF0800 asDouble.
+     'nearest even is upper'.
+     self assert: 16r1FFFFFFFFFFF1F800 asDouble = 16r1FFFFFFFFFFF20000 asDouble.
+     'nearest even is lower'.
+     self assert: 16r1FFFFFFFFFFFF0800 asDouble = 16r1FFFFFFFFFFFF0000 asDouble.
     "
 !
 
@@ -132,6 +167,87 @@
     ^ (self fromInteger:fract numerator) / (self fromInteger:fract denominator)
 !
 
+fromNumerator:numerator denominator:denominator
+    "Create a limited precision real from a Rational.
+     This version will answer the nearest flotaing point value,
+     according to IEEE 754 round to nearest even default mode"
+
+    |a b q r exponent floatExponent n ha hb hq q1 eMin zero precision nSubHq|
+
+    a := numerator abs.
+    b := denominator abs.
+    ha := a highBit.
+    hb := b highBit.    
+
+    zero := self zero.
+    precision := self precision.
+
+    "If both numerator and denominator are represented exactly in floating point number,
+     then fastest thing to do is to use hardwired float division"
+    (ha <= precision and:[hb <= precision]) ifTrue: [
+        ^ (zero coerce:numerator) / (zero coerce:denominator)
+    ].
+
+    "Try and obtain a mantissa with 1 bit excess precision.
+     First guess is rough, we might get one more bit or one less"
+    exponent := ha - hb - precision - 1.
+    exponent > 0
+        ifTrue: [b := b bitShift: exponent]
+        ifFalse: [a := a bitShift: exponent negated].
+    q := a quo: b.
+    r := a - (q * b).
+    hq := q highBit.
+
+    "check for gradual underflow, in which case we should use less bits"
+    floatExponent := exponent + hq - 1.
+    eMin := 2 - (1 bitShift: (self numBitsInExponent - 1)).
+    n := floatExponent >= eMin
+        ifTrue: [precision + 1]
+        ifFalse: [precision + 1 + floatExponent - eMin].
+
+    nSubHq := n - hq.
+
+    nSubHq < 0 ifTrue: [
+        exponent := exponent + hq - n.
+        r := (q bitAnd: (1 bitShift:nSubHq) - 1) * b + r.
+        q := q bitShift:nSubHq
+    ] ifFalse:[nSubHq > 0 ifTrue: [
+        exponent := exponent + hq - n.
+        q1 := (r bitShift:nSubHq) quo: b.
+        q := (q bitShift:nSubHq) bitAnd: q1.
+        r := r - (q1 * b)
+    ]].
+
+    "check if we should round upward.
+     The case of exact half (q bitAnd: 1) isZero not & (r isZero)
+     will be handled by self fromInteger: conversion"
+    ((q bitAnd:1) ~~ 0 and:[r ~~ 0]) ifTrue:[
+        q := q + 1
+    ].
+
+    (numerator negative xor:denominator negative) ifTrue: [
+        q := q negated.
+    ].
+
+    ^ (zero coerce:q) timesTwoPower:exponent
+
+    "
+        Time millisecondsToRun:[
+            1000000  timesRepeat:[
+                Float fromNumerator:12345678901234567890 denominator:987654321
+            ].
+        ]
+
+        |fraction|
+        fraction := 12345678901234567890//987654321.
+        Time millisecondsToRun:[
+            1000000  timesRepeat:[
+                fraction asFloat
+            ].
+        ]
+    "
+!
+
 new:aNumber
     "catch this message - not allowed for floats/doubles"
 
@@ -197,7 +313,7 @@
 epsilon
     "return the maximum relative spacing"
 
-    ^ self radix asFloat raisedTo:(1 - self precision)
+    ^ self radix asFloat raisedToInteger:(1 - self precision)
 
     "
       ShortFloat epsilon  
@@ -392,6 +508,31 @@
 
 floor
     ^ self asTrueFraction floor
+!
+
+timesTwoPower:anInteger
+    "multiply self by a power of two.
+     Implementation takes care of preserving class and avoiding overflow/underflow
+     Thanks to Nicolas Cellier for this code"
+
+    |half two|
+
+    two := self coerce:2.
+
+    anInteger abs highBit >= (self class numBitsInExponent - 2) ifTrue:[
+        half := anInteger // 2.
+        ^ (self * (two raisedToInteger: half)) * (two raisedToInteger: anInteger - half)
+    ] ifFalse: [
+        ^ self * (two raisedToInteger: anInteger).
+    ]
+
+
+
+"
+    (1 asShortFloat timesTwoPower: 3) class = ShortFloat.
+    (1 asLongFloat timesTwoPower: 1024).
+    (1 asFloat timesTwoPower: -1024) timesTwoPower: 1024.
+"
 ! !
 
 !LimitedPrecisionReal methodsFor:'coercing & converting'!
@@ -862,6 +1003,7 @@
    ^ 0
 ! !
 
+
 !LimitedPrecisionReal methodsFor:'testing'!
 
 isFinite
@@ -989,7 +1131,7 @@
 !LimitedPrecisionReal class methodsFor:'documentation'!
 
 version
-    ^ '$Header: /cvs/stx/stx/libbasic/LimitedPrecisionReal.st,v 1.63 2006-02-17 11:50:15 cg Exp $'
+    ^ '$Header: /cvs/stx/stx/libbasic/LimitedPrecisionReal.st,v 1.64 2006-06-30 09:52:52 stefan Exp $'
 ! !
 
 LimitedPrecisionReal initialize!