about Josephus problem

Problem: Josephus problem is a classical problem which can be solved by DP. The problem assume n persons in a circle, then start from the 1st person, we eliminate a person after skipping a fix number of persons. The process stops when there is only one person left. We want to know the last person's index in the original circle. For example, if n = 5 and k = 2, at the beginning, we have 1 2 3 4 5, then we eliminate 2 and 4, then we eliminate 1 and 5, so the last person's index is 3.

Solution: For the general case where k is an arbitrary number, we use F(n, k) to denote the solution for n persons. Then we have F(n, k) = (F(n-1, k) + k ) mod n. Why? we can think the problem as follow:

• The first person we need to eliminate is the kth person, after he is removed, we have n-1 persons now.
• However, we can not directly resort to the subproblem F(n-1, k), since now our process will start at the k+1th person at the original circle. Therefore we need to remap the indices. Since the process is applied on a circle, we can move the first k-1 persons (1st, 2rd, ... k-1th person at the original circle) to the end of the circle. Thus, assume the new index in F(n-1, k) is i, the corresponding index in F(n, k) is (i + k) mod n.

The above approach will take O(N) time. For smaller k, actually O(logN) solution exists. For example, when k = 2, we have F(2n) = 2*F(n) -1 and F(2n+1) = 2*F(n) +1. Actually for k = 2, even an O(1) solution exists: if n = 2^m + l where  0 =< l  < 2^m, the last person will be 2*l + 1.

Calculate Number of Different Ways to Climb a Stairway with Constant Space and Linear Time

Problem: There are three ways to climb a stair, you can climb 1 step, 2 steps or 3 steps each time. Given a n step stairway, how many of ways can you climb up to the top?

Solution: As known, this is a simple DP problem.  If we can maintain an array dp[] with dp[i] representing the different ways to climb up i steps, we can get the final result in O(N). However, we need O(N) space as well. But actually, since dp[i] = dp[i-1] + dp[i-2] + dp[i-3], we don't need O(N) space, instead, we just need a cache with size 3. Then we can have a space complexity O(1) algorithm:

int climb_stair(int n)
{
  int cache[3];
  cache[0] = 1; /*ways to climb 0 step */
  cache[1] = 1; /*ways to climb 1 step */
  cach2[2] = 2; /*ways to climb 2 steps */

  if(n<=2) return cache[n];

  for(int i=3; i<=n; i++)
  {
     cache[i%3] = cache[0]+cache[1]+cache[2];
  }

  return cache[i%3];

}

Similarly, we can have a O(1) space algorithm for Fibonacci numbers.

User Bin Crash problem

Problem: This is a problem from facebook buzz. The problem is that a plane is gong to crash,  you need to throw out W lbs of freight to survive the crash. Each type of freight fi is associated with a weight wi and price pi. If we throw one unit of fi out, we lost pi but the plane becomes lighter by wi. Our mission is to throw out minmum freights measured by their aggregate price and avoid the crash, which means the total weight of those freights should be larger than W.

Solution: The problem is very similar to the 0-1 knapsack. We can also use similar strategy to solve it. The detail is as follow:

• Basically, we need first sort different types of freight by their weight. 
• Then we will need to fill an array dp[n][W+1] where n is the number of different types of freights. dp[i][j] represents the minmum cost to throw out at least j lbs freight and we can throw up to ith type of freights (ranked by their unit weight and from light to heavy). 
• To fill dp[i][j], we need first select the most cost/efficient type of freight among the i types of freights. Cost/efficient is measured by the unit weight of the freight divided by its unit price. Say the most cost/efficient type of freight among the i types of freights is kth type of freights, we have if j <= kth type's unit weight, dp[i][j] = min(dp[i-1][j], kth type's unit price); Otherwise,  dp[i][j] = min(dp[i-1][j], dp[i][j-kth type's unit weight] + kth type's unit price). We don't need to use greedy method, just do regular dp. dp[i][j] = min(dp[i-1][j], dp[i][j-ith type's unit weight] + ith type's unit price).
• The time complexity of this approach is O(nW).

Convert an Array into a Sorted Array with Minimum Cost

Problem: Given an array of positive integers, and you only have two operations : a) decrease an integer by k  with a cost of k; b) delete an integer with a cost of the value of that integer. Try to use these two operations to convert the array into a sorted array with minmum cost. For example, a[] = {5,4,7,3}, the optimal way is to decrease "5" by 1 and delete "3" with the total cost of 4 (1+3). The resulting sorted array is {4, 4, 7}.

