Implemented predicate mediants/2 with Stern-Brocot tree

diff --git a/src/prolog/lib/arithmetic.pl b/src/prolog/lib/arithmetic.pl
index 406e847312e6da6545d753e244d5a0972e39eb27..ccbdc0dd710e27a636d2a209ef6c07e16cbad0ef 100644 (file)
--- a/src/prolog/lib/arithmetic.pl
+++ b/src/prolog/lib/arithmetic.pl
@@ -1,4 +1,4 @@
-:- module(arithmetic, [msb/2, lsb/2]).
+:- module(arithmetic, [lsb/2, msb/2, mediants/2]).
 
 :- use_module(library(error)).
 
@@ -25,3 +25,28 @@ msb_(X, M, N) :-
     X1 is X >> 1,
     M1 is M + 1,
     msb_(X1, M1, N).
+
+mediants_precision(1.0e-15).
+
+mediants(Real, Fraction) :-
+    builtins:must_be_number(Real, mediants),
+    Fraction = P/Q,
+    mediants_precision(Eps),
+    R is abs(Real),
+    stern_brocot(Eps, R, P1/Q),
+    P is sign(Real) * P1.
+
+% If 0 <= Eps <= 1e-16 then the search is for "infinite" precision.
+stern_brocot(Eps, Real, Fraction) :-
+    Rn is Real - Eps,
+    Rp is Real + Eps,
+    stern_brocot_(Rn, Rp, 0/1, 1/0, Fraction).
+
+stern_brocot_(Rn, Rp, A/B, C/D, Fraction) :-
+    M0 is A + C,
+    M1 is B + D,
+    M is M0 / M1,
+    (   M < Rn -> stern_brocot_(Rn, Rp, M0/M1, C/D, Fraction)
+    ;   M > Rp -> stern_brocot_(Rn, Rp, A/B, M0/M1, Fraction)
+    ;   Fraction = M0/M1
+    ).