I am sharing with you the “raising to the power” code that I found earlier this week. This is a mechanical translation from the C++ code found in the book “From Mathematics to Generic Programming” by Stepanov and Rose. (“Mechanical” means I just translated it manually without thinking too much about what it means. It looks like I mostly kept the original names. This code is very procedural.)
EDIT: An earlier book by the same authors has been made freely available online: “Elements of Programming” by Stepanov and Rose. Chapter 3, “Associative Operations” is where this is covered in detail. For me personally the lessons are about the process of writing code, not about any particular algorithm or computational complexity. Anyway…
The idea as I understand it with my limited background is that if an operation is associative, as in:
(AB)C = A(BC)
… and you want to apply an operation to the same element N times in a row, you can “raise to the power” with less than N operations. It isn’t an optimal algorithm for any N (for example N=15), but it is on average quite good – O(log n) instead of O(n).
I wrote it some years ago and forgot about it. It was useful in solving this problem. I don’t want to make a pack/add-on out of it because I don’t have the time. I don’t even feel like documenting it. If anyone wants to do anything at all with the code they have my permission.
There are 3 exported predicates. See the examples on how to use it (in addition to the “bitmask” example in the linked post). The arguments are:
power(Element,
Power,
Operation,
[Identity_element, [Inverse_operation,]]
Result)
You need the identity element if you want to raise to the 0th power, and you need the inverse operation if you want to raise to negative powers. The inverse operation is evaluated only once.
Examples
/* Fibonacci, of course ;-)
* The sequence is 0,1,1,2,3,5,...
* so the 0th Fibonacci number is 0, first is 1, and so on
*/
fib(N, F) :-
( N == 0
-> F = 0
; power(fib(1,0), N, fib_mult, fib(F, _))
).
fib_mult(fib(X1, X2), fib(Y1, Y2), fib(Z1, Z2)) :-
Z1 is X1 * (Y1 + Y2) + X2 * Y1,
Z2 is X1 * Y1 + X2 * Y2.
… and then:
?- fib(3, F).
F = 2.
?- fib(123, F).
F = 22698374052006863956975682.
?- time( fib(1_000_000, _) ).
% 119 inferences, 0.019 CPU in 0.019 seconds (99% CPU, 6333 Lips)
true.
This works because the Fibonacci sequence is a linear recurrence, said the book. The book also suggested that I implement transitive closure and path finding in a directed graph by raising the adjacency matrix to the power of N-1 where N is the number of nodes but I didn’t get as far.
Basic arithmetic:
?- power(3, 4, [X,Y,Z]>>(Z is X + Y), Result). % 3 * 4
Result = 12.
?- power(3, 4, [X,Y,Z]>>(Z is X * Y), Result). % 3^4
Result = 81.
?- power(2, -4, [X,Y,Z]>>(Z is X * Y), 1, [X,Y]>>(Y is 1/X), Result). % 1 / 2^4
Result = 0.0625.
The module
$ cat power.pl
:- module(power, [power/4,
power/5,
power/6]).
:- meta_predicate power(+,+,3,-),
power(+,+,3,+,-),
power(+,+,3,+,2,-).
goal_expansion(odd(N), N /\ 0x1 =:= 1).
goal_expansion(half(N, Half), Half is N >> 1).
power_accumulate_semigroup(R, A, N, Op, P) :-
( N == 0
-> P = R
; power_accumulate_semigroup_1(R, A, N, Op, P)
).
power_accumulate_semigroup_1(R, A, N, Op, P) :-
( odd(N)
-> power_accumulate_semigroup_oddn(R, A, N, Op, P)
; power_accumulate_semigroup_anyn(R, A, N, Op, P)
).
power_accumulate_semigroup_oddn(R, A, N, Op, P) :-
call(Op, R, A, R1),
( N == 1
-> P = R1
; power_accumulate_semigroup_anyn(R1, A, N, Op, P)
).
power_accumulate_semigroup_anyn(R, A, N, Op, P) :-
half(N, N1),
call(Op, A, A, A1),
power_accumulate_semigroup_1(R, A1, N1, Op, P).
power(A, N, Op, P) :-
must_be(positive_integer, N),
power_semigroup_(A, N, Op, P).
power_semigroup_(A, N, Op, P) :-
power_semigroup_evenn(A, N, Op, A0, N0),
( N0 == 1
-> P = A0
; call(Op, A0, A0, A1),
half(N0 - 1, N1),
power_accumulate_semigroup(A0, A1, N1, Op, P)
).
power_semigroup_evenn(A, N, Op, A1, N1) :-
( odd(N)
-> A = A1,
N = N1
; call(Op, A, A, A0),
half(N, N0),
power_semigroup_evenn(A0, N0, Op, A1, N1)
).
power(A, N, Op, Id, P) :-
must_be(nonneg, N),
power_monoid_(A, N, Op, Id, P).
power_monoid_(A, N, Op, Id, P) :-
( N == 0
-> P = Id
; power_semigroup_(A, N, Op, P)
).
power(A, N, Op, Id, Inverse_op, P) :-
must_be(integer, N),
power_group_(A, N, Op, Id, Inverse_op, P).
power_group_(A, N, Op, Id, Inverse_op, P) :-
( N < 0
-> N1 is -N,
call(Inverse_op, A, A1)
; N1 = N,
A1 = A
),
power_monoid_(A1, N1, Op, Id, P).