Running Median

Hey folks,
Today we gonna talk about how to keep track of median, given a stream of integers. We will discuss the most standard approach to solve the problem.
Prerequisite : Heap Data Structure

Initialization
1. Put greater of the first two numbers in a MinHeap and the other one into a MaxHeap.

Algorithm
Insertion 
2. If the incoming number is smaller than the MaxHeaps root then insert it into the MaxHeap otherwise insert it into the MinHeap.
Balancing
3. If the size difference between the MaxHeap and MinHeap is greater than 1, then delete the root of the larger heap and insert it into the smaller heap.

Output
// At any given time
4. if the size difference in two heaps == 1:
       return (root of the larger heap)
   else
       return (avg of roots of both the heaps)

External Links
1. Heap Data Structure
2. http://stackoverflow.com/questions/10657503/find-running-median-from-a-stream-of-integers

return 42;

Advertisements

Diameter of a Tree

Hey folks,
Today we are going to learn how to compute the diameter of a tree, for that we first need to understand what do we mean by the diameter of a given tree?

Diameter of a Tree:
The diameter of a given tree is the maximum distance between two nodes of a given tree.
For Ex:
diameter

An O(E) solution exists for the problem. The Algorithm for which goes as follows:
Step 1 : Take an arbitrary node (say A) and do a DFS/BFS on it.
Step 2 : Now take the node which is farthest from A (say B), now do a DFS/BFS on B.
Step 3 : Now take the node which is farthest from B (say C).
Step 4 : Return the distance between B and C.

Explanation:
The first traversal(BFS/DFS) on A gives me a node B, which is supposed to be farthest from node A. Now let’s assume B is not the node which is one of the answer nodes. Then let us assume that some node D is the part of the solution. Now we know D is certainly not the farthest node from the root, Now let us assume that the lowest common ancestor(say E) of the answer nodes lies below the root, then DE must have been greater than AB but it is not possible so D must have been B. The explanation of this is better if it comes from intuition, so try experimenting with two-three hand written trees run the above algorithms on the tree and each time you will end up with a correct answer.

Implementation :
This is a c++ implementation of the above algorithm

#include <vector>
#include <stack>
#include <cstdio>
 
using namespace std;

#define MAX 100
using namespace std;

vector < int > G[MAX]; /* Graph of vertices MAX */
bool visited[MAX];

pair < int ,int > DFS(int A){
	stack < int > S;
	S.push(A);
	visited[A]=true;
	bool flag ;
	int max_height=0;
	int ans_vertex=A;
	int x;
	while(!S.empty()){
		x = S.top();
		flag = false;
		for(int i=0;i<G[x].size();i++){
			if(!visited[G[x][i]]){
				flag = true;
				S.push(G[x][i]);
				visited[G[x][i]]=true;
				break;
			}
		}
		if(S.size()>max_height){
			max_height=S.size();
			ans_vertex=S.top();
		}
		//printf("%d %d %d\n",S.top(),S.size(),max_height);
		if(!flag)
			S.pop();
	}
	/* the distance of the farthest node will be maximum height of the stack-1.*/
	return make_pair(ans_vertex,max_height-1);
}
int main(){
	int e; /* the number edges*/
	scanf("%d\n",&e);
	int x,y;
	/* set the visited array to false*/
	for(int i=0;i<MAX;i++)
		visited[i]=false;
	for(int i=0;i<e;i++){
		scanf("%d %d",&x,&y);
		G[x].push_back(y);
		G[y].push_back(x);
	}
	pair < int ,int > B,C;
	B = DFS(0);
	for(int i=0;i<MAX;i++)
		visited[i]=false;
	C = DFS(B.first);
	printf("diameter is %d\n",C.second);
	return 0;
}
    Sample Input:
9
0 1
0 2
1 3
1 4
2 5
2 6
3 7
3 8
3 9
    Sample Output:
diameter is 5

External Links:
SPOJ problem on diameter of a tree
My solution for the above problem

return 42;

Range Minimum Query

Hey folks,
Today we are going to learn about methods to find the minimum element in given range for an array, this problem is popularly known as range minimum query. First let us see, what are we looking at :

Given an array say A[0,N-1], we are queried for the minimum element from some index ‘i’ to the index ‘j’ such that j>=i, the notation for rest of the post for this will be RMQ(i,j) s.t. j>=i.

Ex: A = [3 , 2 , -1 , 4 , 2]
RMQ(0,2) = -1
RMQ(3,4) = 2
RMQ(0,4) = -1

