2

I have an Newton-Raphson Square Root Algorithm I am using which computes the single-precision square root of an input value. However using a testbench I input I found that certain input values don't converge to an answer which is closest to the actual square root. When I say actual square root, I mean the result you would get with more precision than 32-bit IEEE-754. As a result, I was wondering what is considered the correct value to be obtained when performing the square root in IEEE-754. Some people on this forum have told me that the closest value is not necessarily the most correct, that is why I am asking.

When computing the square root of the single precision IEEE-754 32-bit value 0x3f7fffff, what is considered the correct result and why?

Furthermore, what is considered the correct result when compute the square root of 0x7F7FFFFF?

animuson
  • 53,861
  • 28
  • 137
  • 147
Veridian
  • 3,531
  • 12
  • 46
  • 80
  • 1
    What result are you getting, and why do you suspect that it's wrong? – John Jul 11 '13 at 17:41
  • I would consider the correct result to be the IEEE-754 32-bit value whose square is the closest to 0x3f7fffff and 0x7F7FFFFF, respectively. – Casey Jul 11 '13 at 17:41
  • 1
    @John, I am getting the result of 1 = 0x3f800000 – Veridian Jul 11 '13 at 17:42
  • @StephenCanon, you know me better than that. This isn't homework. You've answered my questions in the past. – Veridian Jul 11 '13 at 17:42
  • 1
    Basic information any request for help should always include: 1) What did you do? 2) What happened? 3) What did you expect to happen? This is a universal principle, not unique to programming or even to technology – Nemo Jul 11 '13 at 17:50
  • 1
    @Nemo, I want to just know what the accepted value is. I thought extraneous information isn't relevant. – Veridian Jul 11 '13 at 17:53
  • 1
    @starbox: sorry, didn’t notice that you were asking the question. I’ll post an answer. – Stephen Canon Jul 11 '13 at 17:56

1 Answers1

6

0x3f7fffff is 1.0 - u, where u = 2**-24. The Taylor Series for sqrt(1 + x) is:

sqrt(1 + x) = 1 + x/2 - x^2/8 + O(x^3)

If we plug -u in for x, we get:

sqrt(1 - u) = 1 - u/2 - u^2/8 - O(u^3)

The value 1 - u/2 is the exact halfway point between the two closest representable floating-point numbers, 1-u and 1; since the next term in the Taylor series is negative, the value of sqrt(1 - u) is just a tiny bit smaller, and so the result rounds down to 1 - u.

0x7f7fffff is just 2**128*(1-u), so the mathematically exact square root is 2**64*(1 - u/2 - u^2/8 - ...), which rounds down to 2**64 * (1-u), as described above.

Stephen Canon
  • 103,815
  • 19
  • 183
  • 269
  • Hmm, so you're saying the rounding mode should not be round to nearest? it should be round down? – Veridian Jul 11 '13 at 18:07
  • 1
    Sorry, I did the analysis for round-to-nearest only; it rounds down to `1 - u` in round-to-zero, round-to-minus-infinity, round-to-nearest, schoolbook rounding, and round-to-odd. It rounds up to `1` in round-to-plus-infinity. – Stephen Canon Jul 11 '13 at 18:09
  • Hmm, I'm using the algorithm shown here: http://stackoverflow.com/questions/17580859/how-to-fix-this-floating-point-square-root-algorithm, do you know a better one? – Veridian Jul 11 '13 at 18:12
  • 1
    Does IEEE-754 say anything about whether a `sqrt` implementation *must* obtain such a value? – Oliver Charlesworth Jul 11 '13 at 18:24
  • 1
    @OliCharlesworth: yes, of course; `sqrt` is one of the fundamental operations that *shall* be correctly rounded in order to comply with the standard. – Stephen Canon Jul 11 '13 at 19:48
  • 2
    @StephenCanon: The only “to” rounding modes are “to nearest” and “to odd”. The others are all “toward”; e.g. one rounds in the direction toward zero, not all the way to zero. Although that would simplify implementation and improve speed. :-) – Eric Postpischil Jul 11 '13 at 20:26