Python Forum

I'm trying to create a function which calculate root using the Newton method, but with small numbers this doesn't work. Is it problem with Python or with my func? I have tried Decimal, __future__ division, nothing helped.

def square_root(num):
    x1 = (num * 1.0) / 2.0
    x2 = (x1 + (num / x1)) / 2.0
    while abs(x1 - x2) >= num:
        x1 = x2
        x2 = (x1 + (num / x1)) / 2.0
    return x2


print(square_root(4))
print(square_root(4e-128))

Error:
2.0
2.0000000000000003e-64

No it is not a problem with python, nor with your function. It comes from the fact that the computer handles real number with a finite precision. Here is a code to show you which bits are stored in the 64 bits floating point numbers used by python

import struct

def float_to_int64(x):
    """Convert a python float to a 64 bits integer.
    
    Convert a float to an integer representing this
    float in the 64 bits ieee754 double precision format."""
    
    u, v = struct.unpack(">LL", struct.pack(">d", x))
    return (u << 32) | v

def ieee754(x):
    """Returns a representation of the ieee754 of a floating number
    """
    s = '{:0>64}'.format(bin(float_to_int64(x))[2:])
    return ' '.join((s[0], s[1:12], s[12:]))

def square_root(num):
    x1 = (num * 1.0) / 2.0
    x2 = (x1 + (num / x1)) / 2.0
    while abs(x1 - x2) >= num:
        x1 = x2
        x2 = (x1 + (num / x1)) / 2.0
    return x2
 
 
x = 2e-64
y = square_root(4e-128)

print(ieee754(x))
print(ieee754(y))
print('x - y', x - y)

The result is

Output:
0 01100101011 0101000011111111110101000100111101001010011100111101
0 01100101011 0101000011111111110101000100111101001010011100111110
x - y -3.3735033418337674e-80

Your result differs from the best machine approximation of 2e-64 by only -3e-80. As you can see, it differs only on the last bit in machine representation.

You need to read and understand this documentation page