Solution: Here I give a O(N^2) solution by using DP. The key data structure we need is dp[i] which records the minimum cost to covert a subarray a[0, i] into sorted array and the last element of the sorted array should be a[i].  For the example array a[] = {5,4,7,3}, dp[0] = 0, dp[1] = 1 (since we need to decrease "5" by 1), dp[2] = 1, dp[3] = 7. After we had filled up dp[], we need one more data structure: aggr[]. aggr[i] = a[0] + a[1] + ... + a[i]. Then the minmum cost to convert the whole array cost_min = min{dp[i] + aggr[N-1] - aggr[i]}. Here  "aggr[N-1] - aggr[i]" means the cost that we  delete all the elements with a subscript larger than i.

Then the key is to calculate dp[i]. There are several observations here:

1. if a[i] >=  a[j] (j<i),  we can construct a sorted array ended at a[i] by deleting a[j+1], a[j+2], ... , a[i-1] and then append a[i]. The cost of the resulting sorted array is dp[j] + the cost of deleting a[j+1], a[j+2], ... , a[i-1]. 
2. if a[i] < a[j], to have a[i] as the last element, we need to decrease a[j] by a[j] - a[i]. Moreover,  we may not stop at a[j] and we need to keep look backward until we found a a[m] such that a[m] <= a[i], then the cost of the resulting sorted array is dp[m+1] + ∑(a[k] - a[i]) where m<k<j. 
3. To compute dp[i], we need to look at all the j such that  j<i, depending on a[j] and a[i], we need to calculate the cost based on step 1 or 2. Then the minmum cost we encounter is the value of dp[i].

The code is as follow:

int min_convert_cost(int a[], int n)
{
    int dp[n];
    int aggr[n];

    aggr[0] = a[0];
    for(int i=1; i<n; i++)
        aggr[i] = aggr[i-1] + a[i];

    dp[0] = 0;
    for(int i=1; i<n; i++)
    {
       dp[i] = INT_MAX;
       for(int j=i-1; j>=0; j--)
       {
          int cost_i = 0;

          if(a[i] >= a[j])
          {
             cost_i = dp[j] + aggr[i-1] - aggr[j];
          }
          else
          {
              cost_i += aggr[i-1] - aggr[j];
              while(a[j]>a[i] && j >=0)
              {
                cost_i += a[j]-a[i];
                j--;
              }
              
              cost_i += dp[j+1];

          }
          if(cost_i < dp[i])
               dp[i] = cost_i;

       }
    }

   int min = INT_MAX;
   for(int i=0; i<n; i++)
   {
      int cost_i = dp[i] + aggr[n-1] - aggr[i];
      if(cost_i<min) min = cost_i;
   }

   return min;
}

An O(N^2) Algorithm for Optimal BST

Problem: Given a set of ordered keys and their corresponding probabilities to be visited, construct an optimal BST tree with minimum search cost.

Solution: This is the classical optimal BST problem. The well-known solution has an O(N^3) complexity. However, Knuth found a better way that can achieve O(N^2). Basically, when we try to find the root of an optimal BST, previously we need to iterate from i to j. But there is an observation: root[i, j-1] <= root[i,j] <= root[i+1, j]. Given this inequality, we can reduce the range which we look for the root. This new range (root[i, j-1]], root[i+1, j]) is smaller than (i, j). The detail solution and proof can be found here.

Find the Largest Container in a Histogram

Problem: Given a histogram, pick two bars in the histogram as the left and right sides of a container. X axis is consider as the bottom. Then we have a container which can hold water. Find the container that can hold the largest volume of water. The container must be placed horizontally (you can't rotate container).

Solution: We can still use DP to solve this problem in O(n). Basically, there are two factors decides the volume of the container: 1) the minimum of the two sides 2) how wide the container is. If the highest bars are just at the left and right end of the histogram, we know that by picking these two bars we can have the largest container. If the bars at the left and right end are not the two highest, we need to explore further to look at other combinations. Based on these observations, what we need is actually an increasing sequence from left to right and  also an increasing sequence from left to right. The details of the algorithm is as follow:

