Removed predicate mediants/2 and added stern_brocot/3

src/prolog/lib/arithmetic.pl

@@ -1,4 +1,4 @@
-:- module(arithmetic, [lsb/2, msb/2, mediants/2]).
+:- module(arithmetic, [lsb/2, msb/2, stern_brocot/3]).
 
 :- use_module(library(error)).
 
@@ -26,21 +26,65 @@ msb_(X, M, N) :-
     M1 is M + 1,
     msb_(X1, M1, N).
 
-mediants_precision(1.0e-15).
+stern_brocot(E0/E1, Fraction0, Fraction) :-
+    P1/Q1 = Fraction0,
+    !,
+    (   \+ integer(E0) -> type_error(integer, E0, stern_brocot/3)
+    ;   \+ integer(E1) -> type_error(integer, E1, stern_brocot/3)
+    ;   \+ integer(P1) -> type_error(integer, P1, stern_brocot/3)
+    ;   \+ integer(Q1) -> type_error(integer, Q1, stern_brocot/3)
+    ;   S is sign(E0) * sign(E1), S < 0 ->
+            domain_error(not_less_than_zero, E0/E1, stern_brocot/3)
+    ;   P2 is abs(P1),
+        Q2 is abs(Q1),
+        Qn1n is P2 * E1 - Q2 * E0,
+        Qn1d is Q2 * E1,
+        simplify_fraction(Qn1n/Qn1d, Qn1),
+        Qp1n is P2 * E1 + Q2 * E0,
+        Qp1d = Qn1d,
+        simplify_fraction(Qp1n/Qp1d, Qp1),
+        fraction_stern_brocot_(Qn1, Qp1, 0/1, 1/0, P3/Q3),
+        P4 is sign(P1) * sign(Q1) * P3,
+        Fraction = P4/Q3
+    ).
 
-mediants(Real, Fraction) :-
+% If 0 <= Eps0 <= 1e-16 then the search is for "infinite" precision.
-    builtins:must_be_number(Real, mediants),
+stern_brocot(Eps0, Real0, Fraction) :-
-    Fraction = P/Q,
+    (   Real0 = R1/R2 ->
-    mediants_precision(Eps),
+            builtins:must_be_number(R1, stern_brocot/3),
-    R is abs(Real),
+            builtins:must_be_number(R2, stern_brocot/3),
-    stern_brocot(Eps, R, P1/Q),
+            Real is R1/R2
-    P is sign(Real) * P1.
+    ;   builtins:must_be_number(Real0, stern_brocot/3),
+        Real = Real0
+    ),
+    (   Eps0 = E0/E1 ->
+            builtins:must_be_number(E0, stern_brocot/3),
+            builtins:must_be_number(E1, stern_brocot/3),
+            % Eps is Eps0 % doesn't work
+            Eps is E0/E1
+    ;   Eps = Eps0
+    ),
+    S is sign(Eps),
+    (   S < 0 -> domain_error(not_less_than_zero, Eps0, stern_brocot/3)
+    ;   Rn is abs(Real) - Eps,
+        Rp is abs(Real) + Eps,
+        stern_brocot_(Rn, Rp, 0/1, 1/0, P1/Q),
+        P is sign(Real) * P1,
+        Fraction = P/Q
+    ).
 
-% If 0 <= Eps <= 1e-16 then the search is for "infinite" precision.
+fraction_stern_brocot_(Qnn/Qnd, Qpn/Qpd, A/B, C/D, Fraction) :-
-stern_brocot(Eps, Real, Fraction) :-
+    Fn1 is A + C,
-    Rn is Real - Eps,
+    Fd1 is B + D,
-    Rp is Real + Eps,
+    simplify_fraction(Fn1/Fd1, Fn/Fd),
-    stern_brocot_(Rn, Rp, 0/1, 1/0, Fraction).
+    S1 is sign(Fn * Qnd - Fd * Qnn),
+    S2 is sign(Fn * Qpd - Fd * Qpn),
+    (   S1 < 0 ->
+            fraction_stern_brocot_(Qnn/Qnd, Qpn/Qpd, Fn/Fd, C/D, Fraction)
+    ;   S2 > 0 ->
+            fraction_stern_brocot_(Qnn/Qnd, Qpn/Qpd, A/B, Fn/Fd, Fraction)
+    ;   Fraction = Fn/Fd
+    ).
 
 stern_brocot_(Rn, Rp, A/B, C/D, Fraction) :-
     M0 is A + C,
@@ -50,3 +94,8 @@ stern_brocot_(Rn, Rp, A/B, C/D, Fraction) :-
     ;   M > Rp -> stern_brocot_(Rn, Rp, A/B, M0/M1, Fraction)
     ;   Fraction = M0/M1
     ).
+
+simplify_fraction(A0/B0, A/B) :-
+    G is gcd(A0, B0),
+    A is A0 div G,
+    B is B0 div G.