[squeak-dev] The Trunk: Kernel-nice.627.mczcommits at source.squeak.org commits at source.squeak.orgSun Sep 25 10:56:41 UTC 2011
Nicolas Cellier uploaded a new version of Kernel to project The Trunk: http://source.squeak.org/trunk/Kernel-nice.627.mcz ==================== Summary ==================== Name: Kernel-nice.627 Author: nice Time: 25 September 2011, 12:55:51.623 pm UUID: c09ead19-3aa6-482b-b2cc-f6549bc67eb6 Ancestors: Kernel-ul.626 Improve Fraction and Integer asFloat by: - removing unreachable branch - removing un-needed arithmetic - renaming temps - using Float precision explicitly rather than hardcoded magical numbers - clarifying comments The resulting code shall be clearer and a bit faster than previously. =============== Diff against Kernel-ul.626 =============== Item was changed: ----- Method: Fraction>>asFloat (in category 'converting') ----- asFloat "Answer a Float that closely approximates the value of the receiver. + This implementation will answer the closest floating point number to the receiver. + In case of a tie, it will use the IEEE 754 round to nearest even mode. + In case of overflow, it will answer +/- Float infinity." + + | a b mantissa exponent hasTruncatedBits lostBit n ha hb hm | - This implementation will answer the closest floating point number to - the receiver. - It uses the IEEE 754 round to nearest even mode - (can happen in case denominator is a power of two)" - - | a b q r exponent floatExponent n ha hb hq q1 | a := numerator abs. + b := denominator. "denominator is always positive" - b := denominator abs. ha := a highBit. hb := b highBit. + "Number of bits to keep in mantissa plus one to handle rounding." + n := 1 + Float 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 < n and: [hb < n]) ifTrue: [^numerator asFloat / denominator asFloat]. + + "Shift the fraction by a power of two exponent so as to obtain a mantissa with n bits. + First guess is rough, the mantissa might have n+1 bits." + exponent := ha - hb - n. + exponent >= 0 - then fastest thing to do is to use hardwired float division" - (ha < 54 and: [hb < 54]) ifTrue: [^numerator asFloat / denominator asFloat]. - - "Try and obtain a mantissa with 54 bits. - First guess is rough, we might get one more bit or one less" - exponent := ha - hb - 54. - exponent > 0 ifTrue: [b := b bitShift: exponent] ifFalse: [a := a bitShift: exponent negated]. + mantissa := a quo: b. + hasTruncatedBits := a > (mantissa * b). + hm := mantissa highBit. - q := a quo: b. - r := a - (q * b). - hq := q highBit. + "Check for gradual underflow, in which case the mantissa will loose bits. + Keep at least one bit to let underflow preserve the sign of zero." + lostBit := Float emin - (exponent + hm - 1). + lostBit > 0 ifTrue: [n := n - lostBit max: 1]. + + "Remove excess bits in the mantissa." + hm > n + ifTrue: + [exponent := exponent + hm - n. + hasTruncatedBits := hasTruncatedBits or: [mantissa anyBitOfMagnitudeFrom: 1 to: hm - n]. + mantissa := mantissa bitShift: n - hm]. + + "Check if mantissa must be rounded upward. + The case of tie (mantissa odd & hasTruncatedBits not) + will be handled by Integer>>asFloat." + (hasTruncatedBits and: [mantissa odd]) + ifTrue: [mantissa := mantissa + 1]. + - "check for gradual underflow, in which case we should use less bits" - floatExponent := exponent + hq - 1. - n := floatExponent > -1023 - ifTrue: [54] - ifFalse: [54 + floatExponent + 1022]. - - hq > n - ifTrue: [exponent := exponent + hq - n. - r := (q bitAnd: (1 bitShift: hq - n) - 1) * b + r. - q := q bitShift: n - hq]. - hq < n - ifTrue: [exponent := exponent + hq - n. - q1 := (r bitShift: n - hq) quo: b. - q := (q bitShift: n - hq) bitAnd: q1. - r := (r bitShift: n - hq) - (q1 * b)]. - - "check if we should round upward. - The case of exact half (q bitAnd: 1) isZero not & (r isZero) - will be handled by Integer>>asFloat" - ((q bitAnd: 1) isZero or: [r isZero]) - ifFalse: [q := q + 1]. - ^ (self positive + ifTrue: [mantissa asFloat] + ifFalse: [mantissa asFloat negated]) - ifTrue: [q asFloat] - ifFalse: [q = 0 - ifTrue: [Float negativeZero] - ifFalse: [q asFloat negated]]) timesTwoPower: exponent! Item was changed: ----- Method: Integer>>asFloat (in category 'converting') ----- asFloat + "Answer a Float that best approximates the value of the receiver." - "Answer a Float that represents the value of the receiver. - Optimized to process only the significant digits of a LargeInteger. - SqR: 11/30/1998 21:1 + self subclassResponsibility! - This algorithm does honour IEEE 754 round to nearest even mode. - Numbers are first rounded on nearest integer on 53 bits. - In case of exact half difference between two consecutive integers (2r0.1), - there are two possible choices (two integers are as near, 0 and 1) - In this case, the nearest even integer is chosen. - examples (with less than 53bits for clarity) - 2r0.00001 is rounded to 2r0 - 2r1.00001 is rounded to 2.1 - 2r0.1 is rounded to 2r0 (nearest event) - 2r1.1 is rounded to 2.10 (neraest even) - 2r0.10001 is rounded to 2r1 - 2r1.10001 is rounded to 2.10" - - | abs shift sum delta mask trailingBits carry | - self isZero - ifTrue: [^ 0.0]. - abs := self abs. - - "Assume Float is a double precision IEEE 754 number with 53bits mantissa. - We should better use some Float class message for that (Float precision)..." - delta := abs highBit - 53. - delta > 0 - ifTrue: [mask := (1 bitShift: delta) - 1. - trailingBits := abs bitAnd: mask. - "inexact := trailingBits isZero not." - carry := trailingBits bitShift: 1 - delta. - abs := abs bitShift: delta negated. - shift := delta. - (carry isZero - or: [(trailingBits bitAnd: (mask bitShift: -1)) isZero - and: [abs even]]) - ifFalse: [abs := abs + 1]] - ifFalse: [shift := 0]. - - "now, abs has no more than 53 bits, we can do exact floating point arithmetic" - sum := 0.0. - 1 to: abs size do: - [:byteIndex | - sum := ((abs digitAt: byteIndex) asFloat timesTwoPower: shift) + sum. - shift := shift + 8]. - ^ self positive - ifTrue: [sum] - ifFalse: [sum negated]! Item was added: + ----- Method: LargeNegativeInteger>>asFloat (in category 'converting') ----- + asFloat + ^self negated asFloat negated! Item was added: + ----- Method: LargePositiveInteger>>asFloat (in category 'converting') ----- + asFloat + "Answer a Float that best approximates the value of the receiver. + This algorithm is optimized to process only the significant digits of a LargeInteger. + And it does honour IEEE 754 round to nearest even mode in case of excess precision (see details below)." + + "How numbers are rounded in IEEE 754 default rounding mode: + A shift is applied so that the highest 53 bits are placed before the floating point to form a mantissa. + The trailing bits form the fraction part placed after the floating point. + This fractional number must be rounded to the nearest integer. + If fraction part is 2r0.1, exactly between two consecutive integers, there is a tie. + The nearest even integer is chosen in this case. + Examples (First 52bits of mantissa are omitted for brevity): + 2r0.00001 is rounded downward to 2r0 + 2r1.00001 is rounded downward to 2r1 + 2r0.1 is a tie and rounded to 2r0 (nearest even) + 2r1.1 is a tie and rounded to 2r10 (nearest even) + 2r0.10001 is rounded upward to 2r1 + 2r1.10001 is rounded upward to 2r10 + Thus, if the next bit after floating point is 0, the mantissa is left unchanged. + If next bit after floating point is 1, an odd mantissa is always rounded upper. + An even mantissa is rounded upper only if the fraction part is not a tie." + + "Algorihm details: + Floating point hardware will correctly handle the rounding by itself with a single inexact operation if mantissa has one excess bit of precision. + Except in the last case when extra bits are present after an even mantissa, we must round upper by ourselves. + Note 1: the inexact flag in floating point hardware must not be trusted because it won't take into account the bits we truncated by ourselves. + Note 2: the floating point hardware is presumed configured in default rounding mode." + + | mantissa shift sum excess | + + "Check how many bits excess the maximum precision of a Float mantissa." + excess := self highBitOfMagnitude - Float precision. + excess > 1 + ifTrue: + ["Remove the excess bits but one." + mantissa := self bitShift: 1 - excess. + shift := excess - 1. + "Handle the case of extra bits truncated after an even mantissa." + ((mantissa bitAnd: 2r11) = 2r01 and: [self anyBitOfMagnitudeFrom: 1 to: shift]) + ifTrue: [mantissa := mantissa + 1]] + ifFalse: + [mantissa := self. + shift := 0]. + + "Now that mantissa has at most 1 excess bit of precision, let floating point operations perform the final rounding." + sum := 0.0. + 1 to: mantissa digitLength do: + [:byteIndex | + sum := sum + ((mantissa digitAt: byteIndex) asFloat timesTwoPower: shift). + shift := shift + 8]. + ^sum!
More information about the Squeak-dev mailing list |