• First model the histogram as an array h[i]. For example, h[] = {3, 1, 4, 7, 5, 2, 6}. 
• Then find the two increasing sequences. The one from left to right is left[] = {3, 4, 7} and the one from right to left is right[] = {6, 7}.
• begin with the first elements in left[] and right[], calculate the volume of the container formed by these two bars, use max_v to store this volume (which is 3*5=15). 
• Then select the smaller one of the two bars, make one advancement in the array which the smaller bar is from. For our example, between left[0] and right[0], left[0] is the smaller one. Therefore we advance to left[1], which is "4". Then we calculate the volume of the container formed by left[1] and right[0]. The volume is 4*3=12, which is smaller than the previous value 15, so we just keep the old value. 
• Then we repeat the previous step, just advance the array which the smaller bar is from. We will get   left[2] and right[0]. The corresponding volume is 6*2 = 12, which is still smaller than 15.
• Then we need to advance in right[], we will find right[1] and left[2] refer to the same bar. Since at this time, both left[] and right[] have been exhausted, our algorithm just abort.
• One more thing, when two bars  left[i]  == right[j], we need to advance both arrays and inspect   left[i+1] and right[j+1] next.

Make Best Use of the Conference Room

Problem: Given a conference room and a number of presentations with start and end time ( e.g., [4, 9], [5, 10]), try to make an arrangement which allows the conference room to be used for maximum time. Overlapping presentations can't be in the same arrangement.

Solution: We can have a O(nlogn) solution by using DP.  The details are as follows:

• Sort the presentations by their end time. Thus we will have a sorted array end[N].   N is the number of the presentations.
• Have an array Max_arr[N].  Max_arr[i] stands for if we take end[i] as the close time for the conference room, the maximum time the room can be used. It is easy to know that  Max_arr[0] = the duration of the presentation that ends at  end[0]. We need to find out Max_arr[N-1].
• To calculate Max_arr[i], we first get the start time si of the presentation that ends at end[i]. Then we do binary search in end[0] ... end[i-1] for si. Basically we need to find a j such that  end[j] < si &&  end[j+1] > si. Therefore,  Max_arr[i] = max(Max_arr[i-1],  Max_arr[j] + end[i] -  si ).

Zeckendorf's Theorem -- How to Represent a Positive Integer with Fibonacci Numbers?

Zeckendorf's Theorem is about that every positive integer can be represented in the form of the sum of one or more Fibonacci numbers.

For example, 6 = 3 + 2 + 1, 10 = 8 + 2. We can further represent these positive integers in the binary form of Fibonacci numbers. Thus, 6 can be represented as "111" which means fibo(1) +  fibo(2) + fibo(3) -- the sum of the first, second and third Fibonacci numbers. Similarly, 10 can be represented as "10010" which means fibo(2) +  fibo(5)  --  the sum of the fsecond and fifth Fibonacci numbers.

How to find the Fibonacci numbers that sum to a particular positive integer? Basically you can use DP and leverage the binary form.

• 1 can be represented as "1", 2 as "10", 3 as "11", 4 as "101", 5 as "110". Thus, we can observe the binary form of a positive integer i is just the immediate next number of the binary form of a i-1. However, this "immediate next" is slightly different when an integer has a binary form of all "1". For example, 3 is "11". The  "immediate next" of "11" is not "100" but "101". This is due to the nature of Fibonacci numbers. We know "100" = "11" since fibo(n) = fibo(n-1) + fibo(n-2). So here we need to add one more "1" to "100".
• Based on the above analysis, we can get those Fibonacci numbers for a particular positive integer by using DP.

Actually not only positive integers, negative integers can also be represented in the sum of "broadly defined" Fibonacci numbers.

Dynamic Programing Solution to Subset Sum Problem

Problem: Given a set of integers, find a subset of them which sum to a target value s.

Solution: We can definitely first sort the set, then use some recursive approach to solve the problem. The complexity is expected to be exponential. Here we give a DP approach, though the complexity is also not polynomial.

• First we need to define the state. We have a state function Q(i,k) which means if we can find a subset from the set of integers a1,a2, ..., to ai which sum to k. So Q(i,k) is a Boolean function.
• We have  Q(1,k) = ( a1 == k);  Q(i,k) =  Q(i-1,k)  or  ( ai == k) or  Q(i-1,k-ai).
• Then we just need to fill in the matrix of  Q(i,k). We know 0< i <= the size of the set. For k, we have   N=< k <=P, where N is the sum of all the negative integers in the set and P is the sum of all the positive integers in the set.

Enumerate the Ways to Cover a M*N Matrix with Dominoes

Problem: Given a M*N matrix, enumerate all the different ways to cover it with dominoes (1*2 blocks). For example, for a 1*2 matrix, there are only one way; for a 2*2 matrix, there is 2 ways. The original problem is from poj 2411.

