Setting a maximum demoninator for fractions

[Albert Chan]

(10-04-2023 02:33 PM)Albert Chan Wrote:  DOT([p,q], [d1,-n1])                = (d1*q) * (p/q - n1/d1)
DOT([p,q,2*x], [d1,n1,-d1*q]) = (d1*q) * ((p/q-x) + (n1/d1-x))

Perhaps we could use MOD again, to improve accuracy greatly!
Let s = DOT([p,q], [d1,-n1]) = ±1, trivial to get (see previous post)

DOT([p,q,2*x], [d1,n1,-d1*q]) = p*d1 + q*n1 - 2*x*d1*q = s + 2*q*(n1 - d1*x)

We decompose n1 using MOD: n1 = n2 + d2*x:
n1/d1 is a convergent of x    → sign(x) = sign(n1)    → trunc/floor mod, n2 = (n1 MOD x) is exact

x ≈ n1/d1
d1*x ≈ n1 = n2 + d2*x
n2 ≈ (d1-d2)*x

Since |n2| < |x|, (d1-d2) = 0 or 1, no need to use dot product.
Either case, (n2 - (d1-d2)*x) is exact (Sterbenz lemma)

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

CF, version 3:

Code:

EXPORT CF(x, dmax)  // closest fraction
BEGIN
  LOCAL p,q, n1,d1, n2,d2, k, s:=1;
  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; s:=-s;
    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 + d2*x
  IF s*q*(n2-(d1-d2)*x) >= -.5 THEN p:=n1; q:=d1; END;
  RETURN [k*q+p,q];
END;

HOME> CF(PI,3000000000)      → [9253131363,2945363191]

Test is now much more accurate: -0.5890726382 ≥ -0.5 is false --> pick semi-convergent

Confirm in Speedcrunch calculator (v 0.12)

x = 0.14159265359
n1 = 10225711
d1 = 72219220
p = 417041790
q = 2945363191

s = p*d1 - q*n1      → -1
s*q*(n1 - d1*x)      → -0.5890726382