number of primes pi(x)

/*  
    Adrien-Marie Legendre (1752-1833) observed that pi(x), the number of primes p <= x , 
    can be calculated with the inclusion-exclusion principle.

    pi(x) - pi(sqrt(x)) + 1 = x - SIGMA(x/pi) + SIGMA(x/pi/pj) - SIGMA(x/pi/pj/pk) + ...
                                    i            i<j             i<j<k

        x = 26 , sqrt(26) = 5 , primes <= 5 : 2 3 5 , pi(5) = 3
        pi(26) - 3 + 1 = 26 - 26/2 - 26/3 - 26/5 + 26/2/3 + 26/2/5 + 26/3/5 - 26/2/3/5
                       = 26 - 13   - 8    - 5    + 4      + 2      + 1      - 0
                       = 26 - 26 + 7 
        pi(26) = 9 
  
               x         pi(x)  minutes   seconds                          (Athlon 3GHz)
            1e09      50847534               2.01  getPrimes(1e9): 12.5 s
            1e10     455052511              20.90
            1e11    4118054813        4    216.00
            1e12   37607912018       38   2260.00
            1e13  346065536839      399  23900.00
  
    Lagarias, Miller & Odlyzko: pi(1e13) 8 minutes (1985 on an IBM/370 3081 Model K).
*/
using System;
using System.Collections;
using System.Collections.Generic;
using System.Diagnostics;
class pi_x
{
    static Stopwatch sw = new Stopwatch();
    static void Main()
    {
        //test();
        sw.Start();
        ulong r = pi((ulong)1e9);
        sw.Stop();
        Console.Write(sw.Elapsed + " " + r);
        Console.Read();
    }

    /*
    static uint _s32;
    static uint pi32(uint x)
    {
        _p = getPrimes((uint)Math.Sqrt(x));
        _pLen = _p.Length;
        _s32 = 0;
        sxnp32(x, 0, 0);
        return _s32 + (uint)_pLen - 1;
    }
    static void sxnp32(uint q, int n, int pIdx) // sum x divided by n primes
    {
        if ((n & 1) == 0)
            _s32 += q;
        else
            _s32 -= q;
        for (int i = pIdx; i < _pLen; i++)
            if (q < _p[i]) 
                break;
            else
                sxnp32(q / _p[i], n + 1, i + 1);
    }
    */

    static uint[] _p;
    static int _pLen;
    static ulong _s;

    static ulong pi(ulong x)
    {
        if (x == 0) return 0;
        _p = getPrimes((uint)Math.Sqrt(x));
        _pLen = _p.Length;
        _s = 0;
        if (x <= ~0u)
            sxnp((uint)x, 0, 0);
        else
            sxnp(x, 0, 0);
        return _s + (uint)_pLen - 1;
    }

    static void sxnp(uint q, int n, int pIdx)
    {
        if ((n & 1) == 0)
        {
            _s += q;
            for (int i = pIdx; i < _pLen; i++)
                if (q < _p[i]) break;
                else
                    if (q >> 1 < _p[i]) _s--;
                    else sxnp(q / _p[i], n + 1, i + 1);
        }
        else
        {
            _s -= q;
            for (int i = pIdx; i < _pLen; i++)
                if (q < _p[i]) break;
                else
                    if (q >> 1 < _p[i]) _s++;
                    else sxnp(q / _p[i], n + 1, i + 1);
        }
    }

    static void sxnp(ulong q, int n, int pIdx)
    {
        if ((n & 1) == 0)
        {
            _s += q;
            for (int i = pIdx; i < _pLen; i++)
            {
                if (q < _p[i]) break;
                if (q >> 1 < _p[i]) _s--;
                else
                    if (q > ~0u) sxnp(q / _p[i], n + 1, i + 1);
                    else sxnp((uint)q / _p[i], n + 1, i + 1);
            }
        }
        else
        {
            _s -= q;
            for (int i = pIdx; i < _pLen; i++)
            {
                if (q < _p[i]) break;
                if (q >> 1 < _p[i]) _s++;
                else
                    if (q > ~0u) sxnp(q / _p[i], n + 1, i + 1);
                    else sxnp((uint)q / _p[i], n + 1, i + 1);
            }
        }
    }

    static uint[] getPrimes(uint n) // primes <= n (Eratosthenes)   
    {
        int y, d, e, q, r, s, i, j;
        double logN; List<uint> p; BitArray x; uint u;
        if (n < 3) return n == 2 ? new uint[] { 2 } : new uint[] { };
        logN = Math.Log(n);
        p = new List<uint>((int)(n / logN * (1 + 1.2762 / logN)));
        x = new BitArray((int)(n / 3));
        y = x.Length; d = 8; e = 0; q = 3; r = 1; s = 7; i = 0; j = 10;
        while (s < y)
        {
            if (!x[i++])
            {
                int k = s, m = s + q;
                for (; m < y; k += j, m += j) { x[k] = true; x[m] = true; }
                for (; k < y; k += j) x[k] = true;
            }
            e += 8; s += e; r += 8; j += 4;
            if (!x[i++] && s < y)
            {
                int k = s, m = s + r;
                for (; m < y; k += j, m += j) { x[k] = true; x[m] = true; }
                for (; k < y; k += j) x[k] = true;
            }
            d += 16; s += d; q += 4; j += 8;
        }
        p.Add(2); p.Add(3); i = 0; y -= 2; u = 5;
        while (i < y)
        {
            if (!x[i++]) p.Add(u); u += 2;
            if (!x[i++]) p.Add(u); u += 4;
        }
        if (u <= n && !x[i++]) p.Add(u); u += 2;
        if (u <= n && !x[i]) p.Add(u);
        return p.ToArray();
    }

    static void test()
    {
        for (uint x = 0; x <= 50; x++)
            Console.WriteLine("{0,2} {1,2}", x, pi(x));
    }
}