{-# LANGUAGE BangPatterns              #-}
{-# LANGUAGE ExistentialQuantification #-}
-- |
-- Module    : Numeric.Series
-- Copyright : (c) 2016 Alexey Khudyakov
-- License   : BSD3
--
-- Maintainer  : alexey.skladnoy@gmail.com, bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Functions for working with infinite sequences. In particular
-- summation of series and evaluation of continued fractions.
module Numeric.Series (
    -- * Data type for infinite sequences.
    Sequence(..)
    -- * Constructors
  , enumSequenceFrom
  , enumSequenceFromStep
  , scanSequence
    -- * Summation of series
  , sumSeries
  , sumPowerSeries
  , sequenceToList
    -- * Evaluation of continued fractions
  , evalContFractionB
  ) where

import Control.Applicative
import Data.List (unfoldr)

import Numeric.MathFunctions.Constants (m_epsilon)


----------------------------------------------------------------

-- | Infinite series. It's represented as opaque state and step
--   function.
data Sequence a = forall s. Sequence s (s -> (a,s))

instance Functor Sequence where
  fmap :: (a -> b) -> Sequence a -> Sequence b
fmap a -> b
f (Sequence s
s0 s -> (a, s)
step) = s -> (s -> (b, s)) -> Sequence b
forall a s. s -> (s -> (a, s)) -> Sequence a
Sequence s
s0 (\s
s -> let (a
a,s
s') = s -> (a, s)
step s
s in (a -> b
f a
a, s
s'))
  {-# INLINE fmap #-}

instance Applicative Sequence where
  pure :: a -> Sequence a
pure a
a = () -> (() -> (a, ())) -> Sequence a
forall a s. s -> (s -> (a, s)) -> Sequence a
Sequence () (\() -> (a
a,()))
  Sequence s
sA s -> (a -> b, s)
fA <*> :: Sequence (a -> b) -> Sequence a -> Sequence b
<*> Sequence s
sB s -> (a, s)
fB = (s, s) -> ((s, s) -> (b, (s, s))) -> Sequence b
forall a s. s -> (s -> (a, s)) -> Sequence a
Sequence (s
sA,s
sB) (((s, s) -> (b, (s, s))) -> Sequence b)
-> ((s, s) -> (b, (s, s))) -> Sequence b
forall a b. (a -> b) -> a -> b
$ \(!s
sa,!s
sb) ->
    let (a -> b
a,s
sa') = s -> (a -> b, s)
fA s
sa
        (a
b,s
sb') = s -> (a, s)
fB s
sb
    in (a -> b
a a
b, (s
sa',s
sb'))
  {-# INLINE pure  #-}
  {-# INLINE (<*>) #-}

-- | Elementwise operations with sequences
instance Num a => Num (Sequence a) where
  + :: Sequence a -> Sequence a -> Sequence a
(+) = (a -> a -> a) -> Sequence a -> Sequence a -> Sequence a
forall (f :: * -> *) a b c.
Applicative f =>
(a -> b -> c) -> f a -> f b -> f c
liftA2 a -> a -> a
forall a. Num a => a -> a -> a
(+)
  * :: Sequence a -> Sequence a -> Sequence a
(*) = (a -> a -> a) -> Sequence a -> Sequence a -> Sequence a
forall (f :: * -> *) a b c.
Applicative f =>
(a -> b -> c) -> f a -> f b -> f c
liftA2 a -> a -> a
forall a. Num a => a -> a -> a
(*)
  (-) = (a -> a -> a) -> Sequence a -> Sequence a -> Sequence a
forall (f :: * -> *) a b c.
Applicative f =>
(a -> b -> c) -> f a -> f b -> f c
liftA2 (-)
  {-# INLINE (+) #-}
  {-# INLINE (*) #-}
  {-# INLINE (-) #-}
  abs :: Sequence a -> Sequence a
abs         = (a -> a) -> Sequence a -> Sequence a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> a
forall a. Num a => a -> a
abs
  signum :: Sequence a -> Sequence a
signum      = (a -> a) -> Sequence a -> Sequence a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> a
forall a. Num a => a -> a
signum
  fromInteger :: Integer -> Sequence a
fromInteger = a -> Sequence a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (a -> Sequence a) -> (Integer -> a) -> Integer -> Sequence a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> a
forall a. Num a => Integer -> a
fromInteger
  {-# INLINE abs         #-}
  {-# INLINE signum      #-}
  {-# INLINE fromInteger #-}

-- | Elementwise operations with sequences
instance Fractional a => Fractional (Sequence a) where
  / :: Sequence a -> Sequence a -> Sequence a
(/)          = (a -> a -> a) -> Sequence a -> Sequence a -> Sequence a
forall (f :: * -> *) a b c.
Applicative f =>
(a -> b -> c) -> f a -> f b -> f c
liftA2 a -> a -> a
forall a. Fractional a => a -> a -> a
(/)
  recip :: Sequence a -> Sequence a
recip        = (a -> a) -> Sequence a -> Sequence a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> a
forall a. Fractional a => a -> a
recip
  fromRational :: Rational -> Sequence a
fromRational = a -> Sequence a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (a -> Sequence a) -> (Rational -> a) -> Rational -> Sequence a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> a
forall a. Fractional a => Rational -> a
fromRational
  {-# INLINE (/)          #-}
  {-# INLINE recip        #-}
  {-# INLINE fromRational #-}



----------------------------------------------------------------
-- Constructors
----------------------------------------------------------------

-- | @enumSequenceFrom x@ generate sequence:
--
-- \[ a_n = x + n \]
enumSequenceFrom :: Num a => a -> Sequence a
enumSequenceFrom :: a -> Sequence a
enumSequenceFrom a
i = a -> (a -> (a, a)) -> Sequence a
forall a s. s -> (s -> (a, s)) -> Sequence a
Sequence a
i (\a
n -> (a
n,a
na -> a -> a
forall a. Num a => a -> a -> a
+a
1))
{-# INLINE enumSequenceFrom #-}

-- | @enumSequenceFromStep x d@ generate sequence:
--
-- \[ a_n = x + nd \]
enumSequenceFromStep :: Num a => a -> a -> Sequence a
enumSequenceFromStep :: a -> a -> Sequence a
enumSequenceFromStep a
n a
d = a -> (a -> (a, a)) -> Sequence a
forall a s. s -> (s -> (a, s)) -> Sequence a
Sequence a
n (\a
i -> (a
i,a
ia -> a -> a
forall a. Num a => a -> a -> a
+a
d))
{-# INLINE enumSequenceFromStep #-}

-- | Analog of 'scanl' for sequence.
scanSequence :: (b -> a -> b) -> b -> Sequence a -> Sequence b
{-# INLINE scanSequence #-}
scanSequence :: (b -> a -> b) -> b -> Sequence a -> Sequence b
scanSequence b -> a -> b
f b
b0 (Sequence s
s0 s -> (a, s)
step) = (b, s) -> ((b, s) -> (b, (b, s))) -> Sequence b
forall a s. s -> (s -> (a, s)) -> Sequence a
Sequence (b
b0,s
s0) (((b, s) -> (b, (b, s))) -> Sequence b)
-> ((b, s) -> (b, (b, s))) -> Sequence b
forall a b. (a -> b) -> a -> b
$ \(b
b,s
s) ->
  let (a
a,s
s') = s -> (a, s)
step s
s
      b' :: b
b'     = b -> a -> b
f b
b a
a
  in (b
b,(b
b',s
s'))


----------------------------------------------------------------
-- Evaluation of series
----------------------------------------------------------------

-- | Calculate sum of series
--
-- \[ \sum_{i=0}^\infty a_i \]
--
-- Summation is stopped when
--
-- \[ a_{n+1} < \varepsilon\sum_{i=0}^n a_i \]
--
-- where ε is machine precision ('m_epsilon')
sumSeries :: Sequence Double -> Double
{-# INLINE sumSeries #-}
sumSeries :: Sequence Double -> Double
sumSeries (Sequence s
sInit s -> (Double, s)
step)
  = Double -> s -> Double
go Double
x0 s
s0
  where 
    (Double
x0,s
s0) = s -> (Double, s)
step s
sInit
    go :: Double -> s -> Double
go Double
x s
s | Double -> Double
forall a. Num a => a -> a
abs (Double
dDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
x) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
m_epsilon = Double -> s -> Double
go Double
x' s
s'
           | Bool
otherwise              = Double
x'
      where (Double
d,s
s') = s -> (Double, s)
step s
s
            x' :: Double
x'     = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d

-- | Calculate sum of series
--
-- \[ \sum_{i=0}^\infty x^ia_i \]
--
-- Calculation is stopped when next value in series is less than
-- ε·sum.
sumPowerSeries :: Double -> Sequence Double -> Double
sumPowerSeries :: Double -> Sequence Double -> Double
sumPowerSeries Double
x Sequence Double
ser = Sequence Double -> Double
sumSeries (Sequence Double -> Double) -> Sequence Double -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double)
-> Sequence Double -> Sequence Double -> Sequence Double
forall (f :: * -> *) a b c.
Applicative f =>
(a -> b -> c) -> f a -> f b -> f c
liftA2 Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) ((Double -> Double -> Double)
-> Double -> Sequence Double -> Sequence Double
forall b a. (b -> a -> b) -> b -> Sequence a -> Sequence b
scanSequence Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) Double
1 (Double -> Sequence Double
forall (f :: * -> *) a. Applicative f => a -> f a
pure Double
x)) Sequence Double
ser
{-# INLINE sumPowerSeries #-}

-- | Convert series to infinite list
sequenceToList :: Sequence a -> [a]
sequenceToList :: Sequence a -> [a]
sequenceToList (Sequence s
s s -> (a, s)
f) = (s -> Maybe (a, s)) -> s -> [a]
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr ((a, s) -> Maybe (a, s)
forall a. a -> Maybe a
Just ((a, s) -> Maybe (a, s)) -> (s -> (a, s)) -> s -> Maybe (a, s)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. s -> (a, s)
f) s
s



----------------------------------------------------------------
-- Evaluation of continued fractions
----------------------------------------------------------------

-- |
-- Evaluate continued fraction using modified Lentz algorithm.
-- Sequence contain pairs (a[i],b[i]) which form following expression:
--
-- \[
-- b_0 + \frac{a_1}{b_1+\frac{a_2}{b_2+\frac{a_3}{b_3 + \cdots}}}
-- \]
--
-- Modified Lentz algorithm is described in Numerical recipes 5.2
-- "Evaluation of Continued Fractions"
evalContFractionB :: Sequence (Double,Double) -> Double
{-# INLINE evalContFractionB #-}
evalContFractionB :: Sequence (Double, Double) -> Double
evalContFractionB (Sequence s
sInit s -> ((Double, Double), s)
step)
  = let ((Double
_,Double
b0),s
s0) = s -> ((Double, Double), s)
step s
sInit
        f0 :: Double
f0          = Double -> Double
maskZero Double
b0
    in  Double -> Double -> Double -> s -> Double
go Double
f0 Double
f0 Double
0 s
s0
  where
    tiny :: Double
tiny = Double
1e-60
    maskZero :: Double -> Double
maskZero Double
0 = Double
tiny
    maskZero Double
x = Double
x
    
    go :: Double -> Double -> Double -> s -> Double
go Double
f Double
c Double
d s
s
      | Double -> Double
forall a. Num a => a -> a
abs (Double
delta Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
m_epsilon = Double -> Double -> Double -> s -> Double
go Double
f' Double
c' Double
d' s
s'
      | Bool
otherwise                    = Double
f'
      where
          ((Double
a,Double
b),s
s') = s -> ((Double, Double), s)
step s
s
          d' :: Double
d'    = Double -> Double
forall a. Fractional a => a -> a
recip (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double
maskZero (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
d
          c' :: Double
c'    = Double -> Double
maskZero (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
aDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
c 
          delta :: Double
delta = Double
c'Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
d'
          f' :: Double
f'    = Double
fDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
delta