Pell Numbers

//  Pell(1000000) in 207 mS (Athlon X4 640, XP, 2GB)

//  Pell Numbers are defined by:
//  P[0] = 0, P[1] = 1, P[n+2] = 2*P[n+1] + P[n] 
//  The sequence is 0,1,2,5,12,29,70,169,...    

//  Formula's used are:
//     (P[n])^2 = P[n+1]*P[n-1] - (-1)^n
//      P[m+n]  = P[m]*P[n+1] + P[m-1]*P[n]

//  Small values of n use a table.
//  Larger values use a binary powering approach,
//  which scannes n from left to right, 
//  from most to least significant bit.
//  Pells(n) returns P[n] and P[n+1]
//  Each bit requires one square and one multiply.
//  Pell(n) returns P[n]
//  The least significant bit requires one multiply. 

//  Result: mS(Pells(n)) / mS(Pell(n)) ~ 1.41 
//                                     ~ square root of two!?

using Xint = System.Numerics.BigInteger;
using System.Threading.Tasks;
using System.Diagnostics;
using System;

class Pell_Numbers
{
    private static int[] smallPell = { 0, 1, 2, 5, 12, 29, 70, 169,
        408, 985, 2378, 5741, 13860, 33461, 80782, 195025, 470832 };

    public static Xint[] Pells(int n)
    {
        if (n < 16) return new Xint[2] { smallPell[n++], smallPell[n] };
        int i = fL2(n) - 4;
        int m = (n >> i) / 2 & 7 + 8;
        Xint[] P = { smallPell[m++], smallPell[m] };
        for (; i >= 0; i--)
        {
            switch (((3 << i) & n) >> i)
            {
                case 0:
                    P = new Xint[2] { 2 * SQ(P[0]), 2 * MTP(P[1], P[0]) };
                    P = new Xint[2] { P[1] - P[0], P[1] + P[0] + 1 }; break;
                case 2:
                    P = new Xint[2] { 2 * SQ(P[0]), 2 * MTP(P[1], P[0]) };
                    P = new Xint[2] { P[1] - P[0], P[1] + P[0] - 1 }; break;
                case 1:
                    P = new Xint[2] { 2 * MTP(P[1], P[0]), 2 * SQ(P[1]) };
                    P = new Xint[2] { P[1] - P[0] - 1, P[1] + P[0] }; break;
                case 3:
                    P = new Xint[2] { 2 * MTP(P[1], P[0]), 2 * SQ(P[1]) };
                    P = new Xint[2] { P[1] - P[0] + 1, P[1] + P[0] }; break;
            }
        }
        return P;
    }

    public static Xint Pell(int n)
    {
        if (n < 17) return smallPell[n];
        Xint[] P = Pells(n / 2);
        return (n & 1) == 0 ?
            2 * MTP(P[1] - P[0], P[0]) :
            2 * MTP(P[1] + P[0], P[0]) + (1 - (n & 2));
    }

    // 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; for (; n > 0; n /= 2) i++; return i; }

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