BigInteger Algorithms

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

Fibonacci numbers part 2

//  Fibonacci(1000000) in 87 mS (Athlon X4 640, XP, 2GB)
//  Two squares become one multiply.

//                             +------------------------------------+
//                             |   Fibonacci.exe in milliSeconds    |
//                             |------------------------------------|
//                             |    n    | Fib(n)| Fib(n+1)| Fibs(n)|
//                             |---------|-------|---------|--------|
//                             |  156250 |     8 |       8 |     11 |
//                             |  312500 |    18 |      18 |     27 |
//                             |  625000 |    47 |      47 |     68 |
//                             | 1250000 |   127 |     127 |    183 |
//                             | 2500000 |   332 |     334 |    492 |
//                             | 5000000 |   930 |     930 |   1360 |
//                             |10000000 |  2500 |    2510 |   3650 | 
//                             |20000000 |  6770 |    6770 |   9900 |
//                             +------------------------------------+
//                           
//                             mS( Fib(n))      / mS(Fib(n+1))     ~ 1.00
//                             mS(Fibs(n))      / mS(Fib(n))       ~ 1.46
//                    
//  part_2 compared to part_1: mS(Fibs_2(n))    / mS(Fibs_1(n))    ~ 0.74   
//                             mS( Fib_2(2n))   / mS( Fib_1(2n))   ~ 0.83
//                             mS( Fib_2(2n+1)) / mS( Fib_1(2n+1)) ~ 0.65

//  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) };
        int i = fL2(n) - 4;
        int m = n >> i;
        Xint[] F = { fibSmall(m++), fibSmall(m) };
        for (--i; i >= 0; i--)
        {
            Xint G = 2 * F[1] - F[0];
            switch ((3 << i & n) >> i)
            {
                case 0: F = new Xint[2] { MTP(G, F[0]), MTP(G, F[1]) - 1 }; break;
                case 2: F = new Xint[2] { MTP(G, F[0]), MTP(G, F[1]) + 1 }; break;
                case 1: F = new Xint[2] { MTP(G, F[1]) - 1, MTP(G, F[1] + F[0]) - 1 }; break;
                case 3: F = new Xint[2] { MTP(G, F[1]) + 1, MTP(G, F[1] + F[0]) + 1 }; break;
            }
        }
        return F;
    }

    private static int fibSmall(int n)
    {
        int[] f = { 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 };
        return f[n];
    }

    public static Xint Fib(int n)  //  returns F(n)
    {
        if (n < 34) return fibSmall(n);
        int i = fL2(n) - 4;
        int m = n >> i;
        Xint[] F = { fibSmall(m++), fibSmall(m) };
        for (--i; i > 0; i--)
        {
            Xint G = 2 * F[1] - F[0];
            switch ((3 << i & n) >> i)
            {
                case 0: F = new Xint[2] { MTP(G, F[0]), MTP(G, F[1]) - 1 }; break;
                case 2: F = new Xint[2] { MTP(G, F[0]), MTP(G, F[1]) + 1 }; break;
                case 1: F = new Xint[2] { MTP(G, F[1]) - 1, MTP(G, F[1] + F[0]) - 1 }; break;
                case 3: F = new Xint[2] { MTP(G, F[1]) + 1, MTP(G, F[1] + F[0]) + 1 }; break;
            }
        }
        m = (n & 2) - 1; n &= 1; m *= n;
        return MTP(2 * F[1] - F[0], F[n]) + m;
    }

    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        Xint[] FL = { 0, 1 };
        for (int n = 2; n < 20000; n += 2)
        {
            Xint[] FH = Fibs(n);
            if (FH[0] != FL[0] + FL[1]) Console.WriteLine(n);
            if (FH[1] != FH[0] + FL[1]) Console.WriteLine(n);
            FL = FH;
        }
        Console.WriteLine("R E A D Y");
        Console.ReadLine();
    }

    // use the easy to copy faster versions!
    private static Xint MTP(Xint U, Xint V) { return U * V; }
    private static int fL2(int n) { int i = -1; while (n > 0) { i++; n /= 2; } return i; }
}

2011/07/05

Fibonacci numbers part 1

//  Fibonacci(1000000) in 105 mS (Athlon X4 640, XP, 2GB)

//  Two formulas:
//    
//  F[2n-1] = F[n-1]^2 + F[n]^2
//    
//  F[2n]   = (2*F[n-1] + F[n]) * F[n]  
//    
//  http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/fibFormula.html?n=55

