"Power" implementation

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).
5 Likes

Did you know that there is something similar for multiplication in that it solves it in O(n log n)

Integer multiplication in time O(n log n) :wink:

This is too heavy reading for me :slight_smile:

Note that for the code I shared, only the power you are raising to needs to be an integer – the element can be any “type”.

I honestly failed to see how it can be useful to me when I read the book and translated the code to Prolog. However, now that I found one use for it, I know for certain that it isn’t completely useless :slight_smile:

Try this, fun with rational numbers:

fib(X,Y) :- Y is 1 rdiv (1 + X).

Then with (^)/4:

?- call(fib^3, 0, X).
X = 2r3.
?- call(fib^123, 0, X).
X = 22698374052006863956975682r36726740705505779255899443.

But I haven’t figured out some fib_mult based on rational numbers.
Where did you get fib_mult from? What does it do?

P.S.: (^)/4 is a two liner:

^(_, 0, X, Y) :- !, Y = X.
^(F, N, X, Y) :- M is N-1, ^(F, M, X, H), call(F, H, Y).

It does a half-hearted matrix multiplication. I got it from here, look at pp 41-43.

Again, just because I read the book and typed in the code doesn’t mean I understand what is really going on. But it seems that your two-liner (^)/4 is just a applying (calling) the operation N times, is that right?

I personally do not have any use in generating Fibonacci numbers, to me it is more like a parlor trick. I don’t know if you noticed the numbers when you time it.

(skipping the results because they have many digits!)

Here, with your code:

?- time( call(fib^123, 0, X) ).
% 371 inferences, 0.000 CPU in 0.000 seconds (99% CPU, 1114124 Lips)
X = 22698374052006863956975682r36726740705505779255899443.

?- time( call(fib^1_000, 0, _) ).
% 3,001 inferences, 0.006 CPU in 0.006 seconds (99% CPU, 540129 Lips)
true.

?- time( call(fib^10_000, 0, _) ).
% 30,001 inferences, 0.056 CPU in 0.056 seconds (99% CPU, 536734 Lips)
true.

Compare this to using the “power” (which, again, implements an algorithm that was described not later than 3500 BC by the Egyptians, and which I know you know as well):

?- time( fib(123, X) ).
% 70 inferences, 0.000 CPU in 0.000 seconds (96% CPU, 638476 Lips)
X = 22698374052006863956975682.

?- time( fib(1_000, _) ).
% 91 inferences, 0.000 CPU in 0.000 seconds (96% CPU, 774930 Lips)
true.

?- time( fib(10_000, _) ).
% 119 inferences, 0.000 CPU in 0.000 seconds (98% CPU, 465983 Lips)
true.

?- time( fib(1_000_000, _) ).
% 171 inferences, 0.038 CPU in 0.042 seconds (90% CPU, 4517 Lips)
true.

Ok, I found a rational number variant, based on the usual equations.
First lets get invert fib first, just find the inverse of f(x)=1/(1+x):

fib_inv(Y,X) :- X is 1 rdiv Y-1.

Some sanity check whether its really an invert:

?- call(fib_inv^3, 2r3, X).
X = 0.
?- call(fib_inv^123, 22698374052006863956975682r36726740705505779255899443, X).
X = 0.

Now how about this:

fib_mul(X,Y,Z) :-
        J is X*Y,
        Z is (X+Y-J) rdiv (1+J).

Again some sanity check:

?- call(fib^3, 0, X).
X = 2r3.
?- call(fib^7, 0, X).
X = 13r21.
?- fib_mul(2r3, 13r21, X).
X = 55r89.
?- call(fib^10, 0, X).
X = 55r89. 

Now you got everything to plug into your power framework which should allow you to compute in O(log n) time instead of O(n) time.

No, it isn’t mine, really. I literally typed it in and was impressed. Then I forgot about it because I never had use for it. Then, a few years later, I found one use for it, so I thought I might as well share the code in case someone else might want to play with it.

So please go ahead!

I am using a trick now, if I assume that step and mul have the
same name and only differ in arity, I could do this here:

^(_, 0, X, Y) :- !, Y = X.
^(F, M, X, Y) :- M rem 2 =:= 0, !, N is M//2, ^(F, N, X, H), call(F, H, H, Y).
^(F, M, X, Y) :- N is M-1, ^(F, N, X, H), call(F, H, Y).

And then with:

fib(X,Y) :- Y is 1 rdiv (1+X).
fib(X,Y,Z) :- J is X*Y, Z is (X+Y-J) rdiv (1+J).

I get again:

?- call(fib^3, 0, X).
X = 2r3.
?- call(fib^123, 0, X).
X = 22698374052006863956975682r36726740705505779255899443.

But this rather rare that this arity juggling is possible.
Already incorporating fib_inv doesn’t work anymore.

Another solution would be to use Logtalk objects.