Why this algorithm for egyptian fractions doesn't terminate in ~$2$% cases?

I thought up yet another algorithm for egyptian fraction expansion which turned out to be very effective (in terms of the length and the denominator size) - in most cases. However, for some fractions it doesn't terminate at all - it leads to an infinite loop.

Here is the algorithm:

Let $\frac{p}{q}<1$ and $p,q$ coprime. Find the minimal $m$ such that $q/(p+m)$ is an integer. We only need to consider $m \in [1,q-p]$. Represent the fraction as: $$\frac{p}{q}=\frac{p+m}{q}-\frac{m}{q}$$

Now to obtain a positive term instead of a negative one, we split the first fraction in two: $$\frac{p+m}{q}-\frac{m}{q}=\frac{p+m}{2q}+\frac{p+m}{2q}-\frac{m}{q}=\frac{p+m}{2q}+\frac{p-m}{2q}$$

Here is a conditional: if $p<m$, then repeat the previous step (add $\frac{p+m}{4q}$ to the negative fraction). And so on, until both terms are positive. If $p>m$ then $\frac{p}{q} \to \frac{p-m}{2q}$ and we repeat the first step of the algorithm.

The working name is complementary method, so I will use CM to denote it for now.

Despite its simplicity (it's not at all obvious why we are dividing by $2$ instead of using some other way to expand the first term) the algorithm works very well. In a lot of cases it beats every other algorithm I tried.

Since the greedy algorithm and Engel expansion are usually bad in terms of denominator size, I used two other methods to compare: Binary Remainder Method and my own 'Splitting-Joining method' (the details can be found in my Mathematica SE question. I also compared it to a modification of Engel proposed by Daniel Fischer in this answer and CM is mostly better for the examples he provided.

Some examples of the best results (a sequence of denominators is provided in each case):

4/49: CM {14,98}; BR {16,98,196,392,784}; SJ {13,325,925,1813}

3/35: CM {14,70}; BR {16,70,140,560}; SJ {20,28}

47/104: CM {4,8,13}; BR {4,8,16,104,208}; SJ {4, 14, 26, 28, 52, 70, 104, 130, 182}

Some examples of the worst results (but still valid - algorithm terminates):

94/191: CM {4, 8, 16, 32, 64, 256, 512, 1024, 2048, 4096, 8192, 24448, 48896, 97792, 195584, 391168, 782336, 1564672}

65/157: CM {4, 8, 32, 256, 512, 1024, 2048, 4096, 10048, 20096, 40192, 80384, 160768, 643072}

52/139: CM {4, 16, 32, 64, 128, 278, 556, 1112, 2224, 8896, 17792}

However, in these cases both BR and SJ methods also give long expansions with large denominators.

Now the real problem which I'm trying to solve - why in some cases the algorithm doesn't terminate, but leads to loops? From large scale experiments I estimates the proportion of such fractions to be about $1.8$% (for numerators and denominators below $1000$).

The examples of such 'bad' fractions are:

$$\frac{41}{111},\frac{5}{87},\frac{8}{87},\frac{14}{87},\frac{47}{87},\frac{61}{102},\frac{17}{69},\frac{33}{119},\frac{38}{93},\frac{77}{177},\frac{32}{57},\frac{99}{185},\frac{98}{141},\frac{100}{129},\dots$$

The most common denominator is $87$ for some reason. Note that not all of the denominators and/or numerators are prime.

The problem can be solved by using $\frac{p-1}{q}$ instead, but not in every case, for example $7/87$ doesn't work either. However, both $6/87$ and $2/87$ work, and give different denominators, so we can expand $8/87$ after all.

I think the problem might be related to the use of the expansion $1=1/2+1/2$ to divide the fist term. However, when I tried some other schemes, I didn't get good results (for example, I've got repeating fractions when using $1=1/3+2/3$).

The working Mathematica code for the algorithm is:

x=6/87;
p0=Numerator[x];
q0=Denominator[x];
S=0;
Nm=100;
a=Table[1,{k,1,Nm}];
m=Table[1,{k,1,Nm}];
p1=p0;
q1=q0;
j=1;
While[Abs[p0]>1&&j<=Nm&&q0<10^35,M=Catch[Do[If[FractionalPart[q0/(p0+k)]<1/10^55,Throw[k]],{k,0,q0-p0}]];
m[[j]]=M;
a[[j]]=(p0+M)/(2 q0);
p1=Numerator[a[[j]]-M/q0];
q1=Denominator[a[[j]]-M/q0];
While[p1<0,a[[j]]=a[[j]]/2;
p1=Numerator[a[[j]]+p1/q1];
q1=Denominator[a[[j]]+p1/q1]];
If[a[[j]]!=1,S+=a[[j]]];
j++;
p0=p1;q0=q1];
a[[j]]=p1/q1;
S+=a[[j]];
Denominator[Table[a[[k]],{k,1,j}]]

And the second question: how to modify the algorithm so it always terminates?

Update:

Among the first $10000$ fractions with $p \neq 1$ in lexicographic order there are $269$ fractions for which this algorithm doesn't work. (Seems to be more than $2$%). They are:

5/33,5/51,2/55,32/55,4/57,7/57,13/57,23/57,32/57,6/65,43/66,4/69,8/69,11/69,17/69,40/69,50/69,8/85,59/85,4/87,5/87,7/87,8/87,14/87,34/87,47/87,62/87,65/87,76/87,5/93,7/93,10/93,19/93,38/93,50/93,67/93,6/95,8/95,63/95,61/102,9/110,59/110,4/111,7/111,8/111,13/111,16/111,22/111,25/111,31/111,41/111,44/111,59/111,62/111,68/111,82/111,7/114,65/114,71/114,83/114,103/114,6/115,11/115,17/115,63/115,3/119,5/119,10/119,15/119,16/119,33/119,37/119,45/119,61/119,66/119,67/119,71/119,73/119,78/119,96/119,101/119,4/123,5/123,8/123,10/123,11/123,17/123,20/123,26/123,29/123,35/123,46/123,49/123,67/123,70/123,76/123,86/123,92/123,4/129,5/129,10/129,13/129,14/129,19/129,28/129,31/129,37/129,47/129,53/129,71/129,74/129,80/129,91/129,100/129,77/130,53/132,119/132,11/138,31/138,77/138,85/138,91/138,103/138,4/141,5/141,7/141,8/141,14/141,16/141,17/141,23/141,32/141,35/141,41/141,52/141,55/141,74/141,79/141,82/141,88/141,98/141,101/141,110/141,121/141,3/143,7/143,21/143,40/143,42/143,60/143,73/143,80/143,98/143,120/143,138/143,6/145,8/145,13/145,21/145,64/145,79/145,93/145,122/145,6/155,7/155,9/155,12/155,14/155,69/155,99/155,102/155,107/155,131/155,5/159,7/159,10/159,11/159,14/159,19/159,20/159,23/159,32/159,38/159,58/159,64/159,83/159,85/159,91/159,113/159,116/159,125/159,136/159,9/161,11/161,101/161,103/161,16/165,41/165,61/165,116/165,151/165,33/170,101/170,7/174,37/174,43/174,65/174,95/174,97/174,101/174,103/174,115/174,155/174,8/175,11/175,78/175,108/175,111/175,113/175,116/175,148/175,4/177,5/177,8/177,10/177,11/177,13/177,17/177,19/177,20/177,22/177,26/177,29/177,35/177,38/177,44/177,64/177,67/177,70/177,77/177,94/177,95/177,97/177,103/177,122/177,128/177,131/177,137/177,140/177,154/177,4/183,5/183,7/183,10/183,11/183,13/183,14/183,19/183,20/183,22/183,28/183,34/183,37/183,40/183,49/183,65/183,68/183,71/183,74/183

Update 2 (Important)

The question was deleted for a time, because for some of the listed fraction algorithm seems to work just fine if it's done by hand. There is some error in my code, which I hadn't been able to find yet.

But there are fractions which lead to loops by hand as well (such as $41/111$) so the question still stands.

Solution 1:

There's still an error in your code, in the While loop you take if $p<m$:

While[p1<0,a[[j]]=a[[j]]/2;
p1=Numerator[a[[j]]+p1/q1];
q1=Denominator[a[[j]]+p1/q1]];

The problem is that after you update p1, the calculation q1 = Denominator[a[[j]] + p1/q1] is done with the old value of q1 but the new value of p1. A quick fix: change the loop to

While[p1<0,a[[j]]=a[[j]]/2;
{p1,q1}={Numerator[a[[j]]+p1/q1], Denominator[a[[j]]+p1/q1]}];

Here's my own implementation of this algorithm, in convenient function-call form: just run egyptian[41/111] to get back {1/4,1/12,1/48,1/74,1/592}. I've tested it on all fractions with denominator less than $200$, and it works correctly:

split[Rational[p_,q_]]:=Module[{m, t=1},
  m=First@Select[Range[q-p],IntegerQ[q/(p+#)]&,1];
  While[(2^t-1)p < m,t++];
  Return[{(p+m)/(2^t q),((2^t-1)p-m)/(2^t q)}];
];

egyptian[fraction_Rational]:=Module[{first,rest},
  If[Numerator[fraction]==1,Return[{fraction}]];
  {first,rest} = split[fraction];
  Return[Prepend[first] @ egyptian[rest]];
];

(Here, split is a helper function that returns the two parts of the representation of $\frac pq$ as $\frac{p+m}{2q} + \frac{p-m}{2q}$ or as $\frac{p+m}{4q} + \frac{3p-m}{4q}$ or in general as $\frac{p+m}{2^t q} + \frac{(2^t-1)p-m}{2^t q}$ for whichever $t$ is the first to work. Then egyptian just repeatedly applies this to get the Egyptian fraction we wanted.)