Minimal processing of continued fractions

Learn
1

I’m a learner and, as an exercise, I’m trying to write some code to work with finite and infinite continued fractions. I have some code that is OK and seems to be correct: it emits a digit only when it’s certain, and processes another term from the list of coefficients when it needs to. I’m using the mathematics from: https://cdsmithus.medium.com/continued-fractions-haskell-equational-reasoning-property-testing-and-rewrite-rules-in-action-77a16d750e3f

I would like my code to be more elegant. At the moment, I have a special function to “drain” all remaining digits (maybe a finite number, maybe an infinitely repeating decimal) when the coefficients run out. But I feel like it must be possible to have a single emitting function that does both jobs.

matmul (a1, b1, c1, d1) (a2, b2, c2, d2) = (a1*a2 + b1*c2,
                                            a1*b2 + b1*d2,
                                            c1*a2 + d1*c2,
                                            c1*b2 + d1*d2)

cf_to_decimal :: [Integer] -> [Integer]
cf_to_decimal cf = go (1, 0, 0, 1) cf
  where
    emit :: (Integer, Integer, Integer, Integer) -> Maybe Integer
    emit (a,b,c,d)
      | c /= 0 && d /= 0 && y0 == y1 = Just y0
      | otherwise = Nothing
      where
        y0 = fromIntegral (a `div` c)
        y1 = fromIntegral (b `div` d)

    emitAll :: (Integer, Integer, Integer, Integer) -> [Integer]
    emitAll (a,b,c,d)
      | a == 0 = []
      | c /= 0 && d /= 0 && y0 == y1 = y0:go (matmul (10, -10*y0, 0, 1) (a,b,c,d)) []
      where
        y0 = fromIntegral (a `div` c)
        y1 = fromIntegral (b `div` d)

    go :: (Integer, Integer, Integer, Integer) -> [Integer] -> [Integer]
    go (a,b,c,d) (x:xs) =
      case emit (a,b,c,d) of
        Just d' -> d':go (matmul (10, -10*d', 0, 1) (a,b,c,d)) (x:xs)
        Nothing -> go (matmul (a,b,c,d) (x,1,1,0)) xs
    go (a,_,c,_) [] = emitAll (a,a,c,c)

A few experiments seem to have run into laziness problems? Here’s one attempt, inspired by unfoldr, that’s probably a bit of a mess:

matmul (a1, b1, c1, d1) (a2, b2, c2, d2) = (a1*a2 + b1*c2,
                                            a1*b2 + b1*d2,
                                            c1*a2 + d1*c2,
                                            c1*b2 + d1*d2)


cf_to_decimal cf = go (1, 0, 0, 1) cf
  where
    emit (digits, (a,b,c,d))
      | a == 0 = (digits, (a,b,c,d))
      | c /= 0 && d /= 0 && y0 == y1 = emit (digits ++ [y0], (matmul (10, -10*y0, 0, 1) (a,b,c,d)))
      | otherwise = (digits, (a,b,c,d))
      where
        y0 = fromIntegral (a `div` c)
        y1 = fromIntegral (b `div` d)

    go :: (Integer, Integer, Integer, Integer) -> [Integer] -> [Integer]
    go (0,_,_,_) _ = []
    go (a,b,c,d) (x:xs) = case emit ([], (a,b,c,d)) of
       ([], m) ->      go (matmul m (x, 1, 1, 0)) xs
       (d, m') -> d ++ go m' (x:xs)
    go (a,_,c,_) [] = case emit ([], (a,a,c,c)) of
       ([], _) -> []
       (d, m') -> d

Thanks for any comments.

1 Like
2

It would be good to hear what problems you’re having more specifically that you think are laziness related.

For readability it would be very useful to replace the four element tuple with a data type with named fields if possible.

1 Like
3

The four element tuple is just a matrix. I agree that there’s some tidying to be done there.

On the laziness issue: here’s what a session looks like with the broken code:

ghci> take 10 $ cf_to_decimal [1,2]
[1,5]
ghci> take 10 $ cf_to_decimal [1,4]
[1,2,5]
ghci> take 10 $ cf_to_decimal [1,3]
^CInterrupted.
ghci>

So when the output is finite, the code works, but not when it’s a recurring decimal. In contrast, this is what the “working” (if not correct) code does:

ghci> take 10 $ cf_to_decimal [1,2]
[1,5]
ghci> take 10 $ cf_to_decimal [1,4]
[1,2,5]
ghci> take 10 $ cf_to_decimal [1,3]
[1,3,3,3,3,3,3,3,3,3]
ghci>
4

Ok I think I see the loop in your second block. You have an element so you call emit on the empty list and your matrix. That immediately turns the empty list and the matrix again because c == 0 from the initial state. The empty list is then concatenated with go called with the same inputs as it started, meaning you infinitely concatenate the empty list.

My recommendation here would be to make sure that you always reduce your list length if possible, because at least then your algorithm will be finite.

1 Like
5

I see it as a stateful fold with a stateful unfold inside it:

{-# LANGUAGE LambdaCase #-}
import Control.Monad.State.Lazy
import Data.Foldable
import Control.Monad

type Mat = (Integer, Integer, Integer, Integer)

matmul (a1, b1, c1, d1) (a2, b2, c2, d2) = (a1*a2 + b1*c2,
                                            a1*b2 + b1*d2,
                                            c1*a2 + d1*c2,
                                            c1*b2 + d1*d2)

unfoldM :: Monad m => m (Maybe a) -> m [a]
unfoldM m = m >>= \case
  Nothing -> pure []
  Just x -> (x :) <$> unfoldM m

cf_to_decimal :: [Integer] -> [Integer]
cf_to_decimal cf = evalState (foldr feed end cf) (1,0,0,1) where

  emit :: State Mat (Maybe Integer)
  emit = do
    m@(a,b,c,d) <- get
    let y0 = a `div` c
        y1 = b `div` d
    if c /= 0 && d /= 0 && y0 == y1 then do
      put (matmul (10, -10*y0, 0, 1) m)
      pure (Just y0)
    else pure Nothing

  feed :: Integer -> State Mat [Integer] -> State Mat [Integer]
  feed x rest = do
    modify (\m -> matmul m (x,1,1,0))
    (<>) <$> unfoldM emit <*> rest

  end :: State Mat [Integer]
  end = do
   modify (\(a,b,c,d) -> (a,a,c,c))
   unfoldM $ do
     (a,_,_,_) <- get
     if a == 0 then pure Nothing else emit

This also gives you the laziness you want if you use the lazy state monad.

1 Like
6

This formulation also made me wonder whether there’s an opportunity to encode productivity at the type/structure level rather than relying on the y0 == y1 check alone.

For example, you could imagine splitting the state machine into two explicit phases:

  • Accumulate: consume CF terms until the matrix is “stable enough.”
  • Emit: guaranteed to produce at least one digit before needing more input

In other words, make it emit totally once entered, so the only partiality lives in the transition between phases. That might let you model the whole thing as something like a Mealy machine or even a Stream transformer, where non-productivity becomes structurally impossible rather than a runtime condition.

It feels like continued fractions are crying out for that kind of “digit-safe” invariant, similar to how digit-by-digit algorithms for √2 or π enforce progress by construction. I don’t know if it simplifies the code, but it might make the reasoning story even cleaner.

Either way, this fold-inside-unfold pattern is a great lens - it’s much clearer than the recursive versions I was circling around before.

7

Thanks. I think there was a problem with the tail-recursive call to emit. There are two places where I need to extract the first part of emit’s return value. Here’s what seems to be working now:

matmul (a1, b1, c1, d1) (a2, b2, c2, d2) = (a1*a2 + b1*c2,
                                            a1*b2 + b1*d2,
                                            c1*a2 + d1*c2,
                                            c1*b2 + d1*d2)

cf_to_decimal cf = go (1, 0, 0, 1) cf
  where
    emit (digits, m@(a,b,c,d))
      | a == 0 = (digits, m)
      | c /= 0 && d /= 0 && y0 == y1 = 
        let (d', m') = emit ([], (matmul (10, -10*y0, 0, 1) m)) in (y0:d', m')
      | otherwise = (digits, m)
      where
        y0 = fromIntegral (a `div` c)
        y1 = fromIntegral (b `div` d)

    go :: (Integer, Integer, Integer, Integer) -> [Integer] -> [Integer]
    go m (x:xs) = case emit ([], m) of 
       ([], m) ->      go (matmul m (x, 1, 1, 0)) xs
       (d, m') -> d ++ go m' (x:xs)
    go (a,_,c,_) [] = fst $ emit ([], (a,a,c,c))
8

Thanks. This is very clear. I hadn’t encountered the (<>) <$> <*> idiom before.

What’s the practical difference between using the straight State monad and the lazy one, here? I’m curious, for example, about calculating very large numbers of digits in finite memory. I see Kwang's Haskell Blog - Lazy vs Strict State Monad but perhaps memory exhaustion will come about in some other way, not involving the State monad.

9

Since emit may produce more than one digit, it seemed to me that having it return a list of digits made sense. There’s a straightforward equivalence between Nothing and []. On the other hand, while there seems to be a type of non-empty lists in Kmett’s libraries, it seems like a more complex return type than a plain list. Calling emit once per digit seems standard – so would you check the matrix on emit’s return, to see if you needed to call emit again before consuming further terms of the fraction? The matrix will have been modified by emi,t but only if a digit was produced.

I’m keen on an emit function that recurses indefinitely in the case of a rational number, in the “terminal” phase.

10

Yeah, f <$> x <*> y is basically just function application, but where x and y can be applicatives. And you can extend it for as many arguments as you want: f <$> x1 <*> x2 <*> x3 <*> x4.

One gives you the lazy behavior you want and the other has to first compute all the state updates before it can even return one element of the result. I’m not sure about the efficiency. If you want to make sure it is not consuming a linear amount of memory relative to the number of digits, then you need to use either a lazy algorithm or use a manual streaming approach, for example conduit or streamly, or bluefin (from old to new).

11

You might be interested in the quadratic-irrational package, which has some functionality involving continued fractions.

1 Like
12

Interesting, thanks. The description for that package remarks that “Every finite SCF represents a rational number and every infinite, periodic SCF represents a quadratic irrational.” A slightly stronger claim is “The real numbers whose continued fraction eventually repeats are precisely the quadratic irrationals.” (Wikipedia) So repeating continued fractions don’t just cover the rationals, unlike repeating decimals, and quadratic irrationals are part and parcel of the continued fraction picture.

The example of √2 is neat:

Step 1

Integer part: 1
Fractional part: .414…
Inverse of fractional part: 2.414…

Step 2

Integer part: 2
Fractional part: .414…
Inverse of fractional part: 2.414…

and so on, repeating.

It works because it corresponds to rearranging the formula x^2 = 2 as x^2 - 1 = 1 and then (x+1)(x-1) = 1. So the fractional part, .414, equals 1/(x+1), or in other words 1/2.414….

It might be interesting to write some Haskell code which takes any quadratic equation and produces the continued fraction of its real root(s), when they exist, purely symbolically.