{-# LANGUAGE ParallelListComp #-}
{-# LANGUAGE CPP #-}
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
{-# LANGUAGE Safe #-}
#endif
module Math.ContinuedFraction
    ( CF
    , cf, gcf
    , asCF, asGCF
    
    , truncateCF
    
    , equiv
    , setNumerators
    , setDenominators
    
    , partitionCF
    , evenCF
    , oddCF
    
    , convergents
    , steed
    , lentz, lentzWith
    , modifiedLentz, modifiedLentzWith
    
    , sumPartialProducts
    ) where

import Control.Arrow ((***))
import Data.List (tails, mapAccumL)

-- * The 'CF' type and basic operations

-- I think I would like to try refactoring this stuff at some point to use
-- an "Inductive" CF type, something like:
-- 
-- > data CF a
-- >     = CFZero               -- eval CFZero          = 0
-- >     | CFAdd    a (CF a)    -- eval (CFAdd    b x) =      b + eval x
-- >     | CFCont a a (CF a)    -- eval (CFCont a b x) = a / (b + eval x)
-- 
-- Or perhaps Bill Gosper's "∞-centered" representation:
-- 
-- > data CF a
-- >     = CFInfinity           -- eval CFInfinity     = ∞
-- >     | CFCont a a (CF a)    -- eval (CFCont p q x) = p + q / eval x
-- 

-- |A continued fraction.  Constructed by 'cf' or 'gcf'.
data CF a 
    = CF a [a]
    -- ^ Not exported. See 'cf', the public constructor.
    | GCF a [(a,a)]
    -- ^ Not exported. See 'gcf', the public constructor.

-- |Construct a continued fraction from its first term and the 
-- partial denominators in its canonical form, which is the form 
-- where all the partial numerators are 1.
-- 
-- @cf a [b,c,d]@ corresponds to @a + (b \/ (1 + (c \/ (1 + d))))@,
-- or to @GCF a [(1,b),(1,c),(1,d)]@.
cf :: a -> [a] -> CF a
cf :: a -> [a] -> CF a
cf = a -> [a] -> CF a
forall a. a -> [a] -> CF a
CF

-- |Construct a continued fraction from its first term, its partial 
-- numerators and its partial denominators.
--
-- @gcf b0 [(a1,b1), (a2,b2), (a3,b3)]@ corresponds to
-- @b0 + (a1 \/ (b1 + (a2 \/ (b2 + (a3 \/ b3)))))@
gcf :: a -> [(a,a)] -> CF a
gcf :: a -> [(a, a)] -> CF a
gcf = a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
GCF

instance Show a => Show (CF a) where
    showsPrec :: Int -> CF a -> ShowS
showsPrec Int
p (CF a
b0 [a]
ab) = Bool -> ShowS -> ShowS
showParen (Int
pInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
10)
        ( String -> ShowS
showString String
"cf "
        ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> a -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
11 a
b0
        ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Char -> ShowS
showChar Char
' '
        ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [a] -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
11 [a]
ab
        )
    showsPrec Int
p (GCF a
b0 [(a, a)]
ab) = Bool -> ShowS -> ShowS
showParen (Int
pInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>Int
10)
        ( String -> ShowS
showString String
"gcf "
        ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> a -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
11 a
b0
        ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Char -> ShowS
showChar Char
' '
        ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [(a, a)] -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
11 [(a, a)]
ab
        )

instance Functor CF where
    fmap :: (a -> b) -> CF a -> CF b
fmap a -> b
f (CF  a
b0 [a]
cf)  = b -> [b] -> CF b
forall a. a -> [a] -> CF a
CF  (a -> b
f a
b0) ((a -> b) -> [a] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map a -> b
f [a]
cf)
    fmap a -> b
f (GCF a
b0 [(a, a)]
gcf) = b -> [(b, b)] -> CF b
forall a. a -> [(a, a)] -> CF a
GCF (a -> b
f a
b0) (((a, a) -> (b, b)) -> [(a, a)] -> [(b, b)]
forall a b. (a -> b) -> [a] -> [b]
map (a -> b
f (a -> b) -> (a -> b) -> (a, a) -> (b, b)
forall (a :: * -> * -> *) b c b' c'.
Arrow a =>
a b c -> a b' c' -> a (b, b') (c, c')
*** a -> b
f) [(a, a)]
gcf)

-- |Extract the partial denominators of a 'CF', normalizing it if necessary so 
-- that all the partial numerators are 1.
asCF  :: Fractional a => CF a -> (a, [a])
asCF :: CF a -> (a, [a])
asCF (CF  a
b0 [a]
cf) = (a
b0, [a]
cf)
asCF (GCF a
b0 []) = (a
b0, [])
asCF (GCF a
b0 [(a, a)]
cf) = (a
b0, (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
bs [a]
cs)
    where
        (a
a:[a]
as, [a]
bs) = [(a, a)] -> ([a], [a])
forall a b. [(a, b)] -> ([a], [b])
unzip [(a, a)]
cf
        cs :: [a]
cs = a -> a
forall a. Fractional a => a -> a
recip a
a a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a -> a
forall a. Fractional a => a -> a
recip (a
aa -> a -> a
forall a. Num a => a -> a -> a
*a
c) | a
c <- [a]
cs | a
a <- [a]
as]

-- |Extract all the partial numerators and partial denominators of a 'CF'.
asGCF :: (Num a, Eq a) => CF a -> (a,[(a,a)])
asGCF :: CF a -> (a, [(a, a)])
asGCF (CF  a
b0  [a]
cf) = (a
b0, [(a
1, a
b) | a
b <- [a]
cf])
asGCF (GCF a
b0 [(a, a)]
gcf) = (a
b0, ((a, a) -> Bool) -> [(a, a)] -> [(a, a)]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile ((a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/=a
0)(a -> Bool) -> ((a, a) -> a) -> (a, a) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(a, a) -> a
forall a b. (a, b) -> a
fst) [(a, a)]
gcf)

-- |Truncate a 'CF' to the specified number of partial numerators and denominators.
truncateCF :: Int -> CF a -> CF a
truncateCF :: Int -> CF a -> CF a
truncateCF Int
n (CF  a
b0 [a]
ab) = a -> [a] -> CF a
forall a. a -> [a] -> CF a
CF  a
b0 (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
n [a]
ab)
truncateCF Int
n (GCF a
b0 [(a, a)]
ab) = a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
GCF a
b0 (Int -> [(a, a)] -> [(a, a)]
forall a. Int -> [a] -> [a]
take Int
n [(a, a)]
ab)

-- |Apply an equivalence transformation, multiplying each partial denominator 
-- with the corresponding element of the supplied list and transforming 
-- subsequent partial numerators and denominators as necessary.  If the list
-- is too short, the rest of the 'CF' will be unscaled.
equiv :: (Num a, Eq a) => [a] -> CF a -> CF a
equiv :: [a] -> CF a -> CF a
equiv [a]
cs CF a
orig
    = a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
gcf a
b0 ([a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
as' [a]
bs')
    where
        (a
b0, [(a, a)]
terms) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
        ([a]
as,[a]
bs) = [(a, a)] -> ([a], [a])
forall a b. [(a, b)] -> ([a], [b])
unzip [(a, a)]
terms
        
        as' :: [a]
as' = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
cs' (a
1a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
cs')) [a]
as
        bs' :: [a]
bs' = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
cs' [a]
bs
        cs' :: [a]
cs' = [a]
cs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a -> [a]
forall a. a -> [a]
repeat a
1

-- |Apply an equivalence transformation that sets the partial denominators 
-- of a 'CF' to the specfied values.  If the input list is too short, the 
-- rest of the 'CF' will be unscaled.
setDenominators :: (Fractional a, Eq a) => [a] -> CF a -> CF a
setDenominators :: [a] -> CF a -> CF a
setDenominators [a]
denoms CF a
orig
    = a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
gcf a
b0 ([a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
as' [a]
bs')
    where
        (a
b0, [(a, a)]
terms) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
        ([a]
as,[a]
bs) = [(a, a)] -> ([a], [a])
forall a b. [(a, b)] -> ([a], [b])
unzip [(a, a)]
terms
        
        as' :: [a]
as' = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
as ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
cs (a
1a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
cs))
        bs' :: [a]
bs' = ((a -> a) -> a -> a) -> [a -> a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (a -> a) -> a -> a
forall a b. (a -> b) -> a -> b
($) ((a -> a -> a) -> [a] -> [a -> a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a -> a
forall a b. a -> b -> a
const [a]
denoms [a -> a] -> [a -> a] -> [a -> a]
forall a. [a] -> [a] -> [a]
++ (a -> a) -> [a -> a]
forall a. a -> [a]
repeat a -> a
forall a. a -> a
id) [a]
bs
        cs :: [a]
cs = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/) [a]
bs' [a]
bs

-- |Apply an equivalence transformation that sets the partial numerators 
-- of a 'CF' to the specfied values.  If the input list is too short, the 
-- rest of the 'CF' will be unscaled.
setNumerators :: (Fractional a, Eq a) => [a] -> CF a -> CF a
setNumerators :: [a] -> CF a -> CF a
setNumerators [a]
numers CF a
orig
    = a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
gcf a
b0 ([a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
as' [a]
bs')
    where
        (a
b0, [(a, a)]
terms) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
        ([a]
as,[a]
bs) = [(a, a)] -> ([a], [a])
forall a b. [(a, b)] -> ([a], [b])
unzip [(a, a)]
terms
        
        as' :: [a]
as' = ((a -> a) -> a -> a) -> [a -> a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (a -> a) -> a -> a
forall a b. (a -> b) -> a -> b
($) ((a -> a -> a) -> [a] -> [a -> a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a -> a
forall a b. a -> b -> a
const [a]
numers [a -> a] -> [a -> a] -> [a -> a]
forall a. [a] -> [a] -> [a]
++ (a -> a) -> [a -> a]
forall a. a -> [a]
repeat a -> a
forall a. a -> a
id) [a]
as
        bs' :: [a]
bs' = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
bs [a]
cs
        cs :: [a]
cs = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/) [a]
as' ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
as (a
1a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
cs))

-- |Computes the even and odd parts, respectively, of a 'CF'.  These are new
-- 'CF's that have the even-indexed and odd-indexed convergents of the 
-- original, respectively.
partitionCF :: (Fractional a, Eq a) => CF a -> (CF a, CF a)
partitionCF :: CF a -> (CF a, CF a)
partitionCF CF a
orig = case [(a, a)]
terms of
    []          -> (CF a
orig, CF a
orig)
    [(a
a1,a
b1)]   -> 
        let final :: CF a
final = a -> [a] -> CF a
forall a. a -> [a] -> CF a
cf (a
b0 a -> a -> a
forall a. Num a => a -> a -> a
+ a
a1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
b1) []
         in (CF a
final, CF a
final)
    [(a, a)]
_           -> (CF a
evenPart, CF a
oddPart)
    where
        (a
b0, [(a, a)]
terms) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
        ([a]
as, [a]
bs)    = [(a, a)] -> ([a], [a])
forall a b. [(a, b)] -> ([a], [b])
unzip [(a, a)]
terms
        
        pairs :: [b] -> [(b, b)]
pairs (b
a:b
b:[b]
rest) = (b
a,b
b) (b, b) -> [(b, b)] -> [(b, b)]
forall a. a -> [a] -> [a]
: [b] -> [(b, b)]
pairs [b]
rest
        pairs [b
a] = [(b
a,b
0)]
        pairs [b]
_ = []
        
        alphas :: [a]
alphas@(a
alpha1:a
alpha2:[a]
_) = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/) [a]
as ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
bs (a
1a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
bs))
        
        evenPart :: CF a
evenPart = a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
gcf a
b0 ([a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
cs [a]
ds)
            where
                cs :: [a]
cs =     a
alpha1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [(-a
aOdd) a -> a -> a
forall a. Num a => a -> a -> a
* a
aEven  | (a
aEven, a
aOdd)  <- [a] -> [(a, a)]
forall b. Num b => [b] -> [(b, b)]
pairs ([a] -> [a]
forall a. [a] -> [a]
tail [a]
alphas)]
                ds :: [a]
ds = a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
alpha2 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
aOdd a -> a -> a
forall a. Num a => a -> a -> a
+ a
aEven | (a
aOdd,  a
aEven) <- [(a, a)] -> [(a, a)]
forall a. [a] -> [a]
tail ([a] -> [(a, a)]
forall b. Num b => [b] -> [(b, b)]
pairs [a]
alphas)]
        
        oddPart :: CF a
oddPart = a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
gcf (a
b0a -> a -> a
forall a. Num a => a -> a -> a
+a
alpha1) ([a] -> [a] -> [(a, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
cs [a]
ds)
            where
                cs :: [a]
cs = [(-a
aOdd) a -> a -> a
forall a. Num a => a -> a -> a
* a
aEven  | (a
aOdd, a
aEven) <- [a] -> [(a, a)]
forall b. Num b => [b] -> [(b, b)]
pairs [a]
alphas]
                ds :: [a]
ds = [a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
aOdd a -> a -> a
forall a. Num a => a -> a -> a
+ a
aEven | (a
aEven, a
aOdd) <- [a] -> [(a, a)]
forall b. Num b => [b] -> [(b, b)]
pairs ([a] -> [a]
forall a. [a] -> [a]
tail [a]
alphas)]

-- |Computes the even part of a 'CF' (that is, a new 'CF' whose convergents are
-- the even-indexed convergents of the original).
evenCF :: (Fractional a, Eq a) => CF a -> CF a
evenCF :: CF a -> CF a
evenCF = (CF a, CF a) -> CF a
forall a b. (a, b) -> a
fst ((CF a, CF a) -> CF a) -> (CF a -> (CF a, CF a)) -> CF a -> CF a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CF a -> (CF a, CF a)
forall a. (Fractional a, Eq a) => CF a -> (CF a, CF a)
partitionCF

-- |Computes the odd part of a 'CF' (that is, a new 'CF' whose convergents are
-- the odd-indexed convergents of the original).
oddCF :: (Fractional a, Eq a) => CF a -> CF a
oddCF :: CF a -> CF a
oddCF = (CF a, CF a) -> CF a
forall a b. (a, b) -> b
snd ((CF a, CF a) -> CF a) -> (CF a -> (CF a, CF a)) -> CF a -> CF a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CF a -> (CF a, CF a)
forall a. (Fractional a, Eq a) => CF a -> (CF a, CF a)
partitionCF


-- * Evaluating continued fractions

-- |Evaluate the convergents of a continued fraction using the fundamental
-- recurrence formula:
-- 
-- A0 = b0, B0 = 1
--
-- A1 = b1b0 + a1,  B1 = b1
-- 
-- A{n+1} = b{n+1}An + a{n+1}A{n-1}
--
-- B{n+1} = b{n+1}Bn + a{n+1}B{n-1}
--
-- The convergents are then Xn = An/Bn
convergents :: (Fractional a, Eq a) => CF a -> [a]
convergents :: CF a -> [a]
convergents CF a
orig = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
1 ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/) [a]
nums [a]
denoms)
    where
        (a
b0, [(a, a)]
terms) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
        nums :: [a]
nums   = a
1a -> [a] -> [a]
forall a. a -> [a] -> [a]
:a
b0a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
x1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
x0 | a
x0:a
x1:[a]
_ <- [a] -> [[a]]
forall a. [a] -> [[a]]
tails [a]
nums   | (a
a,a
b) <- [(a, a)]
terms]
        denoms :: [a]
denoms = a
0a -> [a] -> [a]
forall a. a -> [a] -> [a]
:a
1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
x1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
x0 | a
x0:a
x1:[a]
_ <- [a] -> [[a]]
forall a. [a] -> [[a]]
tails [a]
denoms | (a
a,a
b) <- [(a, a)]
terms]

-- |Evaluate the convergents of a continued fraction using Steed's method.
-- Only valid if the denominator in the following recurrence for D_i never 
-- goes to zero.  If this method blows up, try 'modifiedLentz'.
--
-- D1 = 1/b1
-- 
-- D{i} = 1 / (b{i} + a{i} * D{i-1})
-- 
-- dx1 = a1 / b1
-- 
-- dx{i} = (b{i} * D{i} - 1) * dx{i-1}
-- 
-- x0 = b0
-- 
-- x{i} = x{i-1} + dx{i}
-- 
-- The convergents are given by @scanl (+) b0 dxs@
steed :: (Fractional a, Eq a) => CF a -> [a]
steed :: CF a -> [a]
steed (CF  a
b0 []) = [a
b0]
steed (GCF a
b0 []) = [a
b0]
steed (CF  a
0 (  a
a  :[a]
rest)) = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/) (CF a -> [a]
forall a. (Fractional a, Eq a) => CF a -> [a]
steed (a -> [a] -> CF a
forall a. a -> [a] -> CF a
CF  a
a [a]
rest))
steed (GCF a
0 ((a
a,a
b):[(a, a)]
rest)) = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
a a -> a -> a
forall a. Fractional a => a -> a -> a
/) (CF a -> [a]
forall a. (Fractional a, Eq a) => CF a -> [a]
steed (a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
GCF a
b [(a, a)]
rest))
steed CF a
orig
    = (a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl a -> a -> a
forall a. Num a => a -> a -> a
(+) a
b0 [a]
dxs
    where
        (a
b0, (a
a1,a
b1):[(a, a)]
gcf) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
        
        dxs :: [a]
dxs = a
a1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
b1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [(a
b a -> a -> a
forall a. Num a => a -> a -> a
* a
d a -> a -> a
forall a. Num a => a -> a -> a
- a
1) a -> a -> a
forall a. Num a => a -> a -> a
* a
dx  | (a
a,a
b) <- [(a, a)]
gcf | a
d <- [a]
ds | a
dx <- [a]
dxs]
        ds :: [a]
ds  =  a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
b1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a -> a
forall a. Fractional a => a -> a
recip (a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
a a -> a -> a
forall a. Num a => a -> a -> a
* a
d) | (a
a,a
b) <- [(a, a)]
gcf | a
d <- [a]
ds]

-- |Evaluate the convergents of a continued fraction using Lentz's method.
-- Only valid if the denominators in the following recurrence never go to
-- zero.  If this method blows up, try 'modifiedLentz'.
--
-- C1 = b1 + a1 / b0
--
-- D1 = 1/b1
-- 
-- C{n} = b{n} + a{n} / C{n-1}
-- 
-- D{n} = 1 / (b{n} + a{n} * D{n-1})
-- 
-- The convergents are given by @scanl (*) b0 (zipWith (*) cs ds)@
lentz :: (Fractional a, Eq a) => CF a -> [a]
lentz :: CF a -> [a]
lentz = (a -> a) -> (a -> a -> a) -> (a -> a) -> CF a -> [a]
forall a b.
(Fractional a, Eq a) =>
(a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b]
lentzWith a -> a
forall a. a -> a
id a -> a -> a
forall a. Num a => a -> a -> a
(*) a -> a
forall a. Fractional a => a -> a
recip

-- |Evaluate the convergents of a continued fraction using Lentz's method,
-- mapping the terms in the final product to a new group before performing
-- the final multiplications.  A useful group, for example, would be logarithms
-- under addition.  In @lentzWith f op inv@, the arguments are:
-- 
-- * @f@, a group homomorphism (eg, 'log') from {@a@,(*),'recip'} to the group
--   in which you want to perform the multiplications.
-- 
-- * @op@, the group operation (eg., (+)).
-- 
-- * @inv@, the group inverse (eg., 'negate').
-- 
-- The 'lentz' function, for example, is given by the identity homomorphism:
-- @lentz@ = @lentzWith id (*) recip@.
-- 
-- The original motivation for this function is to allow computation of 
-- the natural log of very large numbers that would overflow with the naive
-- implementation in 'lentz'.  In this case, the arguments would be 'log', (+),
-- and 'negate', respectively.
-- 
-- In cases where terms of the product can be negative (i.e., the sequence of
-- convergents contains negative values), the following definitions could 
-- be used instead:
-- 
-- > signLog x = (signum x, log (abs x))
-- > addSignLog (xS,xL) (yS,yL) = (xS*yS, xL+yL)
-- > negateSignLog (s,l) = (s, negate l)
{-# INLINE lentzWith #-}
lentzWith :: (Fractional a, Eq a) => (a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b]
lentzWith :: (a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b]
lentzWith a -> b
f b -> b -> b
op b -> b
inv (CF  a
0 (  a
a  :[a]
rest)) = (b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map b -> b
inv              ((a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b]
forall a b.
(Fractional a, Eq a) =>
(a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b]
lentzWith a -> b
f b -> b -> b
op b -> b
inv (a -> [a] -> CF a
forall a. a -> [a] -> CF a
CF  a
a [a]
rest))
lentzWith a -> b
f b -> b -> b
op b -> b
inv (GCF a
0 ((a
a,a
b):[(a, a)]
rest)) = (b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map (b -> b -> b
op (a -> b
f a
a) (b -> b) -> (b -> b) -> b -> b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. b -> b
inv) ((a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b]
forall a b.
(Fractional a, Eq a) =>
(a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b]
lentzWith a -> b
f b -> b -> b
op b -> b
inv (a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
GCF a
b [(a, a)]
rest))
lentzWith a -> b
f b -> b -> b
op b -> b
inv CF a
c = (b -> a -> b) -> b -> [a] -> [b]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl b -> a -> b
opF (a -> b
f a
b0) ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
cs [a]
ds)
   where
       opF :: b -> a -> b
opF b
x a
y = b -> b -> b
op b
x (a -> b
f a
y)
       (a
b0, [a]
cs, [a]
ds) = CF a -> (a, [a], [a])
forall a. (Fractional a, Eq a) => CF a -> (a, [a], [a])
lentzRecurrence CF a
c


-- precondition: b0 /= 0
lentzRecurrence :: (Fractional a, Eq a) => CF a -> (a,[a],[a])
lentzRecurrence :: CF a -> (a, [a], [a])
lentzRecurrence CF a
orig 
    | [(a, a)] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(a, a)]
terms    = (a
b0,[],[])
    | Bool
otherwise = (a
b0, [a]
cs, [a]
ds)
    where
        (a
b0, [(a, a)]
terms) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
        
        cs :: [a]
cs = [   a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
aa -> a -> a
forall a. Fractional a => a -> a -> a
/a
c  | (a
a,a
b) <- [(a, a)]
terms | a
c <- a
b0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
cs]
        ds :: [a]
ds = [a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/(a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
aa -> a -> a
forall a. Num a => a -> a -> a
*a
d) | (a
a,a
b) <- [(a, a)]
terms | a
d <- a
0  a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
ds]


-- |Evaluate the convergents of a continued fraction using Lentz's method,
-- (see 'lentz') with the additional rule that if a denominator ever goes
-- to zero, it will be replaced by a (very small) number of your choosing,
-- typically 1e-30 or so (this modification was proposed by Thompson and 
-- Barnett).  
-- 
-- Additionally splits the resulting list of convergents into sublists, 
-- starting a new list every time the \'modification\' is invoked.  
modifiedLentz :: (Fractional a, Eq a) => a -> CF a -> [[a]]
modifiedLentz :: a -> CF a -> [[a]]
modifiedLentz = (a -> a) -> (a -> a -> a) -> (a -> a) -> a -> CF a -> [[a]]
forall a b.
(Fractional a, Eq a) =>
(a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]]
modifiedLentzWith a -> a
forall a. a -> a
id a -> a -> a
forall a. Num a => a -> a -> a
(*) a -> a
forall a. Fractional a => a -> a
recip

-- |'modifiedLentz' with a group homomorphism (see 'lentzWith', it bears the
-- same relationship to 'lentz' as this function does to 'modifiedLentz',
-- and solves the same problems).  Alternatively, 'lentzWith' with the same
-- modification to the recurrence as 'modifiedLentz'.
{-# INLINE modifiedLentzWith #-}
modifiedLentzWith :: (Fractional a, Eq a) => (a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]]
modifiedLentzWith :: (a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]]
modifiedLentzWith a -> b
f b -> b -> b
op b -> b
inv a
z (CF  a
0 (  a
a  :[a]
rest)) = ([b] -> [b]) -> [[b]] -> [[b]]
forall a b. (a -> b) -> [a] -> [b]
map ((b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map             b -> b
inv ) ((a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]]
forall a b.
(Fractional a, Eq a) =>
(a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]]
modifiedLentzWith a -> b
f b -> b -> b
op b -> b
inv a
z (a -> [a] -> CF a
forall a. a -> [a] -> CF a
CF  a
a [a]
rest))
modifiedLentzWith a -> b
f b -> b -> b
op b -> b
inv a
z (GCF a
0 ((a
a,a
b):[(a, a)]
rest)) = ([b] -> [b]) -> [[b]] -> [[b]]
forall a b. (a -> b) -> [a] -> [b]
map ((b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map (b -> b -> b
op (a -> b
f a
a) (b -> b) -> (b -> b) -> b -> b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. b -> b
inv)) ((a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]]
forall a b.
(Fractional a, Eq a) =>
(a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]]
modifiedLentzWith a -> b
f b -> b -> b
op b -> b
inv a
z (a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
GCF a
b [(a, a)]
rest))
modifiedLentzWith a -> b
f b -> b -> b
op b -> b
inv a
z CF a
orig = [(Bool, b)] -> [[b]]
forall b. [(Bool, b)] -> [[b]]
separate (((Bool, b) -> (Bool, a) -> (Bool, b))
-> (Bool, b) -> [(Bool, a)] -> [(Bool, b)]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (Bool, b) -> (Bool, a) -> (Bool, b)
opF (Bool
False, a -> b
f a
b0) [(Bool, a)]
cds)
    where
        (a
b0, [(Bool, a)]
cs, [(Bool, a)]
ds) = a -> CF a -> (a, [(Bool, a)], [(Bool, a)])
forall a.
(Fractional a, Eq a) =>
a -> CF a -> (a, [(Bool, a)], [(Bool, a)])
modifiedLentzRecurrence a
z CF a
orig
        cds :: [(Bool, a)]
