[[!redirects Square Root in Haskell]]

# 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}

$$