Pagina's

2011/10/11

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();
    }
}

2011/10/06

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();
    }
}