Pagina's

2017/09/25

Direct Search Factorization


// MathWorld: DSF is the simplest (and most simple-minded) prime factorization algorithm.
// It consists of searching for factors of a number by systemetically performing trial divisions,
// usually using a sequence of increasing numbers. Multiples of small primes are commonly excluded
// to reduce the number of trial divisors, but just including them is sometimes faster than the time
// required to exclude them. DSF is very inefficient, and can be used only with fairly small numbers.

using System;                             // i7-4790@3.6GHz
using System.Collections.Generic;
class Program
{
    static void Main()
    {
        findPF(~0uL);                     // warm up
        time(4294967291);                 // largest prime < 2^32          50 µs
        time(4294967291 * 4294967291uL);  //                       6.5 s  
        time(3183958073 * 5793651691uL);  // 2 primes              4.8 s
        time(~0uL - 58);                  // largest prime < 2^64  6.5 s
        time(~0uL);                       // (2^64)-1                     100 µs 
        Console.Read();
    }

    static void time(ulong u)
    {
        Console.Write(" n " + u + "\n" + " t ");
        var sw = System.Diagnostics.Stopwatch.StartNew();
        List<ulong> pf = findPF(u);
        Console.Write(sw.Elapsed + "\n" + " f ");
        foreach (var v in pf) Console.Write(v + " ");
        Console.Write("\n" + "\n");
    }

    static List<ulong> findPF(ulong n)
    {
        var pf = new List<ulong>();
        if (n < 4) { if (n > 0) pf.Add(n); return pf; }
        while (n % 2 == 0) { n /= 2; pf.Add(2); }
        while (n % 3 == 0) { n /= 3; pf.Add(3); }
        while (n % 5 == 0) { n /= 5; pf.Add(5); }
        if (n == 1) return pf;
        ulong d = 1; uint b = 0, rt = sqrt(n);
        if (n <= ~0u) { findPF((uint)n, rt, 1, 0, pf); return pf; }
        for (; ; b = 0x5A28A6)
        {
            while (b > 0)
            {
                d += b & 7; b >>= 3; if (d > rt) { pf.Add(n); return pf; }
                if (n % d == 0)
                {
                    n /= d; pf.Add(d); while (n % d == 0) { n /= d; pf.Add(d); }
                    if (n == 1) return pf;
                    rt = sqrt(n);
                    if (n <= ~0u) { findPF((uint)n, rt, (uint)d, b, pf); return pf; }
                }
            }
            while (d + 30 <= rt && n % (d + 06) > 0 && n % (d + 10) > 0
                                && n % (d + 12) > 0 && n % (d + 16) > 0
                                && n % (d + 18) > 0 && n % (d + 22) > 0
                                && n % (d + 28) > 0 && n % (d + 30) > 0) d += 30;
        }
    }

    static void findPF(uint n, uint rt, uint d, uint b, List<ulong> pf)
    {
        for (; ; b = 0x5A28A6)
        {
            while (b > 0)
            {
                d += b & 7; b >>= 3; if (d > rt) { pf.Add(n); return; }
                if (n % d == 0)
                {
                    n /= d; pf.Add(d); while (n % d == 0) { n /= d; pf.Add(d); }
                    if (n == 1) return;
                    rt = (uint)Math.Sqrt(n);
                }
            }
            while (d + 30 <= rt && n % (d + 06) > 0 && n % (d + 10) > 0
                                && n % (d + 12) > 0 && n % (d + 16) > 0
                                && n % (d + 18) > 0 && n % (d + 22) > 0
                                && n % (d + 28) > 0 && n % (d + 30) > 0) d += 30;
        }
    }

    static uint sqrt(ulong n0)
    {
        if (n0 < 1uL << 52) return (uint)Math.Sqrt(n0);
        ulong n1 = n0 - 1 >> 24; int s = 25;
        if (n1 > 255) { n1 >>= 4; s = 29; }
        if (n1 > 15) { n1 >>= 2; s |= 2; }
        if (n1 > 3) s++;
        ulong r0 = 1uL << s, r1 = r0 + (n0 >> s) >> 1;
        while (r1 < r0) { r0 = r1; r1 = r0 + n0 / r0 >> 1; }
        return (uint)r0;
    }
}



///////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                               //
//                                  Eight threads, four cores.                                   //          
//                                                                                               //
///////////////////////////////////////////////////////////////////////////////////////////////////