cds = ((Bool, a) -> (Bool, a) -> (Bool, a))
-> [(Bool, a)] -> [(Bool, a)] -> [(Bool, a)]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (Bool, a) -> (Bool, a) -> (Bool, a)
forall b. Num b => (Bool, b) -> (Bool, b) -> (Bool, b)
mult [(Bool, a)]
cs [(Bool, a)]
ds
        
        mult :: (Bool, b) -> (Bool, b) -> (Bool, b)
mult (Bool
xa,b
xb) (Bool
ya,b
yb) = (Bool
xa Bool -> Bool -> Bool
|| Bool
ya, b
xb b -> b -> b
forall a. Num a => a -> a -> a
* b
yb)
        opF :: (Bool, b) -> (Bool, a) -> (Bool, b)
opF  (Bool
xa,b
xb) (Bool
ya,a
yb) = (Bool
xa Bool -> Bool -> Bool
|| Bool
ya, b -> b -> b
op b
xb (a -> b
f a
yb))
        
        -- |Takes a list of (Bool,a) and breaks it into sublists, starting
        -- a new one every time it encounters (True,_).
        separate :: [(Bool, b)] -> [[b]]
separate [] = []
        separate ((Bool
_,b
x):[(Bool, b)]
xs) = case ((Bool, b) -> Bool) -> [(Bool, b)] -> ([(Bool, b)], [(Bool, b)])
forall a. (a -> Bool) -> [a] -> ([a], [a])
break (Bool, b) -> Bool
forall a b. (a, b) -> a
fst [(Bool, b)]
xs of
            ([(Bool, b)]
xs, [(Bool, b)]
ys) -> (b
xb -> [b] -> [b]
forall a. a -> [a] -> [a]
:((Bool, b) -> b) -> [(Bool, b)] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map (Bool, b) -> b
forall a b. (a, b) -> b
snd [(Bool, b)]
xs) [b] -> [[b]] -> [[b]]
forall a. a -> [a] -> [a]
: [(Bool, b)] -> [[b]]
separate [(Bool, b)]
ys