//  Fibs(66) 
//  66 = 1000010
//       10000..  = 16 F0 = fibSmall(16) = 987                         = F(16)
//                     F1 = fibSmall(17) = 1597                        = F(17)
//                     F0 = (2*F1-F0)*F0 = (2*1597-987)*987  = 2178309 = F(32)
//                     F1 = F1*F1+F0*F0  = 1597*1597+987*987 = 3524578 = F(33)
//       .....1.       F0 = F1           = 3524578                     = F(33)
//                     F1 = F1+F0        = 3524578+2178309   = 5702887 = F(34)
//                     F0 = (2*F1-F0)*F0 = 27777890035288              = F(66)
//                     F1 = F1*F1+F0*F0  = 44945570212853              = F(67)
//       ......0       no change
//    
//  Fibs(66) = { Fib(66), Fib(67) } = { 27777890035288, 44945570212853 }

//                    +--------------------------------+
//                    | Fibonacci.exe in milliSeconds  |
//                    |--------------------------------|
//                    |    n    | F(n) | F(n+1)|Fibs(n)|
//                    |--------------------------------|
//                    |  156250 |   10 |    13 |    15 |
//                    |  312500 |   23 |    28 |    37 |
//                    |  625000 |   57 |    73 |    93 |
//                    | 1250000 |  148 |   198 |   255 |
//                    | 2500000 |  406 |   515 |   676 |
//                    | 5000000 | 1110 |  1420 |  1860 |
//                    |10000000 | 3010 |  3800 |  4930 | 
//                    |20000000 | 8080 | 10400 | 13500 |
//                    +--------------------------------+

//            !!! Project >> Add Reference >> System.Numerics !!!

using Xint = System.Numerics.BigInteger;
using System.Threading.Tasks;
using System.Diagnostics;
using System;
class program
{
    private static int fibSmall(int n)
    {
        int[] f = { 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 };
        return f[n];
    }
    public static Xint Fib(int n)
    {
        if (n < 34) return fibSmall(n);
        int i = fL2(n) - 4;
        int m = n >> i;
        Xint[] F = { fibSmall(m++), fibSmall(m) };
        for (--i; i > 0; i--)
        {
            F = new Xint[2] { MTP(2 * F[1] - F[0], F[0]), SQ(F[1]) + SQ(F[0]) };
            if (((1 << i) & n) != 0) F = new Xint[2] { F[1], F[1] + F[0] };
        }
        if ((n & 1) == 0) return MTP(2 * F[1] - F[0], F[0]);
        return SQ(F[1]) + SQ(F[0]);
    }
    public static Xint[] Fibs(int n)
    {
        if (n < 33) return new Xint[2] { fibSmall(n++), fibSmall(n) };
        int i = fL2(n) - 4;
        int m = n >> i;
        Xint[] F = { fibSmall(m++), fibSmall(m) };
        for (--i; i >= 0; i--)
        {
            F = new Xint[2] { MTP(2 * F[1] - F[0], F[0]), SQ(F[1]) + SQ(F[0]) };
            if (((1 << i) & n) != 0) F = new Xint[2] { F[1], F[1] + F[0] };
        }
        return F;
    }


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