1. Naive’s Approach
Time Complexity : construction O(1) , Query O(N) <O(1),O(N)>

The idea is simple and straight forward, what we do is nothing but trivially search for the minimum from index ‘i’ to index ‘j’ by actually traversing the array between the indices. The worst case time complexity for the array traversal may be O(N).

#include &lt;cstdio&gt;

int RMQ(int A[],int s,int e){
	int minindex = s;
	for(int i=s;i&lt;=e;i++)
		if(A[i]&lt;A[minindex])
			minindex = i;
	return A[minindex];
}
// driver programme
int main(){
	int A[10] = {3,1,2,-1,5,4,6,0,9,8};
	int s,e;
	scanf(&quot;%d %d&quot;,&amp;s,&amp;e);
	printf(&quot;%d\n&quot;,RMQ(A,s,e));
	return 0;
}

2. Dynamic Programming Approach 1 (Trivial)
Time Complexity : Construction O(N^2) , Query O(1) <O(N^2),O(1)>
Writing recursion for the given problem is no big task.

           j, if A[j] < A[RMQ(i,j-1)]
RMQ(i,j)= 
           RMQ(i,j-1), otherwise

The recursion can easily be memoized by using a two dimensional array of N*N and can be written as:

           j, if A[j] < A[M[i][j-1]]
M[i][j]= 
           M[i][j-1], otherwise
#include &lt;cstdio&gt;
#define MAXN 100
int M[MAXN][MAXN];

int RMQ_DP(int A[],int N){
	for(int i=0;i&lt;N;i++)
		M[i][i]=i;
	for(int i=0;i&lt;N;i++){
		for(int j=i+1;j&lt;N;j++){
			if(A[M[i][j-1]]&gt;A[j])
				M[i][j]=j;
			else
				M[i][j]=M[i][j-1];
		}
	}
}
// driver programme
int main(){
	int A[10] = {3,1,2,-1,5,4,6,0,9,8};
	RMQ_DP(A,10);
	int s,e;
	scanf(&quot;%d %d&quot;,&amp;s,&amp;e);
	printf(&quot;%d\n&quot;,A[M[s][e]]);
	return 0;
}

3. Split and Query
Time Complexity : construction O(N) , Query O(sqrt(N)) <O(N),O(sqrt(N))>
Here we split the array in sqrt(N) equal parts of size sqrt(N), The idea is to store the minimum of each part in the an another array say M(0*sqrt(N)), and then for every query we traverse the array ‘M’ and the remaining elements of the array A.
Let us understand through an example :
A = [2,4,3,1,6,7,8,9,1,7]

RMQ_002
image credits: TC

as it is clear from the image that M[0] is storing the index of the minimum element from A[0] to A[2] then M[1] will store index of the minimum element from A[3] to A[5] … ans so on till it is M[3] which will store A[9] as there are no further elements.

#include &lt;cstdio&gt;
#include &lt;cmath&gt;
#define MAXN 100
int M[MAXN]; 
/* M[i] = the index of the minimum element from A[i*sqrt(N)] to A[i*sqrt(N)+sqrt(N)-1] */
int size_m; // stores the size of the M array
 
void RMQ_SPLIT(int A[],int N){
    size_m =0;
    for(int i=0;i&lt;N;){
        int minindex = i;
        for(int j=0;j&lt;(int)sqrt(N) &amp;&amp; i&lt;N;j++){
            if(A[i]&lt;A[minindex])
                minindex=i;
            i++;
        }
        M[size_m++]=minindex;
    }
}
int RMQ(int A[],int s,int e,int N){
    int start = s/(int)sqrt(N); // this will compute the starting index of the M array 
    int ans = A[s]; 
    int end = e/(int)sqrt(N); // ending index of the M array
    for(int i=s;i&lt;(start+1)*sqrt(N);i++)
      if(A[i]&lt;ans)
        ans = A[i];
    for(int i=start+1;i&lt;end;i++){
        if(A[M[i]]&lt;ans)
            ans=A[M[i]];
    }
    for(int i=end*sqrt(N);i&lt;=e;i++){
        if(A[i]&lt;ans)
            ans = A[i];
    }
    return ans;
}
// driver programme
int main(){
	int A[10] = {3,1,2,-1,5,4,6,0,9,8};
	RMQ_SPLIT(A,10);
	printf(&quot;\n&quot;);
	int s,e;
	scanf(&quot;%d %d&quot;,&amp;s,&amp;e);
	printf(&quot;%d\n&quot;,RMQ(A,s,e,10));
	return 0;
}

