Pagina's

2011/01/28

Eratosthenes, part 3

// The BitArray became 33% smaller in part 2,
// the speed ~20% faster. With a few changes
// this one is ~33% faster. In the IDE primes
// upto 10^9 are found in 15 s. (12.5 in the 
// release version, (Athlon X4 640, XP, 2GB))

using System;
using System.Collections.Generic;
using System.Collections;
using System.Diagnostics;
class Program
{
    private static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        Console.WriteLine(int.MaxValue);
        int n = 1;
        while (n <= 1000000000)     
        {
            sw.Restart();
            int[] primes = getPrimes(n);
            sw.Stop();
            Console.WriteLine("n: " + n + "  primes: " + primes.Length +
                              " in " + sw.ElapsedMilliseconds + " ms");
            n *= 10;
        }
        Console.ReadLine();
    }
    private static int[] getPrimes(int n)
    {
        if (n < 3)
        {
            if (n == 2) { int[] prime2 = { 2 }; return prime2; }
            else { int[] no_prime = { }; return no_prime; }
        }
        BitArray composite = new BitArray(n / 3);
        int len = composite.Length;
        int d1 = 8, d2 = 0, p1 = 3, p2 = 1, s = 7;
        int i = 0;
        int inc = 10;
        while (s < len)
        {
            if (!composite[i++])
            {
                int k = s, m = s + p1;
                for (; m < len; k += inc, m += inc)
                {
                    composite[k] = true;
                    composite[m] = true;
                }
                for (; k < len; k += inc)
                    composite[k] = true;
            }
            d2 += 8; s += d2; p2 += 8; inc += 4;
            if (!composite[i++] && s < len)
            {
                int k = s, m = s + p2;
                for (; m < len; k += inc, m += inc)
                {
                    composite[k] = true;
                    composite[m] = true;
                }
                for (; k < len; k += inc)
                    composite[k] = true;
            }
            d1 += 16; s += d1; p1 += 4; inc += 8;
        }
        List<int> primes = new List<int>();
        double log_n = Math.Log(n);
        primes.Capacity = (int)((n / log_n) * (1 + 1.2762 / log_n));
        primes.Add(2); primes.Add(3);
        int p = 5;
        i = 0;
        while (p <= n)
        {
            if (!composite[i++]) primes.Add(p);
            p += 2;
            if (p <= n && !composite[i++]) primes.Add(p);
            p += 4;
        }
        return primes.ToArray();
    }
}


No comments:

Post a Comment