[Pkg] The Trunk: Kernelnice.1232.mcz
commits at source.squeak.org
commits at source.squeak.org
Mon May 13 15:44:39 UTC 2019
Nicolas Cellier uploaded a new version of Kernel to project The Trunk:
http://source.squeak.org/trunk/Kernelnice.1232.mcz
==================== Summary ====================
Name: Kernelnice.1232
Author: nice
Time: 13 May 2019, 5:43:37.253569 pm
UUID: 880030b991d0244b9198b7b6f9454584
Ancestors: Kernelnice.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 doublerounding 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 inexactrounding 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 Kernelnice.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!
More information about the Packages
mailing list