4. Dynamic Programming Approach 2 (Sparse Table)
Time Complexity : Construction O(NlogN) Query O(1) <O(NlogN),O(1)>
To get an asymptotically faster time bound, we need to think something out of the box then to just look into the trivial comparisons based algorithms. The next algorithm will take use of a special data structure known as Sparse Tables. Sparse tables stores the information from one index ‘i’ to the some index ‘j’ which is at a specific distance from ‘i’.

Here we use Sparse table to store the minimum of the elements between index ‘i’ to i+’2^j’. It can be better understood with the help of an example :
let us say, A = [2,4,3,1,6,7,8,9,1,7]
and the sparse table be two dimensional Array M( N*(log(N)+1) )

RMQ_003

To compute M[i][j] we use dynamic programming

             M[i][j-1]  if(A[M[i][j-1]]<= A[M[i+2^(j-1)-1][j-1]])
M[i][j] = 
             M[i+2^(j-1)-1][j-1]

Now after precomputation of the table ‘M’ The RMQ can be solved in O(1) as follows :
let k = log(j-i+1)

             A[M[i][k]]  if(A[M[i][k]]<=A[M[j-2^k+1][k]])
RMQ(i,j) = 
             A[M[j-2^k+1][k]]
#include &lt;cstdio&gt;
#define MAXN 1000
#define MAXLOGN 10
int M[MAXN][MAXLOGN];
void Compute_ST(int A[],int N){
	int i,j;
	for(i=0;i&lt;N;i++)
		M[i][0]=i;
	for(j=1;1&lt;&lt;j &lt;=N ;j++){
		for(i=0;i+(1&lt;&lt;(j-1))&lt;N;i++){
			if(A[M[i][j-1]]&lt;=A[M[i+(1&lt;&lt;(j-1))][j-1]])
				M[i][j]=M[i][j-1];
			else
				M[i][j]=M[i+(1&lt;&lt;(j-1))][j-1];
		}
	}
}
int RMQ(int A[],int s,int e){
	int k=e-s;
	k=31-__builtin_clz(k+1); // k = log(e-s+1)
	if(A[M[s][k]]&lt;=A[M[e-(1&lt;&lt;k)+1][k]])
		return A[M[s][k]];
	return A[M[e-(1&lt;&lt;k)+1][k]];
}
// driver programme
int main(){
	int A[10] = {3,1,2,-1,5,4,6,0,9,8};
	Compute_ST(A,10);
	int s,e;
	scanf(&quot;%d %d&quot;,&amp;s,&amp;e);
	printf(&quot;%d\n&quot;,RMQ(A,s,e));
	return 0;
}

5. Segment Trees
Time Complexity : Construction O(N) Query O(logN) <O(N),O(logN)>
segment tree are one of the most popular and most powerful data structures for interval update and interval queries. Segment trees are mostly preferred over any other methods described above and are useful in most of the programming competitions. Segment trees are heap like data structures. The segment tree for an array of size 10 would look like this.
RMQ_004

Segment tree can be stored in the form of an array of size ‘2^(logN+1)-1’ say M.The left and the right child of node numbered ‘x’ will be ‘2*x+1’ and ‘2*x+2’ respectively.

#include &lt;cstdio&gt;
#define MAXN 1000
#define MAXSIZE 10000
int M[MAXSIZE];

//node is the index of the segment tree m, start and end are the index of the the array A
void BuildTree(int node,int start,int end,int A[],int N){
	if(start==end)
		M[node]=start;
	else{
		BuildTree(2*node+1,start,(start+end)/2,A,N);
		BuildTree(2*node+2,(start+end)/2+1,end,A,N);
		if(A[M[2*node+1]]&lt;A[M[2*node+2]])
			M[node]=M[2*node+1];
		else
			M[node]=M[2*node+2];
	}
}

int RMQ(int node,int start,int end,int s,int e,int A[]){
	if(s&lt;=start &amp;&amp; e&gt;=end)
		return M[node];
	else if(s&gt;end || e&lt;start)
		return -1;
	int q1 = RMQ(2*node+1,start,(start+end)/2,s,e,A);
	int q2 = RMQ(2*node+2,(start+end)/2+1,end,s,e,A);
	if(q1==-1)
		return q2;
	else if(q2==-1)
		return q1;
	if(A[q1]&lt;A[q2])
		return q1;
	return q2;
}
// driver programme
int main(){
	int A[10] = {3,1,2,-1,5,4,6,0,9,8};
	BuildTree(0,0,10-1,A,10);
	int s,e;
	scanf(&quot;%d %d&quot;,&amp;s,&amp;e);
	printf(&quot;%d\n&quot;,A[RMQ(0,0,10-1,s,e,A)]);
	return 0;
}