Solution: We can definitely use backtracking, but more efficient approaches should be DP. First we need to think what are the states and how to represent them. Obviously, the arrangement of dominoes on the matrix is the state. There might be multiple ways to represent them. Here we just introduce one approach: every cell in the matrix is either 0 or 1; "0" means the domino block is placed horizontally and "1" means the domino is placed vertically. If a cell, say a[i,j] is marked as "1", a[i+1,j] can't be marked as "1". Then we can use the binary sequence to represent the status of one row. For example, "0000" means two dominoes lay horizontally. Besides, there are several key observations:

• For each row, there must be even number of consecutive "zero"s. Otherwise, we cannot fully cover a row.
• The status of current row only depends on the previous row. It means we only need to look back one row.
• Then we build a status matrix d[i][s] with i as the row number and s as a particular state. d[i][s] means the number of ways to arrange a matrix with current state (all i-1 rows are filled and the ith row is with state s). Therefore,  d[i][s] =sum_{all possible d[i-1][k]}. 
• How to decide  d[i-1][k], first, s&k should be zero, which is to conform to the constraint " a[i,j] is marked as "1", a[i+1,j] can't be marked as "1". Second, s|k should be a valid arrangement, which is to conform to the constraint " there must be even number of consecutive zeros". Then the cell d[M][0] stores the answer to our problem.

typedef long long LL;
const int maxn = (1 << 12);

LL h , w , dp[12][maxn];

//to decide if it contains odd number of consecutive "zero"s
bool judge(LL s)
{
     LL cnt = 0;
     for(LL i=0; i < w; i++)
     {
         LL t = s & 1;
         if(t)
         {
             if(cnt & 1)
                 return false;
             else
                 cnt = 0;
         } else
             cnt ++;
         s >>= 1;
     }
     if(cnt & 1)
         return false;
     return true;
}

int main()
{
     while(~scanf("%lld %lld", &h, &w))
     {
         if(h == 0 && w == 0)
             break;
         memset(dp , 0 , sizeof(dp));

         //initialize the state in row 1
         for(LL i=0; i < (1 << w); i++)
         {
             if(judge(i))
             {
                 dp[1][i] = 1;
             }
         }
         for(LL i=2; i <= h; i++)
         {
              for(LL j=0; j < (1 << w); j++)
              {
                 for(LL k=0; k < (1 << w); k++)
                 {
                     LL t = j & k;
                     if(t)
                          continue;
                      //so here "t" is the actual state 
                      //on the current row, since if 
                      //the position in the above 
                      //row is marked as "1", the 
                      //position in current row must 
                      //be "0", otherwise this 
                      //iteration will be skipped by 
                      //previous "continue". But, 
                      //even though it is "0", we 
                      //need to count it as "1", 
                      //since it cannot hold domino.
                      //This is what "t=j|k;" 
                      //really means.

                      t = j | k;
                     if(judge(t))
                      {
                         dp[i][j] += dp[i - 1][k];
                      }
                 }
              }
          }
         printf("%lld\n", dp[h][0]);
     }

}

The code is from xIao.wU.

