[squeakdev] The Trunk: Kernelnice.627.mcz
commits at source.squeak.org
commits at source.squeak.org
Sun Sep 25 10:56:41 UTC 2011
Nicolas Cellier uploaded a new version of Kernel to project The Trunk:
http://source.squeak.org/trunk/Kernelnice.627.mcz
==================== Summary ====================
Name: Kernelnice.627
Author: nice
Time: 25 September 2011, 12:55:51.623 pm
UUID: c09ead193aa6482bb2ccf6549bc67eb6
Ancestors: Kernelul.626
Improve Fraction and Integer asFloat by:
 removing unreachable branch
 removing unneeded 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 Kernelul.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 Squeakdev
mailing list
