The Trunk: Kernel-nice.1231.mcz

Previous Topic Next Topic
classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view

The Trunk: Kernel-nice.1231.mcz

Nicolas Cellier uploaded a new version of Kernel to project The Trunk:

==================== 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') -----
  "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]
  ["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!