Fibonacci number programs

Hi all,

I'm finally preparing some pages to put on the WWW. Better
late than never, eh? :-)  Anyway, it occurred to me that
this list might be interested in some Prolog routines I've
crafted together, which I'll probably be making publicly
available on or via my upcoming homepage.

Below are two programs, the last of which was inspired by
an unfinished Prolog version of a Haskell program posted
in March/April 2002 on comp.lang.prolog by Andrew Bromage.
You'll notice that I've used non-determinism and back-
tracking to minimize memory usage, but that may be easily
changed to a "normal" recursive program, ie. "the Prolog
way", should you so prefer. Feel free to use these programs
-- should you find them worthy of your high standards.

Feedback is welcome. (I expect to get a mouthful from the
LP'ers, but that's alright.)  ;-)

Unless I hear otherwise, I'll post more programs/routines
in the course of my getting organised.

Cheers,
Daniel

PS. I just tried the query fib_range(1, 9999, Fibs) and it
took 74.09s cpu -- but the high 9000 range produces some
increasingly frightening long numbers, besides the fact
that I have too many programs running right now. No savings
to be made, though, by doing it in stages of 1000: in all,
approx. 72.66 cpu (0.62 + 1.36 + 2.60 + 3.44 + 5.88 + 7.37
+ 8.56 + 9.67 + 15.04 + 18.34). (ECLiPSe [TkECLiPSe] 5.2
running on Win2K, 400 MHz PII 256 Mb RAM.) Is that good, or
what?

/********************************************************************

File:    FIB_RANGE.PL
Purpose: Compute a range of Fibonacci numbers using recursion.
Author:  Daniel Dudley <daniel.dudley@chello.no>
Created: 2003.02.02
Notes:   (1) fib_range/3 calls an efficient, tail-recursive
             predicate, fibs/3, which computes the 2nd..N Fibonacci
             numbers in a given range. The returned list is in
             ascending order.
         (2) Include FIB_QMAT_ND.PL
         (3) The program should run with any standard Prolog system.
History: 

********************************************************************/

% fib_range(+From,+To,Fibs) is true if Fibs is the list of the
% range of Fibonacci numbers determined by the indices From and To.
% Examples (ECLiPSe [TkECLiPSe] 5.2 running on Win2K, 400 MHz PII):
% 1) ?- fib_range(1, 72, Fibs).
%    Fibs = [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377,
%            610, 987, 1597, 2584, 4181, ...]
%    Yes (0.01s cpu)
% 2) ?- fib_range(100, 500, Fibs).
%    Fibs = [354224848179261915075, 573147844013817084101, ...]
%    Yes (0.20s cpu)

fib_range(From,To,[Fib|Fibs]):-
  % guards follow:
  integer(From), integer(To),
  From > 0,
  (From =< To
   -> Min is From, Max is To
   ;  Min is To, Max is From
  ),
  % begin computation:
  fib_qmat(Min,[_,Fib,_,_]),
  fibs(Min,Max,Fibs).


% fibs(+Nmin, +Nmax, -Fs) is true if Fs is the list of Nmin..Nmax
% Fibonacci numbers

fibs(N,Max,F1s):-
  (N = Max
   -> F1s = []
   ;  N =< Max,
      N1 is N + 1,
      fib_qmat(N1,[_,F,_,_]),
      F1s = [F|Fs],
      fibs(N1,Max,Fs)
  ).

% end of program


/********************************************************************

File:    FIB_QMAT_ND.PL
Purpose: Compute the nth Fibonacci Q-matrix.
Author:  Daniel Dudley <dudley@online.no>
Created: 2002.04.11
Notes:   1) fib_qmat/2 calls an efficient, non-deterministic
            (iterative) predicate, m4pow_nd/5, which implements
            fast exponentiation of the Q-matrix (via backtracking).
         2) The output variable Qmatrix is a static list of
            the form: [F1,F,F,F0], where
                F1 is the nth+1 Fibonacci number,
                F is the nth Fibonacci number, and
                F0 is the nth-1 Fibonacci number.
            The first (base) Q-matrix is [1,1,1,0].
History: 

********************************************************************/

% fib_qmat(+N,-Q) is true if Q is the Nth Fibonacci Q-matrix
% Eg. ?- fib_qmat(1,Q). => Q = [1,1,1,0]  Yes (0.00s cpu)
%     ?- fib_qmat(2,Q). => Q = [2,1,1,1]  Yes (0.00s cpu)
%     ?- fib_qmat(3,Q). => Q = [3,2,2,1]  Yes (0.00s cpu)
%     ?- fib_qmat(6,Q). => Q = [13,8,8,5]  Yes (0.00s cpu)
%     ?- fib_qmat(49,Q). =>
%     Q = [12586269025, 7778742049, 7778742049, 4807526976]
%     Yes (0.00s cpu)

fib_qmat(N,Qmatrix):-
  integer(N),
  Base = [1,1,1,0],
  (N < 1
   -> Qmatrix = Base
   ;  Exp is N - 1, 
      m4pow_nd(Exp,Base,Base,ExpOut,Qmatrix),
      ExpOut = 0,
      !  % remove choice point on m4pow_nd/5
  ).


% m4pow_nd(+Exp,+Base,+Build,-ExpOut,-Qmat):
% non-deterministic fast exponentiation of a Fibonacci Q-matrix

m4pow_nd(Exp,_,Qmat,Exp,Qmat).
m4pow_nd(Exp,Base,Build,ExpOut,Qmat):-
  Exp > 0,
  (1 is Exp mod 2
   -> m4mult(Build,Base,Build1)
   ;  Build1 = Build
  ),
  m4mult(Base,Base,Base1),
  Exp1 is Exp // 2,
  m4pow_nd(Exp1,Base1,Build1,ExpOut,Qmat).


% m4mult(+M1,+M2,-M3): 2x2 matrix multiply

m4mult([A1,B1,C1,D1],[A2,B2,C2,D2],[A,B,C,D]) :-
  A is A1*A2 + B1*C2,
  B is A1*B2 + B1*D2,
  C is C1*A2 + D1*C2,
  D is C1*B2 + D1*D2.

% end of program