www.pudn.com > NumericalComputingwithMatlabCode.zip > 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', ...
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
```