-- precondition: b0 /= 0
modifiedLentzRecurrence :: (Fractional a, Eq a) => a -> CF a -> (a,[(Bool, a)],[(Bool, a)])
modifiedLentzRecurrence :: a -> CF a -> (a, [(Bool, a)], [(Bool, a)])
modifiedLentzRecurrence a
z CF a
orig
    | [(a, a)] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(a, a)]
terms = (a
b0, [], [])
    | Bool
otherwise  = (a
b0, [(Bool, a)]
cs, [(Bool, a)]
ds)
    where
        (a
b0, [(a, a)]
terms) = CF a -> (a, [(a, a)])
forall a. (Num a, Eq a) => CF a -> (a, [(a, a)])
asGCF CF a
orig
        
        cs :: [(Bool, a)]
cs = [a -> (a -> a) -> (Bool, a)
forall b. a -> (a -> b) -> (Bool, b)
reset (a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
aa -> a -> a
forall a. Fractional a => a -> a -> a
/a
c)    a -> a
forall a. a -> a
id | (a
a,a
b) <- [(a, a)]
terms | a
c <- a
b0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: ((Bool, a) -> a) -> [(Bool, a)] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (Bool, a) -> a
forall a b. (a, b) -> b
snd [(Bool, a)]
cs]
        ds :: [(Bool, a)]
