Setting a maximum demoninator for fractions

[Albert Chan]

(10-04-2023 04:59 PM)Albert Chan Wrote:  We pick smaller denominator convergent, if product of 2 DOT's ≥ 0:

s * (s + 2*q*(n1 - d1*x)) ≥ 0
1 + 2*s*q*(n1 - d1*x) ≥ 0
s*q*(n1 - d1*x) = s*q*(n2 - (d1-d2)*x) ≥ -1/2

But wait! s = sign(x - n1/d1), we don't need to track it!
LHS always negative, we multiply by -1, for equivalent test.

abs(q*(n2 - (d1-d2)*x)) ≤ 1/2

CF, version 4

Code:

EXPORT CF(x, dmax)  // closest fraction
BEGIN
  LOCAL p,q, n1,d1, n2,d2, k;
  p  := x; q  := 1; // x = p / q
  d1 := 1; d2 := 0; // d convergent
  WHILE d2 <= dmax DO
    IF q==0 THEN RETURN [ROUND(x*d2,0),d2]; END;
    k:=p; p:=(p MOD q); k:=(k-p)/q;
    k:=d1+k*d2; d1:=d2; d2:=k;
    k:=p; p:=q; q:=k;       // SWAP p,q
  END;
  k:=ROUND(x,0); x:=x-k;    // reduce numerator
  n1:=ROUND(x*d1,0);        // convergent = [n1,d1]
  q:=d2-CEILING((d2-dmax)/d1)*d1;
  p:=ROUND(x*q,0);          // semi-convergent = [p,q]
  n2:=(n1 MOD x); d2:=(n1-n2)/x;    // n1 = n2 + x*d2
  IF abs(q*(n2-(d1-d2)*x)) <= .5 THEN p:=n1; q:=d1; END;
  RETURN [k*q+p,q];
END;