Odd Greedy Expansions

Last post I mentioned the number theory class I'm in. Dr. Dawsey had us look up open problems in number theory and pick one to explore. I picked the Odd Greedy Expansion problem.

Every positive rational number can be represented as a sum of unique unit fractions. This is called an Egyptian fraction. Fibonacci described a greedy algorithm for obtaining Egyptian fraction expansions for rationals in his Liber Abaci (1202).

It is also true that every rational number with an odd denominator can be represented as a sum of unit fractions with odd denominators, e.g. \[\frac{p}{q}=\sum_{k\in\mathbb{Z}}\frac{1}{2k+1} \] for odd \(q\). The greedy algorithm can be modified to produce only odd denominators; this is the odd greedy algorithm.

I implemented this algorithm in SML for fun. My comments explain most of what's going on, and there's a real explanation below. Just thought I'd post it in case anyone has a use for it. All the code is now on GitHub; the arbitrary-precision rational implementation is originally from here (I modified it somewhat).

fun oge (0, _) = print "0\n" 
  | oge (1, b) = print ("1/" ^ IntInf.toString b ^ "\n")
  | oge (a, b) =
    if b mod 2 = 0 then raise DenominatorNotOdd
    else if a > b then raise ImproperFraction
    else
    (* STEPS:
     * 1. u = min{ odds greater than y div x }
     *    i.e. if floor(y/x) is even then floor(y/x)+1 else floor(y/x)+2
     * 2. include 1/u in the expansion and recurse on x/y - 1/u.
     * 3. terminate when the remainder is a unit fraction
     *    you haven't already included.
     *)
        let
          open rational (* Bring my arbitrary-precision rational
                           datatype into scope *)
          infix == --   (* Allow these functions to be infix operators *)
          val firstr = new (a, b) (* Construct a rational *)
          fun recurse (r : rational, acc) =
            (* If we're given a unit fraction with num = 1,
             * & the denominator is not equal to the previous
             * one, then it must be unique and it is the final
             * term in the sum. We're done; reverse and return the list.
             *)
            if num r = 1 andalso 
               (den (hd acc) handle _ => IntInf.fromInt 1) <> den r 
               then List.rev (r :: acc)
            else
              let
                (* The floor of the inverse, or the nearest
                   integer below this rational *)
                val flr = IntInf.div (den r, num r)
                (* The nearest odd number above this rational *)
                val u = if flr mod 2 = 0 then flr + 1 else flr + 2
                (* Inverse of the nearest odd integer which we
                   haven't used yet *)
                val uinv = if null acc then
                                ratinv u
                           else if u <= den (hd acc) then
                                ratinv (den (hd acc) + 2)
                           else ratinv u
                (* Difference of r and uinv, or p/q - 1/u. *)
                val remaining = r -- uinv

                val intgrstr = rational.show uinv
              in
                print ("\nSize of denom: " ^ Int.toString (size intgrstr));
                print ("\nTerm: " ^ intgrstr);
                (* Stick uinv into the sum and split up the
                   remaining amount *)
                recurse (remaining, uinv :: acc)
              end

          (* Look ma, higher-order functions! *)
          fun printlst R = 
            print (String.concatWith ", " (map rational.show R) ^ "\n")
        in
          (* Start with a/b and an empty list *)
          printlst (recurse (firstr, []))
        end

The algorithm is greedy because at every step it subtracts the largest odd-denominatored (that's definitely a word) unit fraction it can find. If the input is \(\frac{p}{q}\) then it subtracts the inverse of the least odd number greater than \(\frac{q}{p}\), which we haven't already used as a denominator. It repeats this process until it comes up with a unit remainder which hasn't already appeared in the expansion.

For example, \(\frac{2}{3}\) is expanded as \[ \frac{2}{3} = \frac{1}{3} + \frac{1}{5} + \frac{1}{9} + \frac{1}{45}. \]

I found an interesting expansion by iterating over a search space. It seems \(\frac{5}{139}\) has an expansion with at least 17 terms. I terminated the program after the seventeenth term, because that 17th term's denominator has 90700 digits if my code is to be believed. Here were the first few terms: \[ \frac{5}{139} = \frac{1}{29} + \frac{1}{673} + \frac{1}{387\,553} + \frac{1}{131\,422\,274\,281} + \frac{1}{15\,352\,723\,712\,926\,705\,785\,241} + \cdots \]

Sometimes these expansions blow up to enormous sizes for whatever reason, but usually they terminate in roughly four or five terms, from my cursory inspection. This is a reasonably fast and efficient program, but some particular values like \(\frac{5}{139}\) will take a while to compute.

If you play with the files, I made another function called testOverRange (p, q) which expands all fractions with numerators q and valid odd denominators equal to or less than q. Please don't plug in unreduced fractions. :)

As always, email me at ethan.barry@howdytx.technology if you have any comments!

Sources:

MathPages' old file on these problems: https://www.mathpages.com/home/kmath454.htm

The Wikipedia article: https://en.wikipedia.org/wiki/Odd_greedy_expansion

And this implementation of rationals in SML!