[[!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}
$$