Pagina's

2016/04/02

Segmented sieve of Eratosthenes


/* 50,000,000 primes in 1 second.                  70,000,000 ? scroll down
                                                  125,000,000 ? Odd prime sieve
            time (in seconds)     i7@3G6Hz
      n       x      p   isP?     x = composites
    1e9    0.68   0.93   0.83     p = primes <= n < 2^32 ~ 4e9
    2e9    1.39   1.84   1.65     isP? = isPrime[i], i odd
    4e9    2.85   3.72   3.38          = (isP[i >> 6] & 1 << (i >> 1)) != 0
    ~0u    3.06   4.00   3.63
*/
using System;
using System.Threading.Tasks;
using sw = System.Diagnostics.Stopwatch;
namespace primes      
{
    class Program                
    {
        static sw sw = new sw();
        static void Main()
        {
            test((uint)1e9);
            Console.Read();
        }

        static void test(uint n)
        {
            uint pi_n; uint[] p;
            sw.Start();
            p = getPrimes(n);
            sw.Stop();
            pi_n = p[p.Length - 1];
            Console.WriteLine("pi(" + n + ") = " + pi_n);
            if (pi_n > 0)
                Console.Write("largest prime = " + p[pi_n - 1]);
        }

        static uint[] getPrimes(uint n)
        {
            int c, y, i, j; uint u; double logN; int[] x; uint[] p;
            if (n < 3) return n < 2 ? new uint[] { 0 } : new uint[] { 2, 1 };
            logN = Math.Log(n);
            c = (int)(n / logN * (1 + 1.2762 / logN));
            p = new uint[c + 1]; p[0] = 2; p[1] = 3;
            y = (int)(n / 3); x = new int[(y >> 5) + 1];
            if (n < 200000000)
                mark0(y, x);
            else
                mark1(y, x);
            Console.WriteLine(sw.Elapsed + " x");
            y = y <= 2 ? 0 : y - 2;
            j = y < 32 ? 2 : paraPrimes(y >> 5, x, p);
            i = y >> 5 << 5; u = 5 + 3 * (uint)i;
            while (i < y)
            {
                if ((x[i >> 5] & 1 << i++) == 0) p[j++] = u; u += 2;
                if ((x[i >> 5] & 1 << i++) == 0) p[j++] = u; u += 4;
            }
            if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) p[j++] = u; u += 2;
            if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) p[j++] = u;
            p[c] = (uint)j;
            Console.WriteLine(sw.Elapsed + " p");
            return p;
        }

        static void mark0(int y, int[] z)
        {
            int u, v, w, x, s, i, j;
            u = 8; v = 0; w = 3; x = 1; s = 7; i = 0; j = 10;
            while (s < y)
            {
                if ((z[i >> 5] & 1 << i++) == 0)
                {
                    int a = s, b = s + j, c = s + w, d = c + j, k = j << 1;
                    for (; d < y; a += k, b += k, c += k, d += k)
                    {
                        z[a >> 5] |= 1 << a; z[c >> 5] |= 1 << c;
                        z[b >> 5] |= 1 << b; z[d >> 5] |= 1 << d;
                    }
                    if (a < y)
                    {
                        z[a >> 5] |= 1 << a;
                        if (c < y) z[c >> 5] |= 1 << c;
                        if (b < y) z[b >> 5] |= 1 << b;
                    }
                }
                v += 8; s += v; x += 8; j += 4;
                if ((z[i >> 5] & 1 << i++) == 0)
                {
                    int a = s, b = s + j, c = s + x, d = c + j, k = j << 1;
                    for (; d < y; a += k, b += k, c += k, d += k)
                    {
                        z[a >> 5] |= 1 << a; z[c >> 5] |= 1 << c;
                        z[b >> 5] |= 1 << b; z[d >> 5] |= 1 << d;
                    }
                    if (a < y)
                    {
                        z[a >> 5] |= 1 << a;
                        if (c < y) z[c >> 5] |= 1 << c;
                        if (b < y) z[b >> 5] |= 1 << b;
                    }
                }
                u += 16; s += u; w += 4; j += 8;
            }
        }

        static void mark1(int y, int[] z)
        {
            bool bl, bq; int u, v, w, x, s, i, j, m, n; queue q0, q1;
            bq = true; q0 = new queue(); q1 = new queue();
            bl = false; m = 1024; u = 8; v = 0; w = 3; x = 1; s = 7; i = 0; j = 10;
            for (n = m; n < y; n += m)
            {
                if (bl)
                {
                    bl = false;
                    if (bq)
                    {
                        bq = false;
                        while (q0.count > 0)
                        {
                            int a = q0.dequeue(), k = q0.dequeue();
                            for (; a < n; a += k) z[a >> 5] |= 1 << a;
                            if (a < y) q1.enqueue(a, k);
                        }
                        q0.idxIn = q0.idxOut = 0;
                    }
                    else
                    {
                        bq = true;
                        while (q1.count > 0)
                        {
                            int a = q1.dequeue(), k = q1.dequeue();
                            for (; a < n; a += k) z[a >> 5] |= 1 << a;
                            if (a < y) q0.enqueue(a, k);
                        }
                        q1.idxIn = q1.idxOut = 0;
                    }
                }
                while (s < n)
                {
                    if ((z[i >> 5] & 1 << i++) == 0)
                    {
                        bl = true;
                        int a = s, b = s + w;
                        for (; b < n; a += j, b += j) { z[a >> 5] |= 1 << a; z[b >> 5] |= 1 << b; }
                        for (; a < n; a += j) z[a >> 5] |= 1 << a;
                        if (bq)
                        { if (a < y) q0.enqueue(a, j); if (b < y) q0.enqueue(b, j); }
                        else
                        { if (a < y) q1.enqueue(a, j); if (b < y) q1.enqueue(b, j); }
                    }
                    v += 8; s += v; x += 8; j += 4;
                    if ((z[i >> 5] & 1 << i++) == 0)
                    {
                        bl = true;
                        int a = s, b = s + x;
                        for (; b < n; a += j, b += j) { z[a >> 5] |= 1 << a; z[b >> 5] |= 1 << b; }
                        for (; a < n; a += j) z[a >> 5] |= 1 << a;
                        if (bq)
                        { if (a < y) q0.enqueue(a, j); if (b < y) q0.enqueue(b, j); }
                        else
                        { if (a < y) q1.enqueue(a, j); if (b < y) q1.enqueue(b, j); }
                    }
                    u += 16; s += u; w += 4; j += 8;
                }
            }
            while (q0.count > 0)
                for (int a = q0.dequeue(), k = q0.dequeue(); a < y; a += k) z[a >> 5] |= 1 << a;
            while (q1.count > 0)
                for (int a = q1.dequeue(), k = q1.dequeue(); a < y; a += k) z[a >> 5] |= 1 << a;
        }

        static int paraPrimes(int imax, int[] x, uint[] p)
        {
            int j0 = 0, j1 = 0;
            Parallel.For(0, 2, k =>
            {
                int i, j, y; uint u;
                if (k == 0) { i = 0; j = 2; u = 5; }
                else { if (imax < 2) return; i = 1; j = 25; u = 101; }
                for (; ; )
                {
                    y = x[i];
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4;
                    i += 2; if (i >= imax) { if (k == 0) j0 = j; else j1 = j; break; }
                    u += 96; j += bitsSet(~x[i - 1]);
                }
            });
            return j0 > j1 ? j0 : j1;
        }

        static int bitsSet(int i)
        {
            i -= i >> 1 & 0x55555555;
            i = (i & 0x33333333) + (i >> 2 & 0x33333333);
            return (i + (i >> 4) & 0xf0f0f0f) * 0x1010101 >> 24;
        }
    }

    class queue
    {
        public int idxIn, idxOut, count;
        int[] q = new int[26200];
        public void enqueue(int i, int j) { count += 2; q[idxIn++] = i; q[idxIn++] = j; }
        public int dequeue() { count--; return q[idxOut++]; }
    }
}


