BigInteger Algorithms: 2011

2011/12/29

Factorial by binary splitting, part 2

using Xint = System.Numerics.BigInteger;
using System;
class Factorial
{
    public static Xint F(uint n)
    {
        uint a = 0;
        uint s = 0;
        Xint P = 1;
        Xint Q = 1;
        uint b = 1;
        for (int i = fL2((int)(n / 2)); i >= 0; i--)
        {
            a = n >> i;
            s = s + a / 2;
            a = a - 1 | 1;
            P = Q * P;
            Q = OddP(a, b) * Q;
            b = a + 2;
        }
        return Q * P << (int)s;
    }

    private static Xint OddP(uint a, uint b)
    {
        if (a == b) return a;
        uint m = (a + b) / 2;
        m += m & 1;
        return OddP(a, m + 1) * OddP(m - 1, b);
    }

    private static int fL2(int n) { int i = -1; for (; n > 0; n /= 2) i++; return i; }

    // 42!=42*41*..
    //    =41*39*...*3*1 * 42*40*...*4*2 
    //    =41,1? * 42,2?
    //             42,2?=21! *                                                                  2^(42/2)
    //                   21!=21*19*...*3*1 * 20*18*...*4*2                                               
    //                      =21,1? * 20,2?                                                               
    //                               20,2?=10! *                                                2^(20/2)
    //                                     10!=9*7*5*3*1 * 10*8*6*4*2                                
    //                                        =9,1? * 10,2?                                          
    //                                                10,2?=5! *                                2^(10/2)
    //                                                      5!=5*3*1 * 4*2                            
    //                                                        =5,1?  * 4,2?                                  
    //                                                                 4,2?=2! *                2^( 4/2)
    //                                                                      2!=1*1  * 2                
    //                                                                        =1,1? * 2,2?                          
    //                                                                                2,2?=1! * 2^( 2/2)
    //        
    //    = 41,1?  * 21,1?  * 9,1?     * 5,1?   * 2^(42/2+20/2+10/2+4/2+2/2)
    //    = 5,1?   * 9,1?   * 21,1?    * 41,1?  * 2^(21+10+5+2+1)
    //    = 5,1?^4 * 9,7?^3 * 21,11?^2 * 41,23? * 2^39
    //    =    a^4 *    b^3 *      c^2 *    d^1 << 39
    //
    //    = 1                    *  a      
    //    = a^1                  *  a*b    
    //    = a^2 * b^1            *  a*b*c  
    //    = a^3 * b^2 * c^1      *  a*b*c*d
    //    = a^4 * b^3 * c^2* d^1                << 39

    static void Main()
    {
        Console.WriteLine(F(0));
        Console.WriteLine(F(100));
        Console.ReadLine();
    }
}

2011/10/11

Half Companion Pell Numbers

//  Pell(1000000) in 199 mS (Athlon X4 640, XP, 2GB)

//  Half Companion Pell Numbers are defined by:
//  H[0] = 1, H[1] = 1, H[n+2]= 2*H[n+1] + H[n]
//  The sequence is 1,1,3,7,17,41,99,239,...

//  Small values of n use a table.
//  Larger values use a binary powering approach.
//  HCPs(n) returns H[n] and H[n+1]
//  The most important formula used is:
//
//      H[2n] = 2 * H[n]^2 - (-1)^n
//
//  Each bit requires two squares.
//  HCP(n) returns H[n]
//  If n is odd the lowest 1 bit requires a multiply.
//  Trailing zero bits are done by a square each.

//  Pells(n) returns P[n] and P[n+1]
//  It derives a pair of Pell numbers 
//  from a pair of HCP numbers:
//
//      2*P[n]   = H[n+1] - H[n]
//      2*P[n+1] = H[n+1] + H[n]
//
//  Pell(n) returns P[n]
//  The least significant bit requires a multiply.

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

class Half_Companion_Pell_Numbers
{
    private static int[] smallHCP = { 1, 1, 3, 7, 17, 41, 99, 239,
        577, 1393, 3363, 8119, 19601, 47321, 114243, 275807, 665857 };

    public static Xint[] HCPs(int n)
    {
        if (n < 16) return new Xint[2] { smallHCP[n++], smallHCP[n] };
        int i = fL2(n) - 4;
        int m = (n >> i) / 2 & 7 + 8;
        Xint H0 = smallHCP[m++], H1 = smallHCP[m];
        for (; i >= 0; i--)
        {
            H0 = SQ(H0); H1 = SQ(H1);
            switch (((3 << i) & n) >> i)
            {
                case 0: H1 = H1 - H0 + 1; H0 = 2 * H0 - 1; break;
                case 1: H0 = H1 - H0 + 1; H1 = 2 * H1 + 1; break;
                case 2: H1 = H1 - H0 - 1; H0 = 2 * H0 + 1; break;
                case 3: H0 = H1 - H0 - 1; H1 = 2 * H1 - 1; break;
            }
        }
        return new Xint[2] { H0, H1 };
    }

