<div dir="ltr"><div dir="ltr"><br></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Le mer. 1 déc. 2021 à 02:35, <<a href="mailto:commits@source.squeak.org">commits@source.squeak.org</a>> a écrit :<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Nicolas Cellier uploaded a new version of Kernel to project The Trunk:<br>
<a href="http://source.squeak.org/trunk/Kernel-nice.1429.mcz" rel="noreferrer" target="_blank">http://source.squeak.org/trunk/Kernel-nice.1429.mcz</a><br>
<br>
==================== Summary ====================<br>
<br>
Name: Kernel-nice.1429<br>
Author: nice<br>
Time: 1 December 2021, 2:35:08.604257 am<br>
UUID: 1a855b7e-d5db-4907-ba56-ab7c28e4ed38<br>
Ancestors: Kernel-nice.1428<br>
<br>
Fixup SmallInteger>>sqrt so that it answers nearest Float to exact square root.<br>
<br>
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.<br>
<br>
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.<br>
<br>
For example take:<br>
x := SmallInteger maxVal - 64.<br>
<br>
We have:<br>
(2**60)-(2**6)-1 = x.<br>
<br>
or in factored form:<br>
(1-(2** -54) - (2** -60)) * (2**60) = x.<br>
<br>
when taking exact square root:<br>
(1-(2** -54) - (2** -60)) sqrt * (2**30) = x sqrt.<br>
<br>
But (1-(2** -54) - (2** -60)) sqrt is approximated with:<br>
(1-(2** -55) - (2** -61) - eps)<br>
where:<br>
eps<(2** -100)<br>
<br>
This term will be rounded asFloat to 1.0, because much closer to 1.0 than to 1.0 predecessor, that is:<br>
(1-(2** -53))<br>
the boundary between those 2 Floats being:<br>
(1-(2** -54))<br>
<br>
We thus expect:<br>
x sqrt asFloat = (2**30).<br>
<br>
However:<br>
y := x asFloat.<br>
will be rounded down to:<br>
(1-(2** -54)-(2** -60)) asFloat * (2**60) = y.<br>
(1-(2** -53)) * (2**60) = y.<br>
because slightly smaller than the boundary (1-(2** -54)) seen above.<br>
<br>
taking y square root:<br>
(1-(2** -53)) sqrt * (2**30) = y sqrt.<br>
<br>
First term is approximated with:<br>
(1-(2** -54) - eps)<br>
where:<br>
eps<(2** -100)<br>
Once again, slightly below the (1-(2** -54)) boundary.<br>
<br>
This time, the result is rounded down asFloat:<br>
(1-(2** -53)) * (2**30) = y sqrt asFloat.<br>
<br>
That is the predecessor of the expected result.<br>
<br>
We observe a double rounding at work:<br>
(SmallInteger maxVal - 64) asFloat is rounded down<br>
(SmallInteger maxVal - 64) asFloat sqrt is rounded down again because slightly less than boundary between 2r1.0e30 predecessor and 2r1.0e30.<br>
<br>
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.<br>
<br>
We'd thus better resort to more careful algorithm already developped in LargePositiveInteger case.<br>
<br>
=============== Diff against Kernel-nice.1428 ===============<br>
<br>
Item was changed:<br>
----- Method: SmallInteger>>sqrt (in category 'mathematical functions') -----<br>
sqrt<br>
"Answer the square root of the receiver.<br>
If the square root is exact, answer an Integer, else answer a Float approximation"<br>
| floatResult integerResult |<br>
self negative ifTrue: [<br>
^ DomainError signal: 'sqrt undefined for number less than zero.' ].<br>
floatResult := self asFloat sqrt.<br>
integerResult := floatResult truncated.<br>
"Note: truncated works for 60-bit SmallInteger<br>
If self is a square s^2, but asFloat rounds down,<br>
+ f = s^2*(1-e), f^0.5 = s*(1-e)^0.5 = s*(1-0.5*e-O(e^2))<br>
- f = s^2*(1-e), f^0.5 = s*(1-e)^0.5 = s*(1-0.5*e+O(e^2))<br>
since s asFloat is exact, and e <= 0.5*ulp(1),<br>
+ s*(1-0.5*e-O(e^2)) always rounds to s"<br>
- s*(1-0.5*e+O(e^2)) always rounds to s"<br></blockquote><div><br></div><div>Hmm, above comment is not very accurate.</div><div>If ever e=0.5*ulp(1), then (1-0.5*e) is the exact tie, and s*(1-0.5*e-O(e^2)) rounds down to s asFloat predecessor.<div>Having written +O(e^2) did cause this confusion, it's rather a negative correction -O(e^2).</div><div>But once we express the negativity, we see the problem...</div><div><br></div></div><div>Fortunately, there is no such exact square in the range of interest:</div><div> (1<<27+1 to: 1<<30-1) noneSatisfy: [:s | 1.0 predecessor asTrueFraction * s squared = s squared asFloat].</div><div><br></div><div>Indeed, this is impossible, 1.0 predecessor significandAsInteger has 53 bits set.</div><div>Multiplying with an Integer with more than two bits set will result in a span of 54 bits or more of set bits.<br></div><div>The only case when (1.0 predecessor asTrueFraction * s squared) isAnExactFloat is when s squared has a single bit set.</div><div>But that means that s^2 is a power of two and isAnExactFloat, so no rounding problem.</div><div><br></div><div><div>Code is OK, but rationale is not.</div></div><div>The rationale is rather a strict inequality e < 0.5*ulp(1), and (0.5*e+O(e^2)) < 0.25*ulp(1) too.</div><div>Maybe I shall publish again...<br></div><br><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
integerResult * integerResult = self ifTrue: [^integerResult].<br>
+ self isAnExactFloat ifTrue: [^floatResult].<br>
+ "self has more bits than Float precision, so self asFloat and floatResult might be inexact.<br>
+ Use large integer algorithm for a guaranty of returning the nearest float to exact result.<br>
+ Note that shifting by 8 is enough because SmallInteger maxVal highBit - Float precision < 8"<br>
+ ^(self << 8) sqrt timesTwoPower: -4!<br>
- ^floatResult!<br>
<br>
<br>
</blockquote></div></div>