diff --git a/continued-fractions.cabal b/continued-fractions.cabal --- a/continued-fractions.cabal +++ b/continued-fractions.cabal @@ -1,5 +1,5 @@ name: continued-fractions -version: 0.9.0.1 +version: 0.9.1.0 stability: provisional cabal-version: >= 1.6 diff --git a/src/Math/ContinuedFraction.hs b/src/Math/ContinuedFraction.hs --- a/src/Math/ContinuedFraction.hs +++ b/src/Math/ContinuedFraction.hs @@ -16,8 +16,8 @@ , convergents , steed - , lentz - , modifiedLentz + , lentz, lentzWith + , modifiedLentz, modifiedLentzWith , sumPartialProducts ) where @@ -246,17 +246,54 @@ -- -- The convergents are given by @scanl (*) b0 (zipWith (*) cs 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) +lentz = lentzWith id (*) 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) = (negate s, l) +{-# INLINE lentzWith #-} +lentzWith :: Fractional a => (a -> b) -> (b -> b -> b) -> (b -> b) -> CF a -> [b] +lentzWith f op inv (CF 0 ( a :rest)) = map inv (lentzWith f op inv (CF a rest)) +lentzWith f op inv (GCF 0 ((a,b):rest)) = map (op (f a) . inv) (lentzWith f op inv (GCF b rest)) +lentzWith f op inv c = scanl opF (f b0) (zipWith (*) cs ds) + where + opF x y = op x (f y) + (b0, cs, ds) = lentzRecurrence c + + +lentzRecurrence :: Fractional a => CF a -> (a,[a],[a]) +lentzRecurrence orig + | null terms = (b0,[],[]) + | otherwise = (b0, cs, ds) where - (b0, gcf) = asGCF orig + (b0, terms) = asGCF orig - cs = [ b + a/c | (a,b) <- gcf | c <- b0 : cs] - ds = [1/(b + a*d) | (a,b) <- gcf | d <- 0 : ds] + cs = [ b + a/c | (a,b) <- terms | c <- b0 : cs] + ds = [1/(b + a*d) | (a,b) <- terms | d <- 0 : ds] -- |Evaluate the convergents of a continued fraction using Lentz's method, @@ -268,19 +305,37 @@ -- Additionally splits the resulting list of convergents into sublists, -- starting a new list every time the \'modification\' is invoked. 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)) +modifiedLentz = modifiedLentzWith id (*) 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 => (a -> b) -> (b -> b -> b) -> (b -> b) -> a -> CF a -> [[b]] +modifiedLentzWith f op inv z (CF 0 ( a :rest)) = map (map inv ) (modifiedLentzWith f op inv z (CF a rest)) +modifiedLentzWith f op inv z (GCF 0 ((a,b):rest)) = map (map (op (f a) . inv)) (modifiedLentzWith f op inv z (GCF b rest)) +modifiedLentzWith f op inv z orig = separate (scanl opF (False, f b0) cds) where + (b0, cs, ds) = modifiedLentzRecurrence z orig + cds = zipWith mult cs ds + + mult (xa,xb) (ya,yb) = (xa || ya, xb * yb) + opF (xa,xb) (ya,yb) = (xa || ya, op xb (f yb)) + + -- |Takes a list of (Bool,a) and breaks it into sublists, starting + -- a new one every time it encounters (True,_). + separate [] = [] + separate ((_,x):xs) = case break fst xs of + (xs, ys) -> (x:map snd xs) : separate ys + +modifiedLentzRecurrence :: Fractional a => a -> CF a -> (a,[(Bool, a)],[(Bool, a)]) +modifiedLentzRecurrence z orig + | null terms = (b0, [], []) + | otherwise = (b0, cs, ds) + 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] @@ -293,18 +348,9 @@ reset x f | x == 0 = (True, f z) | otherwise = (False, f x) - - -- |Takes a list of (Bool,a) and breaks it into sublists, starting - -- a new one every time it encounters (True,_). - 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 --- |Euler's formula for computing @sum (map product (tail (inits xs)))@. + +-- |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