{- | module: Arithmetic.ContinuedFraction description: Continued fractions license: MIT maintainer: Joe Leslie-Hurd stability: provisional portability: portable -} module Arithmetic.ContinuedFraction where import OpenTheory.Primitive.Natural newtype ContinuedFraction = ContinuedFraction {unContinuedFraction :: (Natural,[Natural])} deriving Eq fromNatural :: Natural -> ContinuedFraction fromNatural n = ContinuedFraction (n,[]) goldenRatio :: ContinuedFraction goldenRatio = ContinuedFraction (1, repeat 1) naturalLogarithmBase :: ContinuedFraction naturalLogarithmBase = ContinuedFraction (2, go 2) where go n = 1 : n : 1 : go (n + 2) convergentsFn :: (Natural -> a) -> (a -> a -> a) -> (a -> a -> a) -> [Natural] -> a -> a -> [a] convergentsFn lift add mult = go where go [] _ _ = [] go (q : qs) y x = z : go qs z y where z = add (mult (lift q) y) x numerators :: (Natural -> a) -> (a -> a -> a) -> (a -> a -> a) -> ContinuedFraction -> [a] numerators lift add mult (ContinuedFraction (q0,qs)) = x : convergentsFn lift add mult qs x one where x = lift q0 one = lift 1 denominators :: (Natural -> a) -> (a -> a -> a) -> (a -> a -> a) -> ContinuedFraction -> [a] denominators lift add mult (ContinuedFraction (_,qs)) = one : convergentsFn lift add mult qs one zero where one = lift 1 zero = lift 0 convergents :: (Natural -> a) -> (a -> a -> a) -> (a -> a -> a) -> (a -> a -> a) -> ContinuedFraction -> [a] convergents lift add mult divf cf = zipWith divf nums dens where nums = numerators lift add mult cf dens = denominators lift add mult cf unstableConvergents :: Eq a => [a] -> [a] unstableConvergents [] = error "empty convergents" unstableConvergents (q0 : qs) = q0 : go q0 qs where go _ [] = [] go x (h : t) = if x == h then [] else h : go h t fractionalConvergents :: Fractional a => ContinuedFraction -> [a] fractionalConvergents = convergents fromIntegral (+) (*) (/) rationalConvergents :: ContinuedFraction -> [Rational] rationalConvergents = convergents fromIntegral (+) (*) (/) toDouble :: ContinuedFraction -> Double toDouble = last . unstableConvergents . fractionalConvergents instance Show ContinuedFraction where show = show . toDouble fromRealFrac :: RealFrac a => a -> ContinuedFraction fromRealFrac x = ContinuedFraction (q0, go y) where go s = if s == 0.0 then [] else let (q,t) = properFraction (1.0 / s) in q : go t (q0,y) = properFraction x invert :: ContinuedFraction -> Maybe ContinuedFraction invert (ContinuedFraction (q0,qs)) = if q0 /= 0 then Just (ContinuedFraction (0, q0 : qs)) else case qs of [] -> Nothing h : t -> Just (ContinuedFraction (h,t))