{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE ConstrainedClassMethods #-}

-- | The Converge type class.
module Data.CReal.Converge
  ( Converge(..)
  ) where

import Control.Arrow ((&&&))
import Data.Coerce (coerce)
import Data.CReal.Internal (CReal(..), atPrecision, crMemoize)
import Data.Function (on)
import Data.Proxy (Proxy)
import GHC.TypeLits (someNatVal, SomeNat(..))

-- $setup
-- >>> :set -XFlexibleContexts
-- >>> import Data.CReal.Internal

-- | If a type is an instance of Converge then it represents a stream of values
-- which are increasingly accurate approximations of a desired value
class Converge a where
  -- | The type of the value the stream converges to.
  type Element a

  -- | 'converge' is a function that returns the value the stream is converging
  -- to. If given a stream which doens't converge to a single value then
  -- 'converge' will not terminate.
  --
  -- If the stream is empty then it should return nothing.
  --
  -- >>> let initialGuess = 1 :: Double
  -- >>> let improve x = (x + 121 / x) / 2
  -- >>> converge (iterate improve initialGuess)
  -- Just 11.0
  --
  -- >>> converge [] :: Maybe [Int]
  -- Nothing
  converge :: a -> Maybe (Element a)

  -- | 'convergeErr' is a function that returns the value the stream is
  -- converging to. It also takes a function err which returns a value which
  -- varies monotonically with the error of the value in the stream. This can
  -- be used to ensure that when 'convergeErr' terminates when given a
  -- non-converging stream or a stream which enters a cycle close to the
  -- solution. See the documentation for the CReal instance for a caveat with
  -- that implementation.
  --
  -- It's often the case that streams generated with approximation
  -- functions such as Newton's method will generate worse approximations for
  -- some number of steps until they find the "zone of convergence". For these
  -- cases it's necessary to drop some values of the stream before handing it
  -- to convergeErr.
  --
  -- For example trying to find the root of the following funciton @f@ with a
  -- poor choice of starting point. Although this doesn't find the root, it
  -- doesn't fail to terminate.
  --
  -- >>> let f  x = x ^ 3 - 2 * x + 2
  -- >>> let f' x = 3 * x ^ 2 - 2
  -- >>> let initialGuess = 0.1 :: Float
  -- >>> let improve x = x - f x / f' x
  -- >>> let err x = abs (f x)
  -- >>> convergeErr err (iterate improve initialGuess)
  -- Just 1.0142132
  convergeErr :: Ord (Element a) => (Element a -> Element a) -> a -> Maybe (Element a)

-- | Every list of equatable values is an instance of 'Converge'. 'converge'
-- returns the first element which is equal to the succeeding element in the
-- list. If the list ends before the sequence converges the last value is
-- returned.
instance {-# OVERLAPPABLE #-} Eq a => Converge [a] where
  type Element [a] = a

  converge :: [a] -> Maybe (Element [a])
converge = [a] -> Maybe a
forall a. [a] -> Maybe a
lastMay ([a] -> Maybe a) -> ([a] -> [a]) -> [a] -> Maybe a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> Bool) -> [a] -> [a]
forall a. (a -> a -> Bool) -> [a] -> [a]
takeWhilePairwise a -> a -> Bool
forall a. Eq a => a -> a -> Bool
(/=)
  {-# INLINE converge #-}

  convergeErr :: (Element [a] -> Element [a]) -> [a] -> Maybe (Element [a])
convergeErr Element [a] -> Element [a]
err [a]
xs = ((a, a) -> a) -> Maybe (a, a) -> Maybe a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (a, a) -> a
forall a b. (a, b) -> b
snd (Maybe (a, a) -> Maybe a)
-> ([(a, a)] -> Maybe (a, a)) -> [(a, a)] -> Maybe a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(a, a)] -> Maybe (a, a)
forall a. [a] -> Maybe a
lastMay ([(a, a)] -> Maybe (a, a))
-> ([(a, a)] -> [(a, a)]) -> [(a, a)] -> Maybe (a, a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((a, a) -> (a, a) -> Bool) -> [(a, a)] -> [(a, a)]
forall a. (a -> a -> Bool) -> [a] -> [a]
takeWhilePairwise (a -> a -> Bool
forall a. Ord a => a -> a -> Bool
(>) (a -> a -> Bool) -> ((a, a) -> a) -> (a, a) -> (a, a) -> Bool
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` (a, a) -> a
forall a b. (a, b) -> a
fst) ([(a, a)] -> Maybe a) -> [(a, a)] -> Maybe a
forall a b. (a -> b) -> a -> b
$ [(a, a)]
es
    where es :: [(a, a)]
es = (a -> a
Element [a] -> Element [a]
err (a -> a) -> (a -> a) -> a -> (a, a)
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& a -> a
forall a. a -> a
id) (a -> (a, a)) -> [a] -> [(a, a)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [a]
xs
  {-# INLINE convergeErr #-}

-- | The overlapping instance for @'CReal' n@ has a slightly different
-- behavior. The instance for 'Eq' will cause 'converge' to return a value when
-- the list converges to within 2^-n (due to the 'Eq' instance for @'CReal' n@)
-- despite the precision the value is requested at by the surrounding
-- computation. This instance will return a value approximated to the correct
-- precision.
--
-- It's important to note when the error function reaches zero this function
-- behaves like 'converge' as it's not possible to determine the precision at
-- which the error function should be evaluated at.
--
-- Find where log x = π using Newton's method
--
-- >>> let initialGuess = 1
-- >>> let improve x = x - x * (log x - pi)
-- >>> let Just y = converge (iterate improve initialGuess)
-- >>> showAtPrecision 10 y
-- "23.1406"
-- >>> showAtPrecision 50 y
-- "23.1406926327792686"
instance {-# OVERLAPPING #-} Converge [CReal n] where
  type Element [CReal n] = CReal n

  converge :: [CReal n] -> Maybe (Element [CReal n])
converge [] = Maybe (Element [CReal n])
forall a. Maybe a
Nothing
  converge [CReal n]
xs =
    CReal n -> Maybe (CReal n)
forall a. a -> Maybe a
Just (CReal n -> Maybe (CReal n)) -> CReal n -> Maybe (CReal n)
forall a b. (a -> b) -> a -> b
$ (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p ->
      case Integer -> Maybe SomeNat
someNatVal (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
p) of
           Maybe SomeNat
Nothing -> [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error [Char]
"Data.CReal.Converge p should be non negative"
           Just (SomeNat (Proxy n
_ :: Proxy p')) ->
             let modifyPrecision :: [CReal n] -> [CReal n]
modifyPrecision = [CReal n] -> [CReal n]
coerce :: [CReal n] -> [CReal p']
             in ([CReal n] -> CReal n
forall a. [a] -> a
last ([CReal n] -> CReal n)
-> ([CReal n] -> [CReal n]) -> [CReal n] -> CReal n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CReal n -> CReal n -> Bool) -> [CReal n] -> [CReal n]
forall a. (a -> a -> Bool) -> [a] -> [a]
takeWhilePairwise CReal n -> CReal n -> Bool
forall a. Eq a => a -> a -> Bool
(/=) ([CReal n] -> [CReal n])
-> ([CReal n] -> [CReal n]) -> [CReal n] -> [CReal n]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [CReal n] -> [CReal n]
modifyPrecision ([CReal n] -> CReal n) -> [CReal n] -> CReal n
forall a b. (a -> b) -> a -> b
$ [CReal n]
xs) CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p)
  {-# INLINE converge #-}

  convergeErr :: (Element [CReal n] -> Element [CReal n])
-> [CReal n] -> Maybe (Element [CReal n])
convergeErr Element [CReal n] -> Element [CReal n]
_ [] = Maybe (Element [CReal n])
forall a. Maybe a
Nothing
  convergeErr Element [CReal n] -> Element [CReal n]
err [CReal n]
xs =
    CReal n -> Maybe (CReal n)
forall a. a -> Maybe a
Just (CReal n -> Maybe (CReal n)) -> CReal n -> Maybe (CReal n)
forall a b. (a -> b) -> a -> b
$ (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p ->
      case Integer -> Maybe SomeNat
someNatVal (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
p) of
           Maybe SomeNat
Nothing -> [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error [Char]
"Data.CReal.Converge p should be non negative"
           Just (SomeNat (Proxy n
_ :: Proxy p')) ->
             let modifyPrecision :: [CReal n] -> [CReal n]
modifyPrecision = [CReal n] -> [CReal n]
coerce :: [CReal n] -> [CReal p']
                 modifyFunPrecision :: (CReal n -> CReal n) -> CReal n -> CReal n
modifyFunPrecision = (CReal n -> CReal n) -> CReal n -> CReal n
coerce :: (CReal n -> CReal n) -> CReal p' -> CReal p'
                 es :: [(CReal n, CReal n)]
es = ((CReal n -> CReal n) -> CReal n -> CReal n
modifyFunPrecision CReal n -> CReal n
Element [CReal n] -> Element [CReal n]
err (CReal n -> CReal n)
-> (CReal n -> CReal n) -> CReal n -> (CReal n, CReal n)
forall (a :: * -> * -> *) b c c'.
Arrow a =>
a b c -> a b c' -> a b (c, c')
&&& CReal n -> CReal n
forall a. a -> a
id) (CReal n -> (CReal n, CReal n))
-> [CReal n] -> [(CReal n, CReal n)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [CReal n] -> [CReal n]
modifyPrecision [CReal n]
xs
                 continue :: (a, a) -> (a, a) -> Bool
continue (a
e1, a
x1) (a
e2, a
x2) = if a
e1 a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then a
x1 a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= a
x2 else a
e1 a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
e2
             in ((CReal n, CReal n) -> CReal n
forall a b. (a, b) -> b
snd ((CReal n, CReal n) -> CReal n)
-> ([(CReal n, CReal n)] -> (CReal n, CReal n))
-> [(CReal n, CReal n)]
-> CReal n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(CReal n, CReal n)] -> (CReal n, CReal n)
forall a. [a] -> a
last ([(CReal n, CReal n)] -> (CReal n, CReal n))
-> ([(CReal n, CReal n)] -> [(CReal n, CReal n)])
-> [(CReal n, CReal n)]
-> (CReal n, CReal n)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((CReal n, CReal n) -> (CReal n, CReal n) -> Bool)
-> [(CReal n, CReal n)] -> [(CReal n, CReal n)]
forall a. (a -> a -> Bool) -> [a] -> [a]
takeWhilePairwise (CReal n, CReal n) -> (CReal n, CReal n) -> Bool
forall a a. (Num a, Eq a, Ord a) => (a, a) -> (a, a) -> Bool
continue ([(CReal n, CReal n)] -> CReal n)
-> [(CReal n, CReal n)] -> CReal n
forall a b. (a -> b) -> a -> b
$ [(CReal n, CReal n)]
es) CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p)
  {-# INLINE convergeErr #-}

takeWhilePairwise :: (a -> a -> Bool) -> [a] -> [a]
takeWhilePairwise :: (a -> a -> Bool) -> [a] -> [a]
takeWhilePairwise a -> a -> Bool
p (a
x1:a
x2:[a]
xs) = if a
x1 a -> a -> Bool
`p` a
x2
                                   then a
x1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> a -> Bool) -> [a] -> [a]
forall a. (a -> a -> Bool) -> [a] -> [a]
takeWhilePairwise a -> a -> Bool
p (a
x2a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
xs)
                                   else [a
x1]
takeWhilePairwise a -> a -> Bool
_ [a]
xs = [a]
xs

lastMay :: [a] -> Maybe a
lastMay :: [a] -> Maybe a
lastMay [] = Maybe a
forall a. Maybe a
Nothing
lastMay [a]
xs = a -> Maybe a
forall a. a -> Maybe a
Just ([a] -> a
forall a. [a] -> a
last [a]
xs)