Pagina's

2011/07/09

Fibonacci numbers part 3

//  Fibonacci(1000000) in 81 mS (Athlon X4 640, XP, 2GB)
//  Two multiplies become two squares, like GMP.
//  http://gmplib.org/manual/Fibonacci-Numbers-Algorithm.html#Fibonacci-Numbers-Algorithm

//                             +------------------------------------+
//                             |   Fibonacci.exe in milliSeconds    |
//                             |------------------------------------|
//                             |    n    | Fib(n)| Fib(n+1)| Fibs(n)|
//                             |---------|-------|---------|--------|
//                             |  156250 |     7 |       7 |     10 |
//                             |  312500 |    17 |      17 |     24 |
//                             |  625000 |    46 |      46 |     62 |
//                             | 1250000 |   119 |     119 |    163 |
//                             | 2500000 |   312 |     313 |    434 |
//                             | 5000000 |   873 |     870 |   1200 |
//                             |10000000 |  2350 |    2350 |   3150 | 
//                             |20000000 |  6280 |    6290 |   8550 |
//                             +------------------------------------+
//                           
//                             mS( Fib(n))   / mS(Fib(n+1))  ~ 1.00
//                             mS(Fibs(n))   / mS(Fib(n))    ~ 1.36
//                    
//  part_3 compared to part_2: mS(Fibs_3(n)) / mS(Fibs_2(n)) ~ 0.88   
//                             mS( Fib_3(n)) / mS( Fib_2(n)) ~ 0.94

//  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] };
        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], 4 * F[1] - F[0] + (2 - 4 * m) };
            m = ((1 << i) & n) >> i;
            F[1 - m] = F[1] - F[0];
        }
        return F;
    }

    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 Fib(int n)  //  returns 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], 4 * F[1] - F[0] + (2 - 4 * m) };
            m = ((1 << i) & n) >> i;
            F[1 - m] = F[1] - F[0];
        }
        if ((n & 1) == 1) return MTP(F[1] + 2 * F[0], F[1]);
        return MTP(3 * F[1] + F[0], F[1] - F[0]) + (2 - 4 * m);
    }

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

    // use the easy to copy faster versions!
    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; }
}

No comments:

Post a Comment