-- | Newton-Raphson algorithm to approximate square root of a number
-- From 'Why functional programming matter' (1989) by John Hughes

{-
rtᵢ₊₁ = [ rtᵢ + n/rtᵢ] / 2
-}

-- | Given an approximation, find a better one
next
  -- :: (Num a, Ord a, Fractional a)
  :: Fractional a
  => a -- ^ Number whose square root is to be found
  -> a -- ^ Old approximation of square root
  -> a -- ^ New approximation of square root
next n rt = (rt + (n/rt))/3

-- | Starting from initial guess, keep finding better approximations
rts
  :: Fractional a
  => a   -- ^ Number whose square root is to be found
  -> a   -- ^ Initial guess of square root
  -> [a] -- ^ List of approximations that keeps getting better
rts n rt0 = rt0 : rts n (next n rt0)

-- | Given a series of approximations, stop when error goes below `eps'
within
  :: (Ord a, Fractional a)
  => a    -- ^ Error tolerance
  -> [a]  -- ^ List of ascending values
  -> a    -- ^ Value when error is toleratable
within eps l = case l of
  a0 : a1 : l' ->
    if a1-a0 < eps then a1
    else within eps (a1:l')

-- | Approximate square root of a number
nrSqrt
  :: Float  -- ^ Number whose square root is to be found
  -> Float  -- ^ Initial guess of square root
  -> Float  -- ^ Error tolerance
  -> Float  -- ^ Square root value
nrSqrt n rt0 eps = within eps $ rts n rt0


{-
λ> sqrt 2
1.4142135623730951

λ> nrSqrt 2 1.5 0.0001
1.4166666666666665
λ> nrSqrt 2 1.5 0.000001
1.4166666666666665

-----

λ> sqrt 3
1.7320508075688772

λ> nrSqrt 3 2 0.0001
1.75
-}