Skip to content

Instantly share code, notes, and snippets.

@jp-diegidio
Last active October 27, 2019 19:13
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save jp-diegidio/eb2d6c1ae432c0e3c75d9f1e5bba602c to your computer and use it in GitHub Desktop.
Save jp-diegidio/eb2d6c1ae432c0e3c75d9f1e5bba602c to your computer and use it in GitHub Desktop.
Numerical Snippets in Prolog
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (SWI-Prolog 7.5.1)
:- module(ilog2,
[ %%
ilog2_n/3, % +X, -Y, -C
ilog2_b/3, % +X, -Y, -C
ilog2_i/3, % +X, -Y, -C
%ilog2_f/3, % +X, -Y, -C
%%
ilog2_print/4, % :G, +L, +H, +X0
ilog2_test/2, % :G, -C
ilog2_case/1 % -X
]).
/** <module> A Collection of Numerical Snippets :: ilog2
This module provides predicates that implement different algorithms for
the computation of the integer binary logarithm in arbitrary precision.
To the purpose of comparing efficiency, predicates in this module return
the count of elementary arithmetic operations performed.
Three variants are offered which are further improvements on the naive
approach that increases the bit mask by one bit at a time:
- ilog2_n/3 : Naive approach.
- ilog2_b/3 : Improved approach that uses plain binary search.
- ilog2_i/3 : Faster approach that reduces the input along the way.
- ilog2_f/3 : Even faster approach that also leverages FP arithmetic. % TODO! #####
*NOTE*: This module is *not* meant to provide production-ready Prolog code.
*/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%! ilog2_n(+X:posint, -Y:nonnegint, -C:nonnegint) is det.
%
% Naive approach.
%
% @arg X Input.
% @arg Y Result.
% @arg C Count of operations.
ilog2_n(X, Y, C) :- % ilog2_n(x, out y) =>
ilog2_n__do(1, 0, X, Y, 0, C). % # _do(1, 0, x, out y)
%C is C1. % (c := c1)
ilog2_n__do(M, K, X, Y, C0, C) :- % _do(m, k, x, out y) =>
X =:= M, !, % if x == m
Y is K, % y := k
C is C0. % (c := c0)
ilog2_n__do(M, K, X, Y, C0, C) :- % else
X < M, !, % if x < m
Y is K - 1, % # y := k - 1
C is C0 + 1. % (c := c0 + 1)
ilog2_n__do(M, K, X, Y, C0, C) :- % else
M1 is M << 1, % # m1 := m * 2
K1 is K + 1, % # k1 := k + 1
C1 is C0 + 2, % (c1 := c0 + 2)
ilog2_n__do(M1, K1, X, Y, C1, C). % # _do(1, 0, x, out y)
%C is C2. % (c := c2)
%! ilog2_b(+X:posint, -Y:nonnegint, -C:nonnegint) is det.
%
% Improved approach that uses plain binary search.
%
% @arg X Input.
% @arg Y Result.
% @arg C Count of operations.
ilog2_b(X, Y, C) :- % ilog2_b(x, out y) =>
ilog2_b__bd(X, 0, 1, L, H, 0, C1), % # _bd(x, 0, 1, out l, out h)
ilog2_b__sc(X, L, H, Y1, C1, C2), % # _sc(x, l, h, out y1)
Y is Y1 - 1, % # y := y1 - 1
C is C2 + 1. % (c := c2 + 1)
ilog2_b__bd(X, L0, H0, L, H, C0, C) :- % _bd(x, l0, h0, out l, out h) =>
B is 1 << H0, % # b := 1 << h0
B >= X, !, % if b >= x
L is L0, % l := l0
H is H0, % h := h0
C is C0 + 1. % (c := c0 + 1)
ilog2_b__bd(X, _, H0, L, H, C0, C) :- % else
L1 is H0, % l1 := h0
H1 is H0 << 1, % # h1 := h0 * 2,
C1 is C0 + 2, % (c1 := c0 + 2)
ilog2_b__bd(X, L1, H1, L, H, C1, C). % # _bd(x, l1, h1, out l, out h)
%C is C2. % (c := c2)
ilog2_b__sc(_, L0, H0, Y, C0, C) :- % _sc(x, l0, h0, out y) =>
L0 >= H0, !, % if l0 >= h0
Y is L0, % y := l0
C is C0. % (c := c0)
ilog2_b__sc(X, L0, H0, Y, C0, C) :- % else
N is L0 + H0, % # n := l0 + h0
M is N >> 1, % # m := n // 2
C1 is C0 + 2, % (c1 := c0 + 2)
ilog2_b__ss(X, L0, H0, M, Y, C1, C). % # _ss(x, l0, h0, m, out y)
%C is C2. % (c := c2)
ilog2_b__ss(X, L0, _, M, Y, C0, C) :- % _ss(x, l0, h0, m, out y) =>
B is 1 << M, % # b := 1 << m
B >= X, !, % if b >= x
L1 is L0, % l1 := l0
H1 is M, % h1 := m
C1 is C0 + 1, % (c1 := c0 + 1)
ilog2_b__sc(X, L1, H1, Y, C1, C). % # _ss(x, l1, h1, out y)
%C is C2. % (c := c2)
ilog2_b__ss(X, _, H0, M, Y, C0, C) :- % else
L1 is M + 1, % # l1 := m + 1
H1 is H0, % h1 := h
C1 is C0 + 2, % (c1 := c0 + 2)
ilog2_b__sc(X, L1, H1, Y, C1, C). % # _sc(x, l1, h1, out y)
%C is C2. % (c := c2)
%! ilog2_i(+X:posint, -Y:nonnegint, -C:nonnegint) is det.
%
% Faster approach that reduces the input along the way.
%
% @arg X Input.
% @arg Y Result.
% @arg C Count of operations.
ilog2_i(X, Y, C) :- % ilog2_i(x, out y) =>
ilog2_i__up(X, 1, M1, X1, 0, C1), % # _up(x, 1, out m1, out x1)
ilog2_i__dn(M1, X1, Y, C1, C). % # _dn(m1, x1, out y)
%C is C2. % (c := c2)
ilog2_i__up(X0, M0, M, X, C0, C) :- % _up(x0, m0, out m, out x) =>
X1 is X0 >> M0, % # x1 := x0 // 2^m0
X1 =\= 0, !, % if x1 =\= 0
M1 is M0 << 1, % # m1 := m0 * 2
C1 is C0 + 2, % (c1 := c0 + 2)
ilog2_i__up(X1, M1, M, X, C1, C). % _up(x1, m1, out m, out x)
%C is C2. % (c := c2)
ilog2_i__up(X0, M0, M, X, C0, C) :- % else
M is M0, % m := m0
X is X0, % x := x0
C is C0 + 1. % (c := c0 + 1)
ilog2_i__dn(M0, 1, Y, C0, C) :- % _dn(m0, x0, out y) =>
M0 =:= 1, !, % if m0 == 1
Y is 0, % y := 0
C is C0. % (c := c0)
ilog2_i__dn(M0, X0, Y, C0, C) :- % else
M1 is M0 >> 1, % # m1 := m0 // 2
Y0 is M0 - 1, % # y0 := m0 - 1
C1 is C0 + 2, % (c1 := c0 + 2)
ilog2_i__do(M1, X0, Y0, Y, C1, C). % _do(m1, x0, y0, out y)
%C is C2. % (c := c2)
ilog2_i__do(M0, X0, Y0, Y, C0, C) :- % _do(m0, x0, y0, out y) =>
X1 is X0 >> M0, % # x1 := x0 // 2^m0
X1 =\= 0, !, % if x1 =\= 0
Y1 is Y0 + M0, % # y1 := y0 + m0
C1 is C0 + 2, % (c1 := c0 + 2)
ilog2_i__do(M0, X1, Y1, Y, C1, C). % _do(m0, x1, y1, out y)
%C is C2. % (c := c2)
ilog2_i__do(M0, X0, Y0, Y, C0, C) :- % else
M0 =\= 1, !, % if m0 =\= 1
M1 is M0 >> 1, % # m1 := m0 // 2
C1 is C0 + 2, % (c1 := c0 + 2)
ilog2_i__do(M1, X0, Y0, Y, C1, C). % _do(m1, x0, y0, out y)
%C is C2. % (c := c2)
ilog2_i__do(1, 1, Y0, Y, C0, C) :- % else
Y is Y0, % y := y0
C is C0 + 1. % (c := c0 + 1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%! ilog2_print(:G:call, +L:posint, +H:posint, +X0:nonnegint) is det.
%
% @arg G Predicate to test.
% @arg L Lower bound for input enumeration.
% @arg H Upper bound for input enumeration.
% @arg X0 Offset for input enumeration.
:- meta_predicate
ilog2_print(3, +, +, +).
ilog2_print(G, L, H, X0) :-
G = _:P,
forall(between(L, H, M),
( X is X0 + 1 << M,
call(G, X, Y, C),
format('~a(~d) = ~d (c=~d)~n', [P, X, Y, C])
)).
%! ilog2_test(:G:call, -C:nonnegint) is det.
%
% @arg G Predicate to test.
% @arg C Total count of operations.
:- meta_predicate
ilog2_test(3, -).
ilog2_test(G, C) :-
Cw = c(0),
forall(ilog2_case(X),
( call(G, X, _, C),
Cw = c(C0),
C1 is C0 + C,
nb_setarg(1, Cw, C1)
)), Cw = c(C).
% ilog2_case(-X:posint) is multi.
%
% @arg X Input enumeration.
ilog2_case(X) :- % ilog2_case(out x) =>
between(0, 100, J), % for j in [0, 100]
X0 is (1 << J), % x0 := 2^j
between(0, 100, I), % for i in [0, 100]
X is X0 + I. % x := x0 + i
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
/*
?- ilog2_print(ilog2_i,100,110,123456789).
ilog2_i(1267650600228229401496826662165) = 100 (c=32)
ilog2_i(2535301200456458802993529867541) = 101 (c=32)
ilog2_i(5070602400912917605986936278293) = 102 (c=34)
ilog2_i(10141204801825835211973749099797) = 103 (c=30)
ilog2_i(20282409603651670423947374742805) = 104 (c=32)
ilog2_i(40564819207303340847894626028821) = 105 (c=32)
ilog2_i(81129638414606681695789128600853) = 106 (c=34)
ilog2_i(162259276829213363391578133744917) = 107 (c=32)
ilog2_i(324518553658426726783156144033045) = 108 (c=34)
ilog2_i(649037107316853453566312164609301) = 109 (c=34)
ilog2_i(1298074214633706907132624205761813) = 110 (c=36)
% _case(out x) =>
% for j in [0, 100]
% for i in [0, 100]
% x := 2^j + i
?- time(ilog2_test(ilog2_n,C)).
?- time(ilog2_test(ilog2_b,C)).
?- time(ilog2_test(ilog2_i,C)).
C = 1053613.
C = 359907.
C = 264235.
% 2,199,366 inferences, 0.219 CPU in 0.233 seconds (94% CPU, 10054245 Lips)
% 902,684 inferences, 0.094 CPU in 0.092 seconds (102% CPU, 9628629 Lips)
% 690,887 inferences, 0.063 CPU in 0.065 seconds (96% CPU, 11054192 Lips)
?- time(forall((N is 2^2^4, between(1,N,X)), ilog2_n(X,_,_))).
?- time(forall((N is 2^2^4, between(1,N,X)), ilog2_b(X,_,_))).
?- time(forall((N is 2^2^4, between(1,N,X)), ilog2_i(X,_,_))).
% 4,390,901 inferences, 0.313 CPU in 0.299 seconds (105% CPU, 14050883 Lips)
% 3,838,521 inferences, 0.281 CPU in 0.282 seconds (100% CPU, 13648075 Lips)
% 3,132,075 inferences, 0.234 CPU in 0.240 seconds (98% CPU, 13363520 Lips)
?- time(forall((M is 2^2^6,N is M+2^2^4, between(M,N,X)), ilog2_n(X,_,_))).
?- time(forall((M is 2^2^6,N is M+2^2^4, between(M,N,X)), ilog2_b(X,_,_))).
?- time(forall((M is 2^2^6,N is M+2^2^4, between(M,N,X)), ilog2_i(X,_,_))).
% 17,498,377 inferences, 2.125 CPU in 2.125 seconds (100% CPU, 8234530 Lips)
% 7,274,583 inferences, 0.750 CPU in 0.745 seconds (101% CPU, 9699444 Lips)
% 4,522,056 inferences, 0.391 CPU in 0.390 seconds (100% CPU, 11576463 Lips)
*/
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment