The Generalized Remainder Method
Here we assume a value p given by the user, and use a close-to-brute-force dynamic programming technique to represent numbers as sums of divisors of p. We then use EgyptPairList because our technique does not necessarily avoid duplicate fractions (and because the chief goal here is to minimize the number of terms in the expansion).
SumDivisors[n_,p_] := DivisorSums[n,p, Table[If[p ~Mod~ i == 0, {{i}}, {}],{i,n}]]; DivisorSums[n_,p_,t_] := If [t[[n]] != {}, t[[n]], DivisorSums[n,p,ExtendDivTable[t,p,n]]]; ExtendTableEntry[e_,d_] := Append[#,d]& /@ Select[e,Last[#] <= d&]; ExtendDivTable[t_,p_,n_] := Table[If[t[[i]] != {}, t[[i]], Join @@ ( ExtendTableEntry[t[[i-#]],#]& /@ Select[Divisors[p], #<i&] ) ], {i,n}] EgyptRemainder[Rational[x_,y_],p_] := Module[{r, s}, q = Quotient[x p, y]; r = Mod[x p, y]; EgyptPairList /@ OuterJoin[SumDivisors[q,p]/p, SumDivisors[r,p]/(p*y)]];
In general, the behavior of this algorithm depends on the number p, which should be chosen to have many divisors (factorials are reasonable choices), and to make q and r both small. With a good choice this seems to be capable of generating very short representations, with reasonably small denominators.
EgyptRemainder[31/311,252]
1 1 1 1 1 1 1 1 {{--, --, ----, ----}, {--, --, ----, -----}, 14 36 2799 8708 14 36 2177 78372 1 1 1 1 1 1 1 1 {--, --, ----, ----}, {--, --, ----, -----}} 12 63 2799 8708 12 63 2177 78372
If y is prime, then any representation of x/y must have some terms 1/ky. In particular, the sum of these terms must be a/by where xb-a is divisible by y. (From this we can quickly see that for some x, the denominators in any representation must be Omega[y log[y]]: each such combination a/b works for only one x so we need at least y different combinations to be available.)
The idea behind the following method is simply to try all small combinations of terms until finding one cancelling the factor of y. The result is a fraction which may well end up with a higher denominator than y, but likely has more divisors.
If y is not prime, we first start by greedily removing fractions d/y where d is a divisor of y; the hope is that this step makes x sufficiently small that a smaller combination a/b will work.
Once we have found a combination that works, we are likely to have a denominator with many more divisors, so we apply some other method better suited to that case (here, the small-denominator reverse greedy method) rather than continuing recursively.
Subsets[s_] := Flatten /@ Distribute[{{},{#}}& /@ s, List,List]; ESMGoodFrac[Rational[a_,b_],x_,y_] := (x*b - a) ~Mod~ y == 0; ESMSubset[x_,y_,i_] := Select[Subsets[Table[1/j,{j,i}]], ESMGoodFrac[Plus@@#,x,y]&]; ESMFind[x_,y_,i_] := If[# == {}, ESMFind[x,y,i+1], #]& @ ESMSubset[x,y,i]; ESMFrac[x_,y_] := Join[EgyptSmallDen[x/y - Plus@@#], #]& /@ (ESMFind[x,y,1] / y); EgyptSmallMult[Rational[x_,y_]] := Reverse[Sort[#]]& /@ ESMFrac[x,y];
It is not clear whether all terms in the representation produced by this method are distinct. However it seems to produce good results in practice for moderate inputs. If y is prime, and this method produces a correct representation in which the largest denominator is divisible by y, it will be the representation with the minimum possible denominator.
EgyptSmallMult[31/311]
1 1 1 1 1 1 1 1 {{--, --, --, ---, ----, ----, ----, ----}, 28 30 36 933 1555 2177 2799 3110 1 1 1 1 1 1 1 {--, --, ---, ----, ----, ----, ----}} 15 36 311 1244 1866 2799 3110
The Short Sequence Method
To make this method less slow (and perhaps even more important, to make it use less memory), we speed up the case for k=2. We wish to solve the equation x/y=1/a + 1/b; this can be rewritten (ax-y)(bx-y)=y^2, and letting the two factors of y^2 be r and s we can solve a=(r+y)/x, b=(s+y)/x. We try all factors r of y^2 for which r<y and collect the ones that make a,b integers. (Note that if a is an integer, b must also be an integer, so we only need one test.)
EgyptTwoTerm[Rational[x_,y_],minden_,preplist_] := Join[preplist,{x/#,x(#-y)/(y #)}]& /@ Select[ y + Take[#,Quotient[Length[#],2]]& @ Divisors[y^2], Mod[#,x]==0 && # >= minden*x & ]
EgyptShort[q_,k_] := EShortRecur[q,k,1,{}]; EShortRecur[q_,k_,b_,a_] := If [q == 0, {a}, If [k <= 2, EgyptTwoTerm[q, b, a], Join @@ Table[EShortRecur[q-1/i,k-1,i+1,Append[a,1/i]], {i,Max[b,Ceiling[1/q]],Floor[k/q]}]]];
This method is quite slow, taking four minutes on my Powerbook 540c to solve the following example. A C++ version that I also coded up performs much more quickly ("solving" the same problem in less than a minute) but because it is less careful than Mathematica about arithmetic overflow it gets the wrong answers.
EgyptShort[31/311,4]
1 1 1 1 1 1 1 1 {{--, ---, -----, --------}, {--, ---, -----, -------}, 11 115 13570 46422970 11 115 13662 1931310 1 1 1 1 1 1 1 1 {--, ---, -----, -------}, {--, ---, ----, --------}, 11 115 13684 1573660 11 116 6728 23016488 1 1 1 1 1 1 1 1 {--, ---, ----, ------}, {--, ---, ----, ------}, 11 116 6842 396836 11 118 3421 403678 1 1 1 1 1 1 1 1 {--, ---, ----, ------}, {--, --, ----, --------}, 11 119 2772 190332 12 62 4628 66927822 1 1 1 1 1 1 1 1 {--, --, ----, ------}, {--, --, ----, ------}, 12 62 4650 964100 12 62 4665 578460 1 1 1 1 1 1 1 1 {--, --, ----, ------}, {--, --, ----, ------}, 12 62 4743 190332 12 63 2124 770658 1 1 1 1 1 1 1 1 {--, --, ----, ------}, {--, --, ----, -----}, 12 63 2142 190332 12 63 2177 78372 1 1 1 1 1 1 1 1 {--, --, ----, ----}, {--, --, ----, ------}, 12 63 2799 8708 12 64 1392 577216 1 1 1 1 1 1 1 1 {--, --, ---, ------}, {--, --, ---, -----}, 12 68 612 190332 12 68 622 31722 1 1 1 1 1 1 1 1 {--, --, ---, -------}, {--, --, ---, ------}, 12 69 540 1931310 12 72 408 190332 1 1 1 1 1 1 1 1 {--, --, ---, ------}, {--, --, ---, -------}, 12 78 284 287053 12 86 212 4252614 1 1 1 1 1 1 1 1 {--, ---, ---, ------}, {--, ---, ---, ------}, 12 102 153 190332 12 119 126 190332 1 1 1 1 1 1 1 1 {--, --, -----, ---------}, {--, --, -----, ---------}, 13 44 35580 791174670 13 44 35581 486890404 1 1 1 1 1 1 1 1 {--, --, -----, --------}, {--, --, -----, --------}, 13 44 35596 71957314 13 44 35607 44295108 1 1 1 1 1 1 1 1 {--, --, -----, --------}, {--, --, -----, --------}, 13 44 35620 30464005 13 44 35646 18760764 1 1 1 1 1 1 1 1 {--, --, -----, -------}, {--, --, -----, -------}, 13 44 35772 6573918 13 44 35893 4059172 1 1 1 1 1 1 1 1 {--, --, -----, -------}, {--, --, -----, -------}, 13 44 36036 2801799 13 44 36076 2579434 1 1 1 1 1 1 1 1 {--, --, -----, -------}, {--, --, -----, -------}, 13 44 36322 1737868 13 44 36387 1601028 1 1 1 1 1 1 1 1 {--, --, -----, ------}, {--, --, -----, ------}, 13 44 40612 287053 13 44 41052 266838 1 1 1 1 1 1 1 1 {--, --, -----, ------}, {--, --, -----, ------}, 13 44 43758 190332 13 44 44473 177892 1 1 1 1 1 1 1 1 {--, --, -----, ------}, {--, --, -----, -----}, 13 44 48516 133419 13 44 56602 95788 1 1 1 1 1 1 1 1 {--, --, ---, ------}, {--, --, ----, ------}, 13 52 284 287053 14 36 2124 770658 1 1 1 1 1 1 1 1 {--, --, ----, ------}, {--, --, ----, -----}, 14 36 2142 190332 14 36 2177 78372 1 1 1 1 1 1 1 1 {--, --, ----, ----}, {--, --, --, -------}, 14 36 2799 8708 15 54 69 1931310 1 1 1 1 1 1 1 1 {--, --, ----, --------}, {--, --, ----, --------}, 16 27 7072 59383584 16 27 7074 17600112 1 1 1 1 1 1 1 1 {--, --, ----, ------}, {--, --, ----, ------}, 16 27 7344 190332 16 27 7464 134352 1 1 1 1 1 1 1 1 {--, --, ----, -----}, {--, --, ---, --------}, 16 27 8397 44784 16 28 683 23790256 1 1 1 1 1 1 1 1 {--, --, ---, -----}, {--, --, ----, ------}, 16 28 688 93611 18 23 1555 643770 1 1 1 1 1 1 1 1 {--, --, ----, ------}, {--, --, ---, ------}, 18 23 1564 190332 18 24 408 190332 1 1 1 1 1 1 1 1 {--, --, ---, ------}, {--, --, --, ------}} 18 28 119 190332 18 34 68 190332
Many of these expansions could have also been obtained by EgyptRemainder, using the following numbers as bases:
Union[((LCM@@ Denominator/@ #) / 311)& /@ %]
{204, 252, 432, 572, 612, 1224, 1276, 1298, 1716, 1860, 2070, 3692, 4004, 4284, 4816, 5060, 5148, 5568, 6210, 7344, 9300, 11076, 12420, 14076, 14868, 16588, 18972, 27348, 36036, 40612, 47124, 56592, 68310, 72644, 74008, 76496, 87516, 142428, 143572, 149270, 190944, 391820, 430404, 462748, 465036, 784212, 1565564, 5087940}
However some (for instance 2070) do not work.
Egyptian Fractions, Number Theory, David Eppstein, ICS, UC Irvine