%!PS
% Time-stamp: <28 Feb 2000 11:13 charpov@cs.unh.edu>

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

% Resolution in dpi
/dpi 24 def

% Colors ?
/color? true def

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

/offset 0 def
% /offset 41 def
% /offset 65 def
% /offset 49921 def
% Change to true to print prime numbers inside the dots
/number? true def

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

% Define some constants
/dot 72 dpi div def
/maxpts papery papery mul dot dup mul div def
% /maxpts 4000 def
% Allocate a string and load a font for printing numbers
number? {
  /buf 10 string def
  /fontsize 9 def % number here is irrelevant
  /Courier findfont fontsize scalefont setfont
} if

% Brute force primality test
% (postscript interpreters are faster than they have memory...)
/prime? {
  /n exch def
   n 0 eq { true } {n 1 eq {true} { % 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
  } ifelse
} bind def

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 setrgbco
+lor
      } {
        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) spira
+l
/markpoint {
  currentpoint % save currentpoint that is lost because of stroke or n
+ewpath
  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
   pri { 0 0 0 setrgbcolor}{ 1 1 1 setrgbcolor} ifelse fill % draw the
+ (previous) dot, always
  number? { % write the (previous) number in white over the dot, if ne
+cessary
    gsave % because we use 'scale'
      pri { 1 setgray } { 0 setgray } ifelse
      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 prin
+t
    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
  %cx cy moveto
}
 bind def
% x_c y_c k r_0 theta_1 theta_2 spiral
/spiral {
  13 dict begin
    /th2 exch def
    /th1 exch def
    /r0 exch def
    /k exch def
    /yc exch def
    /xc exch def
    /precision exch def   % Increase this for better approximation
    /r1 360 k mul r0 add def
    /narcs      % ceiling((precision * r1) / r0)
      r1 precision mul r0 div ceiling cvi def
    /radadjust r1 r0 sub narcs div def  % (r1 - r0) / narcs
    /arcmeasure 360 narcs div def       % Measure of each arc
    /adjust {
      radadjust add
      arcmeasure rotate
      radadjust neg 0 translate
      /th1 th1 arcmeasure sub def
      /th2 th2 arcmeasure sub def
    } def
    matrix currentmatrix
    r0 
    arcmeasure arcmeasure th1 { pop adjust } for
    /centerandrad { dup 0 0 3 -1 roll } def
    centerandrad
    th1 th2 arcmeasure le { th2 } { arcmeasure } ifelse
    arc adjust {
      th2 arcmeasure le { exit } if
      centerandrad 0 arcmeasure arc adjust
    } loop
    th2 0 gt { centerandrad 0 th2 arc } if
    pop
    setmatrix end
} def   

% newpath

/p 1 def
/i 1 def

newpath
/cx paperx 2 div def
/cy papery 2 div def

cx cy moveto % center the picture
spiral? {
  gsave
  newpath
  cx cy translate
  1 dpi div 2 mul dup scale
  % cx cy moveto
  1 0 0 0.1 1 0 360 maxpts sqrt mul spiral stroke % draw better spiral
+ for first 5 revolutions
  grestore
  dot 10 div setlinewidth
} if
% /prime? true def
/x cx def /y cy def
newpath x y dot 2 div 0 360 arc % draw the first dot
cx cy moveto % center the picture
currentpoint
/count 0 def

{
  /radius i sqrt def
  /angle  radius 360 mul def
  /px angle cos radius mul dot mul cx add def
  /py angle sin radius mul dot mul cy add def
  spiral? {
    count 99998 gt {
      count 99999 eq { px py moveto } if
      /r1 i 0.6 sub sqrt def
      /a1 r1 360 mul def
      /r2 i 0.4 sub sqrt def
      /a2 r2 360 mul def
      /x1 a1 cos r1 mul dot mul cx add def
      /y1 a1 sin r1 mul dot mul cy add def
      /x2 a2 cos r2 mul dot mul cx add def
      /y2 a2 sin r2 mul dot mul cy add def
      x1 y1 x2 y2 px py curveto
    } if
  } if
  /count count 1 add def
  /pi i offset add def
  /pri pi prime? def
  % pri { markpoint } if
  markpoint
  /i i 1 add def % increment i
  /x px def /y py def
  x y moveto
  i maxpts ge {exit} if
} loop

% Signature
gsave
  20 20 translate
  1 setlinejoin
  1 setlinewidth
  1 setgray
  newpath
  20 20 20 0 360 arc fill
  0 setgray
  newpath
  10 5 moveto
  0 20 rlineto
  10 -8 rlineto
  10 8 rlineto
  0 -20 rlineto
  10 25 moveto
  4 12 rlineto
  14 5 moveto
  6 12 rlineto
  stroke
  newpath
  22 21 1 0 360 arc fill
grestore
showpage
[download]