-- | Apply algorithms repeatedly until we’re happy with the result.
module Numerics.ConvergentRecursion (
      retryLinearlyUntilPrecision
    , retryExponentiallyUntilPrecision
    , recurseUntilPrecision
) where

import qualified Data.List as L

-- | Search an infinite list for two consecutive values whose relative distance is
-- smaller than the precision parameter.
findCloseConsecutives :: [Double] -> Double -> Double
findCloseConsecutives :: [Double] -> Double -> Double
findCloseConsecutives [Double]
values Double
precision
  = let closeEnoughPair :: (Double, Double) -> Bool
closeEnoughPair (Double
x,Double
y) = (Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
y)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
precision
        Just (Double
_good, Double
evenBetter) = ((Double, Double) -> Bool)
-> [(Double, Double)] -> Maybe (Double, Double)
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Maybe a
L.find (Double, Double) -> Bool
closeEnoughPair ([Double] -> [(Double, Double)]
forall a. [a] -> [(a, a)]
pairsOf [Double]
values)
    in Double
evenBetter

pairsOf :: [a] -> [(a, a)]
pairsOf :: forall a. [a] -> [(a, a)]
pairsOf [a]
xs = [a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
xs ([a] -> [a]
forall a. HasCallStack => [a] -> [a]
tail [a]
xs)

-- | Retry applying a function until two consecutive results are close enough
-- together.
--
-- Each attempt uses twice as many iterations as the last one (hence
-- exponentially).
retryExponentiallyUntilPrecision
    :: (Int -> Double) -- ^ Function of a number of iterations to perform, e.g. integration subdivisions
    -> Double          -- ^ Precision: relative error threshold between two iterations to accept the result
    -> Double          -- ^ Result
retryExponentiallyUntilPrecision :: (Int -> Double) -> Double -> Double
retryExponentiallyUntilPrecision Int -> Double
f Double
precision = [Double] -> Double -> Double
findCloseConsecutives [Int -> Double
f (Int
2Int -> Integer -> Int
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
n) | Integer
n <- [Integer
1..]] Double
precision

-- | Retry applying a function until two consecutive results are close enough
-- together.
--
-- Each attempt uses one more step than the last one (hence linearly).
retryLinearlyUntilPrecision
    :: (Int -> Double) -- ^ Function of a number of iterations to perform, e.g. integration subdivisions
    -> Double          -- ^ Precision: relative error threshold between two iterations to accept the result
    -> Double          -- ^ Result
retryLinearlyUntilPrecision :: (Int -> Double) -> Double -> Double
retryLinearlyUntilPrecision Int -> Double
f Double
precision = [Double] -> Double -> Double
findCloseConsecutives [Int -> Double
f Int
n | Int
n <- [Int
0..]] Double
precision

-- | Recursively apply a function to a value, until the relative distance between
-- two iterations is below the precision parameter.
--
-- Find the a root of \(x^2-1=0\):
--
--
-- >>> let f x = x^2-1
-- >>> recurseUntilPrecision (Numerics.FindRoot.newtonStep 1e-3 f ) 2 1e-10
-- 1.0
recurseUntilPrecision
    :: (Double -> Double) -- ^ Function to iterate
    -> Double             -- ^ Initial value
    -> Double             -- ^ Desired precision
    -> Double
recurseUntilPrecision :: (Double -> Double) -> Double -> Double -> Double
recurseUntilPrecision Double -> Double
f Double
x0 = [Double] -> Double -> Double
findCloseConsecutives ((Double -> Double) -> Double -> [Double]
forall a. (a -> a) -> a -> [a]
iterate Double -> Double
f Double
x0)