Pagina's

2011/09/18

Lucas Numbers

//  Lucas(1000000) in 55 mS (Athlon X4 640, XP, 2GB)

//  http://gmplib.org/manual/Lucas-Numbers-Algorithm.html#Lucas-Numbers-Algorithm
//
//  mpz_lucnum2_ui derives a pair of Lucas numbers from a pair  
//  of Fibonacci numbers with the following simple formulas.
//  
//      L[k]   =   F[k] + 2*F[k-1]
//      L[k-1] = 2*F[k] -   F[k-1]
//  
//  mpz_lucnum_ui is only interested in L[n], and some work can be saved.
//  Trailing zero bits on n can be handled with a single square each.
//  
//      L[2k] = L[k]^2 - 2*(-1)^k
//  
//  And the lowest 1 bit can be handled with one multiply of a pair of
//  Fibonacci numbers, similar to what mpz_fib_ui does.
//  
//      L[2k+1] = 5*F[k-1]*(2*F[k]+F[k-1]) - 4*(-1)^k

//  Results:
//
//      mS(Lucs(n))   / mS(Fibs(n))   ~ 1.00
//      mS(Luc(2n+1)) / mS(Fib(2n+1)) ~ 1.00
//      mS(Luc(2n))   / mS(Fib(2n))   ~ 0.76
//      mS(Luc(n))    / ms(Fib(n))    ~ 0.88
//      mS(Lucs(n))   / ms(Luc(n))    ~ 1.55   

using Xint = System.Numerics.BigInteger;
using System.Threading.Tasks;
using System.Diagnostics;
using System;
class Lucas
{
    private static int[] lucSmall = { 2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, 199,
        322, 521, 843, 1364, 2207, 3571, 5778, 9349, 15127, 24476, 39603, 64079,
        103682, 167761, 271443, 439204, 710647, 1149851, 1860498, 3010349 };

    public static Xint[] Lucs(int n)                        //  { L(n), L(n + 1) }
    {
        if (n < 31) return new Xint[2] { lucSmall[n++], lucSmall[n] };
        Xint[] F = Fibs(n);
        return new Xint[2] { 2 * F[1] - F[0], F[1] + 2 * F[0] };
    }

    public static Xint Luc(int n)                                         //  L(n)
    {
        if (n < 32) return lucSmall[n];
        if ((n & 1) == 1) return Luc_n_is_odd(n);
        int z = fL2(n & -n);                                    //  trailing zeros
        int y = fL2(n) - z - 4;
        if (y < 0)
        {
            z += y;
            Xint L = lucSmall[n >> z];
            for (; z > 0; z--) L = SQ(L) - 2;
            return L;
        }
        else
        {
            Xint L = SQ(Luc_n_is_odd(n >> z)) + 2;
            for (; z > 1; z--) L = SQ(L) - 2;
            return L;
        }
    }

    private static Xint Luc_n_is_odd(int n)
    {
        n = n / 2 - 1;
        Xint[] F = Fibs(n);
        return 4 - 8 * (n & 1) + MTP(5 * F[0], 2 * F[1] + F[0]);
    }

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

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

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

    // 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; while (n > 0) { i++; n /= 2; } return i; }
}

2011/09/05

Inverse Fibonacci filter(1/30.000.000)

// Another simple test for Fibonacci numbers:
// Take the repeated digit sum of a Fibonacci number.
// F11=89  8+9=17 1+7=8
// The repeated digit sum for odd Fibonacci numbers never is 2 or 7.
// Why? I asked Dr Math:
// http://mathforum.org/dr.math/
// The repeated digit sum is the same as taking the remainder after division by 9.
// Next "proof":

//             F%9    
//  F0       0  0------<----
//  F1       1  1--<--      |
//  F2       1  1     |     |
//  F3       2  2     |     |
//  F4       3  3     |     |
//  F5       5  5     |     |
//  F6       8  8     |     |
//  F7      13  4     |     |
//  F8      21  3     |     |
//  F9      34  7     |     |
//  F10     55  1     |     |
//  F11     89  8     |     |   (F(n+2))%9 = ( (F(n+1))%9 + (F(n))%9 ) % 9
//  F12    144  0     |     |
//  F13    233  8     |     |
//  F14    377  8     |     |           377%9 =          8
//  F15    610  7     |     |                   610%9 =    7
//  F16    987  6     |     |   987%9 = 377%9 + 610%9 = (8+7)%9=6
//  F17   1597  4     |     |
//  F18   2584  1     |     |
//  F19   4181  5     |     |
//  F20   6765  6     |     |
//  F21  10946  2     |     |
//  F22  17711  8     |     |
//  F23  28657  1-->--      |
//  F24  46368  0------>----