using System; using System.Threading.Tasks; using sw = System.Diagnostics.Stopwatch; namespace isPrime { class Program { static sw sw = new sw(); static void Main() { test((uint)1e9); Console.Read(); } static void test(uint n) { int[] isP; sw.Start(); isP = buildIsPrime(n); sw.Stop(); // n check int check = 0; // 1e9 1DF99124 for (int i = isP.Length - 1; i >= 0; i--) // 2e9 E9D92C9E check ^= isP[i]; // 4e9 DB1F9720 Console.Write("check: {0:X8}", check); // ~0u D696B791 } static int[] buildIsPrime(uint n) { int y, i, u; int[] x, isP; if (n < 64) return new int[] { 0x64b4cb6e }; isP = new int[(n >> 6) + 1]; isP[0] = 2; y = (int)(n / 3); x = new int[(y >> 5) + 1]; if (n < 200000000) mark0(y, x); else mark1(y, x); Console.WriteLine(sw.Elapsed + " x"); y -= 2; n >>= 1; paraIsPrime(y >> 5, x, isP); u = (y >> 5 << 4) * 3 + 2; i = y >> 5 << 5; while (i < y) { if ((x[i >> 5] & 1 << i++) == 0) isP[u >> 5] |= 1 << u; u += 1; if ((x[i >> 5] & 1 << i++) == 0) isP[u >> 5] |= 1 << u; u += 2; } if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) isP[u >> 5] |= 1 << u; u += 1; if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) isP[u >> 5] |= 1 << u; Console.WriteLine(sw.Elapsed + " p"); return isP; } static void mark0(int y, int[] z) { int u, v, w, x, s, i, j; u = 8; v = 0; w = 3; x = 1; s = 7; i = 0; j = 10; while (s < y) { if ((z[i >> 5] & 1 << i++) == 0) { int a = s, b = s + j, c = s + w, d = c + j, k = j << 1; for (; d < y; a += k, b += k, c += k, d += k) { z[a >> 5] |= 1 << a; z[c >> 5] |= 1 << c; z[b >> 5] |= 1 << b; z[d >> 5] |= 1 << d; } if (a < y) { z[a >> 5] |= 1 << a; if (c < y) z[c >> 5] |= 1 << c; if (b < y) z[b >> 5] |= 1 << b; } } v += 8; s += v; x += 8; j += 4; if ((z[i >> 5] & 1 << i++) == 0) { int a = s, b = s + j, c = s + x, d = c + j, k = j << 1; for (; d < y; a += k, b += k, c += k, d += k) { z[a >> 5] |= 1 << a; z[c >> 5] |= 1 << c; z[b >> 5] |= 1 << b; z[d >> 5] |= 1 << d; } if (a < y) { z[a >> 5] |= 1 << a; if (c < y) z[c >> 5] |= 1 << c; if (b < y) z[b >> 5] |= 1 << b; } } u += 16; s += u; w += 4; j += 8; } } static void mark1(int y, int[] z) { bool bl, bq; int u, v, w, x, s, i, j, m, n; queue q0, q1; bq = true; q0 = new queue(); q1 = new queue(); bl = false; m = 1024; u = 8; v = 0; w = 3; x = 1; s = 7; i = 0; j = 10; for (n = m; n < y; n += m) { if (bl) { bl = false; if (bq) { bq = false; while (q0.count > 0) { int a = q0.dequeue(), k = q0.dequeue(); for (; a < n; a += k) z[a >> 5] |= 1 << a; if (a < y) q1.enqueue(a, k); } q0.idxIn = q0.idxOut = 0; } else { bq = true; while (q1.count > 0) { int a = q1.dequeue(), k = q1.dequeue(); for (; a < n; a += k) z[a >> 5] |= 1 << a; if (a < y) q0.enqueue(a, k); } q1.idxIn = q1.idxOut = 0; } } while (s < n) { if ((z[i >> 5] & 1 << i++) == 0) { bl = true; int a = s, b = s + w; for (; b < n; a += j, b += j) { z[a >> 5] |= 1 << a; z[b >> 5] |= 1 << b; } for (; a < n; a += j) z[a >> 5] |= 1 << a; if (bq) { if (a < y) q0.enqueue(a, j); if (b < y) q0.enqueue(b, j); } else { if (a < y) q1.enqueue(a, j); if (b < y) q1.enqueue(b, j); } } v += 8; s += v; x += 8; j += 4; if ((z[i >> 5] & 1 << i++) == 0) { bl = true; int a = s, b = s + x; for (; b < n; a += j, b += j) { z[a >> 5] |= 1 << a; z[b >> 5] |= 1 << b; } for (; a < n; a += j) z[a >> 5] |= 1 << a; if (bq) { if (a < y) q0.enqueue(a, j); if (b < y) q0.enqueue(b, j); } else { if (a < y) q1.enqueue(a, j); if (b < y) q1.enqueue(b, j); } } u += 16; s += u; w += 4; j += 8; } } while (q0.count > 0) for (int a = q0.dequeue(), k = q0.dequeue(); a < y; a += k) z[a >> 5] |= 1 << a; while (q1.count > 0) for (int a = q1.dequeue(), k = q1.dequeue(); a < y; a += k) z[a >> 5] |= 1 << a; } static void paraIsPrime(int imax, int[] x, int[] z) // z = isP { int pc = Environment.ProcessorCount; for (int n = 0; n < 2; n++) { Parallel.For(0, pc, k => { int m, i, j, u, v, y; if ((k & 1) == n) { j = pc * 2; v = (j - 1) * 48; for (m = 0; m < 2; m++) { i = m + k * 2; u = 2 + i * 48; for (; i < imax; i += j, u += v) { #region y = x[i]; if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2; if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4; if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2; if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4; if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2; if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4; if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2; if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4; if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2; if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4; if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2; if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4; if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2; if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4; if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2; if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1; if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; #endregion } } } }); } } } class queue { public int idxIn, idxOut, count; int[] q = new int[26200]; public void enqueue(int i, int j) { count += 2; q[idxIn++] = i; q[idxIn++] = j; } public int dequeue() { count--; return q[idxOut++]; } } }




/* Slightly faster, pre-sieving 5, 7, 11 & 13.
   New times were measured with seperate versions
   for "isPrime" and "getPrimes".
 
           Old times(seconds)    New times(seconds)  
      n       x      p   isP?       x      p   isP?  
    1e9    0.68   0.93   0.83    0.50   0.66   0.65  
    2e9    1.39   1.84   1.65    1.02   1.31   1.29  
    4e9    2.85   3.72   3.38    2.13   2.65   2.62  
    ~0u    3.06   4.00   3.63    2.28   2.87   2.84
*/
using System;
using System.Threading.Tasks;
using sw = System.Diagnostics.Stopwatch;
namespace primes
{
    class Program
    {
        static void Main()
        {
            uint n = (uint)1e9;
            testIsPrime(n);
            Console.WriteLine();
            testGetPrimes(n);
            Console.Read();
        }

        static void testIsPrime(uint n)
        {
            int[] isP = new soE().buildIsPrime(n);
            int check = 0;                                   // 1e9 1DF99124
            for (int i = isP.Length - 1; i >= 0; i--)        // 2e9 E9D92C9E
                check ^= isP[i];                             // 4e9 DB1F9720
            Console.WriteLine("check: {0:X8}", check);       // ~0u D696B791
        }

        static void testGetPrimes(uint n)
        {
            uint pi_n; uint[] p;
            p = new soE().getPrimes(n);
            pi_n = p[p.Length - 1];
            Console.WriteLine("pi(" + n + ") = " + pi_n);
            if (pi_n > 0)
                Console.Write("largest prime = " + p[pi_n - 1]);
        }
    }

    class soE
    {
        sw sw = new sw();

        public uint[] getPrimes(uint n)
        {
            sw.Restart();
            int c, y, i, j; uint u; double logN; int[] x; uint[] p;
            if (n < 3) return n < 2 ? new uint[] { 0 } : new uint[] { 2, 1 };
            logN = Math.Log(n);
            c = (int)(n / logN * (1 + 1.2762 / logN));
            p = new uint[c + 1]; p[0] = 2; p[1] = 3;
            y = (int)(n / 3); x = new int[(y >> 5) + 1];
            if (n < 50000000) mark0(y, x);
            else { preMark(y >> 5, x); mark1(y, x); }
            Console.WriteLine(sw.Elapsed + " x");
            y = y <= 2 ? 0 : y - 2;
            j = y < 32 ? 2 : paraPrimes(y >> 5, x, p);
            i = y >> 5 << 5; u = 5 + 3 * (uint)i;
            while (i < y)
            {
                if ((x[i >> 5] & 1 << i++) == 0) p[j++] = u; u += 2;
                if ((x[i >> 5] & 1 << i++) == 0) p[j++] = u; u += 4;
            }
            if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) p[j++] = u; u += 2;
            if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) p[j++] = u;
            p[c] = (uint)j;
            sw.Stop(); Console.WriteLine(sw.Elapsed + " p");
            return p;
        }

        public int[] buildIsPrime(uint n)
        {
            sw.Restart();
            int y, i, u; int[] x, isP;
            if (n < 64) return new int[] { 0x64b4cb6e };
            isP = new int[(n >> 6) + 1]; isP[0] = 2;
            y = (int)(n / 3); x = new int[(y >> 5) + 1];
            if (n < 50000000) mark0(y, x);
            else { preMark(y >> 5, x); mark1(y, x); }
            Console.WriteLine(sw.Elapsed + " x");
            y -= 2; n >>= 1;
            paraIsPrime(y >> 5, x, isP);
            u = (y >> 5 << 4) * 3 + 2; i = y >> 5 << 5;
            while (i < y)
            {
                if ((x[i >> 5] & 1 << i++) == 0) isP[u >> 5] |= 1 << u; u += 1;
                if ((x[i >> 5] & 1 << i++) == 0) isP[u >> 5] |= 1 << u; u += 2;
            }
            if (u <= n && ((x[i >> 5] & 1 << i++) == 0)) isP[u >> 5] |= 1 << u; u += 1;
            if (u <= n && ((x[i >> 5] & 1 << i++) == 0)) isP[u >> 5] |= 1 << u;
            sw.Stop(); Console.WriteLine(sw.Elapsed + " isP");
            return isP;
        }

        void mark0(int y, int[] z)
        {
            int u, v, w, x, s, i, j;
            u = 8; v = 0; w = 3; x = 1; s = 7; i = 0; j = 10;
            while (s < y)
            {
                if ((z[i >> 5] & 1 << i++) == 0)
                {
                    int a = s, b = s + j, c = s + w, d = c + j, k = j << 1;
                    for (; d < y; a += k, b += k, c += k, d += k)
                    {
                        z[a >> 5] |= 1 << a; z[c >> 5] |= 1 << c;
                        z[b >> 5] |= 1 << b; z[d >> 5] |= 1 << d;
                    }
                    if (a < y)
                    {
                        z[a >> 5] |= 1 << a;
                        if (c < y) z[c >> 5] |= 1 << c;
                        if (b < y) z[b >> 5] |= 1 << b;
                    }
                }
                v += 8; s += v; x += 8; j += 4;
                if ((z[i >> 5] & 1 << i++) == 0)
                {
                    int a = s, b = s + j, c = s + x, d = c + j, k = j << 1;
                    for (; d < y; a += k, b += k, c += k, d += k)
                    {
                        z[a >> 5] |= 1 << a; z[c >> 5] |= 1 << c;
                        z[b >> 5] |= 1 << b; z[d >> 5] |= 1 << d;
                    }
                    if (a < y)
                    {
                        z[a >> 5] |= 1 << a;
                        if (c < y) z[c >> 5] |= 1 << c;
                        if (b < y) z[b >> 5] |= 1 << b;
                    }
                }
                u += 16; s += u; w += 4; j += 8;
            }
        }

        void preMark(int y, int[] z)
        {
            z[0] = 0x69128480;
            Parallel.For(0, 5, k =>
            {
                if (k == 0) for (int i = 1, j = y; i <= j; i += 5) z[i] = +0x12048120;
                if (k == 1) for (int i = 2, j = y; i <= j; i += 5) z[i] = +0x04812048;
                if (k == 2) for (int i = 3, j = y; i <= j; i += 5) z[i] = -0x7edfb7ee;
                if (k == 3) for (int i = 4, j = y; i <= j; i += 5) z[i] = +0x20481204;
                if (k == 4) for (int i = 5, j = y; i <= j; i += 5) z[i] = +0x48120481;
            });
            Parallel.For(0, 7, k =>
            {
                if (k == 0) for (int i = 1, j = y; i <= j; i += 7) z[i] |= +0x02100840;
                if (k == 1) for (int i = 2, j = y; i <= j; i += 7) z[i] |= +0x40210084;
                if (k == 2) for (int i = 3, j = y; i <= j; i += 7) z[i] |= -0x7bfdeff8;
                if (k == 3) for (int i = 4, j = y; i <= j; i += 7) z[i] |= +0x08402100;
                if (k == 4) for (int i = 5, j = y; i <= j; i += 7) z[i] |= +0x00840210;
                if (k == 5) for (int i = 6, j = y; i <= j; i += 7) z[i] |= +0x10084021;
                if (k == 6) for (int i = 7, j = y; i <= j; i += 7) z[i] |= +0x21008402;
            });
            Parallel.For(0, 11, k =>
            {
                if (k == 00) for (int i = 01, j = y; i <= j; i += 11) z[i] |= +0x20004080;
                if (k == 01) for (int i = 02, j = y; i <= j; i += 11) z[i] |= +0x04080010;
                if (k == 02) for (int i = 03, j = y; i <= j; i += 11) z[i] |= -0x7ffefe00;
                if (k == 03) for (int i = 04, j = y; i <= j; i += 11) z[i] |= +0x10200040;
                if (k == 04) for (int i = 05, j = y; i <= j; i += 11) z[i] |= +0x00040800;
                if (k == 05) for (int i = 06, j = y; i <= j; i += 11) z[i] |= +0x40800102;
                if (k == 06) for (int i = 07, j = y; i <= j; i += 11) z[i] |= +0x00102000;
                if (k == 07) for (int i = 08, j = y; i <= j; i += 11) z[i] |= +0x02000408;
                if (k == 08) for (int i = 09, j = y; i <= j; i += 11) z[i] |= +0x00408001;
                if (k == 09) for (int i = 10, j = y; i <= j; i += 11) z[i] |= +0x08001020;
                if (k == 10) for (int i = 11, j = y; i <= j; i += 11) z[i] |= +0x01020004;
            });
            z[1] |= 0x00800000;
            Parallel.For(0, 13, k =>
            {
                if (k == 00) for (int i = 02, j = y; i <= j; i += 13) z[i] |= +0x00020100;
                if (k == 01) for (int i = 03, j = y; i <= j; i += 13) z[i] |= +0x10000804;
                if (k == 02) for (int i = 04, j = y; i <= j; i += 13) z[i] |= -0x7fbfffe0;
                if (k == 03) for (int i = 05, j = y; i <= j; i += 13) z[i] |= +0x02010000;
                if (k == 04) for (int i = 06, j = y; i <= j; i += 13) z[i] |= +0x00080400;
                if (k == 05) for (int i = 07, j = y; i <= j; i += 13) z[i] |= +0x40002010;
                if (k == 06) for (int i = 08, j = y; i <= j; i += 13) z[i] |= +0x01000080;
                if (k == 07) for (int i = 09, j = y; i <= j; i += 13) z[i] |= +0x08040002;
                if (k == 08) for (int i = 10, j = y; i <= j; i += 13) z[i] |= +0x00201000;
                if (k == 09) for (int i = 11, j = y; i <= j; i += 13) z[i] |= +0x00008040;
                if (k == 10) for (int i = 12, j = y; i <= j; i += 13) z[i] |= +0x04000201;
                if (k == 11) for (int i = 13, j = y; i <= j; i += 13) z[i] |= +0x20100008;
                if (k == 12) for (int i = 14, j = y; i <= j; i += 13) z[i] |= +0x00804000;
            });
        }

        void mark1(int y, int[] z)
        {
            bool bl, bq; int u, v, w, x, s, i, j, n; queue q0, q1;
            bq = true; q0 = new queue(); q1 = new queue();
            bl = false; u = 40; v = 16; w = 11; x = 17; s = 95; i = 4; j = 34;
            for (n = 1024; n < y; n += 1024)
            {
                if (bl)
                {
                    bl = false;
                    if (bq)
                    {
                        bq = false;
                        for (int idxOut = 0; q0.count > 0; q0.count -= 3)
                        {
                            int a = q0.q[idxOut++], b = q0.q[idxOut++], k = q0.q[idxOut++];
                            for (; b < n; a += k, b += k) { z[a >> 5] |= 1 << a; z[b >> 5] |= 1 << b; }
                            for (; a < n; a += k) z[a >> 5] |= 1 << a;
                            if (a < b)
                            { q1.q[q1.idxIn++] = a; q1.q[q1.idxIn++] = b; }
                            else
                            { q1.q[q1.idxIn++] = b; q1.q[q1.idxIn++] = a; }
                            q1.q[q1.idxIn++] = k;
                        }
                        q1.count = q1.idxIn; q0.idxIn = 0;
                    }
                    else
                    {
                        bq = true;
                        for (int idxOut = 0; q1.count > 0; q1.count -= 3)
                        {
                            int a = q1.q[idxOut++], b = q1.q[idxOut++], k = q1.q[idxOut++];
                            for (; b < n; a += k, b += k) { z[a >> 5] |= 1 << a; z[b >> 5] |= 1 << b; }
                            for (; a < n; a += k) z[a >> 5] |= 1 << a;
                            if (a < b)
                            { q0.q[q0.idxIn++] = a; q0.q[q0.idxIn++] = b; }
                            else
                            { q0.q[q0.idxIn++] = b; q0.q[q0.idxIn++] = a; }
                            q0.q[q0.idxIn++] = k;
                        }
                        q0.count = q0.idxIn; q1.idxIn = 0;
                    }
                }
                while (s < n)
                {
                    if ((z[i >> 5] & 1 << i++) == 0)
                    {
                        bl = true; int a = s, b = s + w;
                        for (; b < n; a += j, b += j) { z[a >> 5] |= 1 << a; z[b >> 5] |= 1 << b; }
                        for (; a < n; a += j) z[a >> 5] |= 1 << a;
                        if (bq) q0.enqueue(a, b, j); else q1.enqueue(a, b, j);
                    }
                    v += 8; s += v; x += 8; j += 4;
                    if ((z[i >> 5] & 1 << i++) == 0)
                    {
                        bl = true; int a = s, b = s + x;
                        for (; b < n; a += j, b += j) { z[a >> 5] |= 1 << a; z[b >> 5] |= 1 << b; }
                        for (; a < n; a += j) z[a >> 5] |= 1 << a;
                        if (bq) q0.enqueue(a, b, j); else q1.enqueue(a, b, j);
                    }
                    u += 16; s += u; w += 4; j += 8;
                }
            }
            while (q0.count > 0)
            {
                int a = q0.dequeue(), b = q0.dequeue(), k = q0.dequeue();
                for (; b < y; a += k, b += k) { z[a >> 5] |= 1 << a; z[b >> 5] |= 1 << b; }
                for (; a < y; a += k) z[a >> 5] |= 1 << a;
            }
            while (q1.count > 0)
            {
                int a = q1.dequeue(), b = q1.dequeue(), k = q1.dequeue();
                for (; b < y; a += k, b += k) { z[a >> 5] |= 1 << a; z[b >> 5] |= 1 << b; }
                for (; a < y; a += k) z[a >> 5] |= 1 << a;
            }
        }

        int paraPrimes(int imax, int[] x, uint[] p)
        {
            int j0 = 0, j1 = 0, j2 = 0, j3 = 0;
            Parallel.For(0, 4, k =>
            {
                int i, j, y; uint u; i = 0; j = 2; u = 5;
                if (k == 1) { if (imax < 2) return; i = 1; j = 25; u = 101; }
                if (k == 2) { if (imax < 3) return; i = 2; j = 44; u = 197; }
                if (k == 3) { if (imax < 4) return; i = 3; j = 61; u = 293; }
                for (; ; )
                {
                    y = x[i];
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                    if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4;
                    i += 4;
                    if (i >= imax)
                    {
                        if (k == 0) j0 = j; if (k == 1) j1 = j; if (k == 2) j2 = j; if (k == 3) j3 = j;
                        break;
                    }
                    u += 288; j += bitsSet(i - 3, x) + bitsSet(i - 2, x) + bitsSet(i - 1, x);
                }
            });
            if (j0 < j1) j0 = j1; if (j2 < j3) j2 = j3; return j0 > j2 ? j0 : j2;
        }

        int bitsSet(int i, int[] x)
        {
            if (i < 0) return 0;
            i = ~x[i];
            i -= i >> 1 & 0x55555555;
            i = (i & 0x33333333) + (i >> 2 & 0x33333333);
            return (i + (i >> 4) & 0xf0f0f0f) * 0x1010101 >> 24;
        }

        void paraIsPrime(int imax, int[] x, int[] z)  // z = isP
        {
            int pc = Environment.ProcessorCount;
            for (int n = 0; n < 2; n++)
            {
                Parallel.For(0, pc, k =>
                {
                    int m, i, j, u, v, y;
                    if ((k & 1) == n)
                    {
                        j = pc * 2; v = (j - 1) * 48;
                        for (m = 0; m < 2; m++)
                        {
                            i = m + k * 2; u = 2 + i * 48;
                            for (; i < imax; i += j, u += v)
                            {
                                #region
                                y = x[i];
                                if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                                if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4;
                                if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                                if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4;
                                if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                                if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4;
                                if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                                if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4;
                                if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                                if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4;
                                if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                                if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4;
                                if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                                if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4;
                                if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                                if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                                if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2;
                                #endregion
                            }
                        }
                    }
                });
            }
        }
    }

    class queue
    {
        public int idxIn, idxOut, count;
        public int[] q = new int[20000];
        public void enqueue(int i, int j, int k)
        {
            count += 3;
            if (i < j)
            { q[idxIn++] = i; q[idxIn++] = j; }
            else
            { q[idxIn++] = j; q[idxIn++] = i; }
            q[idxIn++] = k;
        }
        public int dequeue() { count--; return q[idxOut++]; }
    }
}

2016/03/14

isPrime


/*
Another version of the sieve of Eratosthenes, ~30% faster than before.
BitArrays are replaced by int arrays.
"isPrime" was useful for a few record times at Project Euler.
 
Example: isP[0] = 0x64b4cb6e (each byte represents 8 odd numbers) 
    hex:           6        4         b        4        c        b        6        e
    bits:    0 1 1 0  0 1 0 0   1 0 1 1  0 1 0 0  1 1 0 0  1 0 1 1  0 1 1 0  1 1 1 0
    primes:    6 5      5       4   4 4    3      3 2      2   1 1    1 1    0 0 0   
               1 9      3       7   3 1    7      1 9      3   9 7    3 1    7 5 3
 
Results: 
          time (in seconds)    i7@3G6Hz              faster ? scroll down
      n     x      p   isP?    x = composites
    1e9  1.89   2.29   2.48    p = primes
    2e9  4.16   4.94   5.31    isP? = isPrime[i], i odd, 1 <= i <= n < 2^32 ~ 4e9
    4e9  9.02  10.54  11.31         = (isP[i >> 6] & 1 << (i << 26 >> 27)) != 0
                                    = (isP[i >> 6] & 1 << (i >> 1)) != 0
*/
using System;
using sw = System.Diagnostics.Stopwatch;
class isPrime
{
    static sw sw = new sw();
    static void Main()
    {
        test(1000000000);
        Console.Read();
    }

    static void test(uint n)
    {
        int[] isP;
        sw.Start();
        isP = buildIsPrime(n);
        sw.Stop();
        if (n > 0) n -= 1 - (n & 1);
        for (int i = 0; i < 10 && n > 1; n -= 2)
            if ((isP[n >> 6] & 1 << (int)(n >> 1)) != 0)
            { Console.WriteLine(n); i++; }
    }

    static int[] buildIsPrime(uint n)
    {
        int y, d, e, q, r, s, i, j, u; int[] x, isP;
        y = (int)(n / 3); x = new int[(y >> 5) + 1];
        d = 8; e = 0; q = 3; r = 1; s = 7; i = 0; j = 10;
        while (s < y)
        {
            if ((x[i >> 5] & 1 << i++) == 0)
            {
                int k = s, m = s + q;
                for (; m < y; k += j, m += j)
                {
                    x[k >> 5] |= 1 << k;
                    x[m >> 5] |= 1 << m;
                }
                for (; k < y; k += j)
                    x[k >> 5] |= 1 << k;
            }
            e += 8; s += e; r += 8; j += 4;
            if ((x[i >> 5] & 1 << i++) == 0 && s < y)
            {
                int k = s, m = s + r;
                for (; m < y; k += j, m += j)
                {
                    x[k >> 5] |= 1 << k;
                    x[m >> 5] |= 1 << m;
                }
                for (; k < y; k += j)
                    x[k >> 5] |= 1 << k;
            }
            d += 16; s += d; q += 4; j += 8;
        }
        Console.WriteLine(sw.Elapsed + " x");
        i = 0; y -= 2; d = y >> 5 << 5; u = 2; n >>= 1;
        isP = new int[(n >> 5) + 1]; isP[0] = 2;
        while (i < d)
        {
            q = x[i >> 5];
            for (e = 16; e > 0; e--)
            {
                if ((q & 1 << i++) == 0) isP[u >> 5] |= 1 << u; u += 1;
                if ((q & 1 << i++) == 0) isP[u >> 5] |= 1 << u; u += 2;
            }
        }
        while (i < y)
        {
            if ((x[i >> 5] & 1 << i++) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((x[i >> 5] & 1 << i++) == 0) isP[u >> 5] |= 1 << u; u += 2;
        }
        if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) isP[u >> 5] |= 1 << u; u += 1;
        if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) isP[u >> 5] |= 1 << u;
        Console.WriteLine(sw.Elapsed + " p");
        return isP;
    }
}


