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'.

Monday 26 September 2016

Longest Common Subsequence

Code:


#include<bits/stdc++.h>
using namespace std;
void printLCS(string &x,string &y,int**table,int i,int j)
{
	if(i<=0 || j<=0)
		return ;
	if(x[i-1]==y[j-1])
	{
		printLCS(x,y,table,i-1,j-1);
		cout<<x[i-1];
	}
	else
	{
		if(table[i-1][j]>table[i][j-1])
			printLCS(x,y,table,i-1,j);
		else printLCS(x,y,table,i,j-1);
	}
}
int LCS(string &x,string &y) // using a 2-D table for storing solutions to subproblems
{
	int**table;
	table=new int*[x.size()+1];
	for(int i=0;i<=x.size();i++)
		table[i]=new int[y.size()+1];

	for(int i=0;i<=x.size();i++)
		table[i][0]=0;
	for(int i=0;i<=y.size();i++)
		table[0][i]=0;

	for(int i=1;i<=x.size();i++)
	{
		for(int j=1;j<=y.size();j++)
		{
			if(x[i-1]==y[j-1])
				table[i][j]=table[i-1][j-1]+1;
			else
				table[i][j]=max(table[i-1][j],table[i][j-1]);
		}
	}
	printLCS(x,y,table,x.size(),y.size());
	cout<<endl;
	int answer=table[x.size()][y.size()];
	
	for(int i=0;i<=x.size();i++)
		delete table[i];
	delete table;
	return answer;
}

int main()
{
	string x,y;
	cin>>x>>y;
	cout<<"Length of LCS is : "<<LCS(x,y)<<endl;
	return 0;
}


The following are two optimized solution that follow the same approach but are more space efficient. Instead of using a 2-D table, we can use only two 1-D arrays to store the solution to subproblems.

The following is the code:


// just swapping the pointers of the array
int LCSOptimized1(string &x,string &y) // using two 1-D arrays to store the solutions to subproblems
{
	int *table1,*table2;
	table1=new int[y.size()+1];
	table2=new int[y.size()+1];

	for(int i=0;i<=y.size();i++)
		table1[i]=0;

	for(int i=1;i<=x.size();i++)
	{
		for(int j=1;j<=y.size();j++)
		{
			if(x[i-1]==y[j-1])
				table2[j]=1+table1[j-1];
			else table2[j]=max(table2[j-1],table1[j]);
		}
		int *temp;
		temp=table1;
		table1=table2;
		table2=temp;
	}
	return max(table1[y.size()],table2[y.size()]);
}




// a more tricky solution that uses two 1-D arrays.
int LCSOptimized2(string &x,string &y)
{
	int **table;
	table=new int*[2];
	table[0]=new int[y.size()+1];
	table[1]=new int[y.size()+1];

	for(int i=0;i<=y.size();i++)
		table[0][i]=0;

	int index=0;
	for(int i=1;i<=x.size();i++)
	{
		for(int j=1;j<=y.size();j++)
		{
			if(x[i-1]==y[j-1])
				table[index^1][j]=table[index][j-1]+1;
			else
			{
				table[index^1][j]=max(table[index^1][j-1],table[index][j]);
			}
		}
		index=(index^1);
	}
	int answer=table[index][y.size()]; // answer=max(table[0][y.size()],table[1][y.size()]);
	delete table[0];
	delete table[1];
	return answer; 	
}


Largest Sum Contiguous Subarray

Kadane's Algorithm

The following code works even if all the elements of the array are negative.
Code:

#include<bits/stdc++.h>
using namespace std;
int kadane(int*arr,int n)
{
 int sum=0,max=INT_MIN;
 for(int i=0;i<n;i++)
 {
  sum=sum+arr[i];
  if(max<sum)
   max=sum;
  if(sum<0)
   sum=0; 
 }
 return max;
}
int main()
{
 int *arr,n;
 cin>>n;
 arr=new int[n];
 for(int i=0;i<n;i++)
  cin>>arr[i];
 cout<<kadane(arr,n)<<endl;
 return 0;
}


Number of Monotonic Lattice paths - Dynamic Programming