Now, out of the five approaches last two are generally preferred for most of the programming competitions. Segment trees are more popular, because they are more flexible than sparse tables, they can also be used for interval update which is most often clubbed with the interval query in a programming challenges while sparse table can not be used to update the array(as it is based on precomputations).

External Links:
Topcoder Tutorial on RMQ and LCA
Codechef Problem
Sparse Table Solution for the Above Problem
SPOJ Problem BRCKTS
SPOJ Problem GSS3
Spoj Problem HORRIBLE
SPOJ problem GSS2
SPOJ problem SEGSQRSS
SPOJ problem ORDERSET

return 42;

Bridges a.k.a. Cut Edges

Hey folks,
Today we are going to learn that how to compute the edges of importance or essential edges or bridges or cut edges. The bridges are very similar to articulation points that we discussed in the last post.

Bridges :
These are those edges which are essential for the graphs connectivity i.e. If I remove a bridge from the graph the rest of the graph will remain disconnected.

Now how to find the bridges in a graph
Naive’s Approach
remove edges one by one and do a DFS/BFS on the remaining graph if a single DFS/BFS results spanning the full graph than the edge was a waste and it was not a bridge otherwise it is a bridge.

PsuedoCode

FindBridges():
   for e in Edges:
      remove e
      do a DFS on remaining graph
      if the remaining graph is not connected:
            e is a bridge

Time Complexity : O(E*(V+E))

Better Approach
This approach highly resembles with the approach that we discussed to find the articulation points. Here we claim that in a DFS tree an edge is a bridge only when it’s children’s sub-tree doesn’t have a back-edge to one of the parent’s ancestor. Because then the removal of the edge will disconnect the tree.

PsuedoCode

#the 'low' will store the lowest discovery time of the neighbor of the points in it's sub-tree.
#graph(v) is an array of all the vertices attached to vertex 'v'

FindBridges(v):
   time++
   v.discovery_time = time
   v.low = time
   v.used = true
   for i in graph(v):
      child = i
      if(child == v.parent)
         continue
      if(child.used):
         v.low = min(v.low,child.discovery_time)
      else:
         child.parent = v
         FindBridges(child)
         v.low = min(v.low,child.low)
         if(child.low &gt; v.discovery_time &amp;&amp; v.parent!=-1(i.e v is not root)):
            (v,child).is_bridge = true

Time Complexity O(V+E)

External Links
Bridges
bridge searching algo

Articulation Points a.k.a. Cut Vertices

Hey folks,
Today we are going to talk about how to find the number of articulation points in a given connected graph. So let us first write the definition of Cut vertex or Articulation Point.

Cut Vertex / Articulation Point
In simple words cut vertex is the essential vertex for the graph, because if we remove the vertex the remaining graph will not remain connected. So Cut vertices are those vertices in a connected graph whose removal causes the removal of connectivity in the graph.

there are multiple approaches to find the number of articulation points in a graph.

Naive’s Approach
Now you have to find out how essential is the point for the connectivity of the graph then better remove and check it whether the remaining graph is connected or not. So this is how our first algorithm works remove the vertex and do a DFS/BFS on the graph with remaining vertices, if the remaining graph is connected then the vertex is the point of articulation other wise not.

PsuedoCode:

# G[n][n] is the adjacency matrix for the given graph with n vertices
Cut_Vertex(G):
   for i in range(n):
      Remove all edges with one end point in vertex 'i'
      DFS(G)
      if(DFS(G) does not contains all the other n-1 vertices except 'i'):
         i is a cut vertex
  

The above algorithm takes time O(V*(V+E))

Better Approach
So the approach uses a slightly modified DFS on the graph. A point is articulation point in two cases :
1. It is the root of the DFS and have two children.
2. If any one of it’s children have no point in sub-tree of itself which has a back-edge to one of the ancestors of the parent.

The statement 2 may sound somewhat difficult but it is very simple to understand it, It says that if all the children of a point have some point in their sub-tree which has a back-edge to one of the ancestors of the parent then if their parent i.e. the point of focus would not have been there, they could still remain connected with the graph. So the point is not an articulation point or cut vertex.

PsuedoCode

