Wiki2
SquareRoot

Estimating the square root with Haskell

Babylonian method

The first method I am going to try is the Babylonian method. It works as follows.

We want to compute an approximation of the square root of a number SS denoted xx.

x 2=S x^2 = S

The idea here is that, if x 0x_0 is an estimate of xx then Sx 0{S \over x_0} (note: SS=S{S \over \sqrt{S}} = \sqrt{S}) is also an estimate. But the crux to understand here is that if x 0x_0 is an over estimate of xx (x 0=x+ex_0=x+e), then Sx 0{S \over \sqrt{x_0}} is an under estimate. If x 0x_0 is an under estimate then Sx 0{S \over \sqrt{x_0}} is a over-estimate.

One of the following conditions are always true.

x 0xSx 0 x_0 \leq x \leq {S \over x_0}
x 0xSx 0 x_0 \geq x \geq {S \over x_0}

If we use this property we can get a better estimate of xx by computing the mean value of x 0x_0 and Sx 0{S \over x_0} as x 1x_1.

x n+1=x n+Sx n2 x_{n+1} = { x_n + {S \over x_n} \over 2 }

An implementation in Haskell would look like the following…

-- Calculate sqrt(s) in n steps 
babylonianSqrt :: Int -> Double -> Double
babylonianSqrt 0 s = 1
babylonianSqrt n s = (x + s/x) / 2 where x = babylonianSqrt (n-1) s

Newton-Raphson method

Another general method of computing a root, xx, of a real-valued function is the Newton–Raphson method.

f(x)=0 f(x) = 0

The method calculates the value of the function, f(x 0)f(x_0) using an chosen start value, x 0x_0. As an improved estimate of root, xx, we use linear function, g(x)g(x), through the tangent of f(x)f(x) at x=x 0x = x_0 (a linearization of f(x)f(x)) and get the root for that function (which is easier).

g(x)=kx+m g(x) = k x + m
k=f(x 0) k = f'(x_0)
m=f(x 0)kx 0 m = f(x_0) - k x_0
g(x)=f(x 0)x+f(x 0)f(x 0)x 0 g(x) = f'(x_0) x + f(x_0) - f'(x_0) x_0
g(x 1)=0x 1=x 0f(x 0)f(x 0) g(x_1) = 0 \implies x_1 = x_0 - {f(x_0) \over f'(x_0)}
x n=x n1f(x n1)f(x n1) x_n = x_{n-1} - {f(x_{n-1}) \over f'(x_{n-1})}

We can use this general method to calculate square root, xx, of real value ss as well.

f(x)=x 2s=0 f(x) = x^2 - s = 0
f(x)=2x f'(x) = 2 x
x n=x n1x n1 2s2x n1 x_n = x_{n-1} - {x_{n-1}^2 - s \over 2 x_{n-1}}

An implementation in Haskell would look like the following…

-- Calculate sqrt(s) in n steps using Newton-Raphson method
newtonRaphsonSqrt :: Int -> Double -> Double
newtonRaphsonSqrt 0 s = 1
newtonRaphsonSqrt n s = x - ((x*x - s) / (2*x)) where x = newtonRaphsonSqrt (n-1) s

A complete program with both methods presented.

{- Calculate square root according to babylonian method
   See: http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method
   Compile: ghc -O2 SquareRoot.hs
   Run: ./SquareRoot -}

module Main (main) where

-- Calculate sqrt(s) in n steps using babylonian method
babylonianSqrt :: Int -> Double -> Double
babylonianSqrt 0 s = 1
babylonianSqrt n s = (x + s/x) / 2 where x = babylonianSqrt (n-1) s

-- Calculate sqrt(s) in n steps using Newton-Raphson methid
newtonRaphsonSqrt :: Int -> Double -> Double
newtonRaphsonSqrt 0 s = 1
newtonRaphsonSqrt n s = x - ((x*x - s) / (2*x)) where x = newtonRaphsonSqrt (n-1) s

main = print $ "sqrt(2) = " ++ show (babylonianSqrt 3 2) ++ " or " ++ show (newtonRaphsonSqrt 3 2) ++ " or " ++ show (sqrt 2)

Which give the following result.

$ ./SquareRoot 
"sqrt(2) = 1.4142156862745097 or 1.4142156862745099 or 1.4142135623730951"

It might be interesting to note that the Babylonian method and the Newton-Raphson method coincide.

We have:

x n=x n1x n1 2s2x n1=2x n1 2x n1 2+s2x n1=x n1 2+s2x n1=x n1+sx n12 x_n = x_{n-1} - {x_{n-1}^2 - s \over 2 x_{n-1}} = {2x_{n-1}^2 - x_{n-1}^2 + s \over 2 x_{n-1}} ={x_{n-1}^2 + s \over 2 x_{n-1}} ={x_{n-1} + {s \over x_{n-1}} \over 2}