Reverse Greed for Unit Fractions

Here's an interesting algorithm for expanding any fraction into a sum 
of unit fractions.  The method seems to give quite nice expansions 
that look rather "Egyptian".  In fact, it reproduces many of the 
expansions listed in the Rhind Papyrus, etc.  But, in addition, it 
easily generates a nice four-term expansion of 31/311 and in general 
seems to work well even on "difficult" fractions.  Still, it's so 
simple that I could believe the ancient Egyptians used something 
like it.

I suspect the reason I never thought of it before is because I tend 
to think in "greedy" ways, i.e., when trying to split up a fraction 
7/13 my first instinct is to subtract 1/2 and then try to squeeze 
smaller unit fractions into the remainder.  Maybe this is a "modern" 
thought pattern.  It actually makes more sense to choose the smallest
parts first, and work toward the larger parts.  The might, which 
might be called the Reverse Greedy algorithm, is just this:

  Let N/QR be a fraction in lowest terms with Q>1 and (Q,R)=1, 
  and let x,y be the smallest positive solution of Nx-Qy=R.  Then
  1/Qx is a term in the expansion, with remainder y/Rx.  Expand 
  this remainder and iterate until reaching a unit fraction 
  remainder.
   
This is actually even simpler than it sounds.  Take for example the 
"difficult" fraction 31/311.  Here we have N=31, Q=311, R=1. (Notice
that the denominator QR=311, and we require Q>1 and also that Q,R must
have no common factor.  Therefore, we set Q=311 and R=1.)  So for this
first round we need the smallest integer solution of 31x - 311y = 1.
It's easy to see (either by trial and error, or by the Euclidean 
algorithm) that x=301,y=30 is the smallest solution.  Therefore, the 
first term in our expansion is 1/Qx = 1/(311)(301) = 1/93611, with 
the remainder y/Rx = 30/(1)(301).  In other words, we've shown that

                 31         1         30
                 ---  =   -----   +   ---
                 311      93611       301

Now we apply the same procedure to 30/301, which gives

                 30         1         11
                 ---  =   -----   +   ---
                 301       688        112

Then we apply the procedure to 11/112 to find

                 11        1         1
                 ---  =   ---   +   ---
                 112       28        16

and we're done.  We have found

           31        1       1      1       1
           ---  =   ---  +  --- +  ---  + -----
           311       16     28     688    93611

Notice that we found these terms in reverse order, beginning with 
93611'.  Admittedly this denominator is larger than 8708, but it IS
the smallest possible if you allow only one of the denominators to 
be a multiple of 311.  Also, the 2nd and 3rd terms are smaller than
in the 8708 solution.

Interestingly, this method gives 2/91 = 70' + 130' and many others 
from the Rhind 2/n table.  Just for fun I tried it on the fraction 
163/431, which seems like it might be difficult.  The algorithm quite
simply gives the expansion

       163      1       1       1       1        1
       ---  =  ---  +  ---  +  ---  +  ---  +  ------
       431      3       42      55     350     118525

Comparing this method to the 10 methods that David Eppstein describes
on his web page, we notice that in one special circumstance it gives 
the same results as the continued fraction method.  The special 
circumstance is if all the values of R in the process equal 1.  This
rarely happens, because it implies none of the denominators in the 
sequence can be factored into non-trivial coprime factors.  It just 
so happens that one of David's examples, 18/23, is such a case, so 
both methods give 

          18      1       1       1       1       1
          --  =  ---  +  ---  +  ---  +  ---  +  ---
          23      2       6       12      36     207

However, for 7/15 the basic Continued Fraction Method gives

     7      1       1       1       1       1       1       1
    --  =  ---  +  ---  +  ---  +  ---  +  ---  +  ---  +  ---
    15      3       15      35      63      99     143     195

whereas my new method gives

                7       1       1       1
               ---  =  ---  +  ---  +  ---
                15      4       6       20

