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).
LPA Index
Home Page