Pagina's

2016/08/14

Odd prime sieve


/* 105,000,000 Odd primes in 1 second on an i7-4790. It's almost a sieve of Eratosthenes.
Composites are marked by 24 threads running in cache friendly ranges, after they're counted,
a prime array of exactly the right size is filled with ... primes. Bit twiddling hacks used:
"count bits set, in parallel" and "count trailing zero bits with a de Bruijn sequence".
...Code colorized with: C# Syntax Highlighter...                  125,000,000 ? scroll down
 
            Old times      New times    (in seconds).  
      n      x      p       x      p    p(rimes) <= n < 2^32
    1e9   0.50   0.66    0.31   0.42    x: odd composites, marked in 
    2e9   1.02   1.31    0.63   0.84       essentially a bit array.
    4e9   2.13   2.65    1.29   1.70    best of 5
    ~0u   2.28   2.87    1.39   1.82    best of 5
*/
using System;
using System.Threading.Tasks;
using sw = System.Diagnostics.Stopwatch;
class oddPrimes
{
    static sw sw = new sw();
    static void Main()
    {
        uint m = (uint)1e9;
        sw.Start(); uint[] p = getPrimes(m);
        Console.Write(p.Length + " odd primes <= " + m);
        Console.Read();
    }

    static uint[] getPrimes(uint m) { return countPrimes(buildX(m)); }

    static uint[] buildX(uint m)  // odd composites
    {
        uint i, j, t, r; uint[] x, r0, r1;
        t = (uint)Environment.ProcessorCount * 3; r = 32768 * 6;
        m -= m / ~0u; m += m & 1; m /= 2;
        x = new uint[(m >> 5) + 1];
        x[m >> 5] = ~0u << (int)m;
        if (m > 32) x[0] = 0x9b4b3491;
        for (i = j = 0; i < m; i += j += 4) x[i >> 5] |= 1u << (int)i;
        r0 = new uint[t + 1]; r1 = new uint[t + 1];
        for (i = 0; i <= t; i++) { r0[i] = i * r; r1[i] = (i + 1) * r; } r *= t;
        while (r0[0] < m)
        {
            Parallel.For(0, t, k =>
            {
                uint b, c, d, e, f, g, m0, m1;
                if ((m0 = r0[k]) >= m) return;
                if ((m1 = r1[k]) > m) m1 = m;
                for (d = 3, b = e = 4; ; b += e += 4, d += 2)
                {
                    if ((c = b + d) >= m1) break;
                    if ((x[d >> 6] & (1 << (int)(d >> 1))) == 0)
                    {
                        if (c < m0) if ((c += (m0 - c) / d * d) < m0) c += d;
                        for (g = d << 1, f = c + d; f < m1; c += g, f += g)
                        { x[c >> 5] |= 1u << (int)c; x[f >> 5] |= 1u << (int)f; }
                        for (; c < m1; c += g) x[c >> 5] |= 1u << (int)c;
                    }
                }
            });
            for (j = 0; j <= t; j++) { r0[j] += r; r1[j] += r; }
        }
        time("x");
        return x;
    }

    static uint[] countPrimes(uint[] x)
    {
        int t = Environment.ProcessorCount; uint[] c = new uint[t];
        byte[] b = new byte[x.Length]; uint n;
        Parallel.For(0, t, k =>
        {
            uint xi, ck = 0;
            for (int i = x.Length - 1 - k; i >= 0; i -= t)
            {
                xi = x[i]; xi -= (xi >> 1) & 0x55555555;
                xi = (xi & 0x33333333) + (xi >> 2 & 0x33333333);
                ck += b[i] = (byte)(32 - ((xi + (xi >> 4) & 0xf0f0f0f) * 0x1010101 >> 24));
            }
            c[k] = ck;
        });
        for (n = c[0], t--; t > 0; t--) n += c[t];
        time("c");
        return buildPrimes(b, x, n);
    }

    static uint[] buildPrimes(byte[] b, uint[] x, uint n)
    {
        int j = b.Length - 1; uint[] c = new uint[8]; uint[] p = new uint[n];
        sbyte[] bruijn = {01,02,29,03,30,15,25,04,31,23,21,16,26,18,05,09,
                          32,28,14,24,22,20,17,08,27,13,19,07,12,06,11,10};
        for (uint s = 0, k = 1; k < 8 && k <= j; k++) c[k] += s += b[k - 1];
        Parallel.For(0, 8, k =>
        {
            int s, i = k; uint ck = c[k], ci, u, v; long xi;
            if (k > 0 && ck == 0) return;
            for (u = (uint)(k << 6) - 1; ; u += 512)
            {
                v = u; ci = ck; xi = ~x[i];
                while (xi > 0)
                {
                    s = bruijn[((uint)((xi & -xi) * 0x077cb531)) >> 27];
                    p[ci++] = v += (uint)s << 1;
                    xi >>= s;
                }
                i += 8; if (i > j) break;
                ck += (uint)(b[i - 1] + b[i - 2] + b[i - 3] + b[i - 4] +
                             b[i - 5] + b[i - 6] + b[i - 7] + b[i - 8]);
            }
        });
        time("p");
        return p;
    }