/*
With an unrolled loop.
*/
using System;
using sw = System.Diagnostics.Stopwatch;
class Eratosthenes
{
    static sw sw = new sw();
    static void Main()
    {
        test(1000000000);
        Console.Read();
    }

    static void test(uint n)
    {
        uint pi_n; uint[] p;
        sw.Start();
        p = getPrimes(n);
        sw.Stop();
        pi_n = p[p.Length - 1];
        Console.WriteLine("pi(" + n + ") = " + pi_n);
        if (pi_n > 0)
            Console.Write("largest prime = " + p[pi_n - 1]);
    }

    static uint[] getPrimes(uint n)
    {
        int c, y, d, e, q, r, s, i, j; uint u; double logN; int[] x; uint[] p;
        if (n < 3) return n < 2 ? new uint[] { 0 } : new uint[] { 2, 1 };
        logN = Math.Log(n);
        c = (int)(n / logN * (1 + 1.2762 / logN));
        p = new uint[c + 1]; p[0] = 2; p[1] = 3;
        y = (int)(n / 3); x = new int[(y >> 5) + 1];
        d = 8; e = 0; q = 3; r = 1; s = 7; i = 0; j = 10;
        while (s < y)
        {
            if ((x[i >> 5] & 1 << i++) == 0)
            {
                int k = s, m = s + q;
                for (; m < y; k += j, m += j) { x[k >> 5] |= 1 << k; x[m >> 5] |= 1 << m; }
                for (; k < y; k += j) x[k >> 5] |= 1 << k;
            }
            e += 8; s += e; r += 8; j += 4;
            if ((x[i >> 5] & 1 << i++) == 0 && s < y)
            {
                int k = s, m = s + r;
                for (; m < y; k += j, m += j) { x[k >> 5] |= 1 << k; x[m >> 5] |= 1 << m; }
                for (; k < y; k += j) x[k >> 5] |= 1 << k;
            }
            d += 16; s += d; q += 4; j += 8;
        }
        Console.WriteLine(sw.Elapsed + " x");
        j = 2; y -= 2; d = y >> 5; u = 5;
        for (i = 0; i < d; i++)
        {
            q = x[i];
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
            if ((q & 1) == 0) p[j++] = u; u += 2; if ((q & 2) == 0) p[j++] = u; u += 4;
        }
        i <<= 5;
        while (i < y)
        {
            if ((x[i >> 5] & 1 << i++) == 0) p[j++] = u; u += 2;
            if ((x[i >> 5] & 1 << i++) == 0) p[j++] = u; u += 4;
        }
        if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) p[j++] = u; u += 2;
        if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) p[j++] = u;
        p[c] = (uint)j;
        Console.WriteLine(sw.Elapsed + " p");
        return p;
    }
}
//      for (i = 0; i < d; i++)
//      {
//          q = x[i];
//          for (e = 16; e > 0; e--)
//          {
//              if ((q & 1) == 0) p[j++] = u; u += 2;
//              if ((q & 2) == 0) p[j++] = u; u += 4; q >>= 2;
//          }
//      }