#the 'low' will store the lowest discovery time of the neighbor of the points in it's sub-tree.
#graph(v) is an array of all the vertices attached to vertex 'v'

FindArticulation(v):
   time++
   v.discovery_time = time
   v.low = time
   v.used = true
   no_of_children = 0
   for i in graph(v):
      child = i
      if(child == v.parent)
         continue
      if(child.used):
         v.low = min(v.low,child.discovery_time)
      else:
         child.parent = v
         FindArticulation(child)
         v.low = min(v.low,child.low)
         if(child.low >= v.discovery_time && v.parent!=-1(i.e v is not root)):
            v.is_articulation = true
         no_of_children++
   if(no_of_children>1 && v.parent==-1(i.e. v is root)):
      v.is_articulation = true

External Links
Cut vertices
Biconnected Component
Cut Vertices Algo
My Code
Spoj Problem

return 42;

Euclidean Minimum/Maximum Spanning Tree (EMST)

Hey folks,
Today we are going to learn about Euclidean Spanning Tree problem, which is solved by using prim’s algorithm and not by kruskal’s. The idea behind the fact that the problem of euclidean maximum/minimum spanning tree is solved by prim’s is that the complexity of kruskal’s algorithm is (ElogE) where E is the no of edges, kruskal’s algorithm works fine for most of the general programming competition problems, but in case of Euclidean Spanning Tree. To understand that we have to first learn what is the need of this EMST.
Problem
Actually in euclidean MST problem all the ‘n’ points are connected to each other, and becomes the maximally dense graph. By virtue of that the order to compute the MST by using kruskal’s algo becomes (ElogE) where E = V^2 so order is greater than V^2 and in most cases the judges won’t accept it. So we can use prim’s here and can reduce our order to V^2.

Following is the PsuedoCode for computing EMST (Euclidean Minimum Spanning Tree)
PsuedoCode

#Let G[V][V] contains the weights of all the edges.
#visited[V] is initialized to false and keeps track of all the vertices if they are visited or not.

#min_cost[V] is initialized to INF for every vertex , it keeps the minimum distance of a point not currently in the tree from all the points which are currently in the tree.

EMST():
   root = 0 #any random node
   min_cost[root]=0
   for j in range(n):
      selected_vertex = -1
      for i in range(n):
         if( !used[i] &amp;&amp; (selected_vertex==-1 || min_cost[i]&lt;min_cost[selected_vertex]) ):
            selected_vertex = i
     
      if(min_cost[selected_vertex]==INF):
         return &quot;No EMST&quot;
      visited[selected_vertex]=true
   
      #update the min_cost array for rest of the vertices
      for i in range(n):
         if G[selected_vertex][i]&lt;min_cost[i] :
            min_cost[i]=G[selected_vertex][i]

External Links
Euclidean minimum spanning tree
My Code

return 42;

Computing nCr

Hey folks,
Today we are going to talk about how to compute nCr in the programming competitions.

1. Dynamic programming Approach :
From the high school mathematics we know a recursive method to calculate nCr i.e.
C(n,r) + C(n,r-1) = C(n+1,r)
This can also be written as
C(n,r) = C(n-1,r) + C(n-1,r-1)

It’s a very simple DP, as you can see nothing to be afraid of here.

