• Print

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.

tags: , , , ,
  • joem

    I would totally love to see this problem done (minus the graphing) in various other popular languages.

  • http://pyevolve.sourceforge.net/wordpress/ Christian S. Perone

    Hello Mike, good work ! I put a link to your work on the Pyevolve post about this.

  • http://www.spurgeonworld.com Chris Spurgeon

    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…

    Problem 1

    If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.

    Find the sum of all the multiples of 3 or 5 below 1000.

    …but soon become downright diabolical…

    Problem 233

    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

    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

    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

    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

    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

    Isn´t this what Wolfram Alpha is supposed to do?

  • Drake

    Hey, this is really cool. I love taking problems and throwing them into Mathematica and seeing what it does.

  • guardian

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