Fun with Stern-Brocot tree

sbphi.hs

import Control.Monad
import Data.Ratio

data FakeRatio = FakeRatio { numr :: Int, denom :: Int} deriving (Show)

data SBTriple = SBTriple { left :: FakeRatio,
                            mid :: FakeRatio,
                            right :: FakeRatio}
                            
data Directions = TLeft | TRight deriving (Show)
                            
data SBTree = SBTree FakeRatio SBTree SBTree

combine f1 f2 = FakeRatio ((numr f1) + (numr f2)) ((denom f1) + (denom f2))

genSB :: SBTriple -> SBTree
genSB SBTriple {left=l, mid=m, right=r} = SBTree m (genSB (SBTriple l (combine m l) m)) (genSB (SBTriple m (combine m r) r))

sbHead (SBTree m _ _) = m
sbLeft (SBTree _ a _) = a
sbRight (SBTree _ _ a) = a

phiDirections = map f [1..] where
    f n = case (odd n) of
            True -> TLeft
            _    -> TRight

eDirections = TRight : join (zipWith conv [1..100] [2,4..]) where
    conv n r = case (odd n) of
        True -> [TRight,TLeft] ++ (take r $ repeat TRight)
        _    -> [TLeft, TRight] ++ (take r $ repeat TLeft)

findSB [] sb = sbHead sb
findSB (d:ds) sb = findSB ds (f sb) where
    f sb = case d of
        TRight -> sbRight sb
        TLeft  -> sbLeft sb

sbTree = genSB (SBTriple (FakeRatio 0 1) (FakeRatio 1 1) (FakeRatio 1 0))

divideFakeRatio r = (fromIntegral (numr r)) / (fromIntegral (denom r))

genIrrational ds = let r = findSB ds sbTree in divideFakeRatio r

fromReal' n sbTree = let sbh = (sbHead sbTree) in check sbh (divideFakeRatio sbh) where
    check sbh h | h == n = (numr sbh) % (denom sbh)
                | n < h  = fromReal' n (sbLeft sbTree)
                | otherwise = fromReal' n (sbRight sbTree)

sbLogB' b n sbTree = let sbh = (sbHead sbTree) in let sbh' = divideFakeRatio sbh in check sbh' where
    check sh | b**sh == n = sh
             | b**sh < n = sbLogB' b n (sbRight sbTree)
             | otherwise = sbLogB' b n (sbLeft sbTree)

sbLogBase b n = sbLogB' b n sbTree

fromReal n = fromReal' n sbTree

main = do
    print $ genIrrational (take 70 eDirections)
    print $ 1 / (genIrrational (take 20 phiDirections))
    print $ map fromReal [1.5,2..10]
    print $ sbLogBase 2 pi