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

No comments:

Post a Comment