Thursday, October 18, 2012

My approach to counting primes

I like emphasizing that I'm NOT a mathematician as it keeps things simple, while noting my background is in physics, as yes, I did get a lot of mathematical training. But it was mostly in the complex field, and worrying maybe about differential equations, NOT Diophantine equations. Not much call for number theory in the sciences.

But college was over two decades ago, and I'm not even a physicist as I only got my undergraduate degree in physics, but still like challenges, so I found myself doing number theory out of intellectual curiosity.

Which means I've learned most of the number theory I know, late. Which can help explain how I came up with my own way to count prime numbers.

You see for a dedicated math student, who knew a lot of number theory, it'd probably be something they learned early, where counting primes actually involves very easy concepts! So if you can do division and know what a prime number is you should be able to follow along.

For instance, consider the primes up to 16. To count those, what if you go to composites instead? Like 8 of those numbers are even: 2, 4, 6, 8, 10, 12, 14 ,16

And 5 of them are divisible by three: 3, 6, 9, 12, 15

And notice that 16/2 = 8, and 16/3 = 5 1/3, and we don't need that 1/3 to count composites! So we can just throw it away with something called the floor() function, and floor(16/3) = 5. Because typing "floor" a lot can get tedious, another way to do it is to use brackets, so [16/3] = 5.

But we've counted the primes too!

So the even composites up to 16 are really: 8 - 1 = 7. And similarly the composites divisible by 3 are 4, as we need to subtract one for the prime itself in each case.

But notice we have 6 and 12 in both lists! So there is a double-count so we need to subtract 2 back out to correct for counting them twice.

So if we say 7 evens + 4 divisible by three - 2 evens divisible by three, we get: 9 composites

And those 9 composites are: 4, 6, 8, 9, 10, 12, 14, 15, 16

And there is one more number to handle which is 1, itself, as it's not considered prime, and we can now count primes as 16 - 9 - 1 = 6.

And the six primes are: 2, 3, 5, 7, 11, 13

Notice we only needed 2 and 3 to count the primes with this technique, as the first composite that does not have 2 or 3 as a factor is 25. As it is 5(5) and that gives a rule that you only need the primes up to the integer square root. So you only need 2, 3, 5 and 7 using the technique above to count primes up to 100, but also up to 120, as 121 = 11(11) is the start of needing the next prime.

And all of the above was discovered centuries ago, and is the way humanity approached counting primes. Others built on the basic technique figuring out ways to make it faster, as the above gets tedious when you're adding and subtracting for all the corrections. And over time mathematicians built very fast ways to count prime numbers.

And I knew none of the above back in the summer of 2002, when I started thinking about counting primes.

So I came up with my own way.

Hold on to your seats. All of the above is what was known to humanity for centuries.

Now you'll learn my way.

I'm drawing heavily on my post from June 2005 on this blog, as it turns out I don't explain it a lot, unlike many of my other results, and I'm glad I did explain it at least once! Or I might have forgotten how I did it. Turns out I got really bored counting primes, but I digress.

Here we go.

My start is a simple idea: consider a function that gives the count of composites for a particular prime that have that prime as a factor, but no lesser primes as factors.

Why did I start with a function? I don't know, maybe it was my physics background sneaking through.

For instance, up to 16, for 3 that function would give a count of 2, as 9 and 15 are the composites that have 3 as a factor and no lesser primes, while 6 and 12 would be excluded as they have 2 as a factor, and 2 is less than 3.

Now call that function dS(x,pj) where x is where you're counting up to, and pj is the j_th prime, so my earlier result is dS(16,3) = 2.

(Us physics guys like cool looking functions.)

Now generalizing a bit, to get the count of composites that have pj as a factor up to and including some x, it's now time to use that floor() function, which again means to throw away the remainder, as I generalize a bit:

floor(x/pj) - 1

where from now on, I'll use [x] = floor(x), so I have

[x/pj] - 1

for the count of composites with pj as a factor.

And notice that's the idea from before, so it still emerges, where we want the count of composites for a particular prime.

