continued-fractions (empty) → 0.9
raw patch · 3 files changed
+344/−0 lines, 3 filesdep +basesetup-changed
Dependencies added: base
Files
- Setup.lhs +5/−0
- continued-fractions.cabal +29/−0
- src/Math/ContinuedFraction.hs +310/−0
+ Setup.lhs view
@@ -0,0 +1,5 @@+#!/usr/bin/env runhaskell++> import Distribution.Simple+> main = defaultMain+
+ continued-fractions.cabal view
@@ -0,0 +1,29 @@+name: continued-fractions+version: 0.9+stability: provisional++cabal-version: >= 1.6+build-type: Simple++author: James Cook <mokus@deepbondi.net>+maintainer: James Cook <mokus@deepbondi.net>+license: PublicDomain+homepage: /dev/null++category: Math, Numerical+synopsis: Continued fractions.+description: A type and some functions for manipulating and + evaluating continued fractions.++tested-with: GHC == 6.8.3,+ GHC == 6.10.4,+ GHC == 6.12.1, GHC == 6.12.3++source-repository head+ type: darcs+ location: http://code.haskell.org/~mokus/continued-fractions++Library+ hs-source-dirs: src+ exposed-modules: Math.ContinuedFraction+ build-depends: base >= 3 && <5
+ src/Math/ContinuedFraction.hs view
@@ -0,0 +1,310 @@+{-# LANGUAGE ParallelListComp #-}+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)++-- * 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)++-- |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 = 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 = 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)++-- |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 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]++-- |Extract all the partial numerators and partial denominators of a 'CF'.+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)++-- |Truncate a 'CF' to the specified number of partial numerators and denominators.+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)++-- |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 => [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++-- |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 => [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++-- |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 => [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))++-- |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 => 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)]++-- |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 => CF a -> CF a+evenCF = fst . 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 => CF a -> CF a+oddCF = snd . 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 => 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]++-- |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 => 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]++-- |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 => 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]+++-- |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 => a -> CF a -> [[a]]+modifiedLentz z (CF b0 []) = [[b0]]+modifiedLentz z (GCF b0 []) = [[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+ = snd (mapAccumL multSublist b0 (separate cds))+ where+ (b0, gcf) = 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) <- gcf | c <- b0 : map snd cs]+ ds = [reset (b + a*d) recip | (a,b) <- gcf | d <- 0 : map snd 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 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)))@. +-- Successive convergents of the resulting 'CF' are successive partial sums+-- in the series.+sumPartialProducts :: Num a => [a] -> CF a+sumPartialProducts [] = cf 0 []+sumPartialProducts (x:xs) = gcf 0 ((x,1):[(negate x, 1 + x) | x <- xs])