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