// Division by 8 gives 6 different remainders,
// 0%8=0, 1%8=1, 1%8=1, 2%8=2, 3%8=3, 5%8=5, 8%8=0, 5%8=5, 5%8=5,
// 10%8=2,7%8=7, 9%8=1, 8%8=0....6 different remainders(0,1,2,3,5,7)
// Division by 11 gives 7 different remainders.
// Division by 18 gives 11 different remainders.
// Division by F46 gives 66 different remainders.
// Which means the chance a number is a Fibonacci number
// becomes about 1 to 30.000.000(~F46/66).

//        +-----------------------------------+
//        |                      | different  | 
//        |         devisor      | remainders |
//        |----------------------|------------|
//        |          8 | F6      |     6      |
//        |         11 | L5      |     7      |
//        |         18 | L6      |    11      |
//        |         21 | F8      |     9      |
//        |         29 | L7      |    10      |   10/29 < 9/21
//        |         38 | F18/68  |    13      |
//        |         47 | L8      |    15      |
//        |         48 | F12/3   |    15      |
//        |         55 | F10     |    12      |   12*48 < 15*55
//        |      ..... | ....    |    ..      |
//        |   ........ | .....   |    ..      |
//        | 1201881744 | L45/2   |    69      |
//        | 1268860318 | F56/147 |    67      |
//        | 1536404311 | L55/199 |    81      |
//        | 1602508992 | F48/3   |    69      |
//        | 1836311903 | F46     |    66      |
//        +-----------------------------------+

// The devisors found are closely related to the closely related
// Fibonacci and Lucas numbers, see the table below. 
// F12(1,2,3) is F12/1, F12/2, F12/3 is 144, 72, 48 
//
//   +--------------------------------------------------------------------+
//   | F6(1)                   | F32(1,3,7,21,47)                         |
//   | F8(1)                   | F34(1)                                   |
//   | F10(1)                  | F36(1,2,3,4,6,8,9,18,24,48,51)           |
//   | F12(1,2,3)              | F38(1)                                   |
//   | F14(1)                  | F39(2)                                   |
//   | F16(1,3,7)              | F40(3,5,7,11,15,21,35,55)                |
//   | F18(1,2)                | F42(1,2,4,8,26)                          |
//   | F20(1,3)                | F44(1,3)                                 |
//   | F22(1)                  | F46(1,139)                               |
//   | F24(1,2,3,4,7,8,9,48)   | F48(3,4,6,7,8,9,14,18,21,23,24,28,36,47, |
//   | F26(1)                  |     48,56,92,94,96,139,144,161,168,329)  |
//   | F27(2)                  | F50(11,25,55)                            |
//   | F28(1,3)                | F54(436,872)                             |
//   | F30(1,2,4,5,8,10,20,22) | F56(147)                                 |
//   |-------------------------|------------------------------------------|          
//   | L5(1)                   | L24(1)                                   |
//   | L6(1)                   | L25(1,7,11)                              |
//   | L7(1)                   | L27(1,2,4)                               |
//   | L8(1)                   | L28(1)                                   |
//   | L9(1,2)                 | L29(1)                                   |
//   | L10(1)                  | L31(1)                                   |
//   | L11(1)                  | L33(1,2,4)                               |
//   | L12(1)                  | L34(1)                                   |
//   | L13(1)                  | L35(1,29,71)                             |
//   | L14(1)                  | L36(1)                                   |
//   | L15(1)                  | L37(1)                                   |
//   | L16(1)                  | L39(1,2,4)                               |
//   | L17(1)                  | L41(1)                                   |
//   | L18(1)                  | L43(1)                                   |
//   | L19(1)                  | L45(2,4,11,19,31,38,76,124,181,209)      |
//   | L21(1,2,4)              | L49(29)                                  |
//   | L22(1)                  | L55(199)                                 |
//   | L23(1)                  |                                          |
//   +--------------------------------------------------------------------+
//
// F24/48=L12*3, F32/47=L16*21, probably there are more.

