Ulam’s Prime Spiral

To generate a Ulam’s prime spiral the positive integers are arranged in a spiral pathern and prime numbers are highlighted in some way along the spiral. For example we used blue pixels to represent primes and white pixels for composite numbers. The code below uses built in Matlab function to generate the Ulam spiral.

close all;
clc
clear all
sz  = 101; % Size of the NxN square matrix
mat = spiral(sz);
k =1;
for i =1:sz
    for j=1:sz
        if isprime(mat(i,j)) % Check if the number is prime
            % saving indices of primes
            y(k) = i;
            x(k) = j;
            k = k+1;
        end
    end
end
figure('Position',[0 0 800 800]);
hold on;
colormap bone;
scatter(x,y,'.b');
axis([1 sz  1 sz ]);
axis off;
set(gca,'YDir','Reverse');

Now lets construct Varanoi diagram for the prime numbers (higlighted points in blue) in Ulam spiral.

figure('Position',[0 0 800 800]);
hold on;
xy = [x',y'];
[v,c] = voronoin(xy);
for i = 1:length(c)
  patch(v(c{i},1),v(c{i},2),'w');
end
axis([1 sz  1 sz ]); axis off;
set(gca,'YDir','Reverse');
scatter(x,y,'.b');

Cells in Voronoi diagram of the Ulam Prime Spiral

% We observe various poligons with  n-sides such as triangles, quadrangles, pentagons, etc...
% Lets call
% A3 = a sequence of prime numbers which has a Triangular Voronoi cell in Voronoi diagram of the Ulam Prime Spiral
% A4 = a sequence of prime numbers which has a Quadrilateral Voronoi cell in Voronoi diagram of the Ulam Prime Spiral
% A5 = a sequence of prime numbers which has a Pentagonal Voronoi cell in Voronoi diagram of the Ulam Prime Spiral
% A6 = 6 sides
% A7 = 7 sides
% An = n sides
% Q1: What is the maximum value for n, n_max ? In other words is there an uper bound ?
% Q2: If yes, What are the primes whic has a n_max side poligon cells in Voronoi diagram of the Ulam Prime Spiral
% The code below is used to calculate An
k3 = 1; k4 = 1; k5 = 1; k6 = 1; k7 = 1; k8 = 1;
for i = 1:length(c)
  szv = size(v(c{i},1));
  polyN(i) = szv(1);
  switch polyN(i)
    case 3
        A3(k3) = mat(y(i),x(i));
        k3 = k3+1;
        A3xy(k3,:)= [x(i),y(i)];
    case 4
        A4(k4) = mat(y(i),x(i));
        k4 = k4+1;
        A4xy(k4,:)= [x(i),y(i)];
    case 5
        A5(k5) = mat(y(i),x(i));
        k5 = k5+1;
        A5xy(k5,:)= [x(i),y(i)];
    case 6
        A6(k6) = mat(y(i),x(i));
        k6 = k6+1;
        A6xy(k6,:)= [x(i),y(i)];
    case 7
        A7(k7) = mat(y(i),x(i));
        k7 = k7+1;
        A7xy(k7,:)= [x(i),y(i)];
    case 8
        A8(k8) = mat(y(i),x(i));
        k8 = k8+1;
        A8xy(k8,:)= [x(i),y(i)];
  end
end
%
% First 15 terms
A3 = sort(A3);
fprintf('A3 = ');
fprintf('%i, ',A3(1:15));
fprintf('\n');
A4 = sort(A4);
fprintf('A4 = ');
fprintf('%i, ',A4(1:15));
fprintf('\n');
A5 = sort(A5);
fprintf('A5 = ');
fprintf('%i, ',A5(1:15));
fprintf('\n');
A6 = sort(A6);
fprintf('A6 = ');
fprintf('%i, ',A6(1:15));
fprintf('\n');
A7 = sort(A7);
fprintf('A7 = ');
fprintf('%i, ',A7(1:15));
fprintf('\n');
A8 = sort(A8);
fprintf('A8 = ');
fprintf('%i, ',A8(1:15));
fprintf('\n');

