```{-# LANGUAGE CPP #-}
{-# LANGUAGE Safe #-}
#endif
module Math.Sequence.Converge where

import Data.List
import Data.Maybe
import Data.Monoid

-- |Take items from the list until two successive items are equal and
-- return the second of them (or an item is not equal to itself, to handle
-- NaN without a 'RealFloat' context.  In this case, the first item of the
-- pair is returned) .  If the list ends before a match is found,
-- returns the last element of the list.
converge :: Eq a => [a] -> a
converge [] = error "converge: empty list"
converge xs = fromMaybe empty (convergeBy eq Just xs)
where
empty = error "converge: programming error in converge implementation"

eq (a:b:_)
| a == b        = Just b
| b /= b        = Just a
eq _ = Nothing

-- |@convergeTo absEps relEps xs@ takes items from @xs@ until two successive
-- items @x@ and @y@ are within either @absEps@ or @relEps * max (abs x) (abs
-- y)@ of each other, in which case the second of the pair is returned, or
-- until an item is found that does not equal itself (which would typically
-- be a NaN), in which case the preceding item is returned.  If the list ends
-- before a match is found, the last element of the list is returned.
--
-- For example, approximating the golden mean by applying Newton's method to
-- find a root of @x^2 - x - 1@:
--
-- > phi :: Rational
-- > phi = convergeTo 1e-100 0 (iterate (\x -> (x*x + 1) / (2*x-1)) 1)
convergeTo :: (Fractional a, Ord a) => a -> a -> [a] -> a
convergeTo absEps relEps [] = error "convergeTo: empty list"
convergeTo absEps relEps xs = fromMaybe empty (convergeBy eq Just xs)
where
empty = error "convergeTo: programming error in convergeTo implementation"
eq (a:b:_)
| b /= b                            = Just a
| absDiff <= abs absEps             = Just b
| absDiff <= abs relEps * relScale  = Just b
where
absDiff = abs (a-b)
relScale = max (abs a) (abs b)
eq _ = Nothing

-- |@convergeBy f end xs@ looks through @xs@ for the first segment for which
-- @f@ returns a value, and returns that value.  Typically @f@ would be
-- something like:
--
-- > f (a:b:_)
-- >    | abs(a-b) <= eps
-- >    = Just (0.5 * (a + b))
-- > f _ = Nothing
--
-- If no such segment is found, applies @end@ to the last item in the list
-- and returns the result.  If the list was empty, returns 'Nothing'.
--
convergeBy :: ([a] -> Maybe b) -> (a -> Maybe b) -> [a] -> Maybe b
convergeBy f end = listToMaybe . catMaybes . map f' . tails
where
f' xs = case f xs of
Nothing -> end' xs
other   -> other
end' [x] = end x
end'  _  = Nothing
```