Unit Fractions and Fibonacci

Suppose we wish to generate integer solutions of the harmonic equation

                      1       1       1
                     ---  +  ---  =  ---
                      x       y       z

This can be regarded as a special case of a more general expansion 
related to the Fibonacci numbers.  Let s[j], j=0,1,2,... be a sequence 
of integers that satisfy the recurrence s[k] = s[k-1] + s[k-2] with 
arbitrary initial values s[0] and s[1].  It can be shown that for 
any integers m,n with m>n we have

        1                 1            m          1
   -----------   =   -----------  +  SUM    -------------       (1)
   s[n-1] s[n]       s[m] s[m+1]      j=n   s[j-1] s[j+1]


For example, setting s[0]=s[1]=1 and n=5, m=10 gives

       1        1     1     1      1      1      1      1
      ---  =   --- + --- + --- + ---- + ---- + ---- + -----
       40       65   168   442   1155   3026   7920   12816

In general, to expand 1/D into a sum of unit fractions, the method
is to split D into two factors, D = pq.  Then we can set s[0]=p and
s[1]=q and generate the s sequences as follows

           k    s[k]     s[k] s[k-1]      s[k] s[k-2]
          ---  -------   -----------      -----------
           0    p
           1    q        pq
           2    p+q      q(p+q)           p(p+q)
           3    p+2q     (p+q)(p+2q)      q(p+2q)
           4    2p+3q    (p+2q)(2p+3q)    (p+q)(2p+3q)
           5    3p+5q    (2p+3q)(3p+5q)   (p+2q)(3p+5q)
           6    5p+8q    (3p+5q)(5p+8q)   (2p+3q)(5p+8q)
           7    8p+13q   (5p+8q)(8p+13q)  (3p+5q)(8p+13q)
                 etc            etc            etc

We can now express 1/pq as the sum of the inverses of the numbers 
in the third column down to the mth row, plus the inverse of the mth 
number in the second column.  Thus we have

 1/pq  =  1/p(p+q) + 1/q(p+q)

       =  1/p(p+q) + 1/q(p+2q) + 1/(p+q)(p+2q)

       =  1/p(p+q) + 1/q(p+2q) + 1/(p+q)(2p+3q) + 1/(p+2q)(2p+3q)

                              etc.

Of course, we can let m in equation (1) go to infinity, giving the 
infinite unit fraction expansion

                  1              inf         1
             -----------   =    SUM    -------------           (2)
             s[n-1] s[n]         j=n   s[j-1] s[j+1]


This can also be generalized to higher order recurrences.  For
example, if we define the sequence s[j] to satisfy the 3rd order
recurrence s[k] = s[k-2] + s[k-3] with the initial values a,b,c,
then we can generate the following sequences

       k    s[k]      s[k]s[k-1]s[k-2]        s[k]s[k-1]s[k-3]
      ---  ------   --------------------    ----------------------
       0    a
       1    b
       2    c          abc
       3    a+b        bc(a+b)                ac(a+b)
       4    b+c        c(a+b)(b+c)            b(a+b)(b+c)
       5    a+b+c      (a+b)(b+c)(a+b+c)      c(b+c)(a+b+c)
       6    a+2b+c     (b+c)(a+b+c)(a+2b+c)   (a+b)(a+b+c)(a+2b+c)
             etc              etc                     etc

so we have

1/abc = 1/ac(a+b) + 1/bc(a+b)

    = 1/ac(a+b) + 1/b(a+b)(b+c) + 1/c(a+b)(b+c)

  = 1/ac(a+b) + 1/b(a+b)(b+c) + 1/c(b+c)(a+b+c) + 1/(a+b)(b+c)(a+b+c)

and so on.  To illustrate, with a=3,b=7,c=11 this last formula gives

    1/231  =  1/330 + 1/770

           =  1/330 + 1/1260 + 1/1980

           =  1/330 + 1/1260 + 1/3780 + 1/4158

and with a=23,b=c=1 it gives

     1/23  =  1/24 + 1/552

           =  1/48 + 1/50 + 1/552 + 1/1200