[squeakdev] The Trunk: Kernelnice.1429.mcz
Nicolas Cellier
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:
> http://source.squeak.org/trunk/Kernelnice.1429.mcz
>
> ==================== Summary ====================
>
> Name: Kernelnice.1429
> Author: nice
> Time: 1 December 2021, 2:35:08.604257 am
> UUID: 1a855b7ed5db4907ba56ab7c28e4ed38
> Ancestors: Kernelnice.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 Kernelnice.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 60bit SmallInteger
> If self is a square s^2, but asFloat rounds down,
> + f = s^2*(1e), f^0.5 = s*(1e)^0.5 = s*(10.5*eO(e^2))
>  f = s^2*(1e), f^0.5 = s*(1e)^0.5 = s*(10.5*e+O(e^2))
> since s asFloat is exact, and e <= 0.5*ulp(1),
> + s*(10.5*eO(e^2)) always rounds to s"
>  s*(10.5*e+O(e^2)) always rounds to s"
>
Hmm, above comment is not very accurate.
If ever e=0.5*ulp(1), then (10.5*e) is the exact tie, and
s*(10.5*eO(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<<301) 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/squeakdev/attachments/20211201/2e728264/attachment.html>
More information about the Squeakdev
mailing list
