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

[nicolas cellier]

Hello,

here are some enhancements that should answer nearest floating point 
value when converting LargeInteger and Fraction asFloat.

For LargePositiveInteger, i have two versions, one un-optimized and one 
optimized for gst.

All this has already been written for Squeak VW ST/X and Dolphin.

Some possible tests following (to be completed from above dialects 
already written tests)

Nicolas



self assert: ((1 bitShift: 52)+0+(1/4)) asFloatD asExactFraction = ((1 
bitShift: 52)+0).
self assert: ((1 bitShift: 52)+0+(1/2)) asFloatD asExactFraction = ((1 
bitShift: 52)+0).
self assert: ((1 bitShift: 52)+0+(3/4)) asFloatD asExactFraction = ((1 
bitShift: 52)+1).
self assert: ((1 bitShift: 52)+1+(1/4)) asFloatD asExactFraction = ((1 
bitShift: 52)+1).
self assert: ((1 bitShift: 52)+1+(1/2)) asFloatD asExactFraction = ((1 
bitShift: 52)+2).
self assert: ((1 bitShift: 52)+1+(3/4)) asFloatD asExactFraction = ((1 
bitShift: 52)+2).

self assert: ((1 bitShift: 23)+0+(1/4)) asFloatE asExactFraction = ((1 
bitShift: 23)+0).
self assert: ((1 bitShift: 23)+0+(1/2)) asFloatE asExactFraction = ((1 
bitShift: 23)+0).
self assert: ((1 bitShift: 23)+0+(3/4)) asFloatE asExactFraction = ((1 
bitShift: 23)+1).
self assert: ((1 bitShift: 23)+1+(1/4)) asFloatE asExactFraction = ((1 
bitShift: 23)+1).
self assert: ((1 bitShift: 23)+1+(1/2)) asFloatE asExactFraction = ((1 
bitShift: 23)+2).
self assert: ((1 bitShift: 23)+1+(3/4)) asFloatE asExactFraction = ((1 
bitShift: 23)+2).

self assert: ((-1075 to: 1023) allSatisfy: [:i | (1.0 timesTwoPower: i) 
asExactFraction asFloatD = (1.0 timesTwoPower: i)]).

"Note that min denormalized is (1.0 timesTwoPower: -1074)"
self assert: ((1 to: 1075) allSatisfy: [:i | (1 bitShift: i) reciprocal 
negated asFloatD = (1.0 timesTwoPower: i negated) negated]).

"Filed out from GNU Smalltalk version 2.3.1 on 27-Jan-2007  23:29:25"!

!Fraction methodsFor: 'private'!
asFloat: characterization
    "Answer the receiver converted to a Float"

    | 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! !

"Filed out from GNU Smalltalk version 2.3.1 on 28-Jan-2007  0:37:07"!

!LargePositiveInteger methodsFor: 'private'!
asFloat: characterization
    "Answer the receiver converted to a Float"

    | nTruncatedBits result byte |

    "check for number bigger than maximum mantissa"
    nTruncatedBits := self highBit - characterization precision.
    nTruncatedBits > 0 
        ifTrue: 
            [| mantissa exponent mask trailingBits inexact carry | 
            mantissa := self bitShift: nTruncatedBits negated.
            exponent := nTruncatedBits.
            mask := (1 bitShift: nTruncatedBits) - 1.
            trailingBits := self bitAnd: mask.
            inexact := trailingBits ~= 0.
            inexact
                ifTrue:
                    ["Apply IEEE 754 round to nearest even default rounding 
mode"
                    carry := self bitAt: nTruncatedBits.
                    (carry = 0 
                        or: [(trailingBits bitAnd: (mask bitShift: -1)) = 0 
and: [mantissa even]]) 
                            ifFalse: [mantissa := mantissa + 1]].
            ^(characterization coerce: mantissa) timesTwoPower: exponent].

    "conversion can be exact, construct Float by successive mul add operations"
    byte := characterization coerce: 256.
    result := characterization coerce: 0.

    self size to: 1 by: -1 do: [ :index |
        result := result * byte + (self at: index).
    ].
    ^result! !

"Filed out from GNU Smalltalk version 2.3.1 on 28-Jan-2007  1:04:10"!

!LargePositiveInteger methodsFor: 'private'!
asFloat: characterization
    "Answer the receiver converted to a Float"

    | nTruncatedBits result byte |

    "check for number bigger than maximum mantissa"
    nTruncatedBits := self highBit - characterization precision.
    nTruncatedBits > 0 
        ifTrue: 
            [| mantissa exponent carry | 
            mantissa := self bitShift: nTruncatedBits negated.
            exponent := nTruncatedBits.
            "Apply IEEE 754 round to nearest even default rounding mode"
            "inexactFlag := (self checkIfLowBitGreaterThan: nTruncatedBits) not"
            carry := self bitAt: nTruncatedBits.
            (carry = 1 
                and: [mantissa odd or: [(self checkIfLowBitGreaterThan: 
nTruncatedBits - 1) not]]) 
                        ifTrue: [mantissa := mantissa + 1].
            ^(characterization coerce: mantissa) timesTwoPower: exponent].

    "conversion can be exact, construct Float by successive mul add operations"
    byte := characterization coerce: 256.
    result := characterization coerce: 0.

    self size to: 1 by: -1 do: [ :index |
        result := result * byte + (self at: index).
    ].
    ^result! !

"Filed out from GNU Smalltalk version 2.3.1 on 28-Jan-2007  1:12:40"!

!LargePositiveInteger methodsFor: 'private'!
checkIfLowBitGreaterThan: n
    "Answer true if all bit of self lesser or equal to n are zero"

    "fast check byte by byte"
    1 to: n // 8 do: [:i | (self digitAt: i) = 0 ifFalse: [^false]].

    "finish checking bit by bit"
    (n // 8) * 8 + 1 to: n do: [:i | (self bitAt: i) = 0 ifFalse: [^false]].

   "all low bits were zero"
    ^true
! !