From 1d339f74d14d38b8173e43e6983d336c91f8e8d5 Mon Sep 17 00:00:00 2001 From: notoria Date: Wed, 29 Apr 2020 02:44:13 +0200 Subject: [PATCH] Implemented predicate mediants/2 with Stern-Brocot tree --- src/prolog/lib/arithmetic.pl | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/src/prolog/lib/arithmetic.pl b/src/prolog/lib/arithmetic.pl index 406e8473..ccbdc0dd 100644 --- 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 + ). -- 2.54.0