Prolog - 使用累加器查找第 N 个斐波那契数
Prolog - Finding the Nth Fibonacci number using accumulators
我有这段代码可以以相反的顺序生成斐波那契数列列表。
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.
虽然我只需要第一个元素。问题是这段代码还在列表后给出了一个 false. ，所以我获取第一个元素的所有尝试都失败了。有没有什么办法可以得到列表中的第一个元素，或者有什么其他方法可以用累加器计算第 N 个斐波那契数。
提前致谢。
它 给出一个“假” 因为 Prolog 不确定在它提供的第一个解决方案之后是否还有更多解决方案：
?- fib2(4,L).
L = [3,2,1,1,0] ; % maybe more solutions?
false. % no
这不是问题：你可以告诉Prolog在第一个之后确实没有更多的解决方案（或者你不想看到它们）：
?- once(fib2(4,L)).
或
?- fib2(4,L),!.
或者你可以切入每个第一个子句，告诉 Prolog 如果头部匹配，就没有必要尝试另一个子句。这摆脱了流浪的“可能的解决方案”：
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.
可能有问题的是给定的算法存储了所有的fib(i)并在递归调用后进行加法，这意味着Prolog无法优化递归调用进入循环。
对于“基于累加器”（自下而上）的计算方式 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).
然后：
?- fib_bottomup_direct(11,X).
X = 89.
我得到了这个 logarithmic steps O(log n) 解决方案，甚至尾递归。
只是为了好玩，它还可以计算第 n 个卢卡斯数：
<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>
您可以比较：
https://www.wolframalpha.com/input/?i=Fibonacci[100]
https://www.wolframalpha.com/input/?i=LucasN[100]
编辑 28.06.2021：
这是一个非常快速的解释为什么矩阵算法有效。
我们只需要证明斐波那契的一步是线性的。即
这种递归关系导致线性矩阵：
F_{n+2} = F_{n}+F_{n+1}
要查看矩阵，我们必须假设矩阵 M 将向量 b=[Fn,Fn+1] 转换为向量 b'=[F_{n+1}, F_{n+2}]:
b' = M*b
这个矩阵可能是什么？解决一下：
|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}|
此解决方案使用了一种无行李箱，可以随身携带。
公式位于 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>
fib2(120,X), X=[H|_], !. 回答了您的问题，将 H 绑定到该反向列表的头部，因此，第 120 个斐波那契数。
只需将抢眼目标 X=[H|_] 插入 查询。当然如果你真的对榜单不感兴趣，也可以把两个目标合二为一
fib2(120,[H|_]), !.
你的代码执行了 ~ 2N 步，这仍然是 O(N) 就像累加器版本一样，所以，没什么大不了的，它很好。真正的区别是 O(N) space 你的版本，v. 累加器的 O(1)。
但是如果你仔细查看你的代码，
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.
你意识到它在向下递归的最深层次的过程中创建了 N 长的未实例化变量列表，然后计算它们，同时在返回的过程中用计算值填充列表 - - 但只引用最后的 两个 斐波那契数，即该列表中的 前两个 值。所以你最好把它说清楚，最后得到......一个基于累加器的版本，你自己！
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.
除了它仍然不是尾递归的。实现这一点的方法通常是使用附加参数，您可以在另一个 here, by David Tonhofer 中看到这样的代码。但希望你现在能看到它和最后一个之间的清晰路径。
只是为了好玩，下面展示了一个甚至更快的斐波那契版本（即使不使用尾递归）：
% -----------------------------------------------------------------------
% 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.
两个版本（我的和@MostowskiCollapse 的）的时间复杂度是O(lg n), 忽略乘法成本。
使用 SWI-Prolog 版本 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.
基于 Donald Knuth 的 Fibonacci-by-matrix 乘法方法，由
Mostowski Collapse，但更明确。
算法可以在模块文件加上单元测试文件on github:
中找到
该原理基于 Donald Knuth 提供的矩阵恒等式（在 Donald E. Knuth. The Art of Computer Programming. Volume 1. Fundamental
算法，第二版第 80 页）
对于 n >= 1 我们有（对于 n=0，单位矩阵出现在右侧，但不清楚 fib(-1) 是什么):
n
[ fib(n+1) fib(n) ] [ 1 1 ]
[ ] = [ ]
[ fib(n) fib(n-1) ] [ 1 0 ]
但是如果我们使用常量 fib(0) 和 fib(1) 而不假设它们的值分别为 0 和 1 （我们可能正在处理一个特殊的斐波那契数列），那么我们必须规定对于 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 ]
我们将单独计算右侧的“幂矩阵”并显式地与“斐波那契起始矩阵”相乘，因此：
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.
如果您的起始矩阵肯定是 [[1,1],[1,0]]，您可以将主谓词中的两个操作 matrixpow/3 和后面的 matrixmult/3 合并为对 matrixpow/3 的单个调用.
上述算法计算“太多”，因为斐波那契数列矩阵中的两个值可以从另外两个推导出来。我们可以摆脱这种冗余。 Mostowski Collapse 提出了一个紧凑的算法来做到这一点。为了便于理解，下文进行了扩展：
我们的想法是摆脱 matrixmult/3 中的冗余操作，方法是利用我们所有的矩阵都是对称的并且实际上持有
斐波那契数
[ fib(n+1) fib(n) ]
[ ]
[ fib(n) fib(n-1) ]
所以，如果我们将矩阵 A 和 B 相乘得到 C，我们总是有这种形式的东西（即使在初始情况下 B 是
单位矩阵）：
[ A1+A2 A1 ] [ B1+B2 B1 ] [ C1+C2 C1 ]
[ ] * [ ] = [ ]
[ A1 A2 ] [ B1 B2 ] [ C1 C2 ]
我们可以只保留每个矩阵的第二列 w/o 损失
信息。这些向量之间的操作不是一些
像乘法这样的标准运算，我们用⨝来标记：
[ A1 ] [ B1 ] [ C1 ]
[ ] ⨝ [ ] = [ ]
[ A2 ] [ B2 ] [ C2 ]
其中：
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.
这个使用Golden Ratio公式：
<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>
我有这段代码可以以相反的顺序生成斐波那契数列列表。
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.
虽然我只需要第一个元素。问题是这段代码还在列表后给出了一个 false. ，所以我获取第一个元素的所有尝试都失败了。有没有什么办法可以得到列表中的第一个元素，或者有什么其他方法可以用累加器计算第 N 个斐波那契数。
提前致谢。
它 给出一个“假” 因为 Prolog 不确定在它提供的第一个解决方案之后是否还有更多解决方案：
?- fib2(4,L).
L = [3,2,1,1,0] ; % maybe more solutions?
false. % no
这不是问题：你可以告诉Prolog在第一个之后确实没有更多的解决方案（或者你不想看到它们）：
?- once(fib2(4,L)).
或
?- fib2(4,L),!.
或者你可以切入每个第一个子句，告诉 Prolog 如果头部匹配，就没有必要尝试另一个子句。这摆脱了流浪的“可能的解决方案”：
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.
可能有问题的是给定的算法存储了所有的fib(i)并在递归调用后进行加法，这意味着Prolog无法优化递归调用进入循环。
对于“基于累加器”（自下而上）的计算方式 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).
然后：
?- fib_bottomup_direct(11,X).
X = 89.
我得到了这个 logarithmic steps O(log n) 解决方案，甚至尾递归。
只是为了好玩，它还可以计算第 n 个卢卡斯数：
<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>
您可以比较：
https://www.wolframalpha.com/input/?i=Fibonacci[100]
https://www.wolframalpha.com/input/?i=LucasN[100]
编辑 28.06.2021：
这是一个非常快速的解释为什么矩阵算法有效。
我们只需要证明斐波那契的一步是线性的。即
这种递归关系导致线性矩阵：
F_{n+2} = F_{n}+F_{n+1}
要查看矩阵，我们必须假设矩阵 M 将向量 b=[Fn,Fn+1] 转换为向量 b'=[F_{n+1}, F_{n+2}]:
b' = M*b
这个矩阵可能是什么？解决一下：
|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}|
此解决方案使用了一种无行李箱，可以随身携带。
公式位于 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>
fib2(120,X), X=[H|_], !. 回答了您的问题，将 H 绑定到该反向列表的头部，因此，第 120 个斐波那契数。
只需将抢眼目标 X=[H|_] 插入 查询。当然如果你真的对榜单不感兴趣，也可以把两个目标合二为一
fib2(120,[H|_]), !.
你的代码执行了 ~ 2N 步，这仍然是 O(N) 就像累加器版本一样，所以，没什么大不了的，它很好。真正的区别是 O(N) space 你的版本，v. 累加器的 O(1)。
但是如果你仔细查看你的代码，
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.
你意识到它在向下递归的最深层次的过程中创建了 N 长的未实例化变量列表，然后计算它们，同时在返回的过程中用计算值填充列表 - - 但只引用最后的 两个 斐波那契数，即该列表中的 前两个 值。所以你最好把它说清楚，最后得到......一个基于累加器的版本，你自己！
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.
除了它仍然不是尾递归的。实现这一点的方法通常是使用附加参数，您可以在另一个
只是为了好玩，下面展示了一个甚至更快的斐波那契版本（即使不使用尾递归）：
% -----------------------------------------------------------------------
% 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.
两个版本（我的和@MostowskiCollapse 的）的时间复杂度是O(lg n), 忽略乘法成本。
使用 SWI-Prolog 版本 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.
基于 Donald Knuth 的 Fibonacci-by-matrix 乘法方法，由 Mostowski Collapse，但更明确。
算法可以在模块文件加上单元测试文件on github:
中找到该原理基于 Donald Knuth 提供的矩阵恒等式（在 Donald E. Knuth. The Art of Computer Programming. Volume 1. Fundamental 算法，第二版第 80 页）
对于 n >= 1 我们有（对于 n=0，单位矩阵出现在右侧，但不清楚 fib(-1) 是什么):
n
[ fib(n+1) fib(n) ] [ 1 1 ]
[ ] = [ ]
[ fib(n) fib(n-1) ] [ 1 0 ]
但是如果我们使用常量 fib(0) 和 fib(1) 而不假设它们的值分别为 0 和 1 （我们可能正在处理一个特殊的斐波那契数列），那么我们必须规定对于 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 ]
我们将单独计算右侧的“幂矩阵”并显式地与“斐波那契起始矩阵”相乘，因此：
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.
如果您的起始矩阵肯定是 [[1,1],[1,0]]，您可以将主谓词中的两个操作 matrixpow/3 和后面的 matrixmult/3 合并为对 matrixpow/3 的单个调用.
上述算法计算“太多”，因为斐波那契数列矩阵中的两个值可以从另外两个推导出来。我们可以摆脱这种冗余。 Mostowski Collapse 提出了一个紧凑的算法来做到这一点。为了便于理解，下文进行了扩展：
我们的想法是摆脱 matrixmult/3 中的冗余操作，方法是利用我们所有的矩阵都是对称的并且实际上持有
斐波那契数
[ fib(n+1) fib(n) ]
[ ]
[ fib(n) fib(n-1) ]
所以，如果我们将矩阵 A 和 B 相乘得到 C，我们总是有这种形式的东西（即使在初始情况下 B 是 单位矩阵）：
[ A1+A2 A1 ] [ B1+B2 B1 ] [ C1+C2 C1 ]
[ ] * [ ] = [ ]
[ A1 A2 ] [ B1 B2 ] [ C1 C2 ]
我们可以只保留每个矩阵的第二列 w/o 损失 信息。这些向量之间的操作不是一些 像乘法这样的标准运算，我们用⨝来标记：
[ A1 ] [ B1 ] [ C1 ]
[ ] ⨝ [ ] = [ ]
[ A2 ] [ B2 ] [ C2 ]
其中：
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.
这个使用Golden Ratio公式：
<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>