/*****************************************************
"isPrime" as fast as Eratosthenes, ie isP(1e9) 2.29 s.
*/
using System;
using sw = System.Diagnostics.Stopwatch;
class isPrime
{
    static sw sw = new sw();
    static void Main()
    {
        test(1000000000);
        Console.Read();
    }

    static void test(uint n)
    {
        int[] isP;
        sw.Start();
        isP = buildIsPrime(n);
        sw.Stop();
        if (n > 0) n -= 1 - (n & 1);
        for (int i = 0; i < 10 && n > 1; n -= 2)
            if ((isP[n >> 6] & 1 << (int)(n >> 1)) != 0)
            { Console.WriteLine(n); i++; }
    }

    static int[] buildIsPrime(uint n)
    {
        int y, d, e, q, r, s, i, j, u; int[] x, isP;
        isP = new int[(n >> 6) + 1]; isP[0] = 2;
        y = (int)(n / 3); x = new int[(y >> 5) + 1];
        d = 8; e = 0; q = 3; r = 1; s = 7; i = 0; j = 10;
        while (s < y)
        {
            if ((x[i >> 5] & 1 << i++) == 0)
            {
                int k = s, m = s + q;
                for (; m < y; k += j, m += j) { x[k >> 5] |= 1 << k; x[m >> 5] |= 1 << m; }
                for (; k < y; k += j) x[k >> 5] |= 1 << k;
            }
            e += 8; s += e; r += 8; j += 4;
            if ((x[i >> 5] & 1 << i++) == 0 && s < y)
            {
                int k = s, m = s + r;
                for (; m < y; k += j, m += j) { x[k >> 5] |= 1 << k; x[m >> 5] |= 1 << m; }
                for (; k < y; k += j) x[k >> 5] |= 1 << k;
            }
            d += 16; s += d; q += 4; j += 8;
        }
        Console.WriteLine(sw.Elapsed + " x");
        y -= 2; d = y >> 5; u = 2; n >>= 1;
        for (i = 0; i < d; i++)
        {
            q = x[i];
            if ((q & 1) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 2) == 0) isP[u >> 5] |= 1 << u; u += 2;
            if ((q & 4) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 8) == 0) isP[u >> 5] |= 1 << u; u += 2; q >>= 4;
            if ((q & 1) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 2) == 0) isP[u >> 5] |= 1 << u; u += 2;
            if ((q & 4) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 8) == 0) isP[u >> 5] |= 1 << u; u += 2; q >>= 4;
            if ((q & 1) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 2) == 0) isP[u >> 5] |= 1 << u; u += 2;
            if ((q & 4) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 8) == 0) isP[u >> 5] |= 1 << u; u += 2; q >>= 4;
            if ((q & 1) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 2) == 0) isP[u >> 5] |= 1 << u; u += 2;
            if ((q & 4) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 8) == 0) isP[u >> 5] |= 1 << u; u += 2; q >>= 4;
            if ((q & 1) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 2) == 0) isP[u >> 5] |= 1 << u; u += 2;
            if ((q & 4) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 8) == 0) isP[u >> 5] |= 1 << u; u += 2; q >>= 4;
            if ((q & 1) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 2) == 0) isP[u >> 5] |= 1 << u; u += 2;
            if ((q & 4) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 8) == 0) isP[u >> 5] |= 1 << u; u += 2; q >>= 4;
            if ((q & 1) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 2) == 0) isP[u >> 5] |= 1 << u; u += 2;
            if ((q & 4) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 8) == 0) isP[u >> 5] |= 1 << u; u += 2; q >>= 4;
            if ((q & 1) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 2) == 0) isP[u >> 5] |= 1 << u; u += 2;
            if ((q & 4) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((q & 8) == 0) isP[u >> 5] |= 1 << u; u += 2;
        }
        i <<= 5;
        while (i < y)
        {
            if ((x[i >> 5] & 1 << i++) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((x[i >> 5] & 1 << i++) == 0) isP[u >> 5] |= 1 << u; u += 2;
        }
        if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) isP[u >> 5] |= 1 << u; u += 1;
        if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) isP[u >> 5] |= 1 << u;
        Console.WriteLine(sw.Elapsed + " p");
        return isP;
    }
}