using System;
using System.Collections.Generic;
using System.Threading.Tasks;
class Program
{
    static void Main()
    {
        findPF(~0uL);                     // warm up
        time(4294967291);                 // largest prime < 2^32          90 µs
        time(4294967291 * 4294967291uL);  //                              110 µs 
        time(3183958073 * 5793651691uL);  // 2 primes              0.9 s
        time(~0uL - 58);                  // largest prime < 2^64  1.1 s
        time(~0uL);                       // (2^64)-1                      80 µs 
        Console.Read();
    }

    static void time(ulong u)
    {
        Console.Write(" n " + u + "\n" + " t ");
        var sw = System.Diagnostics.Stopwatch.StartNew();
        List<ulong> pf = findPF(u);
        Console.Write(sw.Elapsed + "\n" + " f ");
        foreach (var v in pf) Console.Write(v + " ");
        Console.Write("\n" + "\n");
    }

    static List<ulong> findPF(ulong n)
    {
        var pf = new List<ulong>();
        if (n < 4) { if (n > 0) pf.Add(n); return pf; }
        while (n % 2 == 0) { n /= 2; pf.Add(2); }
        while (n % 3 == 0) { n /= 3; pf.Add(3); }
        while (n % 5 == 0) { n /= 5; pf.Add(5); }
        if (n > 1) { recFindF(n, pf); pf.Sort(); } return pf;
    }

    static void recFindF(ulong n, List<ulong> pf)
    {
        ulong f = findF(n); if (f == n) pf.Add(n);
        else { recFindF(f, pf); recFindF(n / f, pf); }
    }

    static ulong findF(ulong n)
    {
        uint r = sqrt(n); if (n % r == 0) return r;
        var locker = new object();
        Parallel.For(0, 8, k =>
        {
            ulong d0 = k == 0 ? 07u : k == 1 ? 11u : k == 2 ? 13u : k == 3 ? 17u :
                       k == 4 ? 19u : k == 5 ? 23u : k == 6 ? 29u : 31;
            ulong d1 = d0 + 30, d2 = d0 + 60, d3 = d0 + 90;
            while (d0 <= r)
            {
                if (n % d0 == 0 || n % d1 == 0 || n % d2 == 0 || n % d3 == 0)
                    lock (locker)
                    {
                        if (r == 0) return; r = 0;
                        n = n % d0 == 0 ? d0 : n % d1 == 0 ? d1 : n % d2 == 0 ? d2 : d3;
                        return;
                    }
                d0 += 120; d1 += 120; d2 += 120; d3 += 120;
            }
        }); return n;
    }

    static uint sqrt(ulong n0)
    {
        if (n0 < 1uL << 52) return (uint)Math.Sqrt(n0);
        ulong n1 = n0 - 1 >> 24; int s = 25;
        if (n1 > 255) { n1 >>= 4; s = 29; }
        if (n1 > 15) { n1 >>= 2; s |= 2; }
        if (n1 > 3) s++;
        ulong r0 = 1uL << s, r1 = r0 + (n0 >> s) >> 1;
        while (r1 < r0) { r0 = r1; r1 = r0 + n0 / r0 >> 1; }
        return (uint)r0;
    }
}


///////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                               //
//                                               ?                                               //
//                                                                                               //
///////////////////////////////////////////////////////////////////////////////////////////////////

using System;
using System.Collections.Generic;
using System.Threading.Tasks;
class Program
{
    static void Main()
    {
        sqrt(1uL << 52);
        ulong[] n = new ulong[] { ~0u - 4, 0, ~0uL - 172, ~0uL - 58, ~0uL };
        n[1] = n[0] * n[0];
        findPF(n[0]); time(n[0]);  // largest prime < 2^32          50 µs
        findPF(n[1]); time(n[1]);  //                               80 µs 
        findPF(n[2]); time(n[2]);  // 2 primes              0.78 s
        findPF(n[3]); time(n[3]);  // largest prime < 2^64  1.05 s
        findPF(n[4]); time(n[4]);  // (2^64)-1                      70 µs 
        Console.Read();
    }

    static void time(ulong u)
    {
        Console.Write(" n " + u + "\n" + " t ");
        var sw = System.Diagnostics.Stopwatch.StartNew();
        List<ulong> pf = findPF(u);
        Console.Write(sw.Elapsed + "\n" + " f ");
        foreach (var v in pf) Console.Write(v + " ");
        Console.Write("\n" + "\n");
    }

