Finding reciprocal via Newton Raphson: How to determine initial guess?

Question

I'm learning Newton-Raphson to get the reciprocal of any arbitrary value. I need this in order to accurately emulate how the PlayStation 1 does the divide. They used a modified version of the algorithm.

Trying to understand the basic algorithm I came across this video.

Which is all fine and dandy, makes sense. But then as I started writing the code to get the reciprocal, I wasn't sure how to assign the initial guess value. In his example, he set it to 0.1f to find the reciprocal for 7. I found that if I set it to anything else, I would get totally different results. e.g. If I set it to 0.5f, I would get -217.839 in 3 iterations.

Code:

float GetRecip(float Number, float InitialGuess, int Iterations)
{
    float Result;

    Result = InitialGuess;
    for (int i = 0; i < Iterations; i++)
    {
        Result = Result * (2 - Number*Result);
    }

    return (Result);
}

// Example:
float Recip1 = GetRecip(7, 0.1f, 3); // 0.142847776
float Recip2 = GetRecip(7, 0.5f, 3); // -217.839844

Changing the number of iterations doesn't help, it would yield more drastic different results.

How do I go about finding that initial guess of 0.1f?

First step I would take is to write down the iterations: $y_{k+1}=2 y_k - y_k^2 x$. Now one can do the "energy method" I believe (calculate what the difference in energy between steps are, e.g. some kind of norm.). If the energy diff doesn't decrease its bad news. Maybe you can read up on convergence disks and stuff like that too. — Emil, Jan 16 '17 at 20:00

score 1 · Answer 1 · answered Jan 16 '17 at 14:09

1

The first step in the original link, counting leading zeros resp. the number of hex digits of the number, can be replicated with a loop

x0 = 1
while(N*x0>1) x0 = x0/16

where the division by 16 can be implemented as bit manipulation of the exponent in a floating point format or as a bit shift in a fixed point format.

Or directly extract the exponent from the float format of the input and set a corresponding negated exponent in the initial point.

answered Jan 16 '17 at 14:09

Lutz Lehmann

126,666

Thanks. More a partial answer but still helps. Do you know what they're doing with that unr table? What does it represent? – vexe Jan 16 '17 at 15:21

Alan Campbell · Answer 2 · 2019-12-09T22:49:23.233

The closer your initial guess is, the less likely the routine will "break". So, take advantage of the floating point format...

A float is in the general format: 1.xxx * 2^yyy where xxx and yyy are binary numbers.

If you take the negative of yyy you have an approximation for your first guess.

In your case : take 7, round up to 8 = 2^3.

2^(-3) = 0.125 and works quite nicely.

A similar trick is used with Newton's Method when evaluating square roots. For square roots, take yyy and divide by 2 for 1st guess.

EDIT: If your initial value of 7 is rounded up to 10, the reciprocal is 0.1 (the given first guess). Using binary helps you optimise your guess by fiddling with fields in the floating point format. For a reference to (32 bit) float format, see: https://en.wikipedia.org/wiki/Single-precision_floating-point_format

WimC · Accepted Answer · 2021-06-06T14:07:52.440

The algorithm is obfuscated a bit (among other things) because the GTE works exclusively with fixed point numbers. Barring these details, the algorithm to approximate $d^{-1}$ is as follows.

Determine an integer $m\in[0, 255]$ and an integer $k$ such that $(256 + m) {\cdot} 2^k$ is closest to $d$. That is, round $d$ to nine significant bits.
Use a lookup table $\mathrm t[]$ of integers in the range $[0, 255]$ to approximate $$ (256 + m)^{-1} \approx (257 + \mathrm t[m]) {\cdot} 2^{-17}$$ to nine significant bits.
Perform one Newton-Raphson iteration on the approximation $ r = (257 + \mathrm t[m]) {\cdot} 2^{-(17 + k)} \approx d^{-1}$.

For the example $d=7$ this proceeds as follows.

$7 = (256 + 192) {\cdot} 2^{-6}$, so $m=192$ and $k=-6$.
$\mathrm t[192] = 36$ so $(256 + 192)^{-1} \approx (257 + 36) {\cdot} 2^{-17}$.
With $r=(257 + 36) {\cdot} 2^{-11} = 0.143066{\cdots}$ compute the Newton-Raphson iteration $2r - 7 r^2 = 0.142857{\cdots}$. The absolute error is about $3 {\cdot} 10^{-7}$.

Finding reciprocal via Newton Raphson: How to determine initial guess?

3 Answers3