/***********************************************************
Composites(4e9) 8.98 s, pi(4e9) 9.05 s.  Bit Twiddling Hacks  
*/
using System;
using sw = System.Diagnostics.Stopwatch;
class countPrimes
{
    static sw sw = new sw();
    static void Main()
    {
        test(4000000000);
    }

    static void test(uint n)
    {
        sw.Start();
        Console.Write("pi(" + n + ") = " + pi(n));
        sw.Stop();
        Console.Read();
    }

    static uint pi(uint n)
    {
        int y, d, e, q, r, s, i, j; uint c, u; uint[] x;
        if (n < 3) return n >> 1;
        y = (int)(n / 3); x = new uint[(y >> 5) + 1];
        d = 8; e = 0; q = 3; r = 1; s = 7; i = 0; j = 10;
        while (s < y)
        {
            if ((x[i >> 5] & 1 << i++) == 0)
            {
                int k = s, m = s + q;
                for (; m < y; k += j, m += j) { x[k >> 5] |= 1u << k; x[m >> 5] |= 1u << m; }
                for (; k < y; k += j) x[k >> 5] |= 1u << k;
            }
            e += 8; s += e; r += 8; j += 4;
            if ((x[i >> 5] & 1 << i++) == 0 && s < y)
            {
                int k = s, m = s + r;
                for (; m < y; k += j, m += j) { x[k >> 5] |= 1u << k; x[m >> 5] |= 1u << m; }
                for (; k < y; k += j) x[k >> 5] |= 1u << k;
            }
            d += 16; s += d; q += 4; j += 8;
        }
        Console.WriteLine(sw.Elapsed + " x");
        c = 2; y -= 2; d = y >> 5;
        for (i = 0; i < d; i++)
        {                                       // Bit Twiddling Hacks
            u = ~x[i];
            u -= u >> 1 & 0x55555555;
            u = (u & 0x33333333) + (u >> 2 & 0x33333333);
            c += (u + (u >> 4) & 0xf0f0f0f) * 0x1010101 >> 24;
        }
        i <<= 5; u = 5 + 3 * (uint)i;
        while (i < y)
        {
            if ((x[i >> 5] & 1 << i++) == 0) c++;
            if ((x[i >> 5] & 1 << i++) == 0) c++; u += 6;
        }
        if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) c++; u += 2;
        if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) c++;
        Console.WriteLine(sw.Elapsed + " c");
        return c;
    }
}



