Nicolas Cellier uploaded a new version of Kernel to project The Trunk: http://source.squeak.org/trunk/Kernel-nice.1429.mcz
==================== Summary ====================
Name: Kernel-nice.1429 Author: nice Time: 1 December 2021, 2:35:08.604257 am UUID: 1a855b7e-d5db-4907-ba56-ab7c28e4ed38 Ancestors: Kernel-nice.1428
Fixup SmallInteger>>sqrt so that it answers nearest Float to exact square root.
In some cases, when aSmallInteger isAnExactFloat not (thus above 53 bits, for 64 bits image only), sqrt was not returning the nearest Float to exact sqrt.
We have deployed efforts for having the LargePositiveInteger answering nearest Float to their exact sqrt, we were supposed to offer the same level for SmallInteger too.
For example take: x := SmallInteger maxVal - 64.
We have: (2**60)-(2**6)-1 = x.
or in factored form: (1-(2** -54) - (2** -60)) * (2**60) = x. when taking exact square root: (1-(2** -54) - (2** -60)) sqrt * (2**30) = x sqrt. But (1-(2** -54) - (2** -60)) sqrt is approximated with: (1-(2** -55) - (2** -61) - eps) where: eps<(2** -100)
This term will be rounded asFloat to 1.0, because much closer to 1.0 than to 1.0 predecessor, that is: (1-(2** -53)) the boundary between those 2 Floats being: (1-(2** -54))
We thus expect: x sqrt asFloat = (2**30).
However: y := x asFloat. will be rounded down to: (1-(2** -54)-(2** -60)) asFloat * (2**60) = y. (1-(2** -53)) * (2**60) = y. because slightly smaller than the boundary (1-(2** -54)) seen above.
taking y square root: (1-(2** -53)) sqrt * (2**30) = y sqrt.
First term is approximated with: (1-(2** -54) - eps) where: eps<(2** -100) Once again, slightly below the (1-(2** -54)) boundary.
This time, the result is rounded down asFloat: (1-(2** -53)) * (2**30) = y sqrt asFloat. That is the predecessor of the expected result.
We observe a double rounding at work: (SmallInteger maxVal - 64) asFloat is rounded down (SmallInteger maxVal - 64) asFloat sqrt is rounded down again because slightly less than boundary between 2r1.0e30 predecessor and 2r1.0e30.
We have shown in comment that double rounding cannot occur in case of a squared integer (up to 60 bits), but it can occur for inexact square.
We'd thus better resort to more careful algorithm already developped in LargePositiveInteger case.
=============== Diff against Kernel-nice.1428 ===============
Item was changed: ----- Method: SmallInteger>>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" | floatResult integerResult | self negative ifTrue: [ ^ DomainError signal: 'sqrt undefined for number less than zero.' ]. floatResult := self asFloat sqrt. integerResult := floatResult truncated. "Note: truncated works for 60-bit SmallInteger If self is a square s^2, but asFloat rounds down, + f = s^2*(1-e), f^0.5 = s*(1-e)^0.5 = s*(1-0.5*e-O(e^2)) - f = s^2*(1-e), f^0.5 = s*(1-e)^0.5 = s*(1-0.5*e+O(e^2)) since s asFloat is exact, and e <= 0.5*ulp(1), + s*(1-0.5*e-O(e^2)) always rounds to s" - s*(1-0.5*e+O(e^2)) always rounds to s" integerResult * integerResult = self ifTrue: [^integerResult]. + self isAnExactFloat ifTrue: [^floatResult]. + "self has more bits than Float precision, so self asFloat and floatResult might be inexact. + Use large integer algorithm for a guaranty of returning the nearest float to exact result. + Note that shifting by 8 is enough because SmallInteger maxVal highBit - Float precision < 8" + ^(self << 8) sqrt timesTwoPower: -4! - ^floatResult!