algorithm - Newton's Method for finding the reciprocal of a floating point number for division

Question

I am trying to divide two numbers, a numerator N by a divisor D. I am using the Newton–Raphson method which uses Newton's method to find the reciprocal of D (1/D). Then the result of the division can be found by multiplying the numerator N by the reciprocal 1/D to get N/D.

The Newton-Raphson algorithm can be found here

So the first step of the algorithm is to start with an initial guess for 1/D which we call X_0.

X_0 is defined as X_0 = 48/17-39/17*D

However, we must first apply a bit-shift to the divisor D to scale it so that 0.5 ≤ D ≤ 1. The same bit-shift should be applied to the numerator N so that the quotient does not change.

We then find X_(i+1) using the formula X_(i+1) = X_i*(2-D*X_i)

Since both the numerator N, divisor D, and result are all floating point IEEE-754 32-bit format, I am wondering how to properly apply this scaling since my value for 1/D does not converge to a value, it just approaches -Inf or +Inf (depending on D).

What I have found works though is that if I make X_0 less than 1/D, the algorithm seems to always converge. So if I just use a lookup table where I always store a bunch of values of 1/D and I can always ensure I have a stored 1/D value where D > Dmin, then I should be okay. But is that standard practice?

Answer 1

1. To set the sign bit correctly, perform the XOR on the sign of the original dividend and divisor.

2. Make the sign of the divisor and dividend positive now.

3. First set the dividend exponent equal to dividend_exponent- 1 - divisor_exponent - 1 + 127. The +127 is for the bias since we just subtracted it out. This scales the dividend by the same amount we will scale the divisor by.

4. Change the divisor exponent to 126 (biased) or -1 (unbiased). This scales the divisor to between 0.5 and 1.

5. Proceed to find Xo with the new scaled D value from step one. Xo = 48/17-32/17 * D.

6. Proceed to find Xn using the new D until we have iterated enough times so that we have the precision we need. X(i+1) = X(i) * (2-D*X(i)). Also, the number of steps S we need is S = ceil(log_2((P + 1)/log_2(17))). Where P is the number of binary places

7. Multiply Xn * N = 1/D * N = N/D and your result should be correct.

Update: This algorithm works correctly.