/*******************************************************************************
primes(1e9) 2.12 s, primes(2e9) 4.61 s, primes(4e9) 9.90 s, primes(~0u) 10.65 s.
*/
using System;
using System.Threading.Tasks;
using sw = System.Diagnostics.Stopwatch;
class Eratosthenes
{
    static sw sw = new sw();
    static void Main()
    {
        test(1000000000);
        Console.Read();
    }

    static void test(uint n)
    {
        uint pi_n; uint[] p;
        sw.Start();
        p = getPrimes(n);
        sw.Stop();
        pi_n = p[p.Length - 1];
        Console.WriteLine("pi(" + n + ") = " + pi_n);
        if (pi_n > 0)
            Console.Write("largest prime = " + p[pi_n - 1]);
    }

    static uint[] getPrimes(uint n)
    {
        int c, y, d, e, q, r, s, i, j; uint u; double logN; int[] x; uint[] p;
        if (n < 3) return n < 2 ? new uint[] { 0 } : new uint[] { 2, 1 };
        logN = Math.Log(n);
        c = (int)(n / logN * (1 + 1.2762 / logN));
        p = new uint[c + 1]; p[0] = 2; p[1] = 3;
        y = (int)(n / 3); x = new int[(y >> 5) + 1];
        d = 8; e = 0; q = 3; r = 1; s = 7; i = 0; j = 10;
        while (s < y)
        {
            if ((x[i >> 5] & 1 << i++) == 0)
            {
                int k = s, m = s + q;
                for (; m < y; k += j, m += j) { x[k >> 5] |= 1 << k; x[m >> 5] |= 1 << m; }
                for (; k < y; k += j) x[k >> 5] |= 1 << k;
            }
            e += 8; s += e; r += 8; j += 4;
            if ((x[i >> 5] & 1 << i++) == 0 && s < y)
            {
                int k = s, m = s + r;
                for (; m < y; k += j, m += j) { x[k >> 5] |= 1 << k; x[m >> 5] |= 1 << m; }
                for (; k < y; k += j) x[k >> 5] |= 1 << k;
            }
            d += 16; s += d; q += 4; j += 8;
        }
        Console.WriteLine(sw.Elapsed + " x");
        y = y <= 2 ? 0 : y - 2;
        j = y < 32 ? 2 : paraPrimes(y >> 5, x, p);
        i = y >> 5 << 5; u = 5 + 3 * (uint)i;
        while (i < y)
        {
            if ((x[i >> 5] & 1 << i++) == 0) p[j++] = u; u += 2;
            if ((x[i >> 5] & 1 << i++) == 0) p[j++] = u; u += 4;
        }
        if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) p[j++] = u; u += 2;
        if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) p[j++] = u;
        p[c] = (uint)j;
        Console.WriteLine(sw.Elapsed + " p");
        return p;
    }

    static int paraPrimes(int imax, int[] x, uint[] p)
    {
        int j0 = 0, j1 = 0;
        Parallel.For(0, 2, k =>
        {
            int i, j, y; uint u;
            if (k == 0) { i = 0; j = 2; u = 5; }
            else { if (imax < 2) return; i = 1; j = 25; u = 101; }
            for (; ; )
            {
                y = x[i];
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4; y >>= 2;
                if ((y & 1) == 0) p[j++] = u; u += 2; if ((y & 2) == 0) p[j++] = u; u += 4;
                i += 2; if (i >= imax) { if (k == 0) j0 = j; else j1 = j; break; }
                u += 96; j += bitsSet(~x[i - 1]);
            }
        });
        return j0 > j1 ? j0 : j1;
    }

    static int bitsSet(int i)
    {
        i -= i >> 1 & 0x55555555;
        i = (i & 0x33333333) + (i >> 2 & 0x33333333);
        return (i + (i >> 4) & 0xf0f0f0f) * 0x1010101 >> 24;
    }
}


/*******************************************************************
Faster isPrime                                             time (s)
                                                     n       x  isP?
                                                   1e9    1.84  2.00
                                                   2e9    4.08  4.38
                                                   4e9    8.84  9.40
*/
using System;
using System.Threading.Tasks;
using sw = System.Diagnostics.Stopwatch;
class isPrime
{
    static sw sw = new sw();
    static void Main()
    {
        test((uint)1e9);
        Console.Read();
    }

    static void test(uint n)
    {
        int[] isP;
        sw.Start();
        isP = buildIsPrime(n);
        sw.Stop();
        int check = 0;                               //   n    check
        for (int i = isP.Length - 1; i >= 0; i--)    // 1e9 1DF99124
            check ^= isP[i];                         // 2e9 E9D92C9E
        Console.Write("check: {0:X8}", check);       // 4e9 DB1F9720
    }

