# "Power" implementation

**URL:** <https://swi-prolog.discourse.group/t/power-implementation/1937>\
**Category:** Resources\
**Created:** [February 28, 2020, 7:26pm UTC](https://swi-prolog.discourse.group/t/power-implementation/1937 "2020-02-28T19:26:47Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Boris](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/boris/32/7486_2.png) [@Boris](https://swi-prolog.discourse.group/u/Boris)\
**Post date:** [February 28, 2020, 7:26pm UTC](https://swi-prolog.discourse.group/t/power-implementation/1937/1 "2020-02-28T19:26:47Z")

</div>

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](http://elementsofprogramming.com/eop.pdf). 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](https://swi-prolog.discourse.group/t/effective-implementation-of-a-set/1713/8). 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:

```prolog
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

```prolog
/* 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:

```prolog
?- 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:

```prolog
?- 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

```prolog
$ 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).

```

---

<div class="post-metadata">

**Author:** ![EricGT](https://avatars.discourse-cdn.com/v4/letter/e/f1d935/32.png) [@EricGT](https://swi-prolog.discourse.group/u/EricGT)\
**Post date:** [February 28, 2020, 8:02pm UTC](https://swi-prolog.discourse.group/t/power-implementation/1937/2 "2020-02-28T20:02:11Z")

</div>

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)](https://hal.archives-ouvertes.fr/hal-02070778/document) 😉

---

<div class="post-metadata">

**Author:** ![Boris](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/boris/32/7486_2.png) [@Boris](https://swi-prolog.discourse.group/u/Boris)\
**Post date:** [February 29, 2020, 4:38am UTC](https://swi-prolog.discourse.group/t/power-implementation/1937/3 "2020-02-29T04:38:39Z")

</div>

This is too heavy reading for me 🙂

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 🙂

---

<div class="post-metadata">

**Author:** ![Boris](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/boris/32/7486_2.png) [@Boris](https://swi-prolog.discourse.group/u/Boris)\
**Post date:** [February 29, 2020, 9:40pm UTC](https://swi-prolog.discourse.group/t/power-implementation/1937/5 "2020-02-29T21:40:45Z")

</div>

> [@anon95304481](#):
>
> Where did you get `fib_mult` from? What does it do?

It does a half-hearted matrix multiplication. [I got it from here](http://elementsofprogramming.com/eop.pdf), 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:

```prolog
?- 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):

```prolog
?- 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.

```

---

<div class="post-metadata">

**Author:** ![Boris](https://yyz2.discourse-cdn.com/free1/user_avatar/swi-prolog.discourse.group/boris/32/7486_2.png) [@Boris](https://swi-prolog.discourse.group/u/Boris)\
**Post date:** [February 29, 2020, 9:50pm UTC](https://swi-prolog.discourse.group/t/power-implementation/1937/7 "2020-02-29T21:50:07Z")

</div>

> [@anon95304481](#):
>
> your power framework

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!