Lucas Sequences

October 1, 2013

We begin with a function that computes the first n elements of a Lucas sequence, given the first two elements of the sequence:

(define (lucas l2 l1 n)
  (let loop ((n n) (l2 l2) (l1 l1) (ls (list l1 l2)))
    (if (zero? n) (reverse ls)
      (let ((l (+ l1 l2)))
        (loop (- n 1) l1 l (cons l ls))))))

This can be used to calculate the Fibonacci and Lucas numbers:

> (lucas 1 1 20) ; fibonacci numbers
(1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584
4181 6765 10946 17711)
> (lucas 1 3 20) ; lucas numbers
(1 3 4 7 11 18 29 47 76 123 199 322 521 843 1364 2207 3571
5778 9349 15127 24476 39603)

Another approach uses P, Q and the first two elements of the sequence:

(define (lucas p q l2 l1 n)
  (let loop ((n n) (l2 l2) (l1 l1) (ls (list l1 l2)))
    (if (zero? n) (reverse ls)
      (let ((l (- (* p l1) (* q l2))))
        (loop (- n 1) l1 l (cons l ls))))))

And here are the Fibonacci and Lucas numbers again:

> (lucas 1 -1 1 1 20)
(1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584
4181 6765 10946 17711)
> (lucas 1 -1 1 3 20)
(1 3 4 7 11 18 29 47 76 123 199 322 521 843 1364 2207 3571
5778 9349 15127 24476 39603)

Our last functions are mutually-recursive functions that compute the nth element of a U or V sequence in time O(log n):

(define (u n p q)
  (if (< n 1) 1 (if (= n 1) p
    (let ((k (quotient n 2)))
      (if (odd? n)
          (- (* (u (+ k 1) p q) (v k p q)) (expt q k))
          (* (u k p q) (v k p q)))))))

(define (v n p q)
  (if (< n 1) 2 (if (= n 1) p
    (let ((k (quotient n 2)))
      (if (odd? n)
          (- (* (v (+ k 1) p q) (v k p q)) (* p (expt q k)))
          (- (expt (v k p q) 2) (* 2 (expt q k))))))))

Then we can compute particular elements of the Fibonacci or Lucas sequences:

> (u 22 1 -1) ; fibonacci
17711
> (v 22 1 -1) ; lucas
39603

You can run the program at http://programmingpraxis.codepad.org/90DQYRBD.

Pages: 1 2

Posted by programmingpraxis
Filed in Exercises
6 Comments »

6 Responses to “Lucas Sequences”

  1. Graham said
    October 2, 2013 at 1:26 PM

    Here’s a Haskell version that gives two implementations: a linear update a la
    the usual iterative Fibonacci procedure, and a (mutually) recursive one like
    the last solution.

    import Control.Arrow ((&&&))

lucas :: Integer -> Integer -> Integer -> Integer -> [Integer]
lucas x0 x1 p q = x0 : x1 : zipWith (\x y-> p * x - q * y) as bs
    where
        as = lucas x0 x1 p q
        bs = tail as

u1 :: Integer -> Integer -> [Integer]
u1 = lucas 0 1

v1 :: Integer -> Integer -> [Integer]
v1 p = lucas 2 p p

u2 :: Integer -> Integer -> Integer -> Integer
u2 p q n | n <= 0    = 0
         | n == 1    = 1
         | odd n     = u2 p q (k + 1) * v2 p q k - q^k
         | otherwise = u2 p q k * v2 p q k
         where
             k = n `div` 2

v2 :: Integer -> Integer -> Integer -> Integer
v2 p q n | n <= 0    = 2
         | n == 1    = p
         | odd n     = v2 p q (k + 1) * v2 p q k - p * q^k
         | otherwise = v2 p q k ^2  - 2 * q^k
         where
             k = n `div` 2

main :: IO ()
main = do
    let (p,q) = (1,-1)
    mapM_ (print . take 17) [u1 p q, v1 p q]
    print $ map (u2 p q &&& v2 p q) [0..16]
  2. Graham said
    October 2, 2013 at 5:00 PM

    I’m away from my computer at the moment so I’m not sure if this works, but I like the following definition better (if it works):
    [sourecode lang=”css”]
    lucas x0 x1 p q = x0 : scanl (\x y -> p * x – q * y) x1 $ lucas x0 x1 p q
    [/sourcecode]

  3. Graham said
    October 2, 2013 at 5:00 PM

    Apologies for the phone keyboard typo.

  4. Paul said
    October 4, 2013 at 6:56 AM

    In Python.

    # see http://mathworld.wolfram.com/LucasSequence.html
def U(P, Q):
    """generates Lucas U series - the first item is U1"""
    if P ** 2 - 4 * Q <= 0:
        raise ValueError("D is not positive")
    a, b = 0, 1
    while 1:
        yield b
        a, b = b, P * b - Q * a
        
def V(P, Q):
    """generates Lucas V series - the first item is V1"""
    if P ** 2 - 4 * Q <= 0:
        raise ValueError("D is not positive")
    a, b = 2, P
    while 1:
        yield b
        a, b = b, P * b - Q * a
        
def Un(P, Q, n):
    """recursive formula for Lucas U"""
    if n in (0,1):
        return n
    m = n // 2
    if n & 1:
        return Un(P, Q, m+1) * Vn(P, Q, m) - Q ** m
    else:
        return Un(P, Q, m) * Vn(P, Q, m)
  5. Paul said
    October 4, 2013 at 6:58 AM

    Oops, forgot the recursive formula for Lucas V.

    # see mathworld.wolfram.com/LucasSequence.html
def U(P, Q):
    """generates Lucas U series - the first item is U1"""
    if P ** 2 - 4 * Q <= 0:
        raise ValueError("D is not positive")
    a, b = 0, 1
    while 1:
        yield b
        a, b = b, P * b - Q * a
        
def V(P, Q):
    """generates Lucas V series - the first item is V1"""
    if P ** 2 - 4 * Q <= 0:
        raise ValueError("D is not positive")
    a, b = 2, P
    while 1:
        yield b
        a, b = b, P * b - Q * a
        
def Un(P, Q, n):
    """recursive formula for Lucas U"""
    if n in (0,1):
        return n
    m = n // 2
    if n & 1:
        return Un(P, Q, m+1) * Vn(P, Q, m) - Q ** m
    else:
        return Un(P, Q, m) * Vn(P, Q, m)

def Vn(P, Q, n):
    """recursive formula for Lucas V"""
    if n in (0,1):
        return (2, P)[n]
    m = n // 2
    if n & 1:
        return Vn(P, Q, m+1) * Vn(P, Q, m) - P * Q ** m
    else:
        return Vn(P, Q, m) ** 2 - 2 * Q ** m

Leave a comment