fraction between 2 number, minimum denominator

[Albert Chan]

From thread Companion program for Joe Horn’s continued fractions article.
Python code, translated for HP Prime (if micro-Python available, that should work too)

Code:

d_convergent(a,b, d1,d2)
BEGIN
LOCAL ia, ib;
ia, ib := floor(a), floor(b);
IF ia == a  THEN RETURN d2, d1+(ia-1)*d2 END;
IF ia != ib THEN RETURN d2, d1+ia*d2 END;
d_convergent(1/(b-ib), 1/(a-ia), d2, d1+ia*d2);
END;

frac_between(a,b)
BEGIN
LOCAL d, d2;
IF a>b THEN a,b := b,a END;
d, d2 := d_convergent(a,b, 1,0);
REPEAT d += d2 UNTIL (a*d) <= floor(b*d);
RETURN ceiling(a*d)/d;
END;

Once minimum denominator are picked, program pick minimum numerator too.
Code assumed number are *exact* fractions, not float.

Cas> frac_between(18/10,19/10) --> 9/5
Cas> frac_between(18/10,18/10) --> 9/5
Cas> frac_between(18/10,17/10) --> 7/4
Cas> frac_between(18/10,16/10) --> 5/3
Cas> frac_between(18/10,15/10) --> 3/2
Cas> frac_between(18/10,14/10) --> 3/2

Cas> frac_between(999999999977/10^12,999999999978/10^12) --> 43478260869/43478260870

[Albert Chan]

(09-16-2023 07:47 PM)Albert Chan Wrote:  Cas> frac_between(18/10,18/10) --> 9/5

Cas> frac_between(1.8, 1.8) practically hang the emulator.

Cas 48 bits float('1.8') = hexfloat(0x1.cccccccccccc) < 9/5

C:\> spigot -C 0x1.cccccccccccc
1/1
2/1
7/4
9/5
126663739519795/70368744177664

Because of rounding errors, d_convergent(1.8, 1.8, 1,0) only return [1,4], d nowhere close enough.

The trick (for float) is to use mod: n/d = floor_div(n,d) + floor_mod(n,d)/d
First term is CF coef, and we simply flip RHS fraction for others.

If n≥0, d>0, (euclidean_mod = floor_mod = trunc_mod), and trunc_mod is *exact*
see Eli Bendersky's blog, Computing remainders by doubling, recursive steps are exact.

Unfortunately, Cas does not have mod for float. We use Lua to illustrate.
Although we work with float, d convergents are exact. No fudge factor is needed.

Code:

function d_convergent(a,a2, b,b2, d,d2) -- 0 <= (a/a2) <= (b/b2)
    local a0, b0 = a, b
    a, b = a%a2, b%b2 -- map a/a2 to a0+a/a2
    a0, b0 = (a0-a)/a2, (b0-b)/b2 -- CF coef
    d = d + a0*d2
    if a == 0   then return d2, d-d2 end
    if a0 ~= b0 then return d2, d end
    return d_convergent(b2,b, a2,a, d2,d)
end

Code:

function frac_between(a,b)      -- assumed non-negative a,b
    if a>b then a,b = b,a end   -- 0 <= a <= b
    local d, d2 = d_convergent(a,1, b,1, 1,0)
    repeat d = d+d2 until (a*d) <= floor(b*d)
    return ceil(a*d), d
end

lua> x = 0x1.cccccccccccc -- = Cas float("1.8")
lua> d_convergent(x,1, x,1, 1,0)
5      70368744177659
lua> frac_between(x,x) -- = fraction(x)
126663739519795      70368744177664

lua> frac_between(0.999999999977, 0.999999999978)
43478173323      43478173324

Lua only see binary float, not decimal (*) In that sense, above fraction is correct.

To reduce decimal to binary conversion error, instead of (a,b), lets do (1-a, 1-b)
Note: numbers have 10 9's after decimal point.

lua> frac_between(0.23E-10, 0.22E-10)
1      43478260870

Or, we enter (a,b), but input as actual fraction, to remove conversion error.

lua> d_convergent(999999999977,1E12, 999999999978,1E12, 1,0)
1      43478260869

First d semi-convergent = 1 + 43478260869 = 43478260870

(*) I do have Decimal/Fraction version, lua + qmath

C:\>run frac_between.lua 0.999999999977 0.999999999978
43478260869/43478260870

C:\>run frac_between.lua 10/7 13/9
10/7

[Albert Chan]

(09-18-2023 11:28 PM)Albert Chan Wrote:  Lua only see binary float, not decimal (*) In that sense, above fraction is correct.

If machine only see Decimal float, there is no bin↔dec conversion errors to worry about.
Within machine precision limits, we could use fast float, instead of exact integers.

Here is HP71B code, fraction = (numerator, denominator (default=1)).

Quote:10 DESTROY ALL @ COMPLEX A,B
20 INPUT "A, B= ";A,B
30 A=A+(0,IMPT(A)=0) @ B=B+(0,IMPT(B)=0)
40 IF IMPT(A*CONJ(B))<0 THEN VARSWAP A,B
50 D=1 @ D2=0 @ CALL DCONV(A,B,D,D2) @ "DCONV= ";D;D2,
60 REPEAT @ D=D+D2 @ UNTIL REPT(A*(D,IP(D*REPT(B)/IMPT(B))))<=0
70 "FRAC = "; (CEIL(D*REPT(A)/IMPT(A)), D)
80 GOTO 20
100 SUB DCONV(A,B,D,D2)
110 A1=REPT(A) @ A2=IMPT(A) @ B1=REPT(B) @ B2=IMPT(B)
120 A0=A1 @ A1=MOD(A1,A2) @ A0=(A0-A1)/A2
130 B0=B1 @ B1=MOD(B1,B2) @ B0=(B0-B1)/B2
140 D=D+(A0-(A1=0))*D2 @ VARSWAP D,D2 ! NEXT CONVERGENT
150 IF A1 AND A0=B0 THEN VARSWAP A1,B2 @ VARSWAP A2,B1 @ GOTO 120

Again, we assumed non-negative inputs.

>run
A, B= 1.8, 1.9
DCONV=  1  4                        FRAC = (9,5)
A, B= 1.8, 1.8
DCONV=  1  4                        FRAC = (9,5)
A, B= 1.8, 1.7
DCONV=  1  3                        FRAC = (7,4)
A, B= 1.8, 1.6
DCONV=  1  2                        FRAC = (5,3)
A, B= 1.8, 1.5
DCONV=  1  1                        FRAC = (3,2)
A, B= (10,7), (13,9)
DCONV=  2  5                        FRAC = (10,7)
A, B= 3.1415, 3.1416
DCONV=  7  99                      FRAC = (333,106)
A, B= .999999999977, .999999999978
DCONV=  1  43478260869      FRAC = (43478260869,43478260870)

Translated code for other decimal machine welcome ...