    public static Xint HCP(int n)
    {
        if (n < 17) return smallHCP[n];
        if ((n & 1) == 1) return HCP_n_is_odd(n);
        int z = fL2(n & -n);
        int y = fL2(n) - z - 3;
        if (y < 0)
        {
            z += y;
            Xint H = smallHCP[n >> z];
            for (; z > 0; z--) H = 2 * SQ(H) - 1;
            return H;
        }
        else
        {
            Xint H = 2 * SQ(HCP_n_is_odd(n >> z)) + 1;
            for (; z > 1; z--) H = 2 * SQ(H) - 1;
            return H;
        }
    }

    private static Xint HCP_n_is_odd(int n)
    {
        Xint[] H = HCPs(n / 2);
        return (n & 2) - 1 + 2 * MTP(H[1], H[0]);
    }

    private static int[] smallPell = { 0, 1, 2, 5, 12, 29, 70, 169,
        408, 985, 2378, 5741, 13860, 33461, 80782, 195025, 470832 };

    public static Xint[] Pells(int n)
    {
        if (n < 16) return new Xint[2] { smallPell[n++], smallPell[n] };
        Xint[] H = HCPs(n);
        return new Xint[2] { (H[1] - H[0]) / 2, (H[1] + H[0]) / 2 };
    }

    public static Xint Pell(int n)
    {
        if (n < 17) return smallPell[n];
        Xint[] H = HCPs(n / 2);
        return (n & 1) == 0 ?
            MTP(H[1] - H[0], H[0]) :
            MTP(H[1] + H[0], H[0]) + ((n & 2) - 1);
    }

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

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

2011/10/06

Pell Numbers

//  Pell(1000000) in 207 mS (Athlon X4 640, XP, 2GB)

//  Pell Numbers are defined by:
//  P[0] = 0, P[1] = 1, P[n+2] = 2*P[n+1] + P[n] 
//  The sequence is 0,1,2,5,12,29,70,169,...    

//  Formula's used are:
//     (P[n])^2 = P[n+1]*P[n-1] - (-1)^n
//      P[m+n]  = P[m]*P[n+1] + P[m-1]*P[n]

//  Small values of n use a table.
//  Larger values use a binary powering approach,
//  which scannes n from left to right, 
//  from most to least significant bit.
//  Pells(n) returns P[n] and P[n+1]
//  Each bit requires one square and one multiply.
//  Pell(n) returns P[n]
//  The least significant bit requires one multiply. 

//  Result: mS(Pells(n)) / mS(Pell(n)) ~ 1.41 
//                                     ~ square root of two!?

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

class Pell_Numbers
{
    private static int[] smallPell = { 0, 1, 2, 5, 12, 29, 70, 169,
        408, 985, 2378, 5741, 13860, 33461, 80782, 195025, 470832 };

    public static Xint[] Pells(int n)
    {
        if (n < 16) return new Xint[2] { smallPell[n++], smallPell[n] };
        int i = fL2(n) - 4;
        int m = (n >> i) / 2 & 7 + 8;
        Xint[] P = { smallPell[m++], smallPell[m] };
        for (; i >= 0; i--)
        {
            switch (((3 << i) & n) >> i)
            {
                case 0:
                    P = new Xint[2] { 2 * SQ(P[0]), 2 * MTP(P[1], P[0]) };
                    P = new Xint[2] { P[1] - P[0], P[1] + P[0] + 1 }; break;
                case 2:
                    P = new Xint[2] { 2 * SQ(P[0]), 2 * MTP(P[1], P[0]) };
                    P = new Xint[2] { P[1] - P[0], P[1] + P[0] - 1 }; break;
                case 1:
                    P = new Xint[2] { 2 * MTP(P[1], P[0]), 2 * SQ(P[1]) };
                    P = new Xint[2] { P[1] - P[0] - 1, P[1] + P[0] }; break;
                case 3:
                    P = new Xint[2] { 2 * MTP(P[1], P[0]), 2 * SQ(P[1]) };
                    P = new Xint[2] { P[1] - P[0] + 1, P[1] + P[0] }; break;
            }
        }
        return P;
    }

    public static Xint Pell(int n)
    {
        if (n < 17) return smallPell[n];
        Xint[] P = Pells(n / 2);
        return (n & 1) == 0 ?
            2 * MTP(P[1] - P[0], P[0]) :
            2 * MTP(P[1] + P[0], P[0]) + (1 - (n & 2));
    }

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

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

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

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

Pagina's

2011/12/29

2011/10/11

2011/10/06

2011/09/18

2011/09/05

2011/08/10

2011/08/06