which seems much closer to the Egyptian style.  (This is even better 
than what David calls the Grouped Continued Fraction Method, which is 
a quite complicated algorithm, and gives 7/15 = 3' + 9' + 45'.)  All 
in all, this new method seems to consistently produce representations
that are very close to the spirit of the historical Egyptian unit 
fractions.

One interesting aspect of the algorithm described above is that, in
order to yield unique results, must include a rule for factoring the 
denominators.  Interestingly, with Q=1 it's equivalent to the Greedy 
Method, whereas with R=1 it's equivalent to the Continued Fraction 
Method.  Thus it appears that a wide variety of methods can be 
expressed in terms of the above algorithm together with some fixed 
rule for factoring the denominators.  

I would classify these as "local" algorithms, i.e., methods that 
proceed in a step-wise manner by optimising some function of the 
immediate remainder.  This is in contrast to "global" algorithms, 
which optimise some function of the entire completed expansion.  I 
suspect that any sort of general global optimization is NP-complete,
and that no local algorithm can guarantee a globally optimum result 
in every case.

The stipulation that (Q,R)=1 was included as sort of an example of 
one way we could restrict the choices (not meant to imply that this 
condition is required for Nx-Qy=R to be solvable), but of course it 
doesn't completely specify the decomposition of denominators that are
divisible by more than two distinct primes.  My first implementation 
of this method set Q equal to the highest power of the largest prime 
dividing the denominator, the idea being to completely eliminate the 
largest prime from the denominator of the remainder.  (By the way, 
to answer one of David Eppstein's questions, this explains why the 
method gave R=1 at each step of the expansion of 18/23, because the
denominators were all prime powers.)

Of course, we're free to base the algorithm on any fixed factoring 
rule we choose.  David Eppstein has experimented with a couple of 
different rules.  I've tried dropping the (Q,R)=1 requirement and 
just choosing the factors with Q>R such that Q-R is minimized.  This
rule gives an algorithm that is closely related to the "Fibonacci" 
methods.

Overall, the basic algorithm can be regarded as a means of unifying 
and classifying all the various "local" methods, each corresponding 
to a particular rule of factorization.  Given a set of expansions 
from some historical document, you could check to see if some fixed 
factorization rule generates all those expansions.

From a historical standpoint, it seems (from what little I know of 
the surviving documentary evidence) that the ancient Egyptians didn't 
use any single rule for generating their unit fraction expansions.  As
a result, the small quantity of material must be sub-divided into even
smaller sets of examples of any particular method, which makes it 
almost impossible to build conclusive arguments about the thought 
processes they may have used.  If only the 2/n table in the Rhind 
Papyrus had gone up to 2/1001 instead of just 2/101 we would be in 
a much better position to "reverse engineer" their methods!

Recall that we found the expansion

        31         1        1        1         1
        ---  =   -----  +  ---  +  -----  +  -----
        311       16        28      688      93611

Although several local methods yield this four-term expansion, I've 
not found one that gives the expansion with smallest max denominator.
However, we can expand our analysis in a global direction by expanding
into TWO unit fractions at once.  

Given a fraction N/D in lowest terms, we begin with some factorization
of the denominator D=QR and then deduct TWO unit fractions from N/D, 
both of which have denominators divisible by Q.  (Historically this 
comes from the fact that many of the Egyptian expansions seem to have
the final two or even three denominators divisible by a large prime 
factr of D.)  So we have

              N        1   / 1     1 \        n
             ----  -  --- ( --- + --- )  =   ---
              QR       Q   \ x     y /        d

which can be written in the form

                  Nxy - R(x+y)        n
                  ------------   =   ---
                      QRxy            d

We want to eliminate Q from the denominator of n/d, so we have the 
condition

                       Nxy = R(x+y)   (mod Q)

which has the smallest solution x=9,y=28, leaving the remainder 
25/252.  We can then check by the usual method to see if 25/252 has 
a two-term expansion, which it does, 1/14 + 1/36, so we have the 
result

        31        1        1        1         1
        ---  =   ----  +  ---  +  -----  +  -----
        311       14       36      2799      8708

So this algorithm consists of first checking for simple 2-term 
expansions, then (assuming none exist) expanding the first two 
high-order terms, and iterating on the remainder.  (Obviously we 
still need to specify a fixed facotrization rule to split QR at 
each stage.)

I don't imagine the ancient Egyptians practiced anything like this.  
It's more likely that their expansions tended to have pairs (or 
triples) of denominators with common factors because they produced
them by multiplying binomials (or trinomials) of smaller expansions.

One more interesting fact about the ratio 31/311 is that, although 
it has no three-term expansion, the nearest three-term expansion 
differs from the true value by a unit fraction, i.e.,

         31      1       1        1               1
        ---  =  ---  +  ---  +  -----    -   ----------
        311      11     115     13566        5337067890

where the final (negative) denominator is just the least common
multiple of all the other denominators.