The Longest Palindrome Substring (Manacher's algorithm)

Problem: Given a string, find a longest palindrome substring.

Solution: We can use general suffix tree that stores the original string and its reverse, which is an O(N) algorithm. However, here we give a better one with less space overhead while still O(N) complexity. This algorithm is called Manacher's algorithm. If we check a string from left to right, we can leverage the palindrome check we did previously. This is from the symmetry of palindrome. The main idea is as follow:

• Create an array called P[], P[i] stands for the longest palindrome centered at location i. Here i is not the index in the original string. For the original string, the locations we need to check for palindromes contains those characters in string along with the spaces between characters. So if we have a string of length l, we need to have a P[] with length 2*l+1.
• Our goal is to fill in P[]. For a particular position, we check its left and right. If equals, we extend our check further. Otherwise, the longest palindrome centered at location is found.
• However, we need to be smarter. Actually we can leverage previous computed P[i] when we calculate a P[x] where x>i. 
• So here we add two pointers, p1 and p2, which point to the left and right of the current location i such that |i-p1| = |i-p2| and p2>i>p1. We know p1 refers to a palindrome t and i refers to a palindrome s. If the first character of t is strictly on the right of the first character of s, we know P[p2] = P[p1].
• Otherwise, say if the first character of t is not strictly on the right of the first character of s, we have P[p2] >= r - p2. where r is the right bound of the palindrome that centered at i. We then need to check if the palindrome at p2 can be longer than r - p2. The good thing is that we only need to start the characters beyond the length of r - p2.
• When the first character of t is strictly on the right of the first character of s, we don't need to move the current center (i). Only when the first character of t is not strictly on the right of the first character of s, we need to move the current center to p2.
• The total cost is O(N).

The code is as follow:

void manacher(const string &s)
{
    int len = s.size();
    if(len == 0) return;

    int m[2*len+1];
    m[0] = 0;
    m[1] = 1;
    // "cur" is the current center
    // "r" is the right bound of the palindrome
    // that centered at current center
    int cur, r;
    r = 2;
    cur = 1;

    // iterate from 2 to 2*len+1
    for(int p2=2; p2<2*len+1; p2++)
    {
        int p1 = cur- (p2-cur);
        //if p1 is negative, we need to 
        //move "cur" forward

        // re-adjust cur based on p2
        while(p1 < 0)
        {
           cur++;
           r = m[cur] + cur;
           p1 = cur- (p2-cur);

        }

        // If the first character of t is 
        // strictly on the right of the 
        // first character of s

        //

        // Or here, from the symmetry, if

        // the palindrome centered at cur

        // cover the palindrome centered at

        // p1, we know
        if(m[p1] < r - p2)
            m[p2] = m[p1];
        //otherwise
        else
        {

           // we need to explore the length of

           // the palindrome centered at p2

           // if the palindrome centered at cur covers

           // p2, we can start at "k = r-p2"

           // otherwise, we start at "k=0"

           //reset "cur" 
           cur = p2;
           int k = r-p2;
           if(k<0) k = 0;
           while(1) 
           {
              if((p2+k+1)&1)
              {
                if(p2+k+1 < 2*len+1 && p2-k-1 >=0 && s[(p2+k)/2] == s[(p2-k-2)/2])
                  k++;
                else break;
              }
               else
              {
                if(p2+k+1 < 2*len+1 && p2-k-1 >=0)
                  k++;
                else break;
              }

           }
           // set the right boundary to be "p2+k"
           r = p2+k;
           m[p2] = k;
        }


    }

 
}

Find the Longest Sub-sequence that is a Palindrome within a String

Problem: Given a string, you can delete any characters, find the longest sub-sequence (the characters remained after your deletion) that is a palindrome.

Solution: For palindrome problem, one trick often used is to reverse the string. Here we first reverse the string, then find the longest common sub-sequence between the new string and the original one. It is a O(n^2) solution. Remember, there could be multiple longest common sub-sequences, some of them may not be palindrome, you need do some checks.

0-1 Knapsack Problem

Problem:  Given a knapsack that can hold at most W weight items, also given a list of items with their weight wi and value vi (no items share the same weight), try to find a valid assignment which achieve the highest value in the knapsack (can't  be over-weighted at the same time).

Solution: We can use DP to solve this problem. However, one-dimension DP is not enough. If we only record state s[0], s[1], ... s[W], the later state may not be able to reuse previous states. Instead, we need a two-dimension DP here:

• First sort the items based on their weights
• The status we want to calculate is s[i, w], which means if the total weight is w and we can only use up To the ith item (based on weight and non-descending), the optimal maximum value we can get.
•  s[0, w] = 0 and  s[i, 0] = 0. 
• For  s[i, w], if wi > w,  s[i, w] =  s[i-1, w]; otherwise,  s[i, w] = max ( s[i-1, w], s[i-1, w-wi] + vi).

Volume of Water Held by a Histogram

Problem: Pour water on a histogram, calculate the water volume the histogram can hold.

Solution: This is a relatively simple DP problem. Here we only give the main idea.

• For a particular bar bi, if we know the highest bar on its left Li and highest bar on its right Ri.  If the height of bi is smaller than both Li and Ri, the water volume can be held on this bar is min(Li, Ri) - hi; otherwise, it can't hold water.
• To calculate Li and Ri, we just need to record the maximum height we had observed so far from the left (and from the right). Therefore, a O(n) algorithm is straightforward here.

Maximum Rectangle Area within a Histogram

Problem: Given a histogram, find the maximum rectangle area within it.

Solution:  It is easy to find a O(n^2) algorithm by comparing the height of a bar with the rest bars. Here we give a O(n) algorithm.

• For each bar bi, we need to know the number of adjacent bars on the left Li and on the right Ri that are higher than it. Then the maximum rectangle within the histogram with the height hi will be  (Li+Ri+1)*hi. Then we just need to select the largest one.
• When deciding the number of adjacent bars on the left of bi, the key observation is that we don't need to inspect every bar on its left. The key here is to use a stack to track the bars that had been inspected in a smart way. The bars in the stack are in ascending order (from base to top). Besides, before pushing new bar, we need to pop the bars in the stacks that are higher than it. Then we can achieve calculate Li in O(n). Calculating Ri will be the same.

The following only shows how to calculate Li:

//bar[] represents the height of each bar in the histogram
//left[] stores the number of adjacent bars that are taller than bi 
//this stack is used to track inspected bars
stack s; 

for(int i=0; i<N; i++)
{

  int num = 0;
  int idx;
  while(!s.empty())
  {
    idx = s.top();
    if(bar[idx]>bar[i])
    {
       num++;
       s.pop();
       left[i]+= left[idx];  
    }
    else break;
  }

  s.push(i);
  left[i]= num ? num + left[i]:0;

}

Find the Closet Pair of Points

Problem: Given N points on a cartesian ordinate plane, find the two points with the shortest distance.

Solution: We can use the approach that solving P(n) based on P(n-1). We first sort the points based on their X coordinates. Then assume we have solved the problem with the first n-1 points and the shortest distance is d. Then given the nth point, we don't need naively compute the distance between the nth point and all the previous n-1 points. We can leverage the important information we have obtained: "d". If we can further find a closer pair  given the nth point, the other point in the pair should be in a specific area. The area is a rectangle with dimension (d*2d) that sits on the left of the nth point. Moreover, the number of potential points that could be in this area is limited, at most six points (which can be proved using pigeonhole principle). Therefore our algorithm comes as follows:

• two important data structures are used. The first one is an array X_array that contains all the points sorted by their X coordinates. The second is a BST Y_BST that contains the active points sorted by their Y coordinates. Insert the first two points into X_array and Y_BST. Initialize d_min as the distance between the first points.
• using a sweeping line method. Have a vertical line sweep from left to right and start from the third point in X_array.
• when the line touch a point p, given the current d_min,  first we need to remove the points in X_array that have distance larger than d_min to p. We also need to remove these points in Y_BST. Then we add p to Y_BST and adjust the active region within X_array accordingly.
• search the points in the range [p.y-d_min, p.y+d_min] in Y_BST, for each point in the range, calculate its distance to p. If smaller than d_min, update d_min.
• When all the points have been swept, the algorithm finishes.
• The complexity for the whole algorithm is O(nlogn). Why? The initial sorting cost for X_array is O(nlogn).  For each points being removed from the active region in X_array, it takes O(logn). So the total cost for removing points from the active region in X_array  is O(nlogn). Similarly,  the total cost for removing points from  Y_BST is also O(nlogn). Since there are limited number of points in the rectangle with dimension (d*2d). The cost to calculate the distance between them and p is constant. Therefore, the complexity for the whole algorithm is O(nlogn).

typedef pair<double, double> Coordinate;

bool compX(Coordinate p1, Coordinate p2)
{
   if(p1.first != p2.first)
        return p1.first < p2.first;
   else return p1.second < p2.second;
}

bool compY(Coordinate p1, Coordinate p2)
{
  if(p1.second != p2.second)
    return p1.second < p2.second;
  else return p1.first < p2.first;
}

double cal_distance(Coordinate p, Coordinate q)
{
   double a = abs(p.first-q.first);
   double b = abs(p.second-q.second);

   return sqrt(a*a + b*b);
}

double find_closest(vector<Coordinate> &co)
{

   //sort based on X coordinates
   sort(co.begin(), co.end(), compX);

   if(co.size() == 0) return -1;
   if(co.size() == 2) return cal_distance(co[0], co[1]);

   double d_min = cal_distance(co[0], co[1]);
   //inition the y table
   bool(*fn_pt)(Coordinate,Coordinate) = compY;
   set<Coordinate, bool(*)(Coordinate,Coordinate)> y_table(fn_pt);
   y_table.insert(co[0]);
   y_table.insert(co[1]);

   //start from the third point from the left
   int start = 2;

   vector<Coordinate>::iterator tp;
   vector<Coordinate>::iterator it_head = co.begin();
   vector<Coordinate>::iterator it_tail = it_head+1;

   for( ; start<co.size(); start++)
   {
       double l_bound = co[start].first - d_min;
       Coordinate lc(l_bound, 0);
       tp = lower_bound(it_head, it_tail, lc, compX);

      //delete points from y_table
      vector::iterator pp, endp;
      if(tp == it_tail)
      {
         if((*tp).first > (*it_tail).first)
           tp = it_tail+1;
      }

      for( pp = it_head; pp!=tp; pp++)
      {
         y_table.erase(*pp);
      }
      it_head = tp;
      it_tail++;
      y_table.insert(co[start]);

      //select points from y_table
      double y_low = co[start].second - d_min;
      double y_high = co[start].second + d_min;

      //search for these two bounds
      set<Coordinate>::iterator sp, sp_low, sp_high;
      Coordinate y_lc(0, y_low);
      Coordinate y_hc(0, y_high);
      sp_low = lower_bound(y_table.begin(), y_table.end(), y_lc, compY);
      sp_high = lower_bound(y_table.begin(), y_table.end(), y_hc, compY);

      //calculate the distance and compare
      for(sp=sp_low; sp!=sp_high; sp++)
      {
        if(sp->first == co[start].first && sp->second == co[start].second)
             continue;
        double new_d = cal_distance(co[start], *sp);
        if(new_d < d_min) d_min = new_d;
      }

   }

   return d_min;
}

Alternative approach: We can also try to use K-d tree, more precisely, 2-d tree. K-d tree is basically a binary tree. It is constructed by recursively partition the points based on a particular axis (dimension). For example, for a 2-d tree, we first partition by x-coordinate (it is also ok to partition by y-coordinate), given all points, we find the points which x-coordinate is the median of the x-coordinates of all other points. Then we take the points which has this x-coordinate as root. The remaining points are divided into two groups: a) those which x-coordinates are smaller than root.x b) those which x-coordinates are bigger than root.x. Then we recursively partition these two groups, but by y-coordinate. This will help us find the nodes on the second level (assume root is the first level). To find the nodes at the third level, we partition by  x-coordinate. Repeat this until we can't partition any more.

