Prolog - Finding the Nth Fibonacci number using accumulators

Question

I have this code to generate a list of the Fibonacci sequence in reverse order.

fib2(0, [0]).
fib2(1, [1,0]).
fib2(N, [R,X,Y|Zs]) :-
    N > 1,
    N1 is N - 1,
    fib2(N1, [X,Y|Zs]),
    R is X + Y.

I only need the first element though. The problem is that this code also gives out a false. after the list, so all my attempts at getting the first element have failed. Is there any way I can get that first element in the list, or any other way of calculating the Nth Fibonacci number with accumulators.

Thanks in advance.

Answer 1

I got this logarithmic steps O(log n) solution, and even tail recursive.
Just for fun, it can also compute the n-th Lucas number:

<pre id="in">
fib(N, X) :-
   powmat(N, [[0,1],[1,1]], [[1,0],[0,1]], 
             [[_,X],[_,_]]).
luc(N, Z) :-
   powmat(N, [[0,1],[1,1]], [[1,0],[0,1]], 
             [[X,Y],[_,_]]), Z is 2*X+Y.

powmat(0, _, R, R) :- !.
powmat(N, A, R, S) :- N rem 2 =\= 0, !,
   mulmat(A, R, H), M is N//2, mulmat(A, A, B), powmat(M, B, H, S).
powmat(N, A, R, S) :- 
   M is N//2, mulmat(A, A, B), powmat(M, B, R, S).

mulmat([[A11,A12],[A21,A22]], 
       [[B11,B12],[B21,B22]],
       [[C11,C12],[C21,C22]]) :-
   C11 is A11*B11+A12*B21,
   C12 is A11*B12+A12*B22,
   C21 is A21*B11+A22*B21,
   C22 is A21*B12+A22*B22.

?- fib(100,X).
?- luc(100,X).
</pre>
<script src="http://www.dogelog.ch/lib/exchange.js"></script>

You can compare with:

https://www.wolframalpha.com/input/?i=Fibonacci[100]
https://www.wolframalpha.com/input/?i=LucasN[100]

Edit 28.06.2021:
Here is a very quick explanation why the matrix algorithm works. We only need to show that one step of Fibonacci is linear. Namely that this recurrence relation leads to linear matrix:

     F_{n+2} = F_{n}+F_{n+1}

To see the matrix, we have to assume that the matrix M, transforms a vector b=[Fn,Fn+1] into a vector b'=[F_{n+1}, F_{n+2}]:

      b' = M*b

What could this matrix be? Just solve it:

    |F_{n+1}|   |0*F_{n}+1*F_{n+1}|    |0  1|   |F_{n}  |
    |       | = |                 | =  |    | * |       |
    |F_{n+2}|   |1*F_{n}+1*F_{n+1}|    |1  1|   |F_{n+1}|

Answer 2

It gives out a "false" because Prolog is unsure whether there are more solutions after the first one it provides:

?- fib2(4,L).
L = [3,2,1,1,0] ;  % maybe more solutions? 
false.             % no

This is not a problem: You can tell Prolog that there are indeed no more solutions after the first one (or that you are not interested in seeing them):

?- once(fib2(4,L)).

or

?- fib2(4,L),!.

or you can cut in each of the first clauses, telling Prolog that if the head matches, there is no point trying another clause. This gets rid of the stray "possible solution":

fib2(0, [0])   :- !.
fib2(1, [1,0]) :- !.
fib2(N, [R,X,Y|Zs]) :-
    N > 1,
    N1 is N - 1,
    fib2(N1, [X,Y|Zs]),
    R is X + Y.

What may be a problem is that the given algorithm stores all the fib(i) and performs an addition after the recursive call, which means that Prolog cannot optimize the recursive call into a loop.

For the "accumulator-based" (bottom-up) way of computing fib(N):

% -------------------------------------------------------------
% Proceed bottom-up, without using any cache, or rather a cache
% consisting of two additional arguments.
%
% ?- fib_bottomup_direct(10,F).
% F = 55.
% ------------------------------------------------------------

fib_bottomup_direct(N,F) :-
   N>0,
   !,
   const(fib0,FA),
   const(fib1,FB),
   up(1,N,FA,FB,F).
fib_bottomup_direct(0,F0) :-
   const(fib0,F0).

