Haskell script for calculating good rational approximations to real numbers

approximate-ratio.hs

#!/usr/bin/env runhaskell
{-# OPTIONS -Wno-x-partial #-}
import Data.List (inits, transpose, intercalate)
import Data.Ratio
import System.Environment (getArgs)
import System.Exit (die)
import Text.Printf

main :: IO ()
main = do
  args <- getArgs
  (cfg, number) <- parse defaultConfig args
  let approxes = if useFarey cfg then farey number else approximate number
  let filter = case atmost cfg of Just n -> take n; Nothing -> id
  printTable $ filter approxes

data Config = Config { useFarey :: Bool, atmost :: Maybe Int }
defaultConfig = Config { useFarey = False, atmost = Just 10 }

parse :: Config -> [String] -> IO (Config, Double)
parse cfg ("-farey": args) = parse (cfg { useFarey = True }) args
parse cfg ("-continued": args) = parse (cfg { useFarey = False }) args
parse cfg ("-all": args) = parse (cfg { atmost = Nothing }) args
parse cfg ["-n"] = die "-n expects a natural number argument"
parse cfg ("-n": n: args) = parse (cfg { atmost = Just (read n) }) args
parse cfg [n] = return (cfg, read n)
parse cfg [] = die "expected a real number argument"
parse cfg (_:_:_) = die "too many arguments"


---------- DISPLAYING STUFF ----------
-- TODO: also print accuracy/difference from original input
-- NB. assumes nonnegative.
printTable :: [Rational] -> IO ()
printTable [] = error "empty table"
printTable approxes = mapM_ (putStrLn . formatRow) table
  where formatRow [num, denom, int, frac] =
          concat [ padLeft numLen num
                 , " / "
                 , padLeft denomLen denom
                 , "  →  "
                 , padLeft intLen int
                 , frac ]
        [numLen, denomLen, intLen, fracLen] = map (maximum . map length) $ transpose table
        table = [let real :: Double = fromRational r in
                 let (integral, fractional) = span (/= '.') (show real) in
                 [show (numerator r),
                  show (denominator r),
                  integral,
                  fractional]
                | r <- approxes]
        padLeft size s = replicate (size - length s) ' ' ++ s
        padRight size s = s ++ replicate (size - length s) ' '

-- Simpler methods for interactive use.
display :: Rational -> String
display r = show r ++ "  →  " ++ show (fromRational r)

displayN :: Int -> [Rational] -> IO ()
displayN n xs = mapM_ (putStrLn . display) $ take n xs


---------- METHOD 1: CONTINUED FRACTIONS ----------
-- Finds good rational approximations for a target number, from simpler to more complex.
approximate :: RealFrac n => n -> [Rational]
-- The 1st approximation (zero) and 2nd (the integral part of the number) are
-- too boring to list.
approximate r = drop 2 approx ++ take 1 correct
  where (approx, correct) = span ((r /=) . fromRational) $ approximations r

approximations :: RealFrac n => n -> [Rational]
approximations r = map discontinue $ inits $ continued r

-- Finds the continued fraction for a real number.
-- Almost always produces an infinite list.
-- probably runs out of floating-point precision and into useless noise quite fast.
-- NB. may not do the right thing for negative numbers.
continued :: RealFrac n => n -> [Integer]
continued r = i : rest
  where (i, f) = properFraction r
        rest | f /= 0 = continued (recip f)
             | otherwise = []

-- Interprets a (finite!) continued fraction into a rational number.
-- eg. ratio (take 2 (continued pi)) == 22 % 7.
discontinue :: [Integer] -> Rational
discontinue [] = 0 -- not sure if correct/meaningful
discontinue [x] = x % 1
discontinue (x:xs) = (x % 1) + recip (discontinue xs)


---------- METHOD 2: FAREY ----------
-- Finds progressively more accurate bounding intervals via Farey sequences.
-- Often generates more approximations than `approximate`. Good if you want to
-- find all simple, plausibly close ratios. Bad if you want to quickly get good
-- approximations, or find only the 'best' simple approximations.
farey :: RealFrac n => n -> [Rational]
farey n = case loop (round n % 1) (floor n % 1) (ceiling n % 1) of
            [] -> [toRational n]
            xs -> xs
  where distance x = abs (n - fromRational x)
        loop :: Rational -> Rational -> Rational -> [Rational]
        loop best lo hi = if better then middle : rest else rest
          where middle = mediant lo hi
                middleReal = fromRational middle
                better = distance middle < distance best
                best' = if better then middle else best
                rest | n < middleReal  = loop best' lo middle
                     | middleReal < n  = loop best' middle hi
                     | n == middleReal = []
                     | otherwise = error "impossible"

-- Finds bounding intervals. Simpler. Not used by main.
fareyBounds :: RealFrac n => n -> [(Rational, Rational)]
fareyBounds n = loop (floor n % 1) (ceiling n % 1)
  where loop lo hi = (lo, hi) : rest
          where middle = mediant lo hi
                middleReal = fromRational middle
                rest | n == middleReal = [(middle, middle)]
                     | n < middleReal = loop lo middle
                     | middleReal < n = loop middle hi
                     | otherwise = error "impossible"

mediant :: Rational -> Rational -> Rational
mediant n m = (numerator n + numerator m) % (denominator n + denominator m)