    static List<ulong> findPF(ulong n)
    {
        var pf = new List<ulong>();
        if (n < 4) { if (n > 0) pf.Add(n); return pf; }
        while (n % 2 == 0) { n /= 2; pf.Add(2); }
        while (n % 3 == 0) { n /= 3; pf.Add(3); }
        while (n % 5 == 0) { n /= 5; pf.Add(5); }
        if (n > 1)
            if (n <= ~0u) findPF((uint)n, pf);
            else { recFindF(n, pf, 1); pf.Sort(); }
        return pf;
    }

    static void findPF(uint n, List<ulong> pf)
    {
        for (uint d = 1, r = (uint)Math.Sqrt(n), b = 0; ; b = 0x5A28A6)
        {
            while (b > 0)
            {
                d += b & 7; b >>= 3; if (d > r) { pf.Add(n); return; }
                if (n % d == 0)
                {
                    n /= d; pf.Add(d); while (n % d == 0) { n /= d; pf.Add(d); }
                    if (n == 1) return;
                    r = (uint)Math.Sqrt(n);
                }
            }
            while (d + 30 <= r && n % (d + 06) > 0 && n % (d + 10) > 0
                               && n % (d + 12) > 0 && n % (d + 16) > 0
                               && n % (d + 18) > 0 && n % (d + 22) > 0
                               && n % (d + 28) > 0 && n % (d + 30) > 0) d += 30;
        }
    }

    static void recFindF(ulong n, List<ulong> pf, int i)
    {
        ulong f = findF(n);
        if (f == n) { while (i > 0) { pf.Add(n); i--; } }
        else if (n / f == f) { recFindF(f, pf, i << 1); }
        else { recFindF(f, pf, i); recFindF(n / f, pf, i); }
    }

    static ulong findF(ulong n0)
    {
        if (n0 <= ~0u) return findF((uint)n0);
        uint r = sqrt(n0); if (n0 % r == 0) return r;
        var locker = new object();
        Parallel.For(0, 8, k =>
        {
            ulong d0 = k == 0 ? 07u : k == 1 ? 11u : k == 2 ? 13u : k == 3 ? 17u :
                       k == 4 ? 19u : k == 5 ? 23u : k == 6 ? 29u : 31,
                  d1 = d0 + 30, d2 = d0 + 60, d3 = d0 + 90, n1 = n0;
            for (; d0 <= r; d0 += 120, d1 += 120, d2 += 120, d3 += 120)
                if (n1 % d0 == 0 || n1 % d1 == 0 || n1 % d2 == 0 || n1 % d3 == 0)
                    lock (locker)
                    {
                        if (r == 0) return; r = 0;
                        n0 = n0 % d0 == 0 ? d0 : n0 % d1 == 0 ? d1 : n0 % d2 == 0 ? d2 : d3;
                        return;
                    }
        }); return n0;
    }

    static uint findF(uint n0)
    {
        uint r = (uint)Math.Sqrt(n0); if (n0 % r == 0) return r;
        var locker = new object();
        Parallel.For(0, 8, k =>
        {
            uint d0 = k == 0 ? 07u : k == 1 ? 11u : k == 2 ? 13u : k == 3 ? 17u :
                      k == 4 ? 19u : k == 5 ? 23u : k == 6 ? 29u : 31,
                 d1 = d0 + 30, d2 = d0 + 60, d3 = d0 + 90, n1 = n0;
            for (; d0 <= r; d0 += 120, d1 += 120, d2 += 120, d3 += 120)
                if (n1 % d0 == 0 || n1 % d1 == 0 || n1 % d2 == 0 || n1 % d3 == 0)
                    lock (locker)
                    {
                        if (r == 0) return; r = 0;
                        n0 = n0 % d0 == 0 ? d0 : n0 % d1 == 0 ? d1 : n0 % d2 == 0 ? d2 : d3;
                        return;
                    }
        }); return n0;
    }

    static uint sqrt(ulong n0)
    {
        if (n0 < 1uL << 52) return (uint)Math.Sqrt(n0);
        uint n1 = (uint)(n0 >> 48); int s = 25);
        if (n1 > 255u) { n1 >>= 8; s = 29; }
        if (n1 > 15u) { n1 >>= 4; s |= 2; }
        if (n1 > 3u) s++;
        ulong r0 = 1uL << s, r1 = r0 + (n0 >> s) >> 1;
        while (r1 < r0) { r0 = r1; r1 = r0 + n0 / r0 >> 1; }
        return (uint)r0;
    }
}

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