Half Companion Pell Numbers

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

//  Half Companion Pell Numbers are defined by:
//  H[0] = 1, H[1] = 1, H[n+2]= 2*H[n+1] + H[n]
//  The sequence is 1,1,3,7,17,41,99,239,...

//  Small values of n use a table.
//  Larger values use a binary powering approach.
//  HCPs(n) returns H[n] and H[n+1]
//  The most important formula used is:
//
//      H[2n] = 2 * H[n]^2 - (-1)^n
//
//  Each bit requires two squares.
//  HCP(n) returns H[n]
//  If n is odd the lowest 1 bit requires a multiply.
//  Trailing zero bits are done by a square each.

//  Pells(n) returns P[n] and P[n+1]
//  It derives a pair of Pell numbers 
//  from a pair of HCP numbers:
//
//      2*P[n]   = H[n+1] - H[n]
//      2*P[n+1] = H[n+1] + H[n]
//
//  Pell(n) returns P[n]
//  The least significant bit requires a multiply.

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

class Half_Companion_Pell_Numbers
{
    private static int[] smallHCP = { 1, 1, 3, 7, 17, 41, 99, 239,
        577, 1393, 3363, 8119, 19601, 47321, 114243, 275807, 665857 };

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

    public static Xint HCP(int n)
    {
        if (n < 17) return smallHCP[n];
        if ((n & 1) == 1) return HCP_n_is_odd(n);
        int z = fL2(n & -n);
        int y = fL2(n) - z - 3;
        if (y < 0)
        {
            z += y;
            Xint H = smallHCP[n >> z];
            for (; z > 0; z--) H = 2 * SQ(H) - 1;
            return H;
        }
        else
        {
            Xint H = 2 * SQ(HCP_n_is_odd(n >> z)) + 1;
            for (; z > 1; z--) H = 2 * SQ(H) - 1;
            return H;
        }
    }

    private static Xint HCP_n_is_odd(int n)
    {
        Xint[] H = HCPs(n / 2);
        return (n & 2) - 1 + 2 * MTP(H[1], H[0]);
    }

    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] };
        Xint[] H = HCPs(n);
        return new Xint[2] { (H[1] - H[0]) / 2, (H[1] + H[0]) / 2 };
    }

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

    // 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()
    {
        Pell(1000000); sw.Restart();
        for (int n = 0; n < 10; n++) Pell(1000000);
        sw.Stop(); Console.WriteLine(sw.ElapsedMilliseconds);
        Console.ReadLine();
    }
}