CODEWITHSUNDEEP
c#performancealgorithm
ArtVandelay
ArtVandelay

Reputation: 97

Calculating Lucas Sequences efficiently

I'm implementing the p+1 factorization algorithm. For that I need to calculate elements of the lucas sequence which is defined by:

(1) x_0 = 1, x_1 = a
(2) x_n+l  =  2 * a * x_n - x_n-l

I implemented it (C#) recursively but it is inefficient for bigger indexes.

static BigInteger Lucas(BigInteger a, BigInteger Q, BigInteger N)
    {
        if (Q == 0)
            return 1;
        if (Q == 1)
            return a;
        else
            return (2 * a * Lucas(a, Q - 1, N) - Lucas(a, Q - 2, N)) % N;
    }

I also know 

(3) x_2n = 2 * (x_n)^2 - 1
(4) x_2n+1 = 2 * x_n+1 * x_n - a
(5) x_k(n+1) = 2 * x_k * x_kn - x_k(n-1)

(3) and (4) should help to calculate bigger Qs. But I'm unsure how. Somehow with the binary form of Q I think.

Any help is appreciated.

Upvotes: 3

Views: 1316

Answers (4)

Samuel Edwin Ward
Samuel Edwin Ward

Reputation: 6675

You can get some improvements (just a factor of a million...) without resorting to really fancy math.

First let's make the data flow a little more explicit:

    static BigInteger Lucas(BigInteger a, BigInteger Q, BigInteger N)
    {
        if (Q == 0)
        {
            return 1;
        }
        else if (Q == 1)
        {
            return a;
        }
        else
        {
            BigInteger q_1 = Lucas(a, Q - 1, N);
            BigInteger q_2 = Lucas(a, Q - 2, N);
            return (2 * a * q_1 - q_2) % N;
        }
    }

Unsurprisingly, this doesn't really change the performance.

However, it does make it clear that we only need two previous values to compute the next value. This lets us turn the function upside down into an iterative version:

    static BigInteger IterativeLucas(BigInteger a, BigInteger Q, BigInteger N)
    {
        BigInteger[] acc = new BigInteger[2];
        Action<BigInteger> push = (el) => {
                acc[1] = acc[0];
                acc[0] = el;
        };
        for (BigInteger i = 0; i <= Q; i++)
        {
            if (i == 0)
            {
                push(1);
            }
            else if (i == 1)
            {
                push(a);
            }
            else
            {
                BigInteger q_1 = acc[0];
                BigInteger q_2 = acc[1];
                push((2 * a * q_1 - q_2) % N);
            }
        }
        return acc[0];
    }

There might be a clearer way to write this, but it works. It's also much faster. It's so much faster it's kind of impractical to measure. On my system, Lucas(4000000, 47, 4000000) took about 30 minutes, and IterativeLucas(4000000, 47, 4000000) took about 2 milliseconds. I wanted to compare 48, but I didn't have the patience.

You can squeeze a little more out (maybe a factor of two?) using these properties of modular arithmetic:

(a + b) % n = (a%n + b%n) % n
(a * b) % n = ((a%n) * (b%n)) % n

If you apply these, you'll find that a%N occurs a few times so you can win by precomputing it once before the loop. This is particularly helpful when a is a lot bigger than N; I'm not sure if that happens in your application.

There are probably some clever mathematical techniques that would blow this solution out of the water, but I think it's interesting that such an improvement can be achieved just by shuffling a little code around.

Upvotes: 0

Bernhard Barker
Bernhard Barker

Reputation: 55609

(3) x_2n = 2 * (x_n)^2 - 1
(4) x_2n+1 = 2 * x_n+1 * x_n - a

Whenever you see 2n, you should think "that probably indicates an even number", and similarly 2n+1 likely means "that's an odd number".

You can modify the x indices so you have n on the left (as to make it easier to understand how this corresponds to recursive function calls), just be careful regarding rounding.

3) 2n     n
=> n      n/2

4) it is easy to see that if x = 2n+1, then n = floor(x/2)
     and similarly n+1 = ceil(x/2)

So, for #3, we have: (in pseudo-code)

if Q is even
   return 2 * (the function call with Q/2) - 1

And for #4:

else // following from above if
   return 2 * (the function call with floor(Q/2))
            * (the function call with ceil(Q/2)) - a

And then we can also incorporate a bit of memoization to prevent calculating the return value for the same parameters multiple times:

  • Keep a map of Q value to return value.
  • At the beginning of the function, check if Q's value exists in the map. If so, return the corresponding return value.
  • When returning, add Q's value and the return value to the map.

Upvotes: 1

Tyler Durden
Tyler Durden

Reputation: 11532

The n-th Lucas number has the value:

enter image description here

Exponentiation by squaring can be used to evaluate the function. For example, if n=1000000000, then n = 1000 * 1000^2 = 10 * 10^2 * 1000^2 = 10 * 10^2 * (10 * 10^2 )^2. By simplifying in this way you can greatly reduce the number of calculations.

Upvotes: 0

MBo
MBo

Reputation: 80197

Here one can see how to find Nth Fibbonaci number using matrix powering with matrix

      n
(1 1)
(1 0)

You may exploit this approach to calculate Lucas numbers, using matrix (for your case x_n+l = 2 * a * x_n - x_n-l)

        n
(2a -1)
(1   0)

Note that Nth power of matrix could be found with log(N) matrix multiplications by means of exponentiation by squaring

Upvotes: 3

Related Questions