ds = [a -> (a -> a) -> (Bool, a)
forall b. a -> (a -> b) -> (Bool, b)
reset (a
b a -> a -> a
forall a. Num a => a -> a -> a
+ a
aa -> a -> a
forall a. Num a => a -> a -> a
*a
d) a -> a
forall a. Fractional a => a -> a
recip | (a
a,a
b) <- [(a, a)]
terms | a
d <- a
0  a -> [a] -> [a]
forall a. a -> [a] -> [a]
: ((Bool, a) -> a) -> [(Bool, a)] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (Bool, a) -> a
forall a b. (a, b) -> b
snd [(Bool, a)]
ds]
        
        -- The sublist breaking is computed secondarily - initially, 
        -- 'cs' and 'ds' are constructed with this helper function that
        -- adds a marker to the list whenever a term of interest goes to 0,
        -- while also resetting that term to a small nonzero amount.
        -- Then later, 'separate' breaks the list every time it sees one
        -- of these markers.
        reset :: a -> (a -> b) -> (Bool, b)
reset a
x a -> b
f
            | a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0    = (Bool
True,  a -> b
f a
z)
            | Bool
otherwise = (Bool
False, a -> b
f a
x)


-- |Euler's formula for computing @sum (scanl1 (*) xs)@.  
-- Successive convergents of the resulting 'CF' are successive partial sums
-- in the series.
sumPartialProducts :: Num a => [a] -> CF a
sumPartialProducts :: [a] -> CF a
sumPartialProducts [] = a -> [a] -> CF a
forall a. a -> [a] -> CF a
cf a
0 []
sumPartialProducts (a
x:[a]
xs) = a -> [(a, a)] -> CF a
forall a. a -> [(a, a)] -> CF a
gcf a
0 ((a
x,a
1)(a, a) -> [(a, a)] -> [(a, a)]
forall a. a -> [a] -> [a]
:[(a -> a
forall a. Num a => a -> a
negate a
x, a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
x) | a
x <- [a]
xs])