The two main composite counting functions I used I called S and dS, where I probably was thinking about "sum" as S is the sum of composites.

The main workhorse is the dS function though which is where I had to do the most work in figuring out its form:

dS(x,p

_{j}) = [x/p

_{j}] - 1 - (j-1) - S(x/p

_{j}, p

_{j}-1)

where S(x/p

_{j}, p

_{j}-1) is the count of composites that multiply times p

_{j}to give a product less than or equal to x, where notice that p

_{j}must be less than or equal to sqrt(x) or the composite count given by [x/p

_{j}] - 1 will not be correct.

So the dS function is the count of composites for a particular prime excluding composites that are products of lesser primes. So it gets a count of integers with a prime as a factor, subtracts 1 for the prime itself, and then subtracts the count of primes less than it. And finally it subtracts the count of composites multiplied times that prime. Whew!

Here's an example:

dS(16,3) = [16/3] - 1 - 1 - S(16/3,2) = 5 - 1 - 1 - 1 = 2

There are 5 numbers with 3 as a factor below 16, but you subtract 1 for 3, then 1 for 2, to handle 2(3), and then 1 more for the count of composites times 3, where that is for 12.

Here S(16/3,2) = S(5,2) = 1, because 2(2)(3) is the only composite case where 3 is multiplied times a composite up to 16. I don't typically write the floor() inside the S function as I consider it redundant.

The dS function only wanted composites divisible by 3, but not by a lesser prime, so it tossed off the evens, meaning 6 and 12 were not counted! So no need to correct the count.

Now let's do dS(16,2) as a second example:

dS(16,2) = [16/2] - 1 - 0 - S(16/2,1) = 8 -1 - 0 - 0 = 7

Which is where 6 and 12 get counted along with the other evens. And notice if S(x,1) you can't have any composites made up of smaller primes as 2 is the smallest prime, so S(x,1) = 0.

And that's it. The S function then is just a complete count of composites, where of course, if you have the complete count of composites you can count primes.

So with my examples I can now count primes up to 16:

16 - 1 - dS(16,2) - dS(16,3) = 16 - 1 - 7 - 2 = 6

What's weird though is that the simple shift of making dS a function that counts for each prime excluding counts of lesser primes leads to a very compact mathematical form for a prime counting function itself.

I call the prime counting function P for obvious reason:

P(x,p

_{j}) = [x] - S(x,p

_{j}) - 1

which allows me to simplify dS(x, p

_{j}) to

dS(x,p

_{j}) = P(x/p

_{j},p

_{j}-1) - (j-1)

where again p

_{j}has to be less than or equal to sqrt(x).

But notice that P(x,p

_{j}) is the full count of primes up to and including x, if p

_{j}is greater than or equal to the last prime less than sqrt(x), so I can remove all mention of primes with the following

dS(x,y) = (P(x/y,y-1) - P(y-1, sqrt(y-1)))(P(y, sqrt(y)) - P(y-1, sqrt(y-1)))

as if y is not prime then

P(y, sqrt(y)) - P(y-1, sqrt(y-1)) = 0

with the constraint that if y>sqrt(x), then P(x,y) = P(x,sqrt(x)) to keep the correct count.

And of course now I have a difference equation.

One thing of interest then is that the P function acts as a boolean switch, telling the dS function when y is prime.

For instance 6 is not prime and the count of primes up to 6 is 3, which is the same as the count up to 5, which is the third prime, so the switch will give 0 for y = 6.

Notice here that sqrt(x) is the integer square root, so it means to find the largest integer less than or equal to sqrt(x), for instance, the integer square root of 10 is 3.

I also now have

P(x,y) = [x] - S(x,y) - 1

as I can remove mention of primes themselves! So while I started with a focus on primes, determining at each prime, thanks to my boolean switch I can now remove all mention of them.

Now the count of composites for each prime less than or equal to the square root of x when added together just gives the full count of composites, so I have that the S function equals the sum of the dS functions for each prime less than or equal to sqrt(x), where I define S(x,1) = 0, and my prime counting function is complete.

Here's the entire thing.

P(x,y) = [x] - 1 - sum for j=2 to y of {(P([x/j],j-1) - P(j-1, sqrt(j-1)))(P(j, sqrt(j)) - P(j-1, sqrt(j-1)))}

where if y>sqrt(x), then P(x,y) = P(x,sqrt(x)).

It's actually prettier with the summation sign. But Blogger doesn't make that easy for me to do, so I've just used the format above, which is of course correct, just not as pretty looking.

James Harris