    private static Xint MTP(Xint U, Xint V)
    {
        return MTP(U, V, Xint.Max(U.Sign * U, V.Sign * V).ToByteArray().Length << 3);
    }
    private static Xint MTP(Xint U, Xint V, int n)
    {
        if (n <= 3000) return U * V;
        if (n <= 6000) return TC2(U, V, n);
        if (n <= 10000) return TC3(U, V, n);
        if (n <= 40000) return TC4(U, V, n);
        return TC2P(U, V, n);
    }
    private static Xint MTPr(Xint U, Xint V, int n)
    {
        if (n <= 3000) return U * V;
        if (n <= 6000) return TC2(U, V, n);
        if (n <= 10000) return TC3(U, V, n);
        return TC4(U, V, n);
    }
    private static Xint TC2(Xint U1, Xint V1, int n)
    {
        n >>= 1;
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U1 & Mask; U1 >>= n;
        Xint V0 = V1 & Mask; V1 >>= n;
        Xint P0 = MTPr(U0, V0, n);
        Xint P2 = MTPr(U1, V1, n);
        return ((P2 << n) + (MTPr(U0 + U1, V0 + V1, n) - (P0 + P2)) << n) + P0;
    }
    private static Xint TC3(Xint U2, Xint V2, int n)
    {
        n = (int)((long)(n) * 0x55555556 >> 32); // n /= 3;
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U2 & Mask; U2 >>= n;
        Xint U1 = U2 & Mask; U2 >>= n;
        Xint V0 = V2 & Mask; V2 >>= n;
        Xint V1 = V2 & Mask; V2 >>= n;
        Xint W0 = MTPr(U0, V0, n);
        Xint W4 = MTPr(U2, V2, n);
        Xint P3 = MTPr((((U2 << 1) + U1) << 1) + U0, (((V2 << 1) + V1 << 1)) + V0, n);
        U2 += U0;
        V2 += V0;
        Xint P2 = MTPr(U2 - U1, V2 - V1, n);
        Xint P1 = MTPr(U2 + U1, V2 + V1, n);
        Xint W2 = (P1 + P2 >> 1) - (W0 + W4);
        Xint W3 = W0 - P1;
        W3 = ((W3 + P3 - P2 >> 1) + W3) / 3 - (W4 << 1);
        Xint W1 = P1 - (W4 + W3 + W2 + W0);
        return ((((W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0;
    }
    private static Xint TC4(Xint U3, Xint V3, int n)
    {
        n >>= 2;
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U3 & Mask; U3 >>= n;
        Xint U1 = U3 & Mask; U3 >>= n;
        Xint U2 = U3 & Mask; U3 >>= n;
        Xint V0 = V3 & Mask; V3 >>= n;
        Xint V1 = V3 & Mask; V3 >>= n;
        Xint V2 = V3 & Mask; V3 >>= n;

        Xint W0 = MTPr(U0, V0, n);                               //  0
        U0 += U2; U1 += U3;
        V0 += V2; V1 += V3;
        Xint P1 = MTPr(U0 + U1, V0 + V1, n);                     //  1
        Xint P2 = MTPr(U0 - U1, V0 - V1, n);                     // -1
        U0 += 3 * U2; U1 += 3 * U3;
        V0 += 3 * V2; V1 += 3 * V3;
        Xint P3 = MTPr(U0 + (U1 << 1), V0 + (V1 << 1), n);       //  2
        Xint P4 = MTPr(U0 - (U1 << 1), V0 - (V1 << 1), n);       // -2
        Xint P5 = MTPr(U0 + 12 * U2 + ((U1 + 12 * U3) << 2),
                       V0 + 12 * V2 + ((V1 + 12 * V3) << 2), n); //  4
        Xint W6 = MTPr(U3, V3, n);                               //  inf

        Xint W1 = P1 + P2;
        Xint W4 = (((((P3 + P4) >> 1) - (W1 << 1)) / 3 + W0) >> 2) - 5 * W6;
        Xint W2 = (W1 >> 1) - (W6 + W4 + W0);
        P1 = P1 - P2;
        P4 = P4 - P3;
        Xint W5 = ((P1 >> 1) + (5 * P4 + P5 - W0 >> 4) - ((((W6 << 4) + W4) << 4) + W2)) / 45;
        W1 = ((P4 >> 2) + (P1 << 1)) / 3 + (W5 << 2);
        Xint W3 = (P1 >> 1) - (W1 + W5);
        return ((((((W6 << n) + W5 << n) + W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0;
    }
    private static Xint TC2P(Xint A, Xint B, int n)
    {
        n >>= 1;
        Xint Mask = (Xint.One << n) - 1;
        Xint[] U = new Xint[3];
        U[0] = A & Mask; A >>= n; U[2] = A; U[1] = U[0] + A;
        Xint[] V = new Xint[3];
        V[0] = B & Mask; B >>= n; V[2] = B; V[1] = V[0] + B;
        Xint[] P = new Xint[3];
        Parallel.For(0, 3, (int i) => P[i] = MTPr(U[i], V[i], n));
        return ((P[2] << n) + P[1] - (P[0] + P[2]) << n) + P[0];
    }

    private static Xint SQ(Xint U)
    {
        return SQ(U, U.Sign * U.ToByteArray().Length << 3);
    }
    private static Xint SQ(Xint U, int n)
    {
        if (n <= 700) return U * U;
        if (n <= 3000) return Xint.Pow(U, 2);
        if (n <= 6000) return SQ2(U, n);
        if (n <= 10000) return SQ3(U, n);
        if (n <= 40000) return SQ4(U, n);
        return SQ2P(U, n);
    }
    private static Xint SQr(Xint U, int n)
    {
        if (n <= 3000) return Xint.Pow(U, 2);
        if (n <= 6000) return SQ2(U, n);
        if (n <= 10000) return SQ3(U, n);
        return SQ4(U, n);
    }
    private static Xint SQ2(Xint U1, int n)
    {
        n >>= 1;
        Xint U0 = U1 & ((Xint.One << n) - 1); U1 >>= n;
        Xint P0 = SQr(U0, n);
        Xint P2 = SQr(U1, n);
        return ((P2 << n) + (SQr(U0 + U1, n) - (P0 + P2)) << n) + P0;
    }
    private static Xint SQ3(Xint U2, int n)
    {
        n = (int)((long)(n) * 0x55555556 >> 32);
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U2 & Mask; U2 >>= n;
        Xint U1 = U2 & Mask; U2 >>= n;
        Xint W0 = SQr(U0, n);
        Xint W4 = SQr(U2, n);
        Xint P3 = SQr((((U2 << 1) + U1) << 1) + U0, n);
        U2 += U0;
        Xint P2 = SQr(U2 - U1, n);
        Xint P1 = SQr(U2 + U1, n);
        Xint W2 = (P1 + P2 >> 1) - (W0 + W4);
        Xint W3 = W0 - P1;
        W3 = ((W3 + P3 - P2 >> 1) + W3) / 3 - (W4 << 1);
        Xint W1 = P1 - (W4 + W3 + W2 + W0);
        return ((((W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0;
    }
    private static Xint SQ4(Xint U3, int n)
    {
        n >>= 2;
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U3 & Mask; U3 >>= n;
        Xint U1 = U3 & Mask; U3 >>= n;
        Xint U2 = U3 & Mask; U3 >>= n;
        Xint W0 = SQr(U0, n);                                   //  0
        U0 += U2;
        U1 += U3;
        Xint P1 = SQr(U0 + U1, n);                              //  1
        Xint P2 = SQr(U0 - U1, n);                              // -1
        U0 += 3 * U2;
        U1 += 3 * U3;
        Xint P3 = SQr(U0 + (U1 << 1), n);                       //  2
        Xint P4 = SQr(U0 - (U1 << 1), n);                       // -2
        Xint P5 = SQr(U0 + 12 * U2 + ((U1 + 12 * U3) << 2), n); //  4
        Xint W6 = SQr(U3, n);                                   //  inf
        Xint W1 = P1 + P2;
        Xint W4 = (((((P3 + P4) >> 1) - (W1 << 1)) / 3 + W0) >> 2) - 5 * W6;
        Xint W2 = (W1 >> 1) - (W6 + W4 + W0);
        P1 = P1 - P2;
        P4 = P4 - P3;
        Xint W5 = ((P1 >> 1) + (5 * P4 + P5 - W0 >> 4) - ((((W6 << 4) + W4) << 4) + W2)) / 45;
        W1 = ((P4 >> 2) + (P1 << 1)) / 3 + (W5 << 2);
        Xint W3 = (P1 >> 1) - (W1 + W5);
        return ((((((W6 << n) + W5 << n) + W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0;
    }
    private static Xint SQ2P(Xint A, int n)
    {
        n >>= 1;
        Xint[] U = new Xint[3];
        U[0] = A & ((Xint.One << n) - 1); A >>= n; U[2] = A; U[1] = U[0] + A;
        Xint[] P = new Xint[3];
        Parallel.For(0, 3, (int i) => P[i] = SQr(U[i], n));
        return ((P[2] << n) + P[1] - (P[0] + P[2]) << n) + P[0];
    }

    public static int bL(Xint U)
    {
        byte[] bytes = (U.Sign * U).ToByteArray();
        int i = bytes.Length - 1;
        return (i << 3) + bitLengthMostSignificantByte(bytes[i]);
    }
    private static int bitLengthMostSignificantByte(byte b)
    {
        return b < 08 ? b < 02 ? b < 01 ? 0 : 1 :
                                 b < 04 ? 2 : 3 :
                        b < 32 ? b < 16 ? 4 : 5 :
                                 b < 64 ? 6 : 7;
    }

    private static int fL2(int i)
    {
        return
        i < 1 << 15 ? i < 1 << 07 ? i < 1 << 03 ? i < 1 << 01 ? i < 1 << 00 ? -1 : 00 :
                                                                i < 1 << 02 ? 01 : 02 :
                                                  i < 1 << 05 ? i < 1 << 04 ? 03 : 04 :
                                                                i < 1 << 06 ? 05 : 06 :
                                    i < 1 << 11 ? i < 1 << 09 ? i < 1 << 08 ? 07 : 08 :
                                                                i < 1 << 10 ? 09 : 10 :
                                                  i < 1 << 13 ? i < 1 << 12 ? 11 : 12 :
                                                                i < 1 << 14 ? 13 : 14 :
                      i < 1 << 23 ? i < 1 << 19 ? i < 1 << 17 ? i < 1 << 16 ? 15 : 16 :
                                                                i < 1 << 18 ? 17 : 18 :
                                                  i < 1 << 21 ? i < 1 << 20 ? 19 : 20 :
                                                                i < 1 << 22 ? 21 : 22 :
                                    i < 1 << 27 ? i < 1 << 25 ? i < 1 << 24 ? 23 : 24 :
                                                                i < 1 << 26 ? 25 : 26 :
                                                  i < 1 << 29 ? i < 1 << 28 ? 27 : 28 :
                                                                i < 1 << 30 ? 29 : 30;
    }
}

2011/01/28

Eratosthenes, part 3

// The BitArray became 33% smaller in part 2,
// the speed ~20% faster. With a few changes
// this one is ~33% faster. In the IDE primes
// upto 10^9 are found in 15 s. (12.5 in the 
// release version, (Athlon X4 640, XP, 2GB))

using System;
using System.Collections.Generic;
using System.Collections;
using System.Diagnostics;
class Program
{
    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        Console.WriteLine(int.MaxValue);
        int n = 1;
        while (n <= 1000000000)     
        {
            sw.Restart();
            int[] primes = getPrimes(n);
            sw.Stop();
            Console.WriteLine("n: " + n + "  primes: " + primes.Length +
                              " in " + sw.ElapsedMilliseconds + " ms");
            n *= 10;
        }
        Console.ReadLine();
    }
    private static int[] getPrimes(int n)
    {
        if (n < 3)
        {
            if (n == 2) { int[] prime2 = { 2 }; return prime2; }
            else { int[] no_prime = { }; return no_prime; }
        }
        BitArray composite = new BitArray(n / 3);
        int len = composite.Length;
        int d1 = 8, d2 = 0, p1 = 3, p2 = 1, s = 7;
        int i = 0;
        int inc = 10;
        while (s < len)
        {
            if (!composite[i++])
            {
                int k = s, m = s + p1;
                for (; m < len; k += inc, m += inc)
                {
                    composite[k] = true;
                    composite[m] = true;
                }
                for (; k < len; k += inc)
                    composite[k] = true;
            }
            d2 += 8; s += d2; p2 += 8; inc += 4;
            if (!composite[i++] && s < len)
            {
                int k = s, m = s + p2;
                for (; m < len; k += inc, m += inc)
                {
                    composite[k] = true;
                    composite[m] = true;
                }
                for (; k < len; k += inc)
                    composite[k] = true;
            }
            d1 += 16; s += d1; p1 += 4; inc += 8;
        }
        List<int> primes = new List<int>();
        double log_n = Math.Log(n);
        primes.Capacity = (int)((n / log_n) * (1 + 1.2762 / log_n));
        primes.Add(2); primes.Add(3);
        int p = 5;
        i = 0;
        while (p <= n)
        {
            if (!composite[i++]) primes.Add(p);
            p += 2;
            if (p <= n && !composite[i++]) primes.Add(p);
            p += 4;
        }
        return primes.ToArray();
    }
}

2011/01/16

Eratosthenes, part 2

// Primes can be written as 6*N-1 and 6*N+1,
// 5,7,11,13,17,19,23,25?29,.....
// Increase from prime to next prime:
// 2,4,2,4,2,4,2,4.....
// A table of indexes in the BitArray,
// and corresponding primes,
   
//  n prime
//   
//  0 5
//  1 7
//  2 11
//  3 13
//  4 17
//  5 19
//  6 23
//  7 25
//  8 29
//  9 31
// 10 35
// 11 37
// .. ..

// 25,35,55,65..have to be crossed off (in the BitArray).
// I had seen it somewhere, but didn't understand what was going on.
// Somewhere is:
   
// http://code.msdn.microsoft.com/factorialbenchmark/Release/ProjectReleases.aspx?ReleaseId=4516
   
// Result: it's ~20% faster, the BitArray is 33% smaller.

using System;
using System.Collections.Generic;
using System.Collections;
using System.Diagnostics;
class Program
{
    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        int n = 1;
        while (n <= 1000000000)             // primes.exe n=10^9 18 s. 
        {
            sw.Restart();
            int[] primes = getPrimes(n);
            sw.Stop();
            Console.WriteLine("n: " + n + "  primes: " + primes.Length +
                              " in " + sw.ElapsedMilliseconds + " ms");
            n *= 10;
        }
        Console.ReadLine();
    }
    private static int[] getPrimes(int n)
    {
        if (n < 3)
        {
            if (n == 2) { int[] prime2 = { 2 }; return prime2; }
            else { int[] no_prime = { }; return no_prime; }
        }
        BitArray composite = new BitArray(n / 3);
        int d1 = 8, d2 = 8, p1 = 3, p2 = 7, s = 7, s2 = 3;
        int i = 0;
        int len = composite.Length;
        bool toggle = false;
        while (s < len)
        {
            if (!composite[i++])
            {
                int inc = p1 + p2;
                for (int k = s; k < len; k += inc)
                    composite[k] = true;
                for (int k = s + s2; k < len; k += inc)
                    composite[k] = true;
            }
            p1 += 2;
            if (toggle = !toggle)
            { s += d2; d1 += 16; p2 += 2; s2 = p2; }
            else
            { s += d1; d2 += 08; p2 += 6; s2 = p1; }
        }
        List<int> primes = new List<int>();
        double log_n = Math.Log(n);
        primes.Capacity = (int)((n / log_n) * (1 + 1.2762 / log_n));
        primes.Add(2); primes.Add(3);
        int p = 5;
        toggle = false;
        i = 0;
        while (p <= n)
        {
            if (!composite[i++]) primes.Add(p);
            p += (toggle = !toggle) ? 2 : 4;
        }
        return primes.ToArray();
    }
}

2010/12/21

Primes, Sieve of Eratosthenes

// See Knuth's TAOCP section 4.5.4 exercise 8.
// The number of primes less then or equal to n: 
// http://primes.utm.edu/howmany.shtml 
// pi(n) <= n/ln(n)*(1+1.27262/ln(n))
// A fast sieve:
// http://www.primzahlen.de/files/referent/kw/sieb.htm
// (Pentium III 550 Mhz, C/C++, n=10^9, 2.533 s!)

using System;
using System.Collections;
using System.Collections.Generic;
using System.Diagnostics;
class Program
{
    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        int n = 1;
        while (n <= 1000000000)             // primes.exe n=10^9 23 s. 
        {
            sw.Restart();
            int[] primes = getPrimes(n);
            sw.Stop();
            Console.WriteLine("n: " + n + "  primes: " + primes.Length +
                              " in " + sw.ElapsedMilliseconds + " ms");
            n *= 10;
        }
        Console.ReadLine();
    }
    private static int[] getPrimes(int n)
    {
        List<int> primes = new List<int>();
        if (n < 2) return primes.ToArray();
        double log_n = Math.Log(n);
        primes.Capacity = (int)((n / log_n) * (1 + 1.2762 / log_n));
        n = ((n & 1) + n) >> 1;
        BitArray x = new BitArray(n + 1);
        int i = 4, p = 3, q = 4;
        for (; q < n; q += p >> 1 << 2)
        {
            if (!x[q])
            {
                x[q] = true;
                for (i = q + p; i < n; i += p) x[i] = true;
            }
            p += 2;
        }
        primes.Add(2);
        for (i = 1; i < n; i++) if (!x[i]) primes.Add(i * 2 | 1);
        return primes.ToArray();
    }
}

2010/12/19

Binomial Coefficient, n over k, Choose(n,k)

// Another abbreviation

// BNM(5,2) = 5! / 2! / (5-2)! = 10  binomial coefficient

// A sieve of Eratosthenes is used to get the primes,
// GetBnmPrimes(30,15) = {3,3,5,17,19,23,29} , 7 increasing numbers.
// floorLog2(7) = fL2(7) = 2
// 1 << fL2(7) = 4
// BnmSpecialProduct changes the array into
// {17,3*19,3*23,5*29} = {17,57,69,145} , 4 increasing numbers.
// BnmProduct changes it into 
// {17*145,57*69} = {2465,3933}
// BnmProduct finally changes it into
// {2465*3933} = {9694845}
// A small prime, 2 , still has to be handled,
// its power factor is getBnmPower2(30,15) = 4
// 9694845 * 2^4 = 9694845 << 4 = 155117520

// After a few iterations of BnmProduct, 
// the multiplicands are of equal, or almost equal, bitlength,
// a prerequisite to use Karatsuba etc. effectively.
//
using System;
using System.Collections;
using System.Collections.Generic;
using System.Diagnostics;
using System.Threading.Tasks;
using Xint = System.Numerics.BigInteger;
class Program
{
    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        sw.Restart();
        BNM(6400000, 2133333);
        sw.Stop();
        Console.WriteLine(sw.ElapsedMilliseconds);       // BNM.exe 3880 ms
        sw.Restart();
        Console.ReadLine();
    }
    public static Xint BNM(int n, int k)
    {
        if (k > n) return 0;
        if ((k == 0) | (n == k)) return 1;
        if ((k == 1) | (n == k + 1)) return n;
        Xint[] P = GetBnmPrimes(n, k);
        int i = P.Length;
        int j = 1 << fL2(i);
        if (i != j)
        {
            P = BnmSpecialProduct(P, i, j);
            i = j;
        }
        while (i > 1)
        {
            P = BnmProduct(P, i);
            i >>= 1;
        }
        return P[0] << getBnmPower2(n, k);
    }
    private static Xint[] GetBnmPrimes(int n, int k)
    {
        int m = ((n & 1) + n) >> 1;                         // Eratosthenes
        BitArray x = new BitArray(m + 1);
        int i = 4, p = 3, q = 4;
        for (; q < m; q += p >> 1 << 2)
        {
            if (!x[q])
            {
                x[q] = true;
                for (i = q + p; i < m; i += p) x[i] = true;
            }
            p += 2;
        }
        List<Xint> Primes = new List<Xint>();
        int j;                        // prime power factorization binomial 
        for (i = 1; i < m; i++) if (!x[i])
            {
                p = i * 2 | 1;
                q = n / p; j = q; while (q >= p) { q /= p; j += q; }
                q = k / p; j -= q; while (q >= p) { q /= p; j -= q; }
                q = (n - k) / p; j -= q; while (q >= p) { q /= p; j -= q; }
                for (q = 0; q < j; q++) Primes.Add(p);
            }
        return Primes.ToArray();
    }
    private static int getBnmPower2(int n, int k)
    {
        int m = n - k;
        int i = (n >>= 1);
        while (n > 1) i += (n >>= 1);
        while (m > 1) i -= (m >>= 1);
        while (k > 1) i -= (k >>= 1);
        return i;
    }
    private static Xint[] BnmSpecialProduct(Xint[] P, int i, int j)
    {
        Xint[] Q = new Xint[j];
        int m = i - j;
        i = j;
        j -= m;
        int n = 0;
        for (; n < j; n++)
        {
            Q[n] = P[m];
            m++;
        }
        j = 0;
        for (; n < i; n++)
        {
            Q[n] = P[m] * P[j];
            j++;
            m++;
        }
        return Q;
    }
    private static Xint[] BnmProduct(Xint[] P, int i)
    {
        int k = i >> 1;
        Xint[] Q = new Xint[k];
        for (int j = 0; j < k; j++)
        {
            i--;
            Q[j] = MTP(P[j], P[i]);
        }
        return Q;
    }

    private static Xint MTP(Xint U, Xint V)
    {
        return MTP(U, V, Xint.Max(U.Sign * U, V.Sign * V).ToByteArray().Length << 3);
    }
    private static Xint MTP(Xint U, Xint V, int n)
    {
        if (n <= 3000) return U * V;
        if (n <= 6000) return TC2(U, V, n);
        if (n <= 10000) return TC3(U, V, n);
        if (n <= 40000) return TC4(U, V, n);
        return TC2P(U, V, n);
    }
    private static Xint MTPr(Xint U, Xint V, int n)
    {
        if (n <= 3000) return U * V;
        if (n <= 6000) return TC2(U, V, n);
        if (n <= 10000) return TC3(U, V, n);
        return TC4(U, V, n);
    }
    private static Xint TC2(Xint U1, Xint V1, int n)
    {
        n >>= 1;
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U1 & Mask; U1 >>= n;
        Xint V0 = V1 & Mask; V1 >>= n;
        Xint P0 = MTPr(U0, V0, n);
        Xint P2 = MTPr(U1, V1, n);
        return ((P2 << n) + (MTPr(U0 + U1, V0 + V1, n) - (P0 + P2)) << n) + P0;
    }
    private static Xint TC3(Xint U2, Xint V2, int n)
    {
        n = (int)((long)(n) * 0x55555556 >> 32); // n /= 3;
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U2 & Mask; U2 >>= n;
        Xint U1 = U2 & Mask; U2 >>= n;
        Xint V0 = V2 & Mask; V2 >>= n;
        Xint V1 = V2 & Mask; V2 >>= n;
        Xint W0 = MTPr(U0, V0, n);
        Xint W4 = MTPr(U2, V2, n);
        Xint P3 = MTPr((((U2 << 1) + U1) << 1) + U0, (((V2 << 1) + V1 << 1)) + V0, n);
        U2 += U0;
        V2 += V0;
        Xint P2 = MTPr(U2 - U1, V2 - V1, n);
        Xint P1 = MTPr(U2 + U1, V2 + V1, n);
        Xint W2 = (P1 + P2 >> 1) - (W0 + W4);
        Xint W3 = W0 - P1;
        W3 = ((W3 + P3 - P2 >> 1) + W3) / 3 - (W4 << 1);
        Xint W1 = P1 - (W4 + W3 + W2 + W0);
        return ((((W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0;
    }
    private static Xint TC4(Xint U3, Xint V3, int n)
    {
        n >>= 2;
        Xint Mask = (Xint.One << n) - 1;
        Xint U0 = U3 & Mask; U3 >>= n;
        Xint U1 = U3 & Mask; U3 >>= n;
        Xint U2 = U3 & Mask; U3 >>= n;
        Xint V0 = V3 & Mask; V3 >>= n;
        Xint V1 = V3 & Mask; V3 >>= n;
        Xint V2 = V3 & Mask; V3 >>= n;

        Xint W0 = MTPr(U0, V0, n);                               //  0
        U0 += U2; U1 += U3;
        V0 += V2; V1 += V3;
        Xint P1 = MTPr(U0 + U1, V0 + V1, n);                     //  1
        Xint P2 = MTPr(U0 - U1, V0 - V1, n);                     // -1
        U0 += 3 * U2; U1 += 3 * U3;
        V0 += 3 * V2; V1 += 3 * V3;
        Xint P3 = MTPr(U0 + (U1 << 1), V0 + (V1 << 1), n);       //  2
        Xint P4 = MTPr(U0 - (U1 << 1), V0 - (V1 << 1), n);       // -2
        Xint P5 = MTPr(U0 + 12 * U2 + ((U1 + 12 * U3) << 2),
                       V0 + 12 * V2 + ((V1 + 12 * V3) << 2), n); //  4
        Xint W6 = MTPr(U3, V3, n);                               //  inf

        Xint W1 = P1 + P2;
        Xint W4 = (((((P3 + P4) >> 1) - (W1 << 1)) / 3 + W0) >> 2) - 5 * W6;
        Xint W2 = (W1 >> 1) - (W6 + W4 + W0);
        P1 = P1 - P2;
        P4 = P4 - P3;
        Xint W5 = ((P1 >> 1) + (5 * P4 + P5 - W0 >> 4) - ((((W6 << 4) + W4) << 4) + W2)) / 45;
        W1 = ((P4 >> 2) + (P1 << 1)) / 3 + (W5 << 2);
        Xint W3 = (P1 >> 1) - (W1 + W5);
        return ((((((W6 << n) + W5 << n) + W4 << n) + W3 << n) + W2 << n) + W1 << n) + W0;
    }
    private static Xint TC2P(Xint A, Xint B, int n)
    {
        n >>= 1;
        Xint Mask = (Xint.One << n) - 1;
        Xint[] U = new Xint[3];
        U[0] = A & Mask; A >>= n; U[2] = A; U[1] = U[0] + A;
        Xint[] V = new Xint[3];
        V[0] = B & Mask; B >>= n; V[2] = B; V[1] = V[0] + B;
        Xint[] P = new Xint[3];
        Parallel.For(0, 3, (int i) => P[i] = MTPr(U[i], V[i], n));
        return ((P[2] << n) + P[1] - (P[0] + P[2]) << n) + P[0];
    }
    private static int fL2(int i)
    {
        return
        i < 1 << 15 ? i < 1 << 07 ? i < 1 << 03 ? i < 1 << 01 ? i < 1 << 00 ? -1 : 00 :
                                                                i < 1 << 02 ? 01 : 02 :
                                                  i < 1 << 05 ? i < 1 << 04 ? 03 : 04 :
                                                                i < 1 << 06 ? 05 : 06 :
                                    i < 1 << 11 ? i < 1 << 09 ? i < 1 << 08 ? 07 : 08 :
                                                                i < 1 << 10 ? 09 : 10 :
                                                  i < 1 << 13 ? i < 1 << 12 ? 11 : 12 :
                                                                i < 1 << 14 ? 13 : 14 :
                      i < 1 << 23 ? i < 1 << 19 ? i < 1 << 17 ? i < 1 << 16 ? 15 : 16 :
                                                                i < 1 << 18 ? 17 : 18 :
                                                  i < 1 << 21 ? i < 1 << 20 ? 19 : 20 :
                                                                i < 1 << 22 ? 21 : 22 :
                                    i < 1 << 27 ? i < 1 << 25 ? i < 1 << 24 ? 23 : 24 :
                                                                i < 1 << 26 ? 25 : 26 :
                                                  i < 1 << 29 ? i < 1 << 28 ? 27 : 28 :
                                                                i < 1 << 30 ? 29 : 30;
    }
    private static int bL(Xint U)
    {
        byte[] bytes = (U.Sign * U).ToByteArray();
        int i = bytes.Length - 1;
        return i << 3 | bitLengthMostSignificantByte(bytes[i]);
    }
    private static int bitLengthMostSignificantByte(byte b)
    {
        return b < 08 ? b < 02 ? b < 01 ? 0 : 1 :
                                 b < 04 ? 2 : 3 :
                        b < 32 ? b < 16 ? 4 : 5 :
                                 b < 64 ? 6 : 7;
    }
}

Pagina's

2011/07/09

2011/07/05

2011/01/28

2011/01/16

2010/12/21

2010/12/19