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.

About these ads

Pages: 1 2

6 Responses to “Lucas Sequences”

  1. Graham said

    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

    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

    Apologies for the phone keyboard typo.

  4. Paul said

    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

    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 Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 574 other followers

%d bloggers like this: