Hacking Primes in Mathematica: "
If this is too esoteric, skip it. I couldn't figure out anywhere else to put it.
This morning, Tim Bray tweeted about a post on prime numbers and Benford's law. To cut the esoterica short, one of the big problems in prime numbers is that people don't know how they're distributed. This post suggests that Benford's Law describes the distribution of the first digit of prime numbers. One of the comments asked an important question: is this really just an artifact of base 10? Math really doesn't 'know anything' about bases, so if this idea doesn't generalize to bases other than 10, it doesn't mean much.
That challenge was a bit hard to ignore. A bit of futzing with Mathematica later, I got it. Here's a graph of the distribution of prime numbers in base 36, for primes less than 36^5 (about 60 million, 3.5 million primes):
primes.tiff
Here's the code (somewhat improved from the code that generated the graph)--only 5 very Lisp-like lines. For an explanation of the math, see the original article.
With[{base = 36, ades = 5, a = 1.1}, nprimes := PrimePi[base^ades];
alpha := 1/(Log[base^ades] - a);
gbl[ d_] := ((d + 1)^(1 - alpha) - d^(1 - alpha))/(base^(1 - alpha) - 1);
ListPlot[{Table[100*gbl[digit], {digit, 1, base - 1}] ,
(100/nprimes)* Sort[Tally[Table[IntegerDigits[Prime[n], base][[1]], {n, nprimes}]]][[All, 2]]}]]
Sal Mangano, who's writing our Mathematica Cookbook, inspired the last improvement--now, rather than taking about two minutes to compute and sort primes up to 36^5, it takes a few seconds. That's why he's writing the book, not me.
This certainly isn't rigorous math, from the standpoint of proving anything about the distribution of prime numbers. Just some fun hacking--it's great to have a mathematical language in which you can say 'Give me all the primes less than 32^5 and sort them by their first digits in base 36'--in about that much space.
O'Reilly Radar - http://radar.oreilly.com/