Pagina's

2011/08/06

Inverse Fibonacci, code

// X=Fib(n) and n=invFib(X) have the same speed,
// If (pseudo) random numbers with the bitLength of X are used,
// invFib will be very fast.
// X=Fib(23.000.000) takes 7540 mS (X has 15.967.636 bits ~ 16M bits)
//         invFib(X) takes 7640 mS, random numbers 
// with X's bitLength take   78 mS.

// Another way to test Fibonacci numbers is:
// X is a Fibo if 5*X^2 + 4 or 5*X^2 - 4 is a square number.
// squaring a 16M bits number takes 8460 mS,
// taking the square root takes 25000 mS,
// Let's guess it's possible to take the other square root fast
// then it seems to be about 400 times slower when X is not a Fibo,
// 4 times slower when X is (still n will be unknown).

// Did a few tests, 
// Fib(500.000.000) takes 10 minutes (347.120.956 bits)
// Fib(600.000.000) crashes (2GB memory?)
// so n_max is somewhere in between 500M and 600M
// and in the line:
// if (m > 9.07086593448401E-07) return -1;
//         9.07086593448401E-07 should be:
// Math.Pow(Math.E, c2 * Math.Log(500000000) - c3)
// and for my machine I might change:
// if (m > int.MaxValue - 1) return -2;
// into:
// if (m > 500000000) return -2;


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

class inverseFibonacci
{
    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 };

    private static double c0 = Math.Log(Math.Sqrt(5));
    private static double c1 = Math.Log(Math.Sqrt(5) + 1) - Math.Log(2);
    private const double c2 = 1.0157640563849244; // 1?
    private const double c3 = 35.73932166425341;

    private static int invFib(Xint X)
    {
        if (X <= fibSmall[33])
        {
            int x = (int)X;
            if (x < 2) return x == 0 ? 0 : -1;
            int n = 3;
            for (int m = 16; m > 0; )
            {
                for (; (m > 0) & (x > fibSmall[n]); m /= 2) n += m;
                for (; (m > 0) & (x < fibSmall[n]); m /= 2) n -= m;
                if (x == fibSmall[n]) return n;
            }
            return -1;
        }
        else
        {
            double m = (Xint.Log(X) + c0) / c1;
            if (m > int.MaxValue - 1) return -2;
            int n = (int)Math.Round(m);
            m = Math.Abs(n - m);
            if (m > 9.07086593448401E-07) return -1;
            if (n < 69) return m < 7.105427357602E-15 ? n : -1;
            if (n < 129) return m < 1.4210854715203E-14 && isFib(X, n) ? n : -1;
            if (n < 169) return m < 2.8421709430405E-14 && isFib(X, n) ? n : -1;
            if (n < 338) return m < 5.6843418860809E-14 && isFib(X, n) ? n : -1;
            if (n < 1026) return m < 1.13686837721617E-13 && isFib(X, n) ? n : -1;
            if (n < 1087) return m < 2.27373675443233E-13 && isFib(X, n) ? n : -1;
            if (n < 4120) return m < 4.54747350886465E-13 && isFib(X, n) ? n : -1;
            if (n < 7593) return m < 9.09494701772929E-13 && isFib(X, n) ? n : -1;
            if (n < 8522) return m < 1.81898940354587E-12 && isFib(X, n) ? n : -1;
            if (n < 17136) return m < 3.63797880709172E-12 && isFib(X, n) ? n : -1;
            if (n < 35166) return m < 7.27595761418344E-12 && isFib(X, n) ? n : -1;
            if (n < 131083) return m < 1.45519152283670E-11 && isFib(X, n) ? n : -1;
            if (n < 259737) return m < 2.91038304567338E-11 && isFib(X, n) ? n : -1;
            if (n < 272394) return m < 5.82076609134675E-11 && isFib(X, n) ? n : -1;
            if (n < 1048589) return m < 1.16415321826936E-10 && isFib(X, n) ? n : -1;
            if (n < 1823140) return m < 2.32830643653871E-10 && isFib(X, n) ? n : -1;
            if (n < 2179038) return m < 4.65661287307740E-10 && isFib(X, n) ? n : -1;
            if (n < 8388610) return m < 9.31322574615480E-10 && isFib(X, n) ? n : -1;
            if (n < 14011307) return m < 1.86264514923097E-09 && isFib(X, n) ? n : -1;
            if (n < 24678101) return m < 3.72529029846192E-09 && isFib(X, n) ? n : -1;
            return m < Math.Pow(Math.E, c2 * Math.Log(n) - c3) && isFib(X, n) ? n : -1;
        }
    }

    private static bool isFib(Xint X, int n)
    {
        for (int i = 3; i < 34; i++) // Fnk is a multiple of Fn
            if (n % i == 0)
                if (Xint.Remainder(X, fibSmall[i]) > 0)
                    return false;
        return X == Fib(n) ? true : false;
    }

    private static Xint Fib(int 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);
    }

    // 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; }

    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        sw.Restart();
        Xint X = Fib(1000000);
        sw.Stop();
        Console.WriteLine(sw.ElapsedMilliseconds);
        sw.Restart();
        int n = invFib(X);
        sw.Stop();
        Console.WriteLine(sw.ElapsedMilliseconds);
        Console.ReadLine();
    }
}

No comments:

Post a Comment