    static int[] buildIsPrime(uint n)
    {
        int y, i, u; int[] x, isP;
        if (n < 64) return new int[] { 0x64b4cb6e };
        isP = new int[(n >> 6) + 1]; isP[0] = 2;
        y = (int)(n / 3); x = new int[(y >> 5) + 1];
        mark(y, x);
        y -= 2; n >>= 1;
        paraIsPrime(y >> 5, x, isP);
        u = (y >> 5 << 4) * 3 + 2; i = y >> 5 << 5;
        while (i < y)
        {
            if ((x[i >> 5] & 1 << i++) == 0) isP[u >> 5] |= 1 << u; u += 1;
            if ((x[i >> 5] & 1 << i++) == 0) isP[u >> 5] |= 1 << u; u += 2;
        }
        if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) isP[u >> 5] |= 1 << u; u += 1;
        if (u - 3 <= n - 3 && ((x[i >> 5] & 1 << i++) == 0)) isP[u >> 5] |= 1 << u;
        Console.WriteLine(sw.Elapsed + " p");
        return isP;
    }

    static void mark(int y, int[] z)
    {
        int u, v, w, x, s, i, j;
        u = 8; v = 0; w = 3; x = 1; s = 7; i = 0; j = 10;
        while (s < y)
        {
            if ((z[i >> 5] & 1 << i++) == 0)
            {
                int a = s, b = s + j, c = s + w, d = c + j, k = j << 1;
                for (; d < y; a += k, b += k, c += k, d += k)
                {
                    z[a >> 5] |= 1 << a; z[c >> 5] |= 1 << c;
                    z[b >> 5] |= 1 << b; z[d >> 5] |= 1 << d;
                }
                if (a < y)
                {
                    z[a >> 5] |= 1 << a;
                    if (c < y) z[c >> 5] |= 1 << c;
                    if (b < y) z[b >> 5] |= 1 << b;
                }
            }
            v += 8; s += v; x += 8; j += 4;
            if ((z[i >> 5] & 1 << i++) == 0)
            {
                int a = s, b = s + j, c = s + x, d = c + j, k = j << 1;
                for (; d < y; a += k, b += k, c += k, d += k)
                {
                    z[a >> 5] |= 1 << a; z[c >> 5] |= 1 << c;
                    z[b >> 5] |= 1 << b; z[d >> 5] |= 1 << d;
                }
                if (a < y)
                {
                    z[a >> 5] |= 1 << a;
                    if (c < y) z[c >> 5] |= 1 << c;
                    if (b < y) z[b >> 5] |= 1 << b;
                }
            }
            u += 16; s += u; w += 4; j += 8;
        }
        Console.WriteLine(sw.Elapsed + " x");
    }

    static void paraIsPrime(int imax, int[] x, int[] z)  // z = isP
    {
        int pc = Environment.ProcessorCount;
        for (int n = 0; n < 2; n++)
        {
            Parallel.For(0, pc, k =>
            {
                int m, i, j, u, v, y;
                if ((k & 1) == n)
                {
                    j = pc * 2; v = (j - 1) * 48;
                    for (m = 0; m < 2; m++)
                    {
                        i = m + k * 2; u = 2 + i * 48;
                        for (; i < imax; i += j, u += v)
                        {
                            #region
                            y = x[i];
                            if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                            if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4;
                            if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                            if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4;
                            if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                            if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4;
                            if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                            if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4;
                            if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                            if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4;
                            if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                            if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4;
                            if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                            if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2; y >>= 4;
                            if ((y & 1) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 2) == 0) z[u >> 5] |= 1 << u; u += 2;
                            if ((y & 4) == 0) z[u >> 5] |= 1 << u; u += 1;
                            if ((y & 8) == 0) z[u >> 5] |= 1 << u; u += 2;
                            #endregion
                        }
                    }
                }
            });
        }
    }
}

2016/02/10

Partition numbers


/*
                        ms    i7@3G6Hz
    P[0]...P[10000]    100  
    P[0]...P[20000]    305  
    P[0]...P[40000]    970   
    P[0]...P[80000]   3420  
    p[0]..P[160000]  13800 
    P[0]..P[320000]  53000    P[320000] 624 digits
*/
using System;                                                                         
using System.Diagnostics;
using Xint = System.Numerics.BigInteger;
class Program
{
    static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        uint n = 10000;
        sw.Start();
        Xint[] P = A000041(n);
        sw.Stop();
        Console.Write(sw.Elapsed + "\n" + P[n]);
        Console.Read();
    }

    static Xint[] A000041(uint n)
    {
        int i, d, d0, d1; Xint S; Xint[] P;
        P = new Xint[n + 1]; P[0] = 1;
        for (i = 1; i <= n; i++)
        {
            S = d = d1 = 0; d0 = -2;
            do
            {
                d += d0 += 4; if (i >= d) S += P[i - d];
                d -= d1 += 1; if (i >= d) S += P[i - d];
                if (i <= d) break;
                d += d0 += 4; if (i >= d) S -= P[i - d];
                d -= d1 += 1; if (i >= d) S -= P[i - d];
            }
            while (i > d);
            P[i] = S;
        }
        return P;
    }
}

2015/08/27

number of primes pi(x)

/*  
    Adrien-Marie Legendre (1752-1833) observed that pi(x), the number of primes p <= x , 
    can be calculated with the inclusion-exclusion principle.

    pi(x) - pi(sqrt(x)) + 1 = x - SIGMA(x/pi) + SIGMA(x/pi/pj) - SIGMA(x/pi/pj/pk) + ...
                                    i            i<j             i<j<k

        x = 26 , sqrt(26) = 5 , primes <= 5 : 2 3 5 , pi(5) = 3
        pi(26) - 3 + 1 = 26 - 26/2 - 26/3 - 26/5 + 26/2/3 + 26/2/5 + 26/3/5 - 26/2/3/5
                       = 26 - 13   - 8    - 5    + 4      + 2      + 1      - 0
                       = 26 - 26 + 7 
        pi(26) = 9 
  
               x         pi(x)  minutes   seconds                          (Athlon 3GHz)
            1e09      50847534               2.01  getPrimes(1e9): 12.5 s
            1e10     455052511              20.90
            1e11    4118054813        4    216.00
            1e12   37607912018       38   2260.00
            1e13  346065536839      399  23900.00
  
    Lagarias, Miller & Odlyzko: pi(1e13) 8 minutes (1985 on an IBM/370 3081 Model K).
*/
using System;
using System.Collections;
using System.Collections.Generic;
using System.Diagnostics;
class pi_x
{
    static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        //test();
        sw.Start();
        ulong r = pi((ulong)1e9);
        sw.Stop();
        Console.Write(sw.Elapsed + " " + r);
        Console.Read();
    }

    /*
    static uint _s32;
    static uint pi32(uint x)
    {
        _p = getPrimes((uint)Math.Sqrt(x));
        _pLen = _p.Length;
        _s32 = 0;
        sxnp32(x, 0, 0);
        return _s32 + (uint)_pLen - 1;
    }
    static void sxnp32(uint q, int n, int pIdx) // sum x divided by n primes
    {
        if ((n & 1) == 0)
            _s32 += q;
        else
            _s32 -= q;
        for (int i = pIdx; i < _pLen; i++)
            if (q < _p[i]) 
                break;
            else
                sxnp32(q / _p[i], n + 1, i + 1);
    }
    */

    static uint[] _p;
    static int _pLen;
    static ulong _s;

    static ulong pi(ulong x)
    {
        if (x == 0) return 0;
        _p = getPrimes((uint)Math.Sqrt(x));
        _pLen = _p.Length;
        _s = 0;
        if (x <= ~0u)
            sxnp((uint)x, 0, 0);
        else
            sxnp(x, 0, 0);
        return _s + (uint)_pLen - 1;
    }

    static void sxnp(uint q, int n, int pIdx)
    {
        if ((n & 1) == 0)
        {
            _s += q;
            for (int i = pIdx; i < _pLen; i++)
                if (q < _p[i]) break;
                else
                    if (q >> 1 < _p[i]) _s--;
                    else sxnp(q / _p[i], n + 1, i + 1);
        }
        else
        {
            _s -= q;
            for (int i = pIdx; i < _pLen; i++)
                if (q < _p[i]) break;
                else
                    if (q >> 1 < _p[i]) _s++;
                    else sxnp(q / _p[i], n + 1, i + 1);
        }
    }

    static void sxnp(ulong q, int n, int pIdx)
    {
        if ((n & 1) == 0)
        {
            _s += q;
            for (int i = pIdx; i < _pLen; i++)
            {
                if (q < _p[i]) break;
                if (q >> 1 < _p[i]) _s--;
                else
                    if (q > ~0u) sxnp(q / _p[i], n + 1, i + 1);
                    else sxnp((uint)q / _p[i], n + 1, i + 1);
            }
        }
        else
        {
            _s -= q;
            for (int i = pIdx; i < _pLen; i++)
            {
                if (q < _p[i]) break;
                if (q >> 1 < _p[i]) _s++;
                else
                    if (q > ~0u) sxnp(q / _p[i], n + 1, i + 1);
                    else sxnp((uint)q / _p[i], n + 1, i + 1);
            }
        }
    }

    static uint[] getPrimes(uint n) // primes <= n (Eratosthenes)   
    {
        int y, d, e, q, r, s, i, j;
        double logN; List<uint> p; BitArray x; uint u;
        if (n < 3) return n == 2 ? new uint[] { 2 } : new uint[] { };
        logN = Math.Log(n);
        p = new List<uint>((int)(n / logN * (1 + 1.2762 / logN)));
        x = new BitArray((int)(n / 3));
        y = x.Length; d = 8; e = 0; q = 3; r = 1; s = 7; i = 0; j = 10;
        while (s < y)
        {
            if (!x[i++])
            {
                int k = s, m = s + q;
                for (; m < y; k += j, m += j) { x[k] = true; x[m] = true; }
                for (; k < y; k += j) x[k] = true;
            }
            e += 8; s += e; r += 8; j += 4;
            if (!x[i++] && s < y)
            {
                int k = s, m = s + r;
                for (; m < y; k += j, m += j) { x[k] = true; x[m] = true; }
                for (; k < y; k += j) x[k] = true;
            }
            d += 16; s += d; q += 4; j += 8;
        }
        p.Add(2); p.Add(3); i = 0; y -= 2; u = 5;
        while (i < y)
        {
            if (!x[i++]) p.Add(u); u += 2;
            if (!x[i++]) p.Add(u); u += 4;
        }
        if (u <= n && !x[i++]) p.Add(u); u += 2;
        if (u <= n && !x[i]) p.Add(u);
        return p.ToArray();
    }

    static void test()
    {
        for (uint x = 0; x <= 50; x++)
            Console.WriteLine("{0,2} {1,2}", x, pi(x));
    }
}