Then we need to apply the Nearest Neighbor Search (NNS) in 2-d tree. Here I just briefly explain the algorithm.

• First, we need to locate the query point  query point q. This process is similar to a binary search based on  q's coordinates. We stop when we encounter a leaf node. Calculate the distance between this leaf node and  q. This is the current best (nearest neighbor).
• However, we need to move forward since the current best may not be optimal. We just rewind the recursion, go to the parent of the leaf node.  Calculate the distance between this (parent) node and  q. If necessary , update the current best. We also need to see if there remains some nearer neighbors in the other plain divided by this (parent) node. This can be done by checking if the circle centered at  q with the radius current best across the splitting line. If yes, we need to check the other branch of this node. Basically, we repeat previous steps, searching q in that subtree.
• If we are done with the (parent) node, we just keep going up, check the parent of this (parent) node, until we encounter the root node.

The empirical performance for a NNS in 2-d tree is O(logn). Then do the NNS for all the nodes is O(nlogn). Considering the cost for building 2-d tree is O(nlogn), the total cost for this problem is also  O(nlogn).

See Also: The closest pair problem

Find the Sub-Matrix with the Largest Sum in a Matrix

Problem: Given a m*n matrix that is made up of integers (positive and negative), find the sub-matrix with the largest sum.

Solution: There is an O(n) algorithm to find the sub-array with the largest sum in a one-dimension array. So here, we will try to reduce our problem to the one-dimension array problem.

