module Math.ContinuedFraction
( CF
, cf, gcf
, asCF, asGCF
, truncateCF
, equiv
, setNumerators
, setDenominators
, partitionCF
, evenCF
, oddCF
, convergents
, steed
, lentz
, modifiedLentz
, sumPartialProducts
) where
import Control.Arrow ((***))
import Data.List (tails, mapAccumL)
data CF a
= CF a [a]
| GCF a [(a,a)]
cf :: a -> [a] -> CF a
cf = CF
gcf :: a -> [(a,a)] -> CF a
gcf = GCF
instance Show a => Show (CF a) where
showsPrec p (CF b0 ab) = showParen (p>10)
( showString "cf "
. showsPrec 11 b0
. showChar ' '
. showsPrec 11 ab
)
showsPrec p (GCF b0 ab) = showParen (p>10)
( showString "gcf "
. showsPrec 11 b0
. showChar ' '
. showsPrec 11 ab
)
instance Functor CF where
fmap f (CF b0 cf) = CF (f b0) (map f cf)
fmap f (GCF b0 gcf) = GCF (f b0) (map (f *** f) gcf)
asCF :: Fractional a => CF a -> (a, [a])
asCF (CF b0 cf) = (b0, cf)
asCF (GCF b0 []) = (b0, [])
asCF (GCF b0 cf) = (b0, zipWith (*) bs cs)
where
(a:as, bs) = unzip cf
cs = recip a : [recip (a*c) | c <- cs | a <- as]
asGCF :: Num a => CF a -> (a,[(a,a)])
asGCF (CF b0 cf) = (b0, [(1, b) | b <- cf])
asGCF (GCF b0 gcf) = (b0, takeWhile ((/=0).fst) gcf)
truncateCF :: Int -> CF a -> CF a
truncateCF n (CF b0 ab) = CF b0 (take n ab)
truncateCF n (GCF b0 ab) = GCF b0 (take n ab)
equiv :: Num a => [a] -> CF a -> CF a
equiv cs orig
= gcf b0 (zip as' bs')
where
(b0, terms) = asGCF orig
(as,bs) = unzip terms
as' = zipWith (*) (zipWith (*) cs' (1:cs')) as
bs' = zipWith (*) cs' bs
cs' = cs ++ repeat 1
setDenominators :: Fractional a => [a] -> CF a -> CF a
setDenominators denoms orig
= gcf b0 (zip as' bs')
where
(b0, terms) = asGCF orig
(as,bs) = unzip terms
as' = zipWith (*) as (zipWith (*) cs (1:cs))
bs' = zipWith ($) (map const denoms ++ repeat id) bs
cs = zipWith (/) bs' bs
setNumerators :: Fractional a => [a] -> CF a -> CF a
setNumerators numers orig
= gcf b0 (zip as' bs')
where
(b0, terms) = asGCF orig
(as,bs) = unzip terms
as' = zipWith ($) (map const numers ++ repeat id) as
bs' = zipWith (*) bs cs
cs = zipWith (/) as' (zipWith (*) as (1:cs))
partitionCF :: Fractional a => CF a -> (CF a, CF a)
partitionCF orig = case terms of
[] -> (orig, orig)
[(a1,b1)] ->
let final = cf (b0 + a1/b1) []
in (final, final)
_ -> (evenPart, oddPart)
where
(b0, terms) = asGCF orig
(as, bs) = unzip terms
pairs (a:b:rest) = (a,b) : pairs rest
pairs [a] = [(a,0)]
pairs _ = []
alphas@(alpha1:alpha2:_) = zipWith (/) as (zipWith (*) bs (1:bs))
evenPart = gcf b0 (zip cs ds)
where
cs = alpha1 : [(aOdd) * aEven | (aEven, aOdd) <- pairs (tail alphas)]
ds = 1 + alpha2 : [1 + aOdd + aEven | (aOdd, aEven) <- tail (pairs alphas)]
oddPart = gcf (b0+alpha1) (zip cs ds)
where
cs = [(aOdd) * aEven | (aOdd, aEven) <- pairs alphas]
ds = [1 + aOdd + aEven | (aEven, aOdd) <- pairs (tail alphas)]
evenCF :: Fractional a => CF a -> CF a
evenCF = fst . partitionCF
oddCF :: Fractional a => CF a -> CF a
oddCF = snd . partitionCF
convergents :: Fractional a => CF a -> [a]
convergents orig = drop 1 (zipWith (/) nums denoms)
where
(b0, terms) = asGCF orig
nums = 1:b0:[b * x1 + a * x0 | x0:x1:_ <- tails nums | (a,b) <- terms]
denoms = 0:1 :[b * x1 + a * x0 | x0:x1:_ <- tails denoms | (a,b) <- terms]
steed :: Fractional a => CF a -> [a]
steed (CF b0 []) = [b0]
steed (GCF b0 []) = [b0]
steed (CF 0 ( a :rest)) = map (1 /) (steed (CF a rest))
steed (GCF 0 ((a,b):rest)) = map (a /) (steed (GCF b rest))
steed orig
= scanl (+) b0 dxs
where
(b0, (a1,b1):gcf) = asGCF orig
dxs = a1/b1 : [(b * d 1) * dx | (a,b) <- gcf | d <- ds | dx <- dxs]
ds = 1/b1 : [recip (b + a * d) | (a,b) <- gcf | d <- ds]
lentz :: Fractional a => CF a -> [a]
lentz (CF b0 []) = [b0]
lentz (GCF b0 []) = [b0]
lentz (CF 0 ( a :rest)) = map (1 /) (lentz (CF a rest))
lentz (GCF 0 ((a,b):rest)) = map (a /) (lentz (GCF b rest))
lentz orig
= scanl (*) b0 (zipWith (*) cs ds)
where
(b0, gcf) = asGCF orig
cs = [ b + a/c | (a,b) <- gcf | c <- b0 : cs]
ds = [1/(b + a*d) | (a,b) <- gcf | d <- 0 : ds]
modifiedLentz :: Fractional a => a -> CF a -> [[a]]
modifiedLentz z (CF b0 []) = [[b0]]
modifiedLentz z (GCF b0 []) = [[b0]]
modifiedLentz z (GCF b0 ((0,_):_)) = [[b0]]
modifiedLentz z (CF 0 ( a :rest)) = map (map (1 /)) (modifiedLentz z (CF a rest))
modifiedLentz z (GCF 0 ((a,b):rest)) = map (map (a /)) (modifiedLentz z (GCF b rest))
modifiedLentz z orig
| null terms = error "programming error in modifiedLentz implementation"
| otherwise = snd (mapAccumL multSublist b0 (separate cds))
where
(b0, terms) = asGCF orig
multSublist b0 cds = let xs = scanl (*) b0 cds in (last xs, xs)
cds = zipWith (\(xa,xb) (ya,yb) -> (xa || ya, xb * yb)) cs ds
cs = [reset (b + a/c) id | (a,b) <- terms | c <- b0 : map snd cs]
ds = [reset (b + a*d) recip | (a,b) <- terms | d <- 0 : map snd ds]
reset x f
| x == 0 = (True, f z)
| otherwise = (False, f x)
separate :: [(Bool,a)] -> [[a]]
separate [] = []
separate xs = case break fst xs of
([], x:xs) -> case separate xs of
[] -> [[snd x]]
(xs:rest) -> (snd x:xs):rest
(xs, ys) -> map snd xs : separate ys
sumPartialProducts :: Num a => [a] -> CF a
sumPartialProducts [] = cf 0 []
sumPartialProducts (x:xs) = gcf 0 ((x,1):[(negate x, 1 + x) | x <- xs])