# 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 $S$ denoted $x$.

$x^2 = S$

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

One of the following conditions are always true.

$x_0 \leq x \leq {S \over x_0}$
$x_0 \geq x \geq {S \over x_0}$

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

$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, $x$, of a real-valued function is the Newton–Raphson method.

$f(x) = 0$

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

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

We can use this general method to calculate square root, $x$, of real value $s$ as well.

$f(x) = x^2 - s = 0$
$f'(x) = 2 x$
$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_{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}$