Fibonacci Numbers

This code requires the Multi-Precision Integer Arithmetic package.

/* fibonacci(N, F) is true if N is an integer, and F is a string containing  */
/*   the N-th Fibonacci number.  (The first few Fibonacci numbers are        */
/*   0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181,...)    */
/*   This is based on the identities F[2j+1] = F[j]^2+F[j+1]^2 and           */
/*   F[2j+2] = 2F[j]F[j+1]+F[j+1]^2.                                         */
/* e.g. fibonacci(236, `9366947731425726508977331996039353971111632790877`). */
fibonacci(0, `0`):-!.
fibonacci(N, F):-
  bitcount(N, K),
  K1 is K - 1,
  fibonacci_1(K1, N, [1], [1], F0), % F[1]=1, F[2]=1
  bigint_string(F0, F).

fibonacci_1(0, _, F, _, F):-!. % F[n]=X
fibonacci_1(K, N, X, Y, F):-   % X is F[j], Y is F[j+1]
  is_one(N, K), !,
  multiply(X, X, XX),
  multiply(Y, Y, YY),
  multiply(X, Y, XY),
  add(XX, YY, X1),             % F[2j+1]=F[j]^2+F[j+1]^2
  add(XY, XY, XY2),
  add(XY2, YY, Y1),            % F[2j+2]=2F[j]F[j+1]+F[j+1]^2
  K1 is K - 1,
  fibonacci_1(K1, N, X1, Y1, F).
fibonacci_1(K, N, X, Y, F):-   % X is F[j], Y is F[j+1]
  multiply(X, X, XX),
  multiply(Y, Y, YY),
  multiply(X, Y, XY),
  add(XY, XY, XY2),
  subtract(XY2, XX, X1),       % F[2j]=2F[j-1]F[j]+F[j]^2=2F[j]F[j+1]-F[j]^2
  add(XX, YY, Y1),             % F[2j+1]=F[j]^2+F[j+1]^2
  K1 is K - 1,
  fibonacci_1(K1, N, X1, Y1, F).

/* bitcount(N, Bits) is true if Bits is the number of significant binary     */
/*   digits in the positive integer N.                                       */
/* e.g. bitcount(19, 5).                                                     */
bitcount(N, Bits):-
  N >= 0,
  bitcount_1(0, N, Bits).

bitcount_1(I, 0, I):-!.
bitcount_1(I, N, Bits):-
  I1 is I + 1,
  N1 is N // 2,
  bitcount_1(I1, N1, Bits).
  
/* is_one(N, K) is true if the K-th bit of the positive integer N is 1,      */
/*   where the least significant bit corresponds to K=1.                     */
is_one(N, 1):-!,
  N mod 2 =:= 1.
is_one(N, K):-
  N1 is N // 2,
  K1 is K - 1,
  is_one(N1, K1).

LPA Index     Home Page

Valid HTML 4.0 Strict