[squeakdev] The Trunk: Kernelnice.1231.mcz
commits at source.squeak.org
commits at source.squeak.org
Sun May 12 20:58:48 UTC 2019
Nicolas Cellier uploaded a new version of Kernel to project The Trunk:
http://source.squeak.org/trunk/Kernelnice.1231.mcz
==================== Summary ====================
Name: Kernelnice.1231
Author: nice
Time: 12 May 2019, 10:58:37.936831 pm
UUID: fdae6969331b4fa383b94cadb4d84290
Ancestors: Kernelnice.1230
Oups, I forgot one germane Float sqrt edge case  sorry to double commit.
0.0 sqrt should answer 0.0 (don't ask why, read the standard...)
I already corrected original code for inf,nan,gradual underflow (very inaccurate), and incorrectly rounded for a few ETAfloats (ETA is 10^18).
I have maybe 20 copies of this fallback code in my change log, and I still missed one case!
Note entirely my fault, that's one of the 4 original lines of code that I did not change...
=============== Diff against Kernelnice.1230 ===============
Item was changed:
 Method: Float>>sqrt (in category 'mathematical functions') 
sqrt
"Fallback code for absent primitives.
Care to answer a correctly rounded result as mandated by IEEE754."
 guess selfScaled nextGuess exp secator hi lo remainder maxError 
self <= 0.0
ifTrue: [self = 0.0
+ ifTrue: [^ self]
 ifTrue: [^ 0.0]
ifFalse: [^ DomainError signal: 'sqrt undefined for number less than zero.']].
self isFinite ifFalse: [^self].
"scale to avoid loss of precision in case of gradual underflow
(selfScaled between: 1.0 and: 2.0), so it is a good guess by itself"
exp := self exponent // 2.
guess := selfScaled := self timesTwoPower: exp * 2.
"Use NewtonRaphson iteration  it converges quadratically
(twice many correct bits at each loop)"
[nextGuess := selfScaled  (guess * guess) / (2.0 * guess) + guess.
"test if both guess are within 1 ulp"
(nextGuess + guess) / 2.0 = guess]
whileFalse:
["always round odd upper  this avoids infinite loop with alternate flip of last bit"
guess := nextGuess + (nextGuess ulp/2.0)].
"adjust the rounding  the guess can be 2 ulp up or 1 ulp down
Let u = guess ulp.
if (guess+u/2)^2<self, then guess is underestimated
if (guessu/2)^2>self, then guess is overestimated
Note that they can't be equal (because left term has 55 bits).
(guess+u/2)^2=guess^2 + guess*u + u^2/4 < self
==> self  guess^2 > guess*u
(guessu/2)^2=guess^2  guess*u + u^2/4 > self
==> guess^2  self >= guess*u
(guess^2  self) is evaluated with an emulated fusedmultiplyadd"
["Decompose guess in two 26 bits parts hi,lo
the trick is that they do not necessarily have the same sign
If 53 bits are hi,0,lo => (hi,lo) else hi,1,lo=> (hi+1,lo)"
secator := "1<<27+1" 134217729.0.
hi := guess * secator.
hi :=hi + (guess  hi).
lo := guess  hi.
"The operations below are all exact"
remainder := selfScaled  hi squared  (hi * lo * 2.0)  lo squared.
maxError := guess timesTwoPower: 1  Float precision.
remainder > maxError or: [remainder negated >= maxError]]
whileTrue: [guess :=remainder > 0.0
ifTrue: [guess successor]
ifFalse: [guess predecessor]].
"undo the scaling"
^ guess timesTwoPower: exp!
More information about the Squeakdev
mailing list
