a practical prime counting function

Question

Looking at Prime counting functions on Wikipedia, I only found formulas with no hint on how people got there. So, to better understand, I've decided to build one from scratch, starting from a naive kind of Sieves of Eratosthenes. Here is what I've found so far: $$ \text{Let } k = \left\lfloor\frac {x - 1} {6} \right\rfloor \text{ and a function } p1(x) \text{ such that: } $$ $$ p1(x)=k-\sum _{i=1}^{6i^2-2i\leqslant k} \left( \left\lfloor \frac{k+i}{6i-1}\right\rfloor -i+1 \right) - \sum _{i=1}^{6i^2+2i\leqslant k} \left( \left\lfloor \frac{k-i}{6i+1}\right\rfloor -i+1 \right) + \sum _{i=0}^{i\leqslant \frac{k-29}{25}} \left( \left\lfloor \frac{k+5i+6}{30i+35}\right\rfloor -1 \right) + \sum _{i=0}^{i\leqslant \frac{k-29}{35}} \left\lfloor \frac{k-5i-4}{30i+25}\right\rfloor + \sum _{i=0}^{i\leqslant \frac{k-64}{35}} \left( \left\lfloor \frac{k+7i+13}{42i+77}\right\rfloor -1 \right) + \sum _{i=0}^{i\leqslant \frac{k-57}{49}} \left( \left\lfloor \frac{k-7i-8}{42i+49}\right\rfloor -1 \right) $$ $$ \text{and let } k=\left\lfloor \frac{x+1}{6}\right\rfloor \text{ and a function } p2(x) \text{such that}: $$ $$ p2(x)=k-\sum _{i=1}^{i\leqslant \frac{k+1}{7}} \left\lfloor \frac{k-i}{6i-1}\right\rfloor + \sum _{i=0}^{i\leqslant \frac{k-21}{25}} \left( \left\lfloor \frac{k+5i+4}{30i+25}\right\rfloor -1 \right) + \sum _{i=0}^{i\leqslant \frac{k-41}{35}} \left\lfloor \frac{k-5i-6}{30i+35}\right\rfloor + \sum _{i=0}^{i\leqslant \frac{k-41}{35}} \left( \left\lfloor \frac{k+7i+8}{42i+49}\right\rfloor -1 \right) $$ $$ \text{Then for any } x\leqslant 874, \text{we have:} $$ $$ 2+p1(x)+p2(x)= \pi (x) $$ Although the result looks like complicated, this is only very simple sums and programming an algorithm that matches it is really straightforward. Here for example is the corresponding Mathematica code:

Manipulate[
 Module[{panel, p1, p2, pcount},

  p1[x_] :=
   Module[{k, i, sum, sum1, sum2, sum3, sum4, sum5, sum6},
    k = Floor[(x - 1)/6];
    sum1 = (i = 1; sum = 0; 
      While[6*i^2 - 2*i <= k, sum += Floor[(k + i)/(6*i - 1)] - i + 1;
        i++]; sum);
    sum2 = (i = 1; sum = 0; 
      While[6*i^2 + 2*i <= k, sum += Floor[(k - i)/(6*i + 1)] - i + 1;
        i++]; sum);
    sum3 = (i = 0; sum = 0; 
      While[i <= (k - 29)/25, 
       sum += Floor[(k + 5*i + 6)/(30*i + 35)] - 1; i++]; sum);
    sum4 = (i = 0; sum = 0; 
      While[i <= (k - 29)/35, sum += Floor[(k - 5*i - 4)/(30*i + 25)];
        i++]; sum);
    sum5 = (i = 0; sum = 0; 
      While[i <= (k - 64)/35, 
       sum += Floor[(k + 7*i + 13)/(42*i + 77)] - 1; i++]; sum);
    sum6 = (i = 0; sum = 0; 
      While[i <= (k - 57)/49, 
       sum += Floor[(k - 7*i - 8)/(42*i + 49)] - 1; i++]; sum);
    k - sum1 - sum2 + sum3 + sum4 + sum5 + sum6
    ];

  p2[x_] :=
   Module[{km, i, sum, sum1, sum2, sum3, sum4, sum5, sum6},
    km = Floor[(x + 1)/6];
    sum1 = (i = 1; sum = 0; 
      While[i <= (km + 1)/7, sum += Floor[(km - i)/(6*i - 1)]; i++]; 
      sum);
    sum2 = (i = 0; sum = 0; 
      While[i <= (km - 21)/25, 
       sum += Floor[(km + 5*i + 4)/(30*i + 25)] - 1; i++]; sum);
    sum3 = (i = 0; sum = 0; 
      While[i <= (km - 41)/35, 
       sum += Floor[(km - 5*i - 6)/(30*i + 35)]; i++]; sum);
    sum4 = (i = 0; sum = 0; 
      While[i <= (km - 41)/35, 
       sum += Floor[(km + 7*i + 8)/(42*i + 49)] - 1; i++]; sum);
    km - sum1 + sum2 + sum3 + sum4
    ];

  pcount[x_] := 2 + p1[x] + p2[x];

  panel =
   Panel@Column[{
      Grid[{
        {"x", x}, {"\[Pi](x)", PrimePi[x]}, {"pcount(x)", pcount[x]}
        }, Alignment -> {Left, Top}, Frame -> All]
      }];

  panel
  ]
 , {{x, 301}, 1, 1500, 1, Appearance -> "Open"}
 ]

I've found interesting the fact that you don't need to rely on Primes table nor need of trial division.

Although the range $x\leqslant 874$ (in fact, it also works for range $1001\leqslant x\leqslant 1224$) is quite small, I think that the formula can be extended to any range by adding necessary sums, but it then becomes difficult to avoid duplicates counts.

EDIT1: some people misunderstood "by adding necessary sums". I mean in fact that I'm quite sure that $p1(x)$ and $p2(x)$ valid for any $x$ can be written in the form of a closed single double-sum of $k$ (i.e: $pn(x)=\sum \sum f(k)$ ).

I didn't find anything on the net about this or similar kind of formula. Does anybody know something about?

Answer 1

I don't mean to be a pessimist, but the computational complexity of this "formula" is most likely greater than or equal to the complexity of just running the sieve itself. However, there are actually "formulas" for primes, they just have enormous computing complexity. For instance, look here. They also discuss the above issues.

I use "formulas" in quotation because a formula is supposed to vastly simplify the process of finding a solution to a problem. For instance, integration formulas vastly simplify the process of adding an (un)countable amount of things and indefinite amount of times! The formulas mentioned don't do that, so they generally aren't considered formulas.

Basically what you have, translates an algorithm into mathematical formalism.

Answer 2

Quite a few people have sought formulas for primes before; this effort falls into that body of work. You can convert this into a prime test via pcount(x+1)-pcount(x).

Unfortunately, this algorithm does not appear to be simpler than a direct sieve, and does not appear to be more compact than a list of primes (e.g. the storage required to store the algorithm is approximately equal to that required to simply store a list of all primes up to 874). It is not possible to have a single polynomial that yields all primes (this was proven by Legendre). If you allow floor functions (as the OP does), then a single formula suffices (this was proven by Mills).

If the formula were computer-generated, i.e. if 874 were a parameter, which you could turn to 2000 or whatever you like and automatically get the corresponding formula, then THAT part might be interesting. Unfortunately the actual formula presented is of limited interest itself.