continued-fractions 0.9.0.1 → 0.9.1.0

raw patch · 2 files changed

+79/−33 lines, 2 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• continued-fractions.cabal +1/−1
• src/Math/ContinuedFraction.hs +78/−32

@@ -1,5 +1,5 @@ name:                   continued-fractions-version:                0.9.0.1+version:                0.9.1.0 stability:              provisional  cabal-version:          >= 1.6

@@ -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