Pagina's

2011/07/09

Fibonacci numbers part 2

//  Fibonacci(1000000) in 87 mS (Athlon X4 640, XP, 2GB)
//  Two squares become one multiply.

//                             +------------------------------------+
//                             |   Fibonacci.exe in milliSeconds    |
//                             |------------------------------------|
//                             |    n    | Fib(n)| Fib(n+1)| Fibs(n)|
//                             |---------|-------|---------|--------|
//                             |  156250 |     8 |       8 |     11 |
//                             |  312500 |    18 |      18 |     27 |
//                             |  625000 |    47 |      47 |     68 |
//                             | 1250000 |   127 |     127 |    183 |
//                             | 2500000 |   332 |     334 |    492 |
//                             | 5000000 |   930 |     930 |   1360 |
//                             |10000000 |  2500 |    2510 |   3650 | 
//                             |20000000 |  6770 |    6770 |   9900 |
//                             +------------------------------------+
//                           
//                             mS( Fib(n))      / mS(Fib(n+1))     ~ 1.00
//                             mS(Fibs(n))      / mS(Fib(n))       ~ 1.46
//                    
//  part_2 compared to part_1: mS(Fibs_2(n))    / mS(Fibs_1(n))    ~ 0.74   
//                             mS( Fib_2(2n))   / mS( Fib_1(2n))   ~ 0.83
//                             mS( Fib_2(2n+1)) / mS( Fib_1(2n+1)) ~ 0.65

//  Code translated to colorized HTML by http://puzzleware.net/CodeHTMLer/    

//            !!! Project >> Add Reference >> System.Numerics !!!

using Xint = System.Numerics.BigInteger;
using System.Threading.Tasks;
using System.Diagnostics;
using System;
class Fibonacci
{
    public static Xint[] Fibs(int n)  //  returns { F(n), F(n + 1) }
    {
        if (n < 33) return new Xint[2] { fibSmall(n++), fibSmall(n) };
        int i = fL2(n) - 4;
        int m = n >> i;
        Xint[] F = { fibSmall(m++), fibSmall(m) };
        for (--i; i >= 0; i--)
        {
            Xint G = 2 * F[1] - F[0];
            switch ((3 << i & n) >> i)
            {
                case 0: F = new Xint[2] { MTP(G, F[0]), MTP(G, F[1]) - 1 }; break;
                case 2: F = new Xint[2] { MTP(G, F[0]), MTP(G, F[1]) + 1 }; break;
                case 1: F = new Xint[2] { MTP(G, F[1]) - 1, MTP(G, F[1] + F[0]) - 1 }; break;
                case 3: F = new Xint[2] { MTP(G, F[1]) + 1, MTP(G, F[1] + F[0]) + 1 }; break;
            }
        }
        return F;
    }

    private static int fibSmall(int n)
    {
        int[] f = { 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 };
        return f[n];
    }

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

    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        Xint[] FL = { 0, 1 };
        for (int n = 2; n < 20000; n += 2)
        {
            Xint[] FH = Fibs(n);
            if (FH[0] != FL[0] + FL[1]) Console.WriteLine(n);
            if (FH[1] != FH[0] + FL[1]) Console.WriteLine(n);
            FL = FH;
        }
        Console.WriteLine("R E A D Y");
        Console.ReadLine();
    }

    // use the easy to copy faster versions!
    private static Xint MTP(Xint U, Xint V) { return U * V; }
    private static int fL2(int n) { int i = -1; while (n > 0) { i++; n /= 2; } return i; }
}

No comments:

Post a Comment