Find the number of monotonic lattice paths along the edges of a grid with n × n square cells, which do not pass above the diagonal. A monotonic path is one which starts in the lower left corner, finishes in the upper right corner, and consists entirely of edges pointing rightwards or upwards.

Code:

#include<bits/stdc++.h>

using namespace std;
typedef long long int ll;
ll findMonotonicLatticePaths(int n)
{
	// simple formula is Catalan number 
	// 2nCn / (n+1)
	ll**table;
	table=new ll*[n+1];
	for(int i=0;i<=n;i++)
		table[i]=new ll[n+1];
	// table[i][j] denotes the number of shortest paths to reach (0,n) from (n,0)
	table[n][0]=0;
	table[n][1]=1;
	
	for(int j=1;j<=n;j++)
	{
		table[n][j]=1; // base case
		for(int i=n-1;i+j>=n;i--)
		{
			int x=0,y=0;
			if(j-1>=0 && i+j-1>=n)
				x=table[i][j-1];
			if(i+1<=n && i+1+j>=n)
				y=table[i+1][j];
			table[i][j]=x+y;
		}
	}
	return table[0][n];
}

int main()
{
	int n;
	cin>>n;
	cout<<findMonotonicLatticePaths(n)<<endl;
}

Sunday 25 September 2016

Subset Sum Problem - Dynamic Programming

Given a set of n numbers a sum up to M , and any K ≤ M , say whether there is a subset of the numbers such that they sum to K? We assume n might be as big as 1000, but M or K is not too big.

code
#include<bits/stdc++.h>
using namespace std;

// function to print the elements of subset that forms the given sum
void printSubset(int*arr,int n,int S,bool**table,int i,int j)
{
    if(i<=0 || j<=0)
        return;

    if(table[i][j-1]) // the jth element was not taken for sure.
    {
        printSubset(arr,n,S,table,i,j-1);
    }
    else // the jth element was taken for sure.
    {
        cout<<arr[j-1]<<" ";
        printSubset(arr,n,S,table,i-arr[j-1],j-1);
    }
}

bool isSubsetSum(int *arr,int n,int S)
{
    bool **table; // table of size S x n
    table=new bool*[S+1];
    for(int i=0;i<=S;i++)
        table[i]=new bool[n+1];

    // table[i][j] will be true if sum 'i' can be achieved using first j elements of the array.

    //base cases

    // from the first 0 elements you can never achieve any sum.
    for(int i=0;i<=S;i++)
        table[i][0]=false;
    // 0 sum can be achieved by any number of elements.
    for(int i=0;i<=n;i++)
        table[0][i]=true;

    for(int i=1;i<=S;i++)
        for(int j=1;j<=n;j++)
        {
            table[i][j]=table[i][j-1];
            if(i-arr[j-1]>=0)
                table[i][j]=table[i][j] || table[i-arr[j-1]][j-1];
        }

    bool answer=table[S][n];
    if(answer)
        printSubset(arr,n,S,table,S,n);
    
    for(int i=0;i<=S;i++)
        delete table[i];
    delete table;
    
    return answer;
}

int main()
{
    int n;
    cout<<"Enter the number of elements of the array : ";
    cin>>n;
    int *arr;
    arr=new int[n];
    
    cout<<"Enter the elements of the array : ";
    for(int i=0;i<n;i++)
        cin>>arr[i];

    int S;
    cout<<"Enter the sum : ";
    cin>>S;
    cout<<endl<<isSubsetSum(arr,n,S)<<endl;;
    return 0;
}

Remark. The original idea is to use a 2 dimensional array, where each column only depends on the previous column. By a programming trick we just need one column. But we need to write the j-loop in the reversed way to avoid messing things up. 
The following code shows how to achieve this.


bool isSubsetSumOptimized(int*arr,int n,int S)
{
    bool *table;
    table=new bool[S+1];
    table[0]=1;
    for(int i=0;i<n;i++)
        for(int j=S;j>=arr[i];j--)
            table[j]=table[j]||table[j-arr[i]];
    bool answer=table[S];
    delete table;
    return answer;
}