Nicolas Cellier uploaded a new version of Kernel to project The Trunk:
http://source.squeak.org/trunk/Kernel-nice.1231.mcz ==================== Summary ==================== Name: Kernel-nice.1231 Author: nice Time: 12 May 2019, 10:58:37.936831 pm UUID: fdae6969-331b-4fa3-83b9-4cadb4d84290 Ancestors: Kernel-nice.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 ETA-floats (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 Kernel-nice.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 IEEE-754." | 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 Newton-Raphson 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 under-estimated if (guess-u/2)^2>self, then guess is over-estimated 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 (guess-u/2)^2=guess^2 - guess*u + u^2/4 > self ==> guess^2 - self >= guess*u (guess^2 - self) is evaluated with an emulated fused-multiply-add" ["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! |
Free forum by Nabble | Edit this page |