1{-# LANGUAGE BangPatterns              #-}
2{-# LANGUAGE ExistentialQuantification #-}
3-- |
4-- Module    : Numeric.Series
5-- Copyright : (c) 2016 Alexey Khudyakov
6-- License   : BSD3
7--
8-- Maintainer  : alexey.skladnoy@gmail.com, bos@serpentine.com
9-- Stability   : experimental
10-- Portability : portable
11--
12-- Functions for working with infinite sequences. In particular
13-- summation of series and evaluation of continued fractions.
14module Numeric.Series (
15    -- * Data type for infinite sequences.
16    Sequence(..)
17    -- * Constructors
18  , enumSequenceFrom
19  , enumSequenceFromStep
20  , scanSequence
21    -- * Summation of series
22  , sumSeries
23  , sumPowerSeries
24  , sequenceToList
25    -- * Evaluation of continued fractions
26  , evalContFractionB
27  ) where
28
29import Control.Applicative
30import Data.List (unfoldr)
31
32import Numeric.MathFunctions.Constants (m_epsilon)
33
34
35----------------------------------------------------------------
36
37-- | Infinite series. It's represented as opaque state and step
38--   function.
39data Sequence a = forall s. Sequence s (s -> (a,s))
40
41instance Functor Sequence where
42  fmap f (Sequence s0 step) = Sequence s0 (\s -> let (a,s') = step s in (f a, s'))
43  {-# INLINE fmap #-}
44
45instance Applicative Sequence where
46  pure a = Sequence () (\() -> (a,()))
47  Sequence sA fA <*> Sequence sB fB = Sequence (sA,sB) $ \(!sa,!sb) ->
48    let (a,sa') = fA sa
49        (b,sb') = fB sb
50    in (a b, (sa',sb'))
51  {-# INLINE pure  #-}
52  {-# INLINE (<*>) #-}
53
54-- | Elementwise operations with sequences
55instance Num a => Num (Sequence a) where
56  (+) = liftA2 (+)
57  (*) = liftA2 (*)
58  (-) = liftA2 (-)
59  {-# INLINE (+) #-}
60  {-# INLINE (*) #-}
61  {-# INLINE (-) #-}
62  abs         = fmap abs
63  signum      = fmap signum
64  fromInteger = pure . fromInteger
65  {-# INLINE abs         #-}
66  {-# INLINE signum      #-}
67  {-# INLINE fromInteger #-}
68
69-- | Elementwise operations with sequences
70instance Fractional a => Fractional (Sequence a) where
71  (/)          = liftA2 (/)
72  recip        = fmap recip
73  fromRational = pure . fromRational
74  {-# INLINE (/)          #-}
75  {-# INLINE recip        #-}
76  {-# INLINE fromRational #-}
77
78
79
80----------------------------------------------------------------
81-- Constructors
82----------------------------------------------------------------
83
84-- | @enumSequenceFrom x@ generate sequence:
85--
86-- \[ a_n = x + n \]
87enumSequenceFrom :: Num a => a -> Sequence a
88enumSequenceFrom i = Sequence i (\n -> (n,n+1))
89{-# INLINE enumSequenceFrom #-}
90
91-- | @enumSequenceFromStep x d@ generate sequence:
92--
93-- \[ a_n = x + nd \]
94enumSequenceFromStep :: Num a => a -> a -> Sequence a
95enumSequenceFromStep n d = Sequence n (\i -> (i,i+d))
96{-# INLINE enumSequenceFromStep #-}
97
98-- | Analog of 'scanl' for sequence.
99scanSequence :: (b -> a -> b) -> b -> Sequence a -> Sequence b
100{-# INLINE scanSequence #-}
101scanSequence f b0 (Sequence s0 step) = Sequence (b0,s0) $ \(b,s) ->
102  let (a,s') = step s
103      b'     = f b a
104  in (b,(b',s'))
105
106
107----------------------------------------------------------------
108-- Evaluation of series
109----------------------------------------------------------------
110
111-- | Calculate sum of series
112--
113-- \[ \sum_{i=0}^\infty a_i \]
114--
115-- Summation is stopped when
116--
117-- \[ a_{n+1} < \varepsilon\sum_{i=0}^n a_i \]
118--
119-- where ε is machine precision ('m_epsilon')
120sumSeries :: Sequence Double -> Double
121{-# INLINE sumSeries #-}
122sumSeries (Sequence sInit step)
123  = go x0 s0
124  where
125    (x0,s0) = step sInit
126    go x s | abs (d/x) >= m_epsilon = go x' s'
127           | otherwise              = x'
128      where (d,s') = step s
129            x'     = x + d
130
131-- | Calculate sum of series
132--
133-- \[ \sum_{i=0}^\infty x^ia_i \]
134--
135-- Calculation is stopped when next value in series is less than
136-- ε·sum.
137sumPowerSeries :: Double -> Sequence Double -> Double
138sumPowerSeries x ser = sumSeries $ liftA2 (*) (scanSequence (*) 1 (pure x)) ser
139{-# INLINE sumPowerSeries #-}
140
141-- | Convert series to infinite list
142sequenceToList :: Sequence a -> [a]
143sequenceToList (Sequence s f) = unfoldr (Just . f) s
144
145
146
147----------------------------------------------------------------
148-- Evaluation of continued fractions
149----------------------------------------------------------------
150
151-- |
152-- Evaluate continued fraction using modified Lentz algorithm.
153-- Sequence contain pairs (a[i],b[i]) which form following expression:
154--
155-- \[
156-- b_0 + \frac{a_1}{b_1+\frac{a_2}{b_2+\frac{a_3}{b_3 + \cdots}}}
157-- \]
158--
159-- Modified Lentz algorithm is described in Numerical recipes 5.2
160-- "Evaluation of Continued Fractions"
161evalContFractionB :: Sequence (Double,Double) -> Double
162{-# INLINE evalContFractionB #-}
163evalContFractionB (Sequence sInit step)
164  = let ((_,b0),s0) = step sInit
165        f0          = maskZero b0
166    in  go f0 f0 0 s0
167  where
168    tiny = 1e-60
169    maskZero 0 = tiny
170    maskZero x = x
171
172    go f c d s
173      | abs (delta - 1) >= m_epsilon = go f' c' d' s'
174      | otherwise                    = f'
175      where
176          ((a,b),s') = step s
177          d'    = recip $ maskZero $ b + a*d
178          c'    = maskZero $ b + a/c
179          delta = c'*d'
180          f'    = f*delta
181