C[n][r]
/* initialize the two dimensional array by -<strong>INF</strong>
combination(n,r):
   if (C[n][r]!=-<strong>INF</strong>):
      return C[n][r]
   if(r==n || r==0):
      C[n][r]=1
      return
   C[n][r]= combination(n-1,r) + combination(n-1,r-1) 
   return 

Time complexity = O(n*r)
Space complexity = O(n*r)

2. Expand nCr

C(n,k) = n!/((n-k)!k!)
         [n(n-1)...(n-k+1)][(n-k)...(1)]
       = -------------------------------
           [(n-k)...(1)][k(k-1)...(1)]

this can also be written as

C(n,k) = 1,                 if k = 0
       = (n/k)*C(n-1, k-1), otherwise

PuedoCode

combination(n,r):
   if(r==0):
      return 1
   return (n/k)*combination(n-1,k-1)

3. Using Euler’s theorem and modular multiplicative inverse

Modular Multiplicative Inverse :
In modular arithmetic, the modular multiplicative inverse of an integer ‘a’ modulo ‘m’ is an integer ‘x’ such that
ax = 1 (mod m)

Euler’s theorem :
According to Euler’s theorem, if ‘a’ is coprime to ‘m’, that is, gcd(a, m) = 1, then
a^φ(m) = 1 (mod m)

So we can also write it as
a^(φ(m)-1) = a^(-1) (mod m) ; a^(-1) is actually the modular multiplicative inverse of ‘a’

Now in case ‘m’ is a prime number φ(m) will be (m-1).So the equation will turn into
a^(-1) = a^(m-2) (mod m)

the psuedocode below will tell you how to compute nCr for large numbers effectively but there is a constraint that MOD must be a prime number if it is not a prime number, do the prime factorization of MOD and then use Chinese Remainder Theorem to compute the solution.

PsuedoCode

# this function will compute (a^b)%MOD
pow(a,b):
   if(b==0):
      return 1
   if(b%2==1):
      return (a*pow(a,b-1)%MOD)
   temp = pow(a,b/2)
   return temp*temp

Modular_Multiplicative_Inverse(a):
   return (pow(a,b)) 

combination(n,r):
   fact[0]=1
   for i in range(1,n+1):
      fact[i]=(fact[i-1])*i%MOD
   return fact[n] * Modular_Multiplicative_Inverse(fact(n-r)) * Modular_Multiplicative_Inverse(fact(r))

4. Matrix Exponentiation
(useful only for small ‘r’)
Here we have the recurrence relation C(i,k) = C(i-1,k-1) + C(i-1,k).
If we take k=3 we can write,

C(i-1,1) + C(i-1,0) = C(i,1)
C(i-1,2) + C(i-1,1) = C(i,2)
C(i-1,3) + C(i-1,2) = C(i,3)

Now on the left side we have four variables C(i-1,0), C(i-1,1), C(i-1,2) and C(i-1,3).
On the right side we have three variables C(i,1), C(i,2) and C(i,3).
We need those two sets to be the same, except that the right side index numbers should be one higher than the left side index numbers. So we add C(i,0) on the right side. NOw let’s get our all important Matrix.

(.   .   .   .)  ( C(i-1,0) )   ( C(i,0) )
(.   .   .   .)  ( C(i-1,1) ) = ( C(i,1) )
(.   .   .   .)  ( C(i-1,2) )   ( C(i,2) )
(.   .   .   .)  ( C(i-1,3) )   ( C(i,3) )

The last three rows are trivial and can be filled from the recurrence equations above.

(.   .   .   .)  ( C(i-1,0) )   ( C(i,0) )
(1   1   0   0)  ( C(i-1,1) ) = ( C(i,1) )
(0   1   1   0)  ( C(i-1,2) )   ( C(i,2) )
(0   0   1   1)  ( C(i-1,3) )   ( C(i,3) )

The first row, for C(i,0), depends on what is supposed to happen when k = 0. We know that C(i,0) = 1 for all i when k=0. So the matrix reduces to

(.   .   .   .)  ( C(i-1,0) )   ( C(i,0) )
(1   1   0   0)  ( C(i-1,1) ) = ( C(i,1) )
(0   1   1   0)  ( C(i-1,2) )   ( C(i,2) )
(0   0   1   1)  ( C(i-1,3) )   ( C(i,3) )

And this then leads to the general form:

               i
(.   .   .   .)  ( C(0,0) )   ( C(i,0) )
(1   1   0   0)  ( C(0,1) ) = ( C(i,1) )
(0   1   1   0)  ( C(0,2) )   ( C(i,2) )
(0   0   1   1)  ( C(0,3) )   ( C(i,3) )

For example if we want C(4,3) we just raise the above matrix to the 4th power.

              4
(1   0   0   0)  ( 1 )   ( 1 )
(1   1   0   0)  ( 0 ) = ( 4 )
(0   1   1   0)  ( 0 )   ( 6 )
(0   0   1   1)  ( 0 )   ( 4 )

for k=2 matrix will be

(1  0  0)
(1  1  0)
(0  1  1)

after you go through two or three more examples of k, you will realize a pattern in the matrix i.e.
to compute nCk you have to make a [k+1]*[k+1] matrix which will contain a single ‘1’ in first row and all the following k rows will have two ‘1’s shifted toward left by one at every subsequent row.
PsuedoCode :

FormMatrix(n,k):
   for i in range(1,k+1):
      mat[i][i-1]=1
      mat[i][i]=1

the base column matrix varies for different values and so as the main matrix. So here we can see that the complexity is r^3 * log(n) .. So the algorithm is useful only for very small ‘r’.

External Links
Euler’s Theorem
Modular Multiplicative Inverse
My Code for Euler’s theorem method

return 42;