% Carve the constants fib(0) and fib(1) out of the code.

const(fib0,0).
const(fib1,1).

% Tail recursive call moving "bottom up" towards N.
%
% X:  the "current point of progress"
% N:  the N we want to reach
% FA: the value of fib(X-1)
% FB: the value of fib(X)
% F:  The variable that will receive the final result, fib(N)

up(X,N,FA,FB,F) :-
   X<N, % not there yet, compute fib(X+1)
   !,
   FC is FA + FB,
   Xn is X  + 1,
   up(Xn,N,FB,FC,F).
up(N,N,_,F,F).

Then:

?- fib_bottomup_direct(11,X).
X = 89.

Several more algorithms here; a README here.

Answer 3

This solution uses a tick less baggage, that is carried around.
The formulas are found at the end of the wiki fibmat section:

<pre id="in">
fib(N, X) :-
   powvec(N, (1,0), (0,1), (X,_)).
luc(N, Z) :-
   powvec(N, (1,0), (0,1), (X,Y)), Z is X+2*Y.

powvec(0, _, R, R) :- !.
powvec(N, A, R, S) :- N rem 2 =\= 0, !,
   mulvec(A, R, H), M is N//2, mulvec(A, A, B), powvec(M, B, H, S).
powvec(N, A, R, S) :- 
   M is N//2, mulvec(A, A, B), powvec(M, B, R, S).

mulvec((A1,A2), (B1,B2), (C1,C2)) :-
   C1 is A1*(B1+B2)+A2*B1,
   C2 is A1*B1+A2*B2.

?- fib(100,X).
?- luc(100,X).
</pre>
<script src="http://www.dogelog.ch/lib/exchange.js"></script>

Answer 4

fib2(120,X), X=[H|_], !. answers your question, binding H to the head of that reversed list, so, the 120th Fibonacci number.

Just insert the head-taking goal X=[H|_] into the query. Of course if you're really not interested in the list, you can fuse the two goals into one

fib2(120,[H|_]), !.

Your code does ~ 2N steps, which is still O(N) like an accumulator version would, so, not a big deal, it's fine as it is. The real difference is the O(N) space your version takes, v. the O(1) of the accumulator's.

But if you look closely at your code,

fib2(0, [0]).
fib2(1, [1,0]).
fib2(N, [R,X,Y|Zs]) :-
    N > 1,
    N1 is N - 1,
    fib2(N1, [X,Y|Zs]),
    R is X + Y.

you realize that it creates the N-long list of uninstantiated variables on the way down to the deepest level of recursion, then calculates them while populating the list with the calculated values on the way back up -- but only ever referring to the last two Fibonacci numbers, i.e. the first two values in that list. So you might as well make it explicit, and end up with .... an accumulator-based version, yourself!

fib3(0, 0, 0).
fib3(1, 1, 0).
fib3(N, R, X) :-
    N > 1,
    N1 is N - 1,
    fib3(N1, X, Y),
    R is X + Y.

except that it's still not tail-recursive. The way to achieve that is usually with additional argument(s) and you can see such a code in another answer here, by David Tonhofer. But hopefully you now see the clear path between it and this last one right here.

Answer 5

Just for fun, an even faster version of Fibonacci (even without using tail recursion) is presented below:

% -----------------------------------------------------------------------
% FAST FIBONACCI
% -----------------------------------------------------------------------

ffib(N, F) :-
    ff(N, [_, F]).

ff(1, [0, 1]) :- !.
ff(N, R) :-
    M is N // 2,
    ff(M, [A, B]),
    F1 is A^2   + B^2,
    F2 is 2*A*B + B^2,
    (   N mod 2 =:= 0
    ->  R = [F1, F2]
    ;   F3 is F1 + F2,
        R = [F2, F3]   ).

% -----------------------------------------------------------------------
% MOSTOWSKI COLLAPSE VERSION
% -----------------------------------------------------------------------

fib(N, X) :-
   powvec(N, (1,0), (0,1), (X,_)).

powvec(0, _, R, R) :- !.

powvec(N, A, R, S) :-
   N rem 2 =\= 0, !,
   mulvec(A, R, H),
   M is N // 2,
   mulvec(A, A, B),
   powvec(M, B, H, S).

