|
|
|||||
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.
|
|||||
© 2012, O'Reilly Media, Inc. (707) 827-7019(800) 889-8969 All trademarks and registered trademarks appearing on oreilly.com are the property of their respective owners. |
|||||
Comments: 10
joem [ 9 May 2009 06:42 PM]
I would totally love to see this problem done (minus the graphing) in various other popular languages.
Christian S. Perone [ 9 May 2009 09:52 PM]
Hello Mike, good work ! I put a link to your work on the Pyevolve post about this.
Chris Spurgeon [10 May 2009 01:07 AM]
This is as good a place as any to put in a plug for my favorite math diversion, Project Euler (www.projecteuler.net).
Project Euler (named after 18th Century Swiss mathematician Leonhard Euler) is an online series of math puzzles that often depend upon clever bits of programming to find the answer. The problems start out fairly straightforward...
...but soon become downright diabolical...
Let f(N) be the number of points with integer coordinates that are on a circle passing through (0,0), (N,0),(0,N), and (N,N).
It can be shown that f(10000) = 36.
What is the sum of all positive integers N ≤ 10^(11) such that f(N) = 420 ?
A new question is released each week, and is immediately pounced on by Project Euler fanatics. Mathematica is the tool of choice for many of them, though people have solved the problems with everything from Perl to Python to C to Lisp to Haskel.
Falafulu Fisi [10 May 2009 02:38 PM]
This is how Matlab compute primes:
----------------------------------------------
function p = primes(n)
%PRIMES Generate list of prime numbers.
%PRIMES(N) is a row vector of the prime numbers %less than or equal to N. A prime number is one %that has no factors other than 1 and itself.
if length(n)~=1
error('N must be a scalar');
end
if n p = zeros(1,0,class(n));
return,
end
p = 1:2:n;
q = length(p);
p(1) = 2;
for k = 3:2:sqrt(n)
if p((k+1)/2)
p(((k*k+1)/2):k:q) = 0;
end
end
p = p(p>0);
David Feustel [11 May 2009 09:15 AM]
James McCanney claims to have discovered how
to calculate prime numbers. See his web site
at calculateprimes.com where you can buy his
book with dvd lecture on how it's done.
vicaya [11 May 2009 12:04 PM]
James McCanney is a known crackpot. He dares not even publish his findings for peer review. If he's right, he'd proved Riemann Hypothesis, which would make him mainstream famous. Now he's just peddling snake oil to the curious and uninformed.
Rex [11 May 2009 04:19 PM]
Mike
It is fairly easy to show that Benford's Law is independent of base. If it's true for primes in one base it is true in all bases.
Should save you a bit of work!
El Rorro [12 May 2009 12:46 AM]
Isn´t this what Wolfram Alpha is supposed to do?
Drake [24 September 2009 05:25 AM]
Hey, this is really cool. I love taking problems and throwing them into Mathematica and seeing what it does.
guardian [23 November 2009 11:53 AM]
>>This certainly isn't rigorous math...
>>Just some fun hacking
for even more fun you can use the manipulate function to change the base. while operation of the slider is not very smooth (on my machine :) ) it is a neat way to step through different bases.