Setting a maximum demoninator for fractions

[komame]

(10-04-2023 01:03 AM)Albert Chan Wrote:  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;

Wow, this is incredibly fast regardless of the denominator value. Even for a case like PI with dmax=30,000,000,000; the result is returned instantly.

PK