(download)

%!PS

% Postscript program to draw 'prime numbers patterns'. 
% Program options are as follows:
% * paper size (in inches)
% * resolution (in dpi). Choose very low resolution (like 5) to get
%   close-ups without having to zoom in.
% * color: if color is true, colors are defined from primes distribution
% * spiral: if spiral is true, then not only the dots will appear, but the 
%   lines between them are drawn too. This is nice for close-ups.
% * startzero: if true, starts with 0 instead of 1.
% * number: if true, print prime numbers inside the dots.
                         

% Paper size, in dots (letter)
/paperx 8.5    72 mul def
/papery 11     72 mul def

% Resolution in dpi
/dpi 72 def

% Colors ?
/color? false def

% Change to true to make the 'spiral' appear.
/spiral? false def

% Change to true to have the spiral start with 0 instead of 1
/startzero? false def

% Change to true to print prime numbers inside the dots
/number? false def

% ============= CHANGES BELOW AT YOUR OWN RISK ============================

% Define some constants
/dot 72 dpi div def
/maxpts papery papery mul dot dup mul div def

% Allocate a string and load a font for printing numbers
number? {
  /buf 10 string def
  /fontsize 9 def % number here is irrelevant
  /Courrier-Bold findfont fontsize scalefont setfont
} if

% Brute force primality test
% (postscript interpreters are faster than they have memory...)
/prime? {
  /n exch def
  n 1 eq {false} { % 1 is not prime
    n 2 mod 0 eq {n 2 eq} {
      /divide? false def
      3 2 n sqrt {
        n exch div dup ceiling eq {/divide? true def exit} if
      } for
      divide? not
    } ifelse
  } ifelse
} bind def

% A little less stupid test (printers are not that fast, after all...)
/maxarray 65534 def
/primes maxarray 1 add array def
/lprime -1 def
/prime? {
  /n exch def
  /sn n sqrt def
  /d 1 def
  /id 0 def
  n 1 eq {false} { % 1 is still not prime
    n 2 mod 0 eq {n 2 eq} {
      /divide? false def
      {
        id lprime gt {
          /d d 2 add def
        } {
          /d primes id get def
          /id id 1 add def
        } ifelse
        d sn gt {exit} if
        n d mod 0 eq {/divide? true def exit} if
      } loop
      divide? {
        false
      } {
        lprime maxarray lt {
          /lprime lprime 1 add def
          primes lprime n put
        } if
        true
      } ifelse
    } ifelse
  } ifelse
} bind def

% 'random' primes, to compare
% /prime? {
%   dup 2 mod 0 ne exch
%   1 add ln cvi rand exch mod 0 eq and} bind def

% Change colors
% If a prime is close from its predecessor (<ln(x)), it tends to be blue,
% green if is is at a 'normal' distance (ln(x)), red if it is too far
% and black if it is much too far (>2ln(x))
color? {
  /ppi 0 def
  /chgcolor {
    pi ppi sub pi ln lt {
      0 pi ppi sub 2 sub pi ln 2 sub div dup 1 exch sub setrgbcolor
    } {
      pi ppi sub pi ln 2 mul 2 sub lt {
        pi ppi sub pi ln sub pi ln 2 sub div dup 1 exch sub 0 setrgbcolor
      } {
        0 setgray
      } ifelse
    } ifelse
    /ppi pi def
  } bind def
} if

% Mark a point when a number is prime
% We actually draw the PREVIOUS dot, to overwrite the (possible) spiral
/markpoint {
  currentpoint % save currentpoint that is lost because of stroke or newpath
  spiral? {0 setgray stroke} if % draw the spiral, if necessary
  color? {chgcolor} if % change the color, if necessary
  newpath px py dot 2 div 0 360 arc fill % draw the (previous) dot, always
  number? { % write the (previous) number in white over the dot, if necessary
    gsave % because we use 'scale'
      1 setgray
      px py moveto
      pi buf cvs dup stringwidth pop 
          % have the number fit 90% of the dot (50% if it is a single digit)
      dup dot pi 10 lt {.5} {.9} ifelse mul exch div dup scale
      -2 div fontsize -3 div rmoveto show % center in the dot and print
    grestore
  } if
  2 copy /py exch def /px exch def % update previous dot location
  /pi i def % update previous number
  moveto % move back to the point we were
} bind def

% Here is the algorithm for moving along the spiral
/do-p-times {
  1 1 p {
    pop 
    x y rlineto % if spiral? is not on, the path will not be stroked
    i 2 eq {
      currentpoint /py exch def /px exch def /pi 2 def
    } {
      i prime? {markpoint} if
    } ifelse
    /i i 1 add def
  } for
} bind def

/p 1 def
/i startzero? {1} {2} ifelse def

newpath
paperx 2 div papery 2 div moveto % center the picture
spiral? {dot 10 div setlinewidth} if

{
  /x dot def /y 0 def do-p-times
  /x 0 def /y dot def do-p-times
  /p p 1 add def
  /x dot neg def /y 0 def do-p-times
  /x 0 def /y dot neg def do-p-times
  /p p 1 add def
  i maxpts ge {exit} if
} loop

showpage