2014/07/28

Sum of all divisors of all positive integers <= n.

/*
Example:  N  D        SD  SSD  A024916
          1  1         1    1
          2  1,2       3    4
          3  1,3       4    8
          4  1,2,4     7   15
          5  1,5       6   21
          6  1,2,3,6  12   33
  

Sum of all divisors of all positive integers <= 10^n.

                      A072692 as a simple table
 
     n                       a(n)                                  time
                                                          SSD               SSD3
        
     0                                           1  00:00:00.0000061
     1                                          87  00:00:00.0000058
     2                                        8299  00:00:00.0000092
     3                                      823081  00:00:00.0000215
     4                                    82256014  00:00:00.0000555
     5                                  8224740835  00:00:00.0002824
     6                                822468118437  00:00:00.0007623
     7                              82246711794796  00:00:00.0022483
     8                            8224670422194237  00:00:00.0075755
     9                          822467034112360628  00:00:00.0247891
    10                        82246703352400266400  00:00:00.1029980  00:00:00.0536221
    11                      8224670334323560419029  00:00:00.3310292
    12                    822467033425357340138978  00:00:01.0792704  00:00:00.5139608
    13                  82246703342420509396897774  00:00:03.3665267
    14                8224670334241228180927002517  00:00:10.7148327  00:00:05.1645348  
    15              822467033424114009326065894639  00:00:34.2567848
    16            82246703342411333689227187822414  00:01:49.7741595  00:00:53.1228885
    17          8224670334241132270081671519064067  00:05:55.3252935
    18        822467033424113219363487627735401433  00:19:50.2083007
    19      82246703342411321831834750379375676781  01:12:36.6691244
    20    8224670334241132182459113152714495956789  04:38:51.4847697
    21  822467033424113218236966661186847524013409  15:19:44.3197506


 AMD Athlon II X4 640, 3 GHz, 2 GB, XP  
  
    N=10^16     time(s)   CPU Usage(%)
      SSD(N)     110        25 
     SSD3(N)      53        60
     FSSD(N)      92        25
    FSSD3(N)      41        60
    FSSD4(N)      33        80
*/

using System;
using System.Diagnostics;
using System.Threading.Tasks;
using Xint = System.Numerics.BigInteger;
class A024916
{
    static Xint SSD(Xint N)
    {
        Xint D = 1, Q = N, S = 0;
        while (D < Q)
        {
            S += D++ * (Q - D);
            S += Q++ * Q / 2;
            Q = N / D;
        }
        S += Q++ * Q / 2;
        S -= D-- * D-- * D / 6;
        return S;
    }

    static void Main()
    {
        var sw = new Stopwatch();
        int p = 0;
        Xint N = 1, S;
    L0: sw.Restart();
        S = SSD(N);
        sw.Stop();
        Console.WriteLine(p + "  " + S + "  " + sw.Elapsed);
        p++;
        N *= 10;
        goto L0;
    }

    private static Xint[] SA, DA, QA;

    static Xint SSD3(Xint N)
    {
        SA = new Xint[3];
        DA = new Xint[3];
        QA = new Xint[3];
        Parallel.For(0, 3, (int d) => SSDP3(N, d));
        Xint S = SA[0] + SA[1] + SA[2],
             D = DA[0],
             Q = QA[0];
        if (DA[1] < D)
        {
            D = DA[1];
            Q = QA[1];
        }
        if (DA[2] < D)
        {
            D = DA[2];
            Q = QA[2];
        }
        S += Q++ * Q / 2;
        S -= D-- * D-- * D / 6;
        return S;
    }

    static void SSDP3(Xint N, int d)
    {
        Xint D = d + 1, Q = N / D, S = 0;
        while (D < Q)
        {
            S += D * (Q - D - 1);
            S += Q * (Q + 1) / 2;
            D += 3;
            Q = N / D;
        }
        SA[d] = S;
        DA[d] = D;
        QA[d] = Q;
    }

    static Xint FSSD(Xint N)
    {
        Xint D = 1, Q = N, S = 0;
        while (D < Q)
        {
            S += Q * (Q + 1 + 2 * D) / 2;
            D += 1;
            Q = N / D;
        }
        S += Q * (Q + 1) / 2;
        S -= D * D * (D - 1) / 2;
        return S;
    }

    static Xint FSSD3(Xint N)
    {
        SA = new Xint[3];
        DA = new Xint[3];
        QA = new Xint[3];
        Parallel.For(0, 3, (int d) => FSSDP3(N, d));
        Xint S = SA[0] + SA[1] + SA[2],
             D = DA[0],
             Q = QA[0];
        if (DA[1] < D)
        {
            D = DA[1];
            Q = QA[1];
        }
        if (DA[2] < D)
        {
            D = DA[2];
            Q = QA[2];
        }
        S += Q * (Q + 1) / 2;
        S -= D * D * (D - 1) / 2;
        return S;
    }

    static void FSSDP3(Xint N, int d)
    {
        Xint D = d + 1, Q = N / D, S = 0;
        while (D < Q)
        {
            S += Q * (Q + 1 + 2 * D) / 2;
            D += 3;
            Q = N / D;
        }
        SA[d] = S;
        DA[d] = D;
        QA[d] = Q;
    }

    static Xint FSSD4(Xint N)
    {
        SA = new Xint[4];
        DA = new Xint[4];
        QA = new Xint[4];
        Parallel.For(0, 4, (int d) => FSSDP4(N, d));
        Xint S = SA[0] + SA[1] + SA[2] + SA[3],
             D = DA[0],
             Q = QA[0];
        if (DA[1] < D)
        {
            D = DA[1];
            Q = QA[1];
        }
        if (DA[2] < D)
        {
            D = DA[2];
            Q = QA[2];
        }
        if (DA[3] < D)
        {
            D = DA[3];
            Q = QA[3];
        }
        S += Q * (Q + 1) / 2;
        S -= D * D * (D - 1) / 2;
        return S;
    }

    static void FSSDP4(Xint N, int d)
    {
        Xint D = d + 1, Q = N / D, S = 0;
        while (D < Q)
        {
            S += Q * (Q + 1 + 2 * D) / 2;
            D += 4;
            Q = N / D;
        }
        SA[d] = S;
        DA[d] = D;
        QA[d] = Q;
    }
}