    static void time(string s) { Console.WriteLine(sw.ElapsedMilliseconds + " ms " + s); }
}
/* Nothing new under the sun, pre-sieving 3, 5 & 7.
First I did it the old way (scroll down), but a hard
coded table works faster. Actually the table can be
made "on the fly", with a few more small primes?

            Old times      New times (in seconds).
      n      x      p       x      p
    1e9   0.31   0.42    0.24   0.34
    2e9   0.63   0.84    0.48   0.68
    4e9   1.29   1.70    0.99   1.50 (average of 10 runs)
    ~0u   1.39   1.82    1.07   1.60 (average of 10000 ;)
*/
using System;
using System.Threading.Tasks;
using sw = System.Diagnostics.Stopwatch;
class oddPrimes
{
    static sw sw = new sw();
    static void Main()
    {
        uint m = (uint)1e9;
        sw.Start(); uint[] p = getPrimes(m);
        Console.Write(p.Length + " odd primes <= " + m);
        Console.Read();
    }

    static uint[] getPrimes(uint m) { return countPrimes(buildX(m)); }

    static uint[] buildX(uint m)
    {
        uint i, j, t, r; uint[] x, r0, r1;
        t = (uint)Environment.ProcessorCount * 3; r = 32768 * 6;
        m -= m / ~0u; m += m & 1; m >>= 1;
        x = new uint[(m >> 5) + 1]; x[0] = 0x9b4b3491;
        preMark((int)m >> 5, x);
        for (i = j = 0; i < m; i += j += 4) x[i >> 5] |= 1u << (int)i;
        r0 = new uint[t]; r1 = new uint[t];
        for (i = 0; i < t; i++) { r0[i] = i * r; r1[i] = (i + 1) * r; } r *= t;
        while (r0[0] < m)
        {
            Parallel.For(0, t, k =>
            {
                uint b, c, d, e, f, g, m0, m1;
                if ((m0 = r0[k]) >= m) return;
                if ((m1 = r1[k]) > m) m1 = m;
                for (d = 11, b = 60, e = 20; ; b += e += 4, d += 2)
                {
                    if ((c = b + d) >= m1) break;
                    if ((x[d >> 6] & (1 << (int)(d >> 1))) == 0)
                    {
                        if (c < m0) if ((c += (m0 - c) / d * d) < m0) c += d;
                        for (g = d << 1, f = c + d; f < m1; c += g, f += g)
                        { x[c >> 5] |= 1u << (int)c; x[f >> 5] |= 1u << (int)f; }
                        for (; c < m1; c += g) x[c >> 5] |= 1u << (int)c;
                    }
                }
            });
            for (j = 0; j < t; j++) { r0[j] += r; r1[j] += r; }
        }
        x[m >> 5] |= ~0u << (int)m;
        time("x");
        return x;
    }

    static void preMark(int m, uint[] x)
    {
        uint[] a =  // 3*5*7 uints = 420 bytes
        {
            0x6e92ed65,0x59a5b34d,0x96693cf2,0x25dacb36,0x4b669add,0xd279e4b3,0xb5966d2c,  // No
            0xcd35ba4b,0xf3c96696,0x2cda59a4,0x6b74976b,0x92cd2d9a,0xb4b349e7,0xe92ed659,  // co
            0x9a5b34d6,0x6693cf25,0x5dacb369,0xb669add2,0x279e4b34,0x5966d2cd,0xd35ba4bb,  // de
            0x3c96696c,0xcda59a4f,0xb74976b2,0x2cd2d9a6,0x4b349e79,0x92ed659b,0xa5b34d6e,  // an
            0x693cf259,0xdacb3696,0x669add25,0x79e4b34b,0x966d2cd2,0x35ba4bb5,0xc96696cd,  // aly
            0xda59a4f3,0x74976b2c,0xcd2d9a6b,0xb349e792,0x2ed659b4,0x5b34d6e9,0x93cf259a,  // sis
            0xacb36966,0x69add25d,0x9e4b34b6,0x66d2cd27,0x5ba4bb59,0x96696cd3,0xa59a4f3c,  // iss
            0x4976b2cd,0xd2d9a6b7,0x349e792c,0xed659b4b,0xb34d6e92,0x3cf259a5,0xcb369669,  // ue
            0x9add25da,0xe4b34b66,0x6d2cd279,0xba4bb596,0x6696cd35,0x59a4f3c9,0x976b2cda,  // s
            0x2d9a6b74,0x49e792cd,0xd659b4b3,0x34d6e92e,0xcf259a5b,0xb3696693,0xadd25dac,  // we
            0x4b34b669,0xd2cd279e,0xa4bb5966,0x696cd35b,0x9a4f3c96,0x76b2cda5,0xd9a6b749,  // re
            0x9e792cd2,0x659b4b34,0x4d6e92ed,0xf259a5b3,0x3696693c,0xdd25dacb,0xb34b669a,  // de
            0x2cd279e4,0x4bb5966d,0x96cd35ba,0xa4f3c966,0x6b2cda59,0x9a6b7497,0xe792cd2d,  // tec
            0x59b4b349,0xd6e92ed6,0x259a5b34,0x696693cf,0xd25dacb3,0x34b669ad,0xcd279e4b,  // te
            0xbb5966d2,0x6cd35ba4,0x4f3c9669,0xb2cda59a,0xa6b74976,0x792cd2d9,0x9b4b349e,  // d.
        };
        Parallel.For(0, 8, k =>
        {
            for (int i, j, m0 = 1 + k * 105, m1 = m0 + 104; m0 <= m; m0 += 840, m1 += 840)
            { if (m1 > m) m1 = m; j = 0; i = m0; while (i <= m1) x[i++] = a[j++]; }
        });
    }

