function [Sout,Pout] = primespiral(n,c)
% PRIMESPIRAL  Ulam's prime number spiral.
%   PRIMESPIRAL(n,c) plots the prime numbers in the n-by-n matrix
%   generated by storing c:c+n^2 in a spiral pattern starting in
%   the center.  The concentration of primes on some of the diagonals
%   is remarkable.  This phenomenon was discovered by Stanislaw Ulam
%   in 1963, and featured on the cover of Scientific American in
%   March, 1964.  For n <= 42, primespiral(n,c) shows the spiral 
%   numbering scheme.  The default value of c is 1.
%
%   With one or two output arguments,
%      S = primespiral(n,c), or
%      [S,P] = primespiral(n,c)
%   does no plotting, but returns S = prime spiral and P = isprime(S).
%
%   Examples:
%      primespiral(7)
%      primespiral(17,17)
%      primespiral(41,41)
%      primespiral(250)

if nargout > 0
   [Sout,Pout] = spiralprimes(n,c);
   return
end

left = findobj('string','<');
right = findobj('string','>');
if isempty(left) | isempty(right)
   clf
   shg
   set(gcf,'doublebuffer','on','name','Prime spiral', ...
      'menu','none','numbertitle','off');
   left = uicontrol('style','toggle','pos',[20 20 30 20], ...
      'string','<','fontweight','bold','callback','primespiral(''<'')');
   right = uicontrol('style','toggle','pos',[60 20 30 20], ...
      'string','>','fontweight','bold','callback','primespiral(''>'')');
end

if nargin < 1, n = 250; end
if nargin < 2, c = 1; end
if isstr(n)
   if n == '<', set(right,'value',0), end
   if n == '>', set(left,'value',0), end
   n = get(left,'userdata'); c = get(right,'userdata');
end
set(left,'userdata',n)

while 1
   [S,P] = spiralprimes(n,c);
   spydetail(S,P)
   title(['c = ' int2str(c)])
   xlabel([int2str(nnz(P)) ' primes'])
   drawnow
   if get(left,'value'), c = max(c-1,1); end
   if get(right,'value'), c = c+1; end
   if c == 1
      set(left,'enable','off')
   else
      set(left,'enable','on')
   end
   set(right,'userdata',c)
   if get(left,'value') == 0 & get(right,'value') == 0
      break
   end
end


% ------------------------

function [S,P] = spiralprimes(n,c)
% SPIRALPRIMES
%   [S,P] = spiralprimes(n,c) returns two n-by-n matrices.
%   S is the spiral numbering of c:c+n^2-1.
%   P is true for the primes in S.

m = n^2-1+c;
P = zeros(m,1);
P(primes(m)) = 1;
S = spiral(n)+c-1;
P = reshape(P(S(:)),n,n);


% ------------------------

function S = spiral(n)
%SPIRAL SPIRAL(n) is an n-by-n matrix with elements
%   1:n^2 arranged in a rectangular spiral pattern.

S = [];
for m = 1:n
   S = rot90(S,2);
   S(m,m) = 0;
   p = m^2-m+1;
   v = (m-1:-1:0);
   S(:,m) = p-v';
   S(m,:) = p+v;
end
if mod(n,2)==1
   S = rot90(S,2);
end


% ------------------------

function spydetail(S,P)
% SPYDETAIL
%   SPYDETAIL(S,P) is like SPY(P) with the element values from S.

[n,n] = size(P);
if n <= 42
   delete(gca)
   axis([0 n+1 0 n+1])
   axis square
   axis ij
   for i = 1:n
      for j = 1:n
         if P(i,j)
            color = [0 0 1];
            fs = max(1,floor(16-2*log2(n)));
         else
            color = [1/3 1/3 1/3];
            fs = max(1,floor(16-3*log2(n)));
         end
         text(j,i,int2str(S(i,j)),'fontsize',fs,'color',color)
      end
   end
else
   if n <= 100, ms = 6; else, ms = 1; end
   [i,j] = find(P);
   plot(j,i,'.','markersize',ms)
   axis([0 n+1 0 n+1])
   axis square
   axis ij
end