1. Assume the sub-matrix is between row a and b where 0<=a<=b<=m,  then we discard all the rows that are out of the range [a, b]. We then squeeze the sub-matrix into a one-dimension array. How to squeeze? Just use the sum of the elements remaining in each row to replace that row! Therefore we get a one dimensional array with n elements (each element is a sum). 
2. Then apply the O(n) algorithm we mentioned previously, we will get the sub-array with the largest sum in this one-dimension array, which is actually the sub-matrix between row a and b with the largest sum.
3. Until now, we have examined one pair a and b. If we try all the combination of a and b, we will be able to find the the sub-matrix in the with the largest sum in the m*n matrix. Assume n<m, we need repeat n^2 such operation. Then the overall time complexity is O(n^3). Naive way may take O(n^4). 
4. In order to make the "squeeze" operation (summing the column elements) with a constant time complexity, we need to do some pre-processing. We just need to calculate the prefix-sum. For example, c[i][j] = a[0][j] + a[1][j] +...+ a[i][j] . It is easy to calculate this for the whole matrix with time complexity O(n^2).

Find the Kth Positive Integer That Can Be Divided Only by 3, 5 or 7

Problem: find the kth positive integer that only have factors as 3, 5, or 7. For example, the first integer is 3, then 5, 7, 9, 15, 21......

