Friday, 30 September 2016

Ulam Spiral

The following post shows the distribution of primes when numbers are arranged in a spiral format as depicted below:



The above shown grid is a 7x7 grid.

If similarly we have a grid of 251x251 elements and highlight the primes withe red pixel, the image will look like the following:


And similarly if we randomly put a 0 or 1 on every pixel of 251x251 grid, we will get an image that will be somewhat like this:


We can clearly see that there is a pattern visible in case of Ulam Spiral where many prime numbers are present along different diagonals, where as random numbers follow no specific patterns.

If you are interested into the code of generating this ulam spiral and getting the above plots, have a look at below mentioned code. Basically I have used C++ to get a 2-D matrix where 1 represents that number is prime and 0 represents that it is composite. This all data is written into a file called ulam.txt and random data is written into random.txt.

Now these two file can be fed to the scilab code to get the above plots.

C++ code that will generate random.txt and ulam.txt :

#include<bits/stdc++.h>
#include<bits/stdc++.h>
#include<stdlib.h>
#include<stdlib.h>
using namespace std;
#define MAX 10000000
bitset<MAX+1>isPrime;

void getPrimes()
{
    isPrime.reset();
    isPrime.flip();
    isPrime[0]=isPrime[1]=0;
    for(int i=2;i*i<=MAX;i++)
        if(isPrime[i])
            for(int j=i*i;j<=MAX;j+=i)
                isPrime[j]=0;
}

int**getUlamSpiral(int N)
{
    int **M;
    M=new int*[2*N+1];
    for(int i=0;i<2*N+1;i++)
        M[i]=new int[2*N+1];

    int r,c,current,counter;
    
    M[N][N]=1;
    r=N;
    c=N+1;
    counter=2;
    current=2;

    for(int j=0;j<N;j++,counter+=2)
    {
        for(int i=0;i<counter;i++)
            M[r--][c]=current++;
        r++;
        c--;
        for(int i=0;i<counter;i++)
            M[r][c--]=current++;
        c++;
        r++;
        for(int i=0;i<counter;i++)
            M[r++][c]=current++;
        r--;
        c++;
        for(int i=0;i<counter;i++)
            M[r][c++]=current++;
    }

    return M;
}
void printUlamSpiral(int **M,int N)
{
    for(int i=0;i<2*N+1;i++)
    {
        for(int j=0;j<2*N+1;j++)
            cout<<M[i][j]<<"\t";
        cout<<endl;
    }
}
int**filterPrimes(int**M,int N)
{
    int **M01;
    M01=new int*[2*N+1];
    for(int i=0;i<2*N+1;i++)
        M01[i]=new int[2*N+1];

    for(int i=0;i<1+2*N;i++)
    {
        for(int j=0;j<1+2*N;j++)
        {
            if(isPrime[M[i][j]])
                M01[i][j]=1;
            else M01[i][j]=0;
        }
    }
    return M01;
}
void writeToFile(string fileName,int**M,int N)
{
    FILE*fp;
    fp=fopen(fileName.c_str(),"w");
    if(fp)
    {
        fprintf(fp,"%d\n",N);
        for(int i=0;i<1+2*N;i++)
        {
            for(int j=0;j<1+2*N;j++)
            {
                fprintf(fp,"%d ",M[i][j]);
            }
            fprintf(fp,"\n");
        }
    }
    else cout<<"there is an error in opening the file\n";
}
int**fillRandom(int N)
{
    int **m;
    m=(int**)malloc(sizeof(int*)*(2*N+1));
    for(int i=0;i<2*N+1;i++)
        m[i]=(int*)malloc(sizeof(int)*(2*N+1));
    srand(time(NULL));
    for(int i=0;i<2*N+1;i++)
        for(int j=0;j<2*N+1;j++)
            m[i][j]=rand()%2;
    return m;
}
int main()
{
    int N;
    getPrimes();
    cin>>N;
    int **M=getUlamSpiral(N);
    //printUlamSpiral(M,N);
    int **M01=filterPrimes(M,N);
    int**R01=fillRandom(N);
    writeToFile("ulam.txt",M01,N);
    writeToFile("random.txt",R01,N);
    return 0;
}

Scilab Code
function DrawSpiral(filename)
    fid=mopen(filename,"r");
    [num_read,N]=mfscanf(fid, "%d");
    N
    M = hypermat([2*N+1 2*N+1]);
    for i=1:1+2*N
        for j=1:1+2*N
            [num_read,M(i,j)]=mfscanf(fid, "%d");
        end
    end
    mclose(fid);
    
    clf();
    ax = gca();//get current axes handle
    ax.data_bounds = [0, 0; 100, 100]; //set the data_bounds
    ax.box = 'on'; //draw a box
    Matplot1(M*20, [0,0,100,100])
endfunction


DrawSpiral('ulam.txt')



Steps of Execution :
First of all save the C++ code and compile it and then run it. The code will ask you to input a number N and it will create 2*N+1 by 2*N+1 matrix.

Now run the Scilab code to get the required plots.
You need to call DrawSpiral function 2 times, once with 'ulam.txt' and then with 'random.txt'.

No comments:

Post a Comment