{-# LANGUAGE CPP #-}
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
{-# 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 :: [a] -> a
converge [] = [Char] -> a
forall a. HasCallStack => [Char] -> a
error [Char]
"converge: empty list"
converge [a]
xs = a -> Maybe a -> a
forall a. a -> Maybe a -> a
fromMaybe a
forall a. a
empty (([a] -> Maybe a) -> (a -> Maybe a) -> [a] -> Maybe a
forall a b. ([a] -> Maybe b) -> (a -> Maybe b) -> [a] -> Maybe b
convergeBy [a] -> Maybe a
forall a. Eq a => [a] -> Maybe a
eq a -> Maybe a
forall a. a -> Maybe a
Just [a]
xs)
    where
        empty :: a
empty = [Char] -> a
forall a. HasCallStack => [Char] -> a
error [Char]
"converge: programming error in converge implementation"
        
        eq :: [a] -> Maybe a
eq (a
a:a
b:[a]
_)
            | a
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
b        = a -> Maybe a
forall a. a -> Maybe a
Just a
b
            | a
b a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
b        = a -> Maybe a
forall a. a -> Maybe a
Just a
a
        eq [a]
_ = Maybe a
forall a. Maybe a
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 :: a -> a -> [a] -> a
convergeTo a
absEps a
relEps [] = [Char] -> a
forall a. HasCallStack => [Char] -> a
error [Char]
"convergeTo: empty list"
convergeTo a
absEps a
relEps [a]
xs = a -> Maybe a -> a
forall a. a -> Maybe a -> a
fromMaybe a
forall a. a
empty (([a] -> Maybe a) -> (a -> Maybe a) -> [a] -> Maybe a
forall a b. ([a] -> Maybe b) -> (a -> Maybe b) -> [a] -> Maybe b
convergeBy [a] -> Maybe a
eq a -> Maybe a
forall a. a -> Maybe a
Just [a]
xs)
    where
        empty :: a
empty = [Char] -> a
forall a. HasCallStack => [Char] -> a
error [Char]
"convergeTo: programming error in convergeTo implementation"
        eq :: [a] -> Maybe a
eq (a
a:a
b:[a]
_)
            | a
b a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
b                            = a -> Maybe a
forall a. a -> Maybe a
Just a
a
            | a
absDiff a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a -> a
forall a. Num a => a -> a
abs a
absEps             = a -> Maybe a
forall a. a -> Maybe a
Just a
b
            | a
absDiff a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a -> a
forall a. Num a => a -> a
abs a
relEps a -> a -> a
forall a. Num a => a -> a -> a
* a
relScale  = a -> Maybe a
forall a. a -> Maybe a
Just a
b
            where
                absDiff :: a
absDiff = a -> a
forall a. Num a => a -> a
abs (a
aa -> a -> a
forall a. Num a => a -> a -> a
-a
b)
                relScale :: a
relScale = a -> a -> a
forall a. Ord a => a -> a -> a
max (a -> a
forall a. Num a => a -> a
abs a
a) (a -> a
forall a. Num a => a -> a
abs a
b)
        eq [a]
_ = Maybe a
forall a. Maybe a
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 :: ([a] -> Maybe b) -> (a -> Maybe b) -> [a] -> Maybe b
convergeBy [a] -> Maybe b
f a -> Maybe b
end = [b] -> Maybe b
forall a. [a] -> Maybe a
listToMaybe ([b] -> Maybe b) -> ([a] -> [b]) -> [a] -> Maybe b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Maybe b] -> [b]
forall a. [Maybe a] -> [a]
catMaybes ([Maybe b] -> [b]) -> ([a] -> [Maybe b]) -> [a] -> [b]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([a] -> Maybe b) -> [[a]] -> [Maybe b]
forall a b. (a -> b) -> [a] -> [b]
map [a] -> Maybe b
f' ([[a]] -> [Maybe b]) -> ([a] -> [[a]]) -> [a] -> [Maybe b]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> [[a]]
forall a. [a] -> [[a]]
tails
    where 
        f' :: [a] -> Maybe b
f' [a]
xs = case [a] -> Maybe b
f [a]
xs of
            Maybe b
Nothing -> [a] -> Maybe b
end' [a]
xs
            Maybe b
other   -> Maybe b
other
        end' :: [a] -> Maybe b
end' [a
x] = a -> Maybe b
end a
x
        end'  [a]
_  = Maybe b
forall a. Maybe a
Nothing