Pagina's

2011/08/10

Faster inverse Fibonacci

A simple test to see if a number X is not a Fibonacci number.
Even Fibonacci numbers end binary with 000 or 010

The "proof" starts with F0=000 and F1=001
F2=F1+F0, F3=F2+F1, etc.

 F0  000------<----
 F1  001--<--      |
 F2  001     |     |
 F3  010     |     |
 F4  011     |     |
 F5  101     |     |
 F6  000     |     |
 F7  101     |     |
 F8  101     |     |
 F9  010     |     |
 F10 111     |     |
 F11 001-->--      |
 F12 000------>----

After 11 additions we're back at the starting point.
Even numbers end with 000 or 010 or 100 or 110
So 50% of the even numbers will be recognized as not
Fibonacci numbers. 50% Of the numbers are even.
So 25% will be recognized as not Fibonacci numbers,
just looking at the 3 least significant bits.

Let's look at the 5 least significant bits,
again a pattern can be found,
and a formula:
if (((X & 1)) == 0) && ((X & 7) != 0) && (((X - 2) & 31) != 0))
then X is not a Fibonacci number.
11 Of 32 numbers are excluded, that's 34.375%.
The average computing time decreases by about one third.

I tried to look at more than 5 least significant bits,
sorry, to hard for me today. I expect it's possible to exclude
nearly 3 of 8 numbers, an optimum close to 37.5%.
Is there a trick for odd Numbers? Stay tuned!

Another improvement for small X values <= Fib(46)(= 1836311093),
the largest integer Fibonacci number:

switch ((int)X))
{
    case 0: return 0;
    case 2: return 3;
    case 3: return 4;
    case 5: return 5;
    case 8: return 6;
    case 13: return 7;
//  ..................
//  ......................
//  ...........................
    case 1134903170: return 45;
    case 1836311093: return 46;
    default: return -1;
}

It's more than two times faster.
Isn't it a waste to hard-code all those Fibonacci numbers?
Isn't it a waste to use an array of Fibonacci numbers,
int[] fibSmall = { 0, 1 , 1, 2, 3, 5,....}
It's sufficient to use:
int[] fibSmall = { 0, 1 }
The code length is nothing compared to the magnitude of the
numbers computed. I will go for the fast way.

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

Inverse Fibonacci, invFib(Fib(n))=n

// Detecting when a number is a Fibonacci number, and which one if so.
// It's based on Binet's formula:
// Fn = (((1 + 5^0.5)/2)^n - ((1 - 5^0.5)/2)^n) / 5^0.5
// For larger n:
// Fn ~ (1 + 5^0.5)/2)^n  / 5^0.5
// For small numbers , <= F33 , a binary search will be used,
// for larger numbers the code below was used to get error value's.

    private static void error()
    {
        double c4 = Math.Log(Math.Sqrt(5));
        double c5 = Math.Log(Math.Sqrt(5) + 1) - Math.Log(2);
        Xint[] F = { 0, 1 };
        double e0 = 0;
        for (int n = 1; n < int.MaxValue - 1; n++)
        {
            F = new Xint[] { F[1], F[1] + F[0] };
            double m = (Xint.Log(F[0]) + c4) / c5;
            double e1 = Math.Abs(n - m);
            if ((n > 33) & (e1 > e0))
            {
                Console.WriteLine(n + "  " + e0);
                e0 = e1;
            }
        }
    }

// Below is a graph of log10(e) against log10(n),
// it's quite linear.