Compute reciprocal in Python using Newton's method

Question

I am attempting to write a function to compute the reciprocal of a number using Newton's method, however it gives me wrong results despite having directly translated the formula from Wikipedia to code, the result produced seems to go far into negatives when I am expecting a positive result, any ideas as to why this won't work?

Formula:

Code:

def rcp(a: float) -> float:
    a = float(a)
    x0 = a
    for _ in range(5):
        x1 = x0 * (2 - (a * x0))
        x0 = x1

    return x0


print(rcp(5.0))

Note that the wikipedia page doesn't recommend starting with `x0 = a`; for instance, to compute the inverse of 17, they start with `x0 = 0.1`. — Stef, Aug 29 '21 at 11:38
I will try that adjustment and see if that helps Edit: That will make no difference for inputs significantly far away from 17 (e.g 100) — txk2048, Aug 29 '21 at 11:52
Just to bound the iteration so we don't loop forever, however it shouldn't make any difference to the fact that results spiral completely the wrong way — txk2048, Aug 29 '21 at 12:12
@Toothless204 The point is, the starting point `x0` shouldn't be too far from `1/a`. If you start with `a=100`, you want `x0` to be relatively close to `0.01`. The wikipedia page even includes a method to find a good starting point (which boils down to "it's easy for computer to find the reciprocals of powers of 2, so find a power of 2 close to `a`, and start with `x0 =` reciprocal of that power of 2".) — Stef, Aug 30 '21 at 09:03
@Stef How would you find the reciprocal of a power of 2? I did see mentioned it about it being possible using bit shifts but it never went into detail about how that would work — txk2048, Aug 30 '21 at 10:52
@Toothless204 First notice that in our decimal system, finding the reciprocal of a power of 10 is easy. For instance, 1/1000 = 0.001; this is a lot easier to calculate than, say, 1/897≈0.0011148272. Second, remember that numbers are stored using binary in computers; for this reason, the equivalent to "it's easy to calculate the reciprocal of a power of 10 in decimal" is "it's easy to compute the reciprocal of a power of 2 in binary". However, it still requires understanding how floating-point numbers are stored; — Stef, Aug 30 '21 at 11:10
In memory, a floating-point number uses one bit to store its `sign` (positive or negative), a few bits to store its `exponent` in binary, an integer number; and a few bits to store its `mantissa` in binary, again an integer (actually, there is a trick in the representation of the mantissa, but that's not very important). Then these bits represent the number `(-1)**sign * mantissa * 2**(exponent - 1023)`. If you want the reciprocal of a power of 2, you just have to set `sign = 0` (positive), `mantissa = 1`, and choose `exponent` so that `exponent-1023` is minus the exponent of your power of 2. — Stef, Aug 30 '21 at 11:18

score 1 · Answer 1 · answered Aug 29 '21 at 14:09

This is not a coding but a mathematical issue. The Wikipedia page on Newton’s method you referred to contains a section on failure analysis i.e. various ways in which Newton’s method may fail to converge to a root of the function. In general, for this method to work the starting value x0 must be sufficiently close to the root (which in your case is 1/a). What is “sufficiently close” depends on the function and the root location, but you are selecting a starting point which is too far.

Compute reciprocal in Python using Newton's method

1 Answers1