Algorithms for Egyptian Fractions

Brute Force Methods

The Generalized Remainder Method

This is closely related to the binary remainder method, and to similar methods cited in the discussion there. We find a number p with many divisors, compute xp=qy+r, represent q and r as sums of divisors of p, and use these representations to expand x/y=q/p + r/yp.

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

The Small Multiple Method

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

This method, proposed by Stewart [Ste92], is based on the following idea. If we are looking for any Egyptian fraction representation for q, 1/u + 1/v + ..., then each term must of course be no larger than q itself; so u >= 1/q. However if the terms are in sorted order u > v > ..., and there are a total of k terms, then the largest term 1/u must be at least q/k. So this gives upper and lower bounds for u; we simply try all possibilities and continue recursively. (In the recursive call we have a further bound on the fractions from the assumption that they are generated in sorted order.) Incidentally, since we have a k-level recursion, and each call generates a finite number of recursive calls, the whole call tree must be finite and there are only finitely many k-term representations of q.

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

Formatted by nb2html and filter. Last update: .