[squeak-dev] The Trunk: Kernel-nice.1429.mcz
nicolas.cellier.aka.nice at gmail.com
Wed Dec 1 02:30:25 UTC 2021
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:
> ==================== 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
> 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
> 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)
> 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).
> 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)
> 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')
> "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
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
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
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...
More information about the Squeak-dev