www.pudn.com > NumericalComputing.rar > primespiral.m, change:2003-12-23,size:3485b


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