From that count I need to subtract the number of primes less than pj, which is j-1, so I have:

[x/pj] - 1 - (j-1).

The reason is that for each of those smaller primes I know I'm multiplying times pj.

For instance, for 16 and 3, j = 2, as 3 is the second prime, and j-1 = 1, where that composite is 6, as 2(3) is the one composite.

Now I need the count of composites up to and including x that have primes less than pj as a factor, which is another function.

That is, I now need a function that is the count of composites up to and including some x for all primes less than pj, and I'll call it S(x,pj-1), as the count of composites up to and including x that have primes less than pj as factors.

Here's where things are a little tricky as I now have my full dS(x,pj) function as

dS(x,pj) = [x/pj] - 1 - (j-1) - S(x/pj, pj-1)

where S(x/pj, pj-1) is the count of composites that multiply times pj to give a product less than or equal to x, where notice that pj must be less than or equal to sqrt(x) or the composite count given by [x/pj] - 1 will not be correct.

So sticking in numbers, you have:

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

And again those two composites are 9 and 15. So yeah! Our functions worked!

And should I have used the floor() function inside S? Looks ok to me--study the definition of S to see why--and actually one trick is to just declare everything to be natural numbers so it's understood without using floor() all over the place.

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.

Now let's do dS(16,2).

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 now I can count primes up to 16:

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

Which of course is the same answer as before, as the prime count can't change. And remember 1 is being subtracted for the number 1, which is not a prime number. Then the evens are subtracted which is dS(16,2) followed by composites divisible by 3, which are NOT even, which is dS(16,3).

My approach simply removes the correction piece because 6 and 12 are not double-counted.

Ok, it's time to get more mathematical.

Now I know that in general the count of primes up to and including x is equal to x minus the count of composites minus 1, and that gives an idea for one more function, which I'll call P(x,pj), where

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

which allows me to simplify dS(x, pj) to

dS(x,pj) = P(x/pj,pj-1) - (j-1)

where again pj has to be less than or equal to sqrt(x).

But notice that P(x,pj) is the full count of primes up to and including x, if pj 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.

That is, if a number y is NOT prime then the prime count will not move. 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 our clever function can call itself to see if a number is prime or not!

And don't worry, we didn't wander away from integers! That sqrt(x) here 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.

It's still number theory, and integers rule.

I also now have

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

and notice that dS(x,y) is the count of composites that have y as a factor that do not have primes less than y as a factor, which means that if y is composite dS(x,y) = 0, which is how I can remove mention of specific primes.

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)).

And notice, no need to mention primes at all! And I have my own way to count prime numbers.

And I guess there is no way I'd have figured that thing out if I'd been taught the traditional way first.

Turns out you can connect it back to the approaches used traditionally by mathematicians only with a form where you DO keep in primes, which looks simpler.

P(x,n) = [x] - 1 - sum for j=1 to n of {P([x/p_j],j-1) - (j-1)}

where if n is greater than the count of primes up to and including sqrt(x) then n is reset to that count.

Which is really nice and compact. I do like it.

And that form is faster if you actually program it. The original fully mathematicized form keeps figuring out numbers are prime, as it recurses to itself to determine primeness. That takes a lot of unnecessary computing time.

And if you followed through all of the above, good for you! I'll admit that it's easier for me to rely mostly on the copy and paste from my original post and trust that I got it right, as, um, did I mention counting primes is boring to me now? It's been years since I sat down and thought through this derivation again.

Oh yeah, so there is a differential form as well, so I did get to use calculus, and as a physics guy that was just way cool too. But then you're getting away from exactness, as using it requires approximate methods for the numerical integration.

If you're bored, you can figure out the partial differential equation that follows from the fully mathematicized form, or just find it on my blog as I have it in several places. Seemed like a big deal to me a long time ago.

But now I've had this information for over a decade. So it's just I guess familiar now.

I guess it is kind of cool though, as I refresh on how things work for this post.

James Harris

No comments: