#### Lucas Sequences

/* lucas(I, P, Q, U, V) is true if U is the (I+1)-th element and V is the    */
/*   I-th element respectively of the two Lucas sequences characterized by   */
/*   P and Q.                                                                */
/* e.g. lucas(71,1,-1,F,L) gives F=498454011879264 (the 72-nd Fibonacci      */
/*      number) and L=688846502588399 (the 71-st Lucas number)               */
/*                                                                           */
/* The elements of the sequences are defined by:                             */
/*   U[0]=0, U[1]=1, U[j]=PU[j-1]-QU[j-2]                                    */
/*   V[0]=2, V[1]=P, V[j]=PV[j-1]-QV[j-2]                                    */
/* It can be shown that:                                                     */
/*   U[2j]=U[j]V[j]                                                          */
/*   U[2j+1]=U[j+1]V[j]-Q^j                                                  */
/*   V[2j]=V[j]^2-2Q^j                                                       */
/*   V[2j+1]=V[j+1]V[j]-PQ^j                                                 */
/* so that U[i] and V[i] can be calculated in O(log i) steps.                */
lucas(0, _, _, 1, 2):-!.
lucas(I, P, Q, U, V):-
Qj is Q,       % Q^1
Ub is P,       % U[2]
Va is P,       % V[1]
Vb is P*P-2*Q, % V[2]
bitcount(I, K),
K1 is K - 1,
lucas_1(K1, I, P, Q, Qj, Ub, Va, Vb, U, V).

lucas_1(0, _, _, _, _, U, V, _, U, V):-!.
lucas_1(K, I, P, Q, Qj, Ub, Va, Vb, U, V):- % Q^j, U[j+1], V[j], V[j+1]
is_one(I, K), !,
QjQ is Qj*Q,
Qj1 is QjQ*Qj,                            % Q^(2j+1)=Q^j*Q^j*Q
Ub1 is Ub*Vb,                             % U[2j+2]=U[j+1]V[j+1]
Va1 is Va*Vb-P*Qj,                        % V[2j+1]=V[j]V[j+1]-PQ^j
Vb1 is Vb*Vb-2*QjQ,                       % V[2j+2]=V[j+1]^2-2Q^(j+1)
K1 is K - 1,
lucas_1(K1, I, P, Q, Qj1, Ub1, Va1, Vb1, U, V).
lucas_1(K, I, P, Q, Qj, Ub, Va, Vb, U, V):- % Q^j, U[j+1], V[j], V[j+1]
Qj1 is Qj*Qj,                             % Q^(2j)=Q^j*Q^j
Ub1 is Ub*Va-Qj,                          % U[2j+1]=U[j+1]V[j]-Q^j
Vb1 is Va*Vb-P*Qj,                        % V[2j+1]=V[j]V[j+1]-PQ^j
Va1 is Va*Va-2*Qj,                        % V[2j]=V[j]^2-2Q^j
K1 is K - 1,
lucas_1(K1, I, P, Q, Qj1, Ub1, Va1, Vb1, U, V).

/* 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(N, 0, Bits).

bitcount_1(0, Bits, Bits):-!.
bitcount_1(N, Bits0, Bits):-
Bits1 is Bits0 + 1,
N1 is N // 2,
bitcount_1(N1, Bits1, 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):-!,
1 is N mod 2.
is_one(N, K):-
N1 is N // 2,
K1 is K - 1,
is_one(N1, K1).