[Pkg] The Trunk: Kernel-nice.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/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!



More information about the Packages mailing list