Sieve of Eratostheness

sieve.cs

using System;
using System.Collections.Generic;
using System.IO;

unsafe class Program
{
    static void Main()
    {
        var start = DateTime.Now;
        var fs = new FileStream("prime.txt", FileMode.Create, FileAccess.Write, FileShare.Read);
        var w = new BinaryWriter(fs);

        var sieve = stackalloc bool[65536];
        {    // 0 ~ 65536
            ulong* pt = (ulong*)sieve;
            *(pt++) = 0x0100010001010000;
            for (int i = 1; i < 65536 / 8; i++)
                *(pt++) = 0x0100010001000100;

            for (int i = 3; i < 256; i += 2)
                if (sieve[i])
                    for (int j = i * i; j < 65536; j += 2 * i)
                        sieve[j] = false;
        }
        for (int i = 0, j = 0; i < 65536 / 8; i++)
        {
            byte value = 0;
            for (int k = 0; k < 8; j++, k++)
                if (sieve[j])
                    value |= (byte)(1 << k);
            w.Write(value);
        }

        int[] smallPrimes = new int[0], smallPrimeIndexes = new int[0];
        {
            List<int> smallPrimesList = new List<int>(), smallPrimeIndexesList = new List<int>();
            for (int i = 3; i < 65536; i += 2)
                if (sieve[i])
                {
                    smallPrimesList.Add(i);
                    var temp = (i * i) - 65537;
                    smallPrimeIndexesList.Add(temp < 0 ?
                        ((((int)Math.Ceiling(((65537 + i) >> 1) / (double)i) << 1) - 1) * i - 65537) >> 1 :
                        temp >> 1);
                }
            smallPrimes = smallPrimesList.ToArray();
            smallPrimeIndexes = smallPrimeIndexesList.ToArray();
        }

        for (int i = 1; i < 65536 / 2; i++)
        {
            ulong* pt = (ulong*)sieve;
            for (int j = 0; j < 65536 / 8; j++)
                *(pt++) = 0x0101010101010101;

            for (int j = 0; j < smallPrimes.Length; smallPrimeIndexes[j++] -= 65536)
                while (smallPrimeIndexes[j] < 65536)
                {
                    sieve[smallPrimeIndexes[j]] = false;
                    smallPrimeIndexes[j] += smallPrimes[j];
                }

            for (int j = 0, k = 0; j < 65536 / 4; j++)
            {
                byte value = 0;
                for (int l = 0; l < 4; k++, l++)
                    if (sieve[k])
                        value |= (byte)(1 << (l * 2 + 1));
                w.Write(value);
            }

            if ((i & 511) == 0) Console.WriteLine(i);    // for debug purpose
        }

        {    // for 4294901759 ~ 4294967295
            ulong* pt = (ulong*)sieve;
            for (int j = 0; j < 65536 / 2 / 8; j++)
                *(pt++) = 0x0101010101010101;

            for (int j = 0; j < smallPrimes.Length; j++)
                while (smallPrimeIndexes[j] < 65536 / 2)
                {
                    sieve[smallPrimeIndexes[j]] = false;
                    smallPrimeIndexes[j] += smallPrimes[j];
                }

            for (int j = 0, k = 0; j < 65536 / 2 / 4; j++)
            {
                byte value = 0;
                for (int l = 0; l < 4; k++, l++)
                    if (sieve[k])
                        value |= (byte)(1 << (l * 2 + 1));
                w.Write(value);
            }
        }

        w.Close();
        var end = DateTime.Now;
        Console.WriteLine(end - start);
    }
}