A3 = 313, 389, 1283, 1399, 1669, 1787, 2087, 2143, 2713, 2801, 3469, 4091, 4787, 4789, 4903, 
A4 = 23, 31, 47, 59, 71, 73, 79, 131, 139, 167, 173, 181, 229, 239, 251, 
A5 = 2, 3, 11, 13, 17, 19, 29, 37, 53, 83, 97, 101, 103, 107, 109, 
A6 = 5, 7, 41, 43, 89, 127, 179, 193, 233, 263, 283, 317, 379, 383, 397, 
A7 = 61, 157, 199, 311, 349, 409, 463, 509, 557, 601, 641, 691, 727, 757, 823, 
A8 = 67, 491, 613, 1013, 1117, 1201, 1249, 1301, 1373, 1543, 1753, 1907, 2017, 2339, 2411, 

% Plotting
% The color code used
% n = 3, Triangle yellow
% n = 4, Tetragon green
% n = 5, Pentagon magenta
% n = 6, cyan
% n = 7, red
% n = 8, yellow

figure('Position',[0 0 800 800]);
hold on;
xy = [x',y'];
[v,c] = voronoin(xy);
for i = 1:length(c)
  patch(v(c{i},1),v(c{i},2),'w');
end
axis([1 sz  1 sz ]); axis off;
set(gca,'YDir','Reverse');
scatter(A3xy(:,1),A3xy(:,2),'.k');
scatter(A4xy(:,1),A4xy(:,2),'.g');
scatter(A5xy(:,1),A5xy(:,2),'.m');
scatter(A6xy(:,1),A6xy(:,2),'.c');
scatter(A7xy(:,1),A7xy(:,2),'.r');
scatter(A8xy(:,1),A8xy(:,2),'.y');
% Note that the last terms can be wrong. They corespond to the points on the outer
% edges of the spiral which might be altered when considering a larger spiral.

Triangular Cells with Integer Area

k3 = 1;
k3i = 1;
fprintf('Prime, Cell Area, Perimeter \n');
for i = 1:length(c)
  szv = size(v(c{i},1));
  polyN(i) = szv(1);
    if(polyN(i) == 3)
       % Sides of a triangle
         A = v(c{i}(1, 1),:) - v(c{i}(1, 2),:);
         B = v(c{i}(1, 1),:) - v(c{i}(1, 3),:);
         C = v(c{i}(1, 2),:) - v(c{i}(1, 3),:);
        % Perimeter of a triangle
         pv = norm(A)+ norm(B)+ norm(C);
         P3(k3) = pv;
        % Area of a triangle
        AB = cross([A,0], [B,0], 2);
        S3(k3) =  1/2 * sum(sqrt(sum(AB.^2, 2)));
        A3(k3) = mat(y(i),x(i));
        fprintf('%i %4.4f %4.4f \n',A3(k3), S3(k3),P3(k3))
        %  Primes with triangle Voronoi cells with integer area
        if S3(k3)== round(S3(k3)) % check if area is integer
            A3i(k3i) = mat(y(i),x(i));
            k3i = k3i +1;
        end
            k3 = k3+1;
    end
end
% Primes with triangular Voronoi cells with integer area
fprintf('\n A3i = ');
A3i = sort(A3i);
fprintf('%i, ',A3i);

Prime, Cell Area, Perimeter 
8999 9.0000 14.4853 
9001 9.0000 14.4853 
8627 4.5000 10.2426 
10093  NaN  Inf 
6869 8.0000 13.6569 
9419 6.0000 11.4049 
5867 6.0000 11.4049 
3469 6.0000 11.4049 
2801 6.0000 11.4049 
1669 4.5000 10.2426 
9439 4.5000 10.2426 
4787 6.0000 11.4049 
4789 6.0000 11.4049 
2143 4.5000 10.2426 
1787 4.5000 10.2426 
4933 6.0000 11.4049 
313 6.0000 11.4049 
389 6.0000 11.4049 
1399 6.0000 11.4049 
2713 4.5000 10.2426 
1283 6.0000 11.4049 
2087 4.0000 9.6569 
4091 6.0000 11.4049 
4903 6.2500 12.0711 
8111 6.0000 11.4049 
6037 6.2500 12.0711 
10007  NaN  Inf 
9929  NaN  Inf 

 A3i = 313, 389, 1283, 1399, 2087, 2801, 3469, 4091, 4787, 4789, 4933, 5867, 6869, 8111, 8999, 9001, 9419,