Setting a maximum demoninator for fractions

[Albert Chan]

(10-03-2023 02:44 PM)Albert Chan Wrote:  if abs((n1-d1*x)*q) <= abs((p-q*x)*d1) THEN p:=n1; q:=d1; END;

This test has trouble if numerator is huge, with possibly catastrophic cancellation issue.
We like numerator small, so my next version remove rounded integer part of x, before the test.

I also squared both sides to remove abs, with equivalent test as product of 2 dot products.
HP Prime HOME dot product work with higher precisions, thus more accurate.

Cas> factor(((p-q*x)*d1)^2 - ((n1-d1*x)*q)^2)      → (q*n1-p*d1)*(2*x*q*d1-q*n1-p*d1)

Code:

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