powvec(N, A, R, S) :-
   M is N // 2,
   mulvec(A, A, B),
   powvec(M, B, R, S).

mulvec((A1,A2), (B1,B2), (C1,C2)) :-
   C1 is A1*(B1 + B2) + A2*B1,
   C2 is A1*B1 + A2*B2.

% -----------------------------------------------------------------------
% COMPARISON
% -----------------------------------------------------------------------

comparison :-
   format('n     fib   ffib  speed~n'),
   forall( between(21, 29, E),
      (  N is 2^E,
         cputime(fib( N, F1), T1),
         cputime(ffib(N, F2), T2),
         F1 = F2,        % confirm that both versions compute same answer!
         catch(R is T1/T2, _, R = 1),
         format('2^~w~|~t~2f~6+~|~t~2f~6+~|~t~2f~6+~n', [E, T1, T2, R]))).

cputime(Goal, Time) :-
   T0 is cputime,
   call(Goal),
   Time is cputime - T0.

The time complexity of both versions (mine and @MostowskiCollapse's) is O(lg n), ignoring multiplication cost.

Some simple empirical results (time in seconds) obtained with SWI-Prolog, version 8.2.4:

?- comparison.
n     fib   ffib  speed
2^21  0.05  0.02  3.00
2^22  0.09  0.05  2.00
2^23  0.22  0.09  2.33
2^24  0.47  0.20  2.31
2^25  1.14  0.45  2.52
2^26  2.63  1.02  2.58
2^27  5.89  2.34  2.51
2^28 12.78  5.28  2.42
2^29 28.97 12.25  2.36
true.

Answer 6

This one uses the Golden Ratio formula:

<pre id="in">
fib(N, S) :-
   powrad(N,(1,1),(1,0),(_,X)),
   powrad(N,(1,-1),(1,0),(_,Y)),
   S is (X-Y)//2^N.
luc(N, S) :-
   powrad(N,(1,1),(1,0),(X,_)),
   powrad(N,(1,-1),(1,0),(Y,_)),
   S is (X+Y)//2^N.

powrad(0, _, R, R) :- !.
powrad(N, A, R, S) :- N rem 2 =\= 0, !,
   mulrad(A, R, H), M is N//2, mulrad(A, A, B), powrad(M, B, H, S).
powrad(N, A, R, S) :-
   M is N//2, mulrad(A, A, B), powrad(M, B, R, S).

mulrad((A,B),(C,D),(E,F)) :-
   E is A*C+B*D*5,
   F is A*D+B*C.

?- fib(100,X).
?- luc(100,X).
</pre>
<script src="http://www.dogelog.ch/lib/exchange.js"></script>

Answer 7

The Donald Knuth based Fibonacci-by-matrix multiplication approach, as provided by Mostowski Collapse, but more explicit.

Algorithms can be found in the a module file plus a unit tests file on github:

The principle is based on a matrix identity provided by Donald Knuth (in Donald E. Knuth. The Art of Computer Programming. Volume 1. Fundamental Algorithms, p.80 of the second edition)

For n >= 1 we have (for n=0, the identity matrix appears on the right-hand side, but it is unclear what fib(-1) is):

                                  n
  [ fib(n+1) fib(n)   ]   [ 1  1 ]
  [                   ] = [      ]
  [ fib(n)   fib(n-1) ]   [ 1  0 ]

But if we work with constants fib(0) and fib(1) without assuming their value to be 0 and 1 respectively (we might be working with a special Fibonacci sequence), then we must stipulate that for n >= 1:

                                                       n-1
  [ fib(n+1) fib(n)   ]   [ fib(2) fib(1) ]   [ 1  1 ]
  [                   ] = [               ] * [      ]
  [ fib(n)   fib(n-1) ]   [ fib(1) fib(0) ]   [ 1  0 ]

We will separately compute the the "power matrix" on the right and explicitly multiply with the "fibonacci starter matrix", thus:

const(fib0,0).
const(fib1,1).

fib_matrixmult(N,F) :-
   N>=1,
   !,
   Pow is N-1,
   const(fib0,Fib0),
   const(fib1,Fib1),
   Fib2 is Fib0+Fib1,
   matrixpow(
      Pow,
      [[1,1],[1,0]],
      PowMx),
   matrixmult(
      [[Fib2,Fib1],[Fib1,Fib0]],
      PowMx,
      [[_,F],[F,_]]).
fib_matrixmult(0,Fib0) :-
   const(fib0,Fib0).

matrixpow(Pow, Mx, Result) :-
   matrixpow_2(Pow, Mx, [[1,0],[0,1]], Result).

matrixpow_2(Pow, Mx, Accum, Result) :-
   Pow > 0,
   Pow mod 2 =:= 1,
   !,
   matrixmult(Mx, Accum, NewAccum),
   Powm is Pow-1,
   matrixpow_2(Powm, Mx, NewAccum, Result).
matrixpow_2(Pow, Mx, Accum, Result) :-
   Pow > 0,
   Pow mod 2 =:= 0,
   !,
   HalfPow is Pow div 2,
   matrixmult(Mx, Mx, MxSq),
   matrixpow_2(HalfPow, MxSq, Accum, Result).
matrixpow_2(0, _, Accum, Accum).

matrixmult([[A11,A12],[A21,A22]],
           [[B11,B12],[B21,B22]],
           [[C11,C12],[C21,C22]]) :-
   C11 is A11*B11+A12*B21,
   C12 is A11*B12+A12*B22,
   C21 is A21*B11+A22*B21,
   C22 is A21*B12+A22*B22.

If your starter matrix is sure to be [[1,1],[1,0]] you can collapse the two operations matrixpow/3 followed by matrixmult/3 in the main predicate into a single call to matrixpow/3.

The above algorithm computes "too much" because two of the values in the matrix of Fibonacci numbers can be deduced from the other two. We can get rid of that redundancy. Mostowski Collapse presented a compact algorithm to do just that. Hereunder expanded for comprehensibility:

The idea is to get rid of redundant operations in matrixmult/3, by using the fact that all our matrices are symmetric and actually hold Fibonacci numbers

[ fib(n+1)  fib(n)   ]
[                    ]
[ fib(n)    fib(n-1) ]

So, if we multiply matrices A and B to yield C, we always have something of this form (even in the starter case where B is the identity matrix):

[ A1+A2  A1 ]   [ B1+B2  B1 ]   [ C1+C2  C1 ]
[           ] * [           ] = [           ]
[ A1     A2 ]   [ B1     B2 ]   [ C1     C2 ]

We can just retain the second columns of each matrix w/o loss of information. The operation between these vectors is not some standard operation like multiplication, let's mark it with ⨝:

[ A1 ]    [ B1 ]   [ C1 ]
[    ] ⨝ [    ] = [    ]
[ A2 ]    [ B2 ]   [ C2 ]

where:

C1 = B1*(A1+A2) + B2*A1 or A1*(B1+B2) + A2*B1
C2 = A1*B1 + A2*B2

fib_matrixmult_streamlined(N,F) :-
   N>=1,
   !,
   Pow is N-1,
   const(fib0,Fib0),
   const(fib1,Fib1),
   matrixpow_streamlined(
      Pow,
      v(1,0),
      PowVec),
   matrixmult_streamlined(
      v(Fib1,Fib0),
      PowVec,
      v(F,_)).
fib_matrixmult_streamlined(0,Fib0) :-
   const(fib0,Fib0).

matrixpow_streamlined(Pow, Vec, Result) :-
   matrixpow_streamlined_2(Pow, Vec, v(0,1), Result).

matrixpow_streamlined_2(Pow, Vec, Accum, Result) :-
   Pow > 0,
   Pow mod 2 =:= 1,
   !,
   matrixmult_streamlined(Vec, Accum, NewAccum),
   Powm is Pow-1,
   matrixpow_streamlined_2(Powm, Vec, NewAccum, Result).
matrixpow_streamlined_2(Pow, Vec, Accum, Result) :-
   Pow > 0,
   Pow mod 2 =:= 0,
   !,
   HalfPow is Pow div 2,
   matrixmult_streamlined(Vec, Vec, VecVec),
   matrixpow_streamlined_2(HalfPow, VecVec, Accum, Result).
matrixpow_streamlined_2(0, _, Accum, Accum).

matrixmult_streamlined(v(A1,A2),v(B1,B2),v(C1,C2)) :-
   C1 is A1*(B1+B2) + A2*B1,
   C2 is A1*B1 + A2*B2.