Le mer. 1 déc. 2021 à 02:35, <commits at source.squeak.org> a écrit :

> 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"
>

Hmm, above comment is not very accurate.
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.
Having written +O(e^2) did cause this confusion, it's rather a negative
correction -O(e^2).
But once we express the negativity, we see the problem...

Fortunately, there is no such exact square in the range of interest:
    (1<<27+1 to: 1<<30-1) noneSatisfy: [:s | 1.0 predecessor asTrueFraction
* s squared = s squared asFloat].

Indeed, this is impossible, 1.0 predecessor significandAsInteger has 53
bits set.
Multiplying with an Integer with more than two bits set will result in a
span of 54 bits or more of set bits.
The only case when (1.0 predecessor asTrueFraction * s squared)
isAnExactFloat is when s squared has a single bit set.
But that means that s^2 is a power of two and isAnExactFloat, so no
rounding problem.

Code is OK, but rationale is not.
The rationale is rather a strict inequality e < 0.5*ulp(1), and
(0.5*e+O(e^2)) < 0.25*ulp(1) too.
Maybe I shall publish again...


        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!
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.squeakfoundation.org/pipermail/squeak-dev/attachments/20211201/2e728264/attachment.html>