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 commenters 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

Sal Mangano, who's writing our Mathematica Cookbook, says my approach is very inefficient, because it's using character matching to sort the primes. He's probably right; that's why he's writing the book. I'm still happy with it. Here's the code--only 7 very Lisp-like lines:

With[{base = 36, ades = 5, a = 1.1}, 
   nprimes := PrimePi[base^ades];  
   fcharnum[n_, b_] :=  First[Characters[ToString[BaseForm[n, b]]]]; 
   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[ fcharnum[Prime[n], base], {n, nprimes}]]][[All, 2]]}] ]

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.