Solution: It is a DP problem where we should leverage on previous status. Then we can avoid redundant computation. The main ideas are stated as follow:

1.  Have three queues, let us say Q3, Q5 and Q7. Put 3, 5 and 7 in these queues, receptively.
2. Then select the smallest element among the three queues by just comparing the heads of the queues. If the smallest element s is from Q3, Q3.enque(s*3), Q5.enque(s*5) and Q7.enque(s*7); If the smallest element s is from Q5, Q5.enque(s*5) and Q7.enque(s*7);   If the smallest element s is from Q7, Q7.enque(s*7). Since we strictly follow the ascending order, we guarantee that there is not redundant computation.
3. Repeat the selection (step 2) for k times, what we get is the kth positive integer that only have factors as 3, 5, or 7.
4. The time complexity and space complexity are both O(k).

Longest Common Subsequence

Problem: Longest common subsequence is the subsequence that appear in multiple sequences (usually two) simultaneously. For example, the longest common sequence of "13486" and "23861' is "386". The diff program is based on this algorithm. Besides, this algorithm is frequently used in bioinformatic research.

Solution: Still we will use DP, but different from longest Increasing subsequence, the solution we have is O(n*m), where n and m are the length of two sequence. If we have more than two sequence, the time complexity will be O(n1*n2*n3*...). The following gives the algorithm for two sequences:

1. look at the last elements A[n], B[m] of two sequences, if A[n] == B[m], we know the longest common subsequence (LCS) between  A[n] and B[m],  LCS(A[1:n], B[1:m]) = LCS(A[1:n-1], B[1:m-1]) +1. 
2. if A[n] <> B[m], then  LCS(A[1:n], B[1:m]) = MAX(LCS(A[1:n-1], B[1:m]), LCS(A[1:n], B[1:m-1])).
3. then we will have a matrix for LCS, after filling up the matrix, we will get the answer.

We can further do some optimization.  One is to optimize the space, the LCS matrix. If most diffs center on the middle of sequence, we can trim the common part first.  Also to be fast, we don't need to compare each character, we can hash the string or the line. Then we just compare the hash value sequences.


More: The solution to LCS problem can be easily applied to Edit Distance problem. Given two string Xi and Yj, and a number of operations such as copy(), delete(), insert(), etc. Find the optimal way to transform Xi to Yj. Each operation has a cost and our goal is to minimize the total cost of the transformation. Here only gives the recursive relations:

1. Define T(i,j) as the total cost to transform Xi to Yj (one direction). There are several recursive relations.
2. T(i,j)  = T(i-1, j-1) + cost(copy), if X[i] == Y[j]
3. T(i,j)  = T(i-1, j-1) + cost(replace), if X[i] <> Y[j]
4. T(i,j)  = T(i-1, j) + cost(delete)
5. T(i,j)  = T(i, j-1) + cost(insert)
6. find T(i,j)  as the minimum T(i,j) in step 2,3,4,5

Longest Increasing Subsequence

Problem: Longest increasing sequence is the subsequence within a sequence that the elements in the subsequence is in an increasing order. We are trying to find the longest subsequence that has this property. For example, one of the longest increasing subsequences in "17385" is "138".

Solution: This is a classical DP problem. Naive DP solution will give you a O(n^2) algorithm by keep the longest increasing sequence ended at i. A better solution is to use DP alongside with binary search, which will give you a O(nlogn) algorithm. The solution is as follow:

1. having an array M, M[j] stores the ending position of the increasing sequence with length j in the original sequence X. Thus, X[M[j]] should be an increasing sequence as well.
2. iterate the original sequence X, when visiting X[i], do binary search in X[M[j]] to look for the closest position to X[i]. If X[i] is bigger than all the existing X[M[j]], it means if we append X[i] to the end of the current longest increasing sequence, we get a new longest increasing sequence which is longer by 1. Then we update M[L+1] = i (L is the length of current longest increasing sequence).
3. when searching in  in X[M[j]] for X[i], we can also encounter the situation that  X[M[j]]<X[i]< X[M[j+1]]. Then we need to update M[j+1] = i. The reason is that after we have scanned the ith element in X, we need to assure for M[j], we store the increasing sequence of length  j that has the smallest ending element. This will maximize our chance to be able to append further elements to existing increasing sequence.