This is the code I wrote to compute the newton's square root of a number.
I aimed for extreme clarity over efficiency and I have put an abundant docstring at the start of the newton_sqrt
function where I describe the main formula and explain all the parameters.
This function is not meant to be fast (standard library is 4x times faster and more accurate), I wrote just to learn some good math coding techniques in Python.
- This was written in Python 3 but is fully compatible with Python 2.
- I created constants for the error messages. Is that a good practice? Or should I just write the errors down verbatim when need occurs?
- I chose
n / 2
as my initial guess. Do you know a better initial guess? - Are all the option parameters that I am using necessary or should I delete some of them?
from __future__ import division
import time
def newton_sqrt(n,x0=None,tolerance=10 ** (-3),max_iter=15,epsilon=10 ** (-5)):
"""
This function computes the sqare root of a number
using Newton's method.
After chosing a starting approximate value for x runs:
x = x - f(x) / derivative_of_f(x)
After a given number of iterations or if the result is good enough
the programme will stop and return x.
Parameters:
Compulsory:
@ n: the number of which the square root must be taken.
Optional:
@ x0: the initial approximation of n.
@ tolerance: the maximum error permitted.
@ max_iter: the maximum number of iterations.
@ epsilon: a number so small that dividing by it would break
@ the programme, giving vastly inapproximate results.
"""
DIVISION_BY_TINY_DERIVATIVE_ERROR = """The derivative of f is {} and therefore under the {} limit chosen.
Dividing by it would give a result with a vast error."""
NOT_ENOUGH_ACCURACY_IN_ITER_ERROR = """Solution of {} accuracy for {} was not found in {} iterations"""
if x0 is None:
x0 = n / 2
f = lambda x: x ** 2 - n
fprime = lambda x: 2 * x # The derivative of f(x)
for i in range(max_iter):
y = f(x0)
yprime = fprime(x0)
if(abs(yprime) < epsilon):
raise Exception(DIVISION_BY_TINY_DERIVATIVE_ERROR.format(
yprime, epsilon))
x1 = x0 - y / yprime
if(abs(x1 - x0) / abs(x1) < tolerance):
return x1
x0 = x1
raise Exception(NOT_ENOUGH_ACCURACY_IN_ITER_ERROR.format(
tolerance,n,max_iter))
def test_and_timing(max_):
start = time.time()
n_sqrt = [newton_sqrt(i) for i in range(1, max_)]
print("Newton's sqrt took {} seconds to execute for i 1:{}".format(
time.time() - start, max_))
start = time.time()
b_sqrt = [i ** 0.5 for i in range(1, max_)]
print("Built-in sqrt took {} seconds to execute for i 1:{}".format(
time.time() - start, max_))
print("""Sum of all newton sqrts is: {}
Sum of all built-in sqrts is: {}""".format(sum(n_sqrt), sum(b_sqrt)))
if __name__ == "__main__":
test_and_timing(100000)