Lucas Numbers

//  Lucas(1000000) in 55 mS (Athlon X4 640, XP, 2GB)

//  http://gmplib.org/manual/Lucas-Numbers-Algorithm.html#Lucas-Numbers-Algorithm
//
//  mpz_lucnum2_ui derives a pair of Lucas numbers from a pair  
//  of Fibonacci numbers with the following simple formulas.
//  
//      L[k]   =   F[k] + 2*F[k-1]
//      L[k-1] = 2*F[k] -   F[k-1]
//  
//  mpz_lucnum_ui is only interested in L[n], and some work can be saved.
//  Trailing zero bits on n can be handled with a single square each.
//  
//      L[2k] = L[k]^2 - 2*(-1)^k
//  
//  And the lowest 1 bit can be handled with one multiply of a pair of
//  Fibonacci numbers, similar to what mpz_fib_ui does.
//  
//      L[2k+1] = 5*F[k-1]*(2*F[k]+F[k-1]) - 4*(-1)^k

//  Results:
//
//      mS(Lucs(n))   / mS(Fibs(n))   ~ 1.00
//      mS(Luc(2n+1)) / mS(Fib(2n+1)) ~ 1.00
//      mS(Luc(2n))   / mS(Fib(2n))   ~ 0.76
//      mS(Luc(n))    / ms(Fib(n))    ~ 0.88
//      mS(Lucs(n))   / ms(Luc(n))    ~ 1.55   

using Xint = System.Numerics.BigInteger;
using System.Threading.Tasks;
using System.Diagnostics;
using System;
class Lucas
{
    private static int[] lucSmall = { 2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, 199,
        322, 521, 843, 1364, 2207, 3571, 5778, 9349, 15127, 24476, 39603, 64079,
        103682, 167761, 271443, 439204, 710647, 1149851, 1860498, 3010349 };

    public static Xint[] Lucs(int n)                        //  { L(n), L(n + 1) }
    {
        if (n < 31) return new Xint[2] { lucSmall[n++], lucSmall[n] };
        Xint[] F = Fibs(n);
        return new Xint[2] { 2 * F[1] - F[0], F[1] + 2 * F[0] };
    }

    public static Xint Luc(int n)                                         //  L(n)
    {
        if (n < 32) return lucSmall[n];
        if ((n & 1) == 1) return Luc_n_is_odd(n);
        int z = fL2(n & -n);                                    //  trailing zeros
        int y = fL2(n) - z - 4;
        if (y < 0)
        {
            z += y;
            Xint L = lucSmall[n >> z];
            for (; z > 0; z--) L = SQ(L) - 2;
            return L;
        }
        else
        {
            Xint L = SQ(Luc_n_is_odd(n >> z)) + 2;
            for (; z > 1; z--) L = SQ(L) - 2;
            return L;
        }
    }

    private static Xint Luc_n_is_odd(int n)
    {
        n = n / 2 - 1;
        Xint[] F = Fibs(n);
        return 4 - 8 * (n & 1) + MTP(5 * F[0], 2 * F[1] + F[0]);
    }

    private static int[] fibSmall = { 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144,
        233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 
        75025, 121393, 196418, 317811, 514229, 832040, 1346269, 2178309, 3524578 };

    public static Xint[] Fibs(int n)                         //  { F(n), F(n + 1) }
    {
        if (n < 33) return new Xint[2] { fibSmall[n++], fibSmall[n] };
        n++;
        int i = fL2(n) - 4;
        int m = n >> i;
        Xint[] F = { fibSmall[m - 1], fibSmall[m] };
        m = ((1 << i) & n) >> i;
        for (--i; i >= 0; i--)
        {
            F = new Xint[2] { SQ(F[0]), SQ(F[1]) };
            F = new Xint[2] { F[1] + F[0], 2 - 4 * m + 4 * F[1] - F[0] };
            m = (1 << i & n) >> i;
            F[1 - m] = F[1] - F[0];
        }
        return F;
    }

    public static Xint Fib(int n)                                          //  F(n)
    {
        if (n < 34) return fibSmall[n];
        n++;
        int i = fL2(n) - 4;
        int m = n >> i;
        Xint[] F = { fibSmall[m - 1], fibSmall[m] };
        m = ((1 << i) & n) >> i;
        for (--i; i > 0; i--)
        {
            F = new Xint[2] { SQ(F[0]), SQ(F[1]) };
            F = new Xint[2] { F[1] + F[0], 2 - 4 * m + 4 * F[1] - F[0] };
            m = ((1 << i) & n) >> i;
            F[1 - m] = F[1] - F[0];
        }
        return (n & 1) == 0 ?
            2 - 4 * m + MTP(3 * F[1] + F[0], F[1] - F[0]) :
            MTP(F[1] + 2 * F[0], F[1]);
    }

    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        Luc(1000000);
        for (int i = 0; i < 10; i++)
        {
            sw.Restart(); Luc(1000000); sw.Stop();
            Console.WriteLine(sw.ElapsedMilliseconds);
        }
        Console.ReadLine();
    }

    // use the easy to copy faster versions:
    // http://bigintegers.blogspot.com/2011/fibonacci-numbers-part-1.html
    private static Xint MTP(Xint U, Xint V) { return U * V; }
    private static Xint SQ(Xint U) { return U * U; }
    private static int fL2(int n) { int i = -1; while (n > 0) { i++; n /= 2; } return i; }
}