{-# 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 [] = 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