Nicolas Cellier uploaded a new version of Kernel to project The Trunk:
http://source.squeak.org/trunk/Kernel-nice.1232.mcz ==================== Summary ==================== Name: Kernel-nice.1232 Author: nice Time: 13 May 2019, 5:43:37.253569 pm UUID: 880030b9-91d0-244b-9198-b7b6f9454584 Ancestors: Kernel-nice.1231 Make sure that LargePositiveInteger sqrt are always correctly rounded to nearest Float if inexact. Note that I accidentally introduced a rounding bug in LargePositiveInteger sqrt which was not in original Integer>>#sqrt: (self bitAnd: 1) should have been (integerResult bitAnd: 1). Apologizes. But this formulation was subject to double-rounding error anyway (asFloat did perform a second rounding, which may lead to incorrect rounding), so that did not really matter... Concerning performance, here are the bench results, after correcting the original version (I send twice self asFloat, which was a slip, and the original did not test mightBeASquare either when self asFloat = Float infinity) - When (self mighBeASquare), new code is equivalent, and becomes more performant for huge values This is due to sqrtRem which avoids evaluating a #squared for no additional price Note that 7 large ints out of 256 mightBeASquare - else it depends on bit range, new code is : * as performant in range (SmallInteger maxVal highBit + 1 to: Float precision * 2) that's the case of vast majority of large int in a "normal" image * less performant in range (Float precision * 2 + 1 to: Float emax) by 3x to 5x that's because asFloat sqrt is cheaper than sqrtRem * more performant in range (Float emax + 1 to: infinity) by large this is due to evaluating sqrtFloor on a much smaller LargePositiveInteger (but this trick could have been added in inexact-rounding variant too) So the exact rounding: - does not require additional code (both are equivalent) - does not cost any performance penalty in most common cases (small LargePositiveInteger); the penalty is only in medium range (107 to 1023 bits). In a word, it's worth. =============== Diff against Kernel-nice.1231 =============== Item was changed: ----- Method: LargePositiveInteger>>sqrt (in category 'mathematical functions') ----- sqrt "Answer the square root of the receiver. If the square root is exact, answer an Integer, else answer a Float approximation. + Make sure the result is correctly rounded (i.e. the nearest Float to the exact square root)" - Handle some tricky case of Float overflow or inaccuracy" + | floatResult integerResult guardBit highBit sr | + (highBit := self highBit) < (Float precision * 2) - | selfAsFloat floatResult integerResult | - selfAsFloat := self asFloat. - floatResult := self asFloat sqrt. - floatResult isFinite ifTrue: + ["the sqrt of self asFloat is correctly rounded, so use it" + floatResult := self asFloat sqrt. + self mightBeASquare ifFalse: [^floatResult]. - [self mightBeASquare ifFalse: ["we know for sure no exact solution exists, - just answer the cheap float approximation without wasting time." - ^floatResult]. "Answer integerResult in case of perfect square" integerResult := floatResult truncated. + integerResult squared = self ifTrue: [^integerResult]. + ^floatResult]. - integerResult squared = self ifTrue: [^integerResult]]. - - "In this case, maybe it failed because we are such a big integer that the Float method becomes - inexact, even if we are a whole square number. So, try the slower but more general method" - selfAsFloat >= Float maxExactInteger asFloat squared - ifTrue: [ - integerResult := self sqrtFloor. - integerResult squared = self ifTrue: [ - ^integerResult ]. - - "Nothing else can be done. No exact answer means answer must be a Float. - Answer the best we have which is the rounded sqrt." - integerResult := (self bitShift: 2) sqrtFloor. - ^((integerResult bitShift: -1) + (self bitAnd: 1)) asFloat]. + "Eventually use a guard bit for handling correct rounding direction" + guardBit := highBit <= (Float precision + 1 * 2) + ifTrue: + ["Add one guard bit for rounding correctly" + 1] + ifFalse: + [self mightBeASquare + ifTrue: + ["Keep all the bits in case we are a perfect square" + 0] + ifFalse: + ["Remove superfluous bit that won't change the Float approximation" + Float precision + 1 - (highBit // 2)]]. + + "Get truncated sqrt and remainder for the same price" + sr := (self bitShift: guardBit * 2) sqrtRem. + + "Handle case of perfect square" + integerResult := sr first. + sr last isZero ifTrue: [^integerResult bitShift: guardBit negated]. + + "Answer the best we have which is the sqrt correctly rounded at Float precision." + ^((integerResult bitShift: Float precision - integerResult highBit) + + (integerResult bitAt: integerResult highBit - Float precision)) asFloat + timesTwoPower: integerResult highBit - Float precision - guardBit! - "We need an approximate result" - ^floatResult! Item was added: + ----- Method: LargePositiveInteger>>sqrt2 (in category 'mathematical functions') ----- + sqrt2 + "Answer the square root of the receiver. + If the square root is exact, answer an Integer, else answer a Float approximation. + Handle some tricky case of Float overflow or inaccuracy" + + | selfAsFloat floatResult integerResult | + selfAsFloat := self asFloat. + floatResult := selfAsFloat sqrt. + floatResult isFinite + ifTrue: + [self mightBeASquare ifFalse: ["we know for sure no exact solution exists, + just answer the cheap float approximation without wasting time." + ^floatResult]. + "Answer integerResult in case of perfect square" + integerResult := floatResult truncated. + integerResult squared = self ifTrue: [^integerResult]]. + + "In this case, maybe it failed because we are such a big integer that the Float method becomes + inexact, even if we are a whole square number. So, try the slower but more general method" + selfAsFloat >= Float maxExactInteger asFloat squared + ifTrue: [ + integerResult := self sqrtFloor. + integerResult squared = self ifTrue: [ + ^integerResult ]. + + "Nothing else can be done. No exact answer means answer must be a Float. + Answer the best we have which is the rounded sqrt." + integerResult := (self bitShift: 2) sqrtFloor. + ^((integerResult bitShift: -1) + (self bitAnd: 1)) asFloat]. + + "We need an approximate result" + ^floatResult! |
Free forum by Nabble | Edit this page |