# HG changeset patch # User Stefan Vogel # Date 1151661172 -7200 # Node ID afabe44196a448f228b88a0610ba5b5c583d5259 # Parent cdeae08708cd117e44f30d703b62e28a6a70fa0b Fixes: bug 388: Integer asFloat didn't honor the IEEE default "round-to-nearest even" rounding mode bug 390: #timesTwoPower: generated an overflow and did't keep the class diff -r cdeae08708cd -r afabe44196a4 LimitedPrecisionReal.st --- 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!