    static uint[] countPrimes(uint[] x)
    {
        uint[] c = new uint[8]; byte[] b = new byte[x.Length];
        Parallel.For(0, 8, k =>
        {
            uint xi, ck = 0;
            for (int i = x.Length - 1 - k; i >= 0; i -= 8)
            {
                xi = x[i]; xi -= xi >> 1 & 0x55555555;
                xi = (xi & 0x33333333) + (xi >> 2 & 0x33333333);
                ck += b[i] = (byte)(32 - ((xi + (xi >> 4) & 0xf0f0f0f) * 0x1010101 >> 24));
            }
            c[k] = ck;
        });
        time("c");
        return buildPrimes(b, x, c[0] + c[1] + c[2] + c[3] + c[4] + c[5] + c[6] + c[7]);
    }

    static uint[] buildPrimes(byte[] b, uint[] x, uint n)
    {
        int j = b.Length - 1; byte[] c = new byte[8]; uint[] p = new uint[n];
        sbyte[] bruijn = {01,02,29,03,30,15,25,04,31,23,21,16,26,18,05,09,
                          32,28,14,24,22,20,17,08,27,13,19,07,12,06,11,10};
        for (byte s = 0, k = 1; k < 8 && k <= j; k++) c[k] += s += b[k - 1];
        Parallel.For(0, 8, k =>
        {
            int s, i = k; uint ck = c[k], ci, u, v; long xi;
            if (k > 0 && ck == 0) return;
            for (u = (uint)(k << 6) - 1; ; u += 512)
            {
                v = u; ci = ck; xi = ~x[i];
                while (xi > 0)
                {
                    s = bruijn[((uint)((xi & -xi) * 0x077cb531)) >> 27];
                    p[ci++] = v += (uint)s << 1;
                    xi >>= s;
                }
                i += 8; if (i > j) break;
                ck += (uint)(b[i - 1] + b[i - 2] + b[i - 3] + b[i - 4] +
                             b[i - 5] + b[i - 6] + b[i - 7] + b[i - 8]);
            }
        });
        time("p");
        return p;
    }

    static void time(string s) { Console.WriteLine(sw.ElapsedMilliseconds + " ms " + s); }
}

/*
    static void preMark(int m, uint[] x)  // the old way
    {
        Parallel.For(0, 3, k =>
        {
            if (k == 0) for (int i = 1; i <= m; i += 3) x[i] = 0x24924924;
            if (k == 1) for (int i = 2; i <= m; i += 3) x[i] = 0x49249249;
            if (k == 2) for (int i = 3; i <= m; i += 3) x[i] = 0x92492492;
        });
        Parallel.For(0, 5, k =>
        {
            if (k == 0) for (int i = 1; i <= m; i += 5) x[i] |= 0x42108421;
            if (k == 1) for (int i = 2; i <= m; i += 5) x[i] |= 0x10842108;
            if (k == 2) for (int i = 3; i <= m; i += 5) x[i] |= 0x84210842;
            if (k == 3) for (int i = 4; i <= m; i += 5) x[i] |= 0x21084210;
            if (k == 4) for (int i = 5; i <= m; i += 5) x[i] |= 0x08421084;
        });                                                                
        Parallel.For(0, 7, k =>                                            
        {                                                                  
            if (k == 0) for (int i = 1; i <= m; i += 7) x[i] |= 0x08102040;
            if (k == 1) for (int i = 2; i <= m; i += 7) x[i] |= 0x40810204;
            if (k == 2) for (int i = 3; i <= m; i += 7) x[i] |= 0x04081020;
            if (k == 3) for (int i = 4; i <= m; i += 7) x[i] |= 0x20408102;
            if (k == 4) for (int i = 5; i <= m; i += 7) x[i] |= 0x02040810;
            if (k == 5) for (int i = 6; i <= m; i += 7) x[i] |= 0x10204081;
            if (k == 6) for (int i = 7; i <= m; i += 7) x[i] |= 0x81020408;
        });
    }
*/

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