// For X being a (pseudo) random number with 16.000.000 bits,
// invFib(X) took 78 mS, now it takes less than 6 mS. 

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

class inverseFibonacci
{
    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;
    private const double c3 = 35.73932166425341;

    private static int invFib(Xint X)
    {
        if (X <= int.MaxValue) 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 21: return 8;
                case 34: return 9;
                case 55: return 10;
                case 89: return 11;
                case 144: return 12;
                case 233: return 13;
                case 377: return 14;
                case 610: return 15;
                case 987: return 16;
                case 1597: return 17;
                case 2584: return 18;
                case 4181: return 19;
                case 6765: return 20;
                case 10946: return 21;
                case 17711: return 22;
                case 28657: return 23;
                case 46368: return 24;
                case 75025: return 25;
                case 121393: return 26;
                case 196418: return 27;
                case 317811: return 28;
                case 514229: return 29;
                case 832040: return 30;
                case 1346269: return 31;
                case 2178309: return 32;
                case 3524578: return 33;
                case 5702887: return 34;
                case 9227465: return 35;
                case 14930352: return 36;
                case 24157817: return 37;
                case 39088169: return 38;
                case 63245986: return 39;
                case 102334155: return 40;
                case 165580141: return 41;
                case 267914296: return 42;
                case 433494437: return 43;
                case 701408733: return 44;
                case 1134903170: return 45;
                case 1836311903: return 46;
                default: return -1;
            }

        if (X.IsEven)
        {
            int x = (int)(X & 31);
            if ((x != 2) && ((x & 7) != 0)) return -1;
        }

        int r = (int)Xint.Remainder(X, 1836311903);
        if ((r != 0) && (r != 1) && (r != 1134903170) &&
            (r != 2) && (r != 5) && (r != 13) && (r != 34) && (r != 89) && (r != 233) &&
            (r != 610) && (r != 1597) && (r != 4181) && (r != 10946) && (r != 28657) &&
            (r != 75025) && (r != 196418) && (r != 514229) && (r != 1346269) &&
            (r != 3524578) && (r != 9227465) && (r != 24157817) &&
            (r != 63245986) && (r != 165580141) && (r != 433494437) &&
            (r != 3) && (r != 8) && (r != 21) && (r != 55) && (r != 144) &&
            (r != 377) && (r != 987) && (r != 2584) && (r != 6765) && (r != 17711) &&
            (r != 46368) && (r != 121393) && (r != 317811) && (r != 832040) &&
            (r != 2178309) && (r != 5702887) && (r != 14930352) && (r != 39088169) &&
            (r != 102334155) && (r != 267914296) && (r != 701408733) &&
            (r != 1568397607) && (r != 1733977748) && (r != 1797223734) &&
            (r != 1821381551) && (r != 1830609016) && (r != 1834133594) &&
            (r != 1835479863) && (r != 1835994092) && (r != 1836190510) &&
            (r != 1836265535) && (r != 1836294192) && (r != 1836305138) &&
            (r != 1836309319) && (r != 1836310916) && (r != 1836311526) &&
            (r != 1836311759) && (r != 1836311848) && (r != 1836311882) &&
            (r != 1836311895) && (r != 1836311900) && (r != 1836311902))
            return -1;

        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++)
            if (n % i == 0)
                if (Xint.Remainder(X, fibSmall[i]) > 0)
                    return false;
        return X == Fib(n) ? true : false;
    }

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

    // 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; while (n > 0) { i++; n /= 2; } return i; }

    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        Xint X = RND(16000000);
        invFib(X);
        sw.Restart();
        for (int i = 0; i < 100; i++)
        {
            if (invFib(X) != -1) Console.WriteLine(i);
        }
        sw.Stop();
        Console.WriteLine(sw.ElapsedMilliseconds);
        Console.ReadLine();
    }

    private static int seed;
    public static Xint RND(int n)
    {
        if (n < 2) return n;
        if (seed == int.MaxValue) seed = 0; else seed++;
        Random rand = new Random(seed);
        byte[] bytes = new byte[(n + 15) / 8];
        rand.NextBytes(bytes);
        int i = bytes.Length - 1;
        bytes[i] = 0;
        n = i * 8 - n;
        i--;
        bytes[i] >>= n;
        bytes[i] |= (byte)(128 >> n);
        return new Xint(bytes);
    }
}