Pages

Showing posts with label topcoder. Show all posts
Showing posts with label topcoder. Show all posts

Thursday, August 12, 2010

TC SRM 487 D1-500 KiwiJuice

Category: Dynamic Programming
URL: http://www.topcoder.com/stat?c=problem_statement&pm=11019


class KiwiJuice {
public:
   int theProfit(int, vector <int>, vector <int>);
};

int dp[51][1<<16];
vector<int> data;
vector<int> price;
int _C;

// f(current volume of the special bottle, bitmask of the used bottles)
int f(int i, int mask) {
   // caching mechanism
   if (dp[i][mask] != -1) return dp[i][mask];
   // base case
   if (mask+1 == (1 << (data.size()))-1) {
      return price[i];
   }
   int res = 0;
   for (int j = 1; j < data.size(); j++) {
      if (mask & (1 << j)) continue;   // already used
      // the next bottle to process
      int next = data[j];
      int d1 = 0;
      if (i + next >= _C) {
         // one of the bottles is full, the other is i+next-_C volume
         d1 = price[_C] + f(i+next-_C, mask | (1 << j));
      } else {
         // change to i + next (and used up the empty for bottle[j])
         d1 = price[0] + f(i+next, mask | (1 << j));
      }
      // don't add more to our current bottle
      // we swap with bottle j and "use" up our bottle
      int d2 = f(next, mask | (1 << j)) + price[i];
      
      // take the maximum of the decision branches
      res >?= max(d1,d2);
   }
   return dp[i][mask] = res;
}

int KiwiJuice::theProfit(int C, vector <int> bottles, vector <int> prices) {
   price = prices;
   int res = 0;
   _C = C;
   memset(dp,-1,sizeof(dp));
   data = bottles;
   return res = f(data[0],0);   // arbitrarily assign bottle 0 to be our special bottle
}

TC SRM 487 D1-250 CarrotJumping

Category: Brute Force, Mathematics
URL: http://www.topcoder.com/stat?c=problem_statement&pm=11022


class CarrotJumping {
public:
   int theJump(int);
};

long long MOD = 1000000007;

// compute 2^n (mod MOD)
long long pow(int v) {
   if (v == 0) return 1;
   if (v == 1) return 2;
   long long x = pow(v/2);
   if (v % 2 != 0) return (x * x * 2) % MOD;
   return (x * x) % MOD;
}

int CarrotJumping::theJump(int init) {
   int res = 100001;
   for (int a = 0; a <= 100000; a++) {
      for (int b = 0; b < 3; b++) {
         long long v = pow(a*3 + b*2);
         long long c = (v * init + (v - 1 + MOD) % MOD) % MOD;
         if (c == 0) {
            res <?= (a + b);
         }
      }
   }
   return res >= 100001 ? -1 : res;
}

Thursday, May 13, 2010

TC TCO10 Qualification Round 2

Overview:

This was the second qualification round for TCO10 with around 1300 competitors. As 600 competitors advance through to the next round (provided that at least 600 score a positive score), it provided ample opportunity to qualify for Round 1. The problems were rather simple with two greedy problems and a simple DP problem as the 1000 pointer. Having said that, the system test claimed a decent portion of the 500 pointer solutions due to a corner case which wasn't covered in the examples.

Problem 1 (250 Points): JingleRingle
Category: Greedy
URL: http://www.topcoder.com/stat?c=problem_statement&pm=10896

This problem required us to calculate the maximum profit we can make by buying jingles and selling them to various buyers. We are assumed to have infinite wealth so we can essentially buy all the jingles but we also want to maximise our profits.

The first key observation is to note that we can only use each buyer or seller once. The second key observation is that we only buy jingles from sellers if we can find a suitable buyer in which we can make a net profit (after tax). Now the question that remains is which sellers do we buy from and whom do we sell it to. Since each seller only gives us 1 jingle and we can only sell back 1 jingle at a time, this constraint essentially means that we can match a seller to a buyer (i.e. we are buying at a lower price and re-selling it at a higher price for a net profit). A simple greedy strategy would be then to buy at the lowest price and selling it to the highest price for maximum profit.

Does this greedy strategy work? Let's informally prove it. As we need to match a seller and buyer, if we buy from N people then we also have to sell it to N people. To maximise our profit we need to maximise our revenue, so the optimal choice is to choose the N highest buyers. Then to minimise our expenses we need to choose the N lowest sellers of jingles. It turns out it doesn't matter what order we pair the buyers and sellers as the total profit is the same if we choose the N highest buyers and N lowest sellers. Now all we need is to determine the value of N. Our greedy algorithm keeps pairing the highest available buyer to the lowest available seller until we reach a situation where we no longer profit, this is the value of N. If we assume this is not the value of N, then we have to somehow match a lower buyer to a higher seller or vice versa. But our value of N is at least as good as this on all cases hence our deduced value of N is correct.

A straightforward implementation:

class JingleRingle {
public:
   int profit(vector <int>, vector <int>, int);
};

int calcTax(int p, int tax) {
   return (p * tax) / 100;
}

int JingleRingle::profit(vector <int> buyOffers, 
      vector <int> sellOffers, int tax) {
   int res = 0;
   int buyJ = buyOffers.size()-1;
   sort(buyOffers.begin(),buyOffers.end());
   sort(sellOffers.begin(),sellOffers.end());
   for (int i = 0; i < sellOffers.size(); i++) {
      if (buyJ >= 0 && buyOffers[buyJ] - calcTax(buyOffers[buyJ],tax) 
         - sellOffers[i] >= 0) {
         res += buyOffers[buyJ] - calcTax(buyOffers[buyJ],tax) - sellOffers[i];
      }
      buyJ--;
   }
   return res;
}

Problem 2 (500 Points): FuzzyLife
Category: Greedy, Brute Force
URL: http://www.topcoder.com/stat?c=problem_statement&pm=10911

This problem has us calculate the maximum number of live cells on an infinite 2D plane after one iteration of time. The 2D plane is consisted of cells which are considered either dead or alive at a given slice in time. However, there are unknown cells denoted as '?' and our job is to determine whether each of these cells should be considered alive or dead for the next evolution in time to have the maximum number of alive cells.

The rules are:
- Each live cell with less than 2 or more than 3 live neighbours are dead
- Each dead cell with exactly 3 live neighbours becomes alive
- All other cell remain the same

Each cell has exactly 8 neighbours (even the ones on the "border" of the input because of the infinite plane) - 2 horizontal, 2 vertical and 4 diagonal.

The main observation required for this problem is the constraints of the problem allow us to adopt a greedy strategy and consider the neighbouring cells of a '?' as its own problem subinstance. As no cell will have more than 1 neighbour of type '?' this means that every cell is at most affected by 1 unknown cell. All we need to do is try each '?' cell by placing a '1' and a '0' in its place and simulate how many alive cells we are left with on the next iteration. We then simply choose the maximum of these two choices.

One tricky case is that we also need to take into account the '?' cell, if we place a '1' in it's position does it become alive or dead in the next iteration - vice versa if we place a '0'. We need to add this to our alive count before we compare which one is better. The other tricky case was to add a padding of dead cells around the grid as some of these can turn into alive cells (which need to be included in the answer) - this case however, was covered in the example cases. After we have filled in each of the '?' cells we can then do a single iteration pass and work out the number of alive cells.

The implementation using the ideas above:

class FuzzyLife {
public:
   int survivingCells(vector <string>);
};

// movement vectors
int dx[] = {-1, -1, -1, 0, 0, 1, 1, 1};
int dy[] = {-1, 0, 1, -1, 1, -1, 0, 1};

int FuzzyLife::survivingCells(vector <string> grid) {
   int res = 0;
   // pad the grid with extra 0's around the perimeter
   string str(grid[0].size()+2,'0');
   for (int i = 0; i < grid.size(); i++) {
      grid[i] = '0' + grid[i] + '0';
   }
   grid.insert(grid.begin(),str);
   grid.push_back(str);
   for (int i = 0; i < grid.size(); i++) {
      for (int j = 0; j < grid[i].size(); j++) {
         if (grid[i][j] == '?') {
            // compute options ('?' = 1)
            grid[i][j] = '1';
            int c1 = 0;
            for (int m = 0; m < 8; m++) {
               int neighX = i + dx[m];
               int neighY = j + dy[m];
               if (neighX < 0 || neighY < 0 || neighX >= grid.size() 
               || neighY >= grid[0].size()) continue;
               int alive = 0;
               for (int n = 0; n < 8; n++) {
                  int mX = neighX + dx[n];
                  int mY = neighY + dy[n];
                  if (mX < 0 || mY < 0 || mX >= grid.size() 
                  || mY >= grid[0].size()) continue;
                  if (grid[mX][mY] == '1') alive++;
               }
               if (grid[neighX][neighY] == '1' && (alive == 2 || alive == 3)) c1++;
               else if (grid[neighX][neighY] == '0' && (alive == 3)) c1++;
            }
            // consider the square itself under '?' = 1
            int selfAlive = 0;
            for (int n = 0; n < 8; n++) {
               int selfX = i + dx[n];
               int selfY = j + dy[n];
               if (selfX < 0 || selfY < 0 || selfX >= grid.size() 
               || selfY >= grid[0].size()) continue;
               if (grid[selfX][selfY] == '1') selfAlive++; 
            }
            if (selfAlive == 2 || selfAlive == 3) c1++;
            // compute options ('?' = 0)
            grid[i][j] = '0';
            int c2 = 0;
            for (int m = 0; m < 8; m++) {
               int neighX = i + dx[m];
               int neighY = j + dy[m];
               if (neighX < 0 || neighY < 0 || neighX >= grid.size() 
               || neighY >= grid[0].size()) continue;
               int alive = 0;
               for (int n = 0; n < 8; n++) {
                  int mX = neighX + dx[n];
                  int mY = neighY + dy[n];
                  if (mX < 0 || mY < 0 || mX >= grid.size() 
                  || mY >= grid[0].size()) continue;
                  if (grid[mX][mY] == '1') alive++;
               }
               if (grid[neighX][neighY] == '1' && (alive == 2 || alive == 3)) c2++;
               else if (grid[neighX][neighY] == '0' && (alive == 3)) c2++;
            }
            // consider the square itself under '?' = 0
            selfAlive = 0;
            for (int n = 0; n < 8; n++) {
               int selfX = i + dx[n];
               int selfY = j + dy[n];
               if (selfX < 0 || selfY < 0 || selfX >= grid.size() 
               || selfY >= grid[0].size()) continue;
               if (grid[selfX][selfY] == '1') selfAlive++; 
            }
            if (selfAlive == 3) c2++;
            // choose the one with better live cells
            if (c1 >= c2) grid[i][j] = '1';
         }
      }
   }   
   // iterate one period and determine number of alive
   for (int i = 0; i < grid.size(); i++) {
      for (int j = 0; j < grid[i].size(); j++) {
         int alive = 0;
         for (int n = 0; n < 8; n++) {
            int mX = i + dx[n];
            int mY = j + dy[n];
            if (mX < 0 || mY < 0 || mX >= grid.size() 
            || mY >= grid[0].size()) continue;
            if (grid[mX][mY] == '1') alive++;
         }
         if (grid[i][j] == '1' && (alive == 2 || alive == 3)) res++;
         else if (grid[i][j] == '0' && alive == 3) res++;
      }
   }
   return res;
}

Problem 3 (1000 Points): HandlesSpelling
Category: Dynamic Programming, String Manipulation
URL: http://www.topcoder.com/stat?c=problem_statement&pm=10864

This problem required us to play a game and work out the maximum score based on a set of choices we are forced to make. The game is that we have a string of letters in which we can "stamp" a set of pre-determined strings onto it. We need to maximise the score based on A^2 - B where A is the longest contiguous sequence of letters covered by the stampings and B is the number of uncovered letters.

For example, if the string we are given is HELLO and we are given 3 badges namely: E, HE, L. We can stamp them as follows: {HE}{L}{L}O. We stamped HE with the second badge and the 2 L's with the third badge. The value of A is 4 because we have a contiguous sequence of stampings (HELL) and B is 1 because we didn't have a way to stamp the letter O. Hence the score given this set of operations is 4^2 - 1 = 15, which turns out to be the maximum score we can get.

One strategy is to determine the maximum value of A. Why? A is where our score gets maximised, no other measure increases our score except for the value of A so it's only natural that we want to get A as high as possible. We also note that the maximum contiguous sequence (i.e. maximum value A) is in the optimal solution. Why? Consider two configurations:

xxxxx....----n----....xxxxxxxx
xxxxx--m--.......--m--xxxxxxxx

Here the x's denote similar configurations, so we only need to focus on maximising the centre part of m's, n's and dots. We let n be the maximum contiguous length possible in our string, let that be denoted as length x. Let's assign m to be length x-1. Then our goal is to prove that it is always optimal to choose n over the 2 m's.

Then the score for our first configuration is:
F = (x * x) - (x-1) - (x-1) = x^2 - 2x + 2 (length x ^ 2 minus two sections of x-1 uncovered)

The score for our second configuration is:
G = (x-1) * (x-1) - x = x^2 - 3x + 1 (length (x-1) ^ 2 minus one section of x uncovered)

Hence for all positive values of length x, the score for the first is better than the second. In fact, this is a stricter argument than what is required, as if the sections of n and m were disjoint they'll both be covered, hence there must be some conflict between n and m which forces us to choose between one or the other. Our argument will hold for such a scenario.

Our first task is then to find the maximum contiguous sequence of our string - this can be done easily through basic dynamic programming. Then we are tasked with trying all configurations which have the same maximum contiguous sequence to determine the minimum number of uncovered letters. Again, this can be done via a similar dynamic programming algorithm where we divide our main string into a left and right sub-string and count the maximum number of covered cells we can get. Then we simply subtract this by the size of the left and right sub-string respectively to get the number of uncovered letters.

The implementation which uses the strategy above can be seen as:

class HandlesSpelling {
public:
   int spellIt(vector <string>, vector <string>);
};

#define UNI -12345678
#define INF 12345678

vector<string> _badges;
string _S;

// returns true if str is a prefix to S starting at idx
bool isPrefix(int idx, const string& str, string S) {
   if (idx + str.size() > S.size()) return false;
   for (int i = idx; i < idx + str.size(); i++) {
      if (S[i] != str[i-idx]) return false;
   }
   return true;
}

int prefix[1001][51];   // true if badge j is a prefix at index i
int A[1001];   // keeps track of the longest sequence ending with index i
int B[1001];   // left hand side
int C[1001];   // right hand side
int maxSize;

// maximise coverage using ptr as the memoisation table
int func(int idx, int* ptr) {
   int res = 0;
   if (idx >= maxSize) return 0;
   if (ptr[idx] != -1) return ptr[idx];
   for (int i = 0; i < _badges.size(); i++) {
      if (prefix[idx][i] && idx + _badges[i].size() <= maxSize) {
         res >?= func(idx + _badges[i].size(), ptr) + _badges[i].size();
      }
   }
   res >?= func(idx+1, ptr);
   return ptr[idx] = res;
}

int HandlesSpelling::spellIt(vector <string> parts, vector <string> badges) {
   string data;
   for (int i = 0; i < parts.size(); i++) data += parts[i];
   _badges = badges;
   _S = data;
   // pre-compute the prefixes for speed
   for (int i = 0; i < data.size(); i++) {
      for (int j = 0; j < badges.size(); j++) {
         if (isPrefix(i, badges[j], _S)) { prefix[i][j] = 1; }
      }
   }
   // determine the longest contiguous sequence in our string
   for (int i = _S.size()-1; i >= 0; i--) {
      for (int j = 0; j < _badges.size(); j++) {
         if (isPrefix(i, _badges[j], _S)) {
            A[i] >?= A[i+_badges[j].size()] + _badges[j].size();
         }
      }
   }
   int longest = 0;
   for (int i = 0; i < _S.size(); i++) {
      longest >?= A[i];
   }
   // choose the best score out of all combinations which use A[i] == longest
   int res = UNI;
   for (int i = 0; i < _S.size(); i++) {
      if (A[i] == longest) {
         // branch off to left and right
         string left = _S.substr(0, i);
         string right = _S.substr(i+A[i]);
         // reset memo table for each instance
         memset(B,-1,sizeof(B));
         memset(C,-1,sizeof(C));
         int lOffset = 0, rOffset = 0;
         maxSize = i;
         if (left.size() > 0) lOffset = left.size() - func(0,B);
         maxSize = _S.size();
         if (right.size() > 0) rOffset = right.size() - func(i+A[i],C);
         res >?= (A[i] * A[i]) - (lOffset) - (rOffset);
      }
   }
   return res;
}

Monday, January 25, 2010

TC SRM 459 D2-1000 (ParkAmusement)

Category: Dynamic Programming, Graph Theory
URL: http://www.topcoder.com/stat?c=problem_statement&pm=10723

The problem describes having N landings in the amusement park. The landings can be terminating (i.e. the ride ends at this landing) which are either proper exits or crocodile ponds. The landings can also be non-terminating which means that you will slide down further or until you are stalled with no valid adjacent landings. All of the N landings have different heights and we can only slide down from a taller landing to a shorter one (by the laws of gravity). We start at a random landing and we wish to calculate the probability that we started at a specific landing given that we only went through exactly K different pipes to get safely home (i.e. at a proper exit and not at either a crocodile pond or a dead-end landing).

This problem is best split into two parts. The first is to calculate the probability we reach home from a specific starting point. The second is to calculate the probability that out of all the viable starting points we chose the specific starting landing. We are already given the adjacency matrix in terms of a vector/array of strings which we proceed to convert into a more usable format. We first keep track of exits - i.e. whether or not they are proper exits (home) or crocodile ponds. This is done by specifically checking whether or not landings[i][i] is '0' or not (based on the problem definition). If it's not '0', then the current vertex is an exit and we proceed to determine which type. Otherwise, we add the edges into our adjacency list. This is shown below:

vector<vector<int> > adjlist;   // adjacency list
int table[51];   // the exit table for each vertex, 0 if it's not an exit, 1 if it's home exit, 2 if it's a pond exit

...
double ParkAmusement::getProbability(vector <string> landings, int startLanding, 
      int K) {
   adjlist = vector<vector<int> >(landings.size(), vector<int>());
   for (int i = 0; i < landings.size(); i++) {
      if (landings[i][i] != '0') {
         // exit
         if (landings[i][i] == 'P') table[i] = 2;   // pond
         else table[i] = 1; // exit
         continue;
      }
      for (int j = 0; j < landings[i].size(); j++) {
         if (landings[i][j] == '1') adjlist[i].push_back(j);
      }
   }
...


Once we have parsed the input into a better format, we need to find out the probability of reaching home from a fixed starting vertex/landing. We can accomplish this through a DFS-like recursive function by keeping track of how many steps we have left to reach home. This involves keeping track of which vertex/landing we are currently at and how many hops left we must use. The base case is when the number of steps left is zero, if we are at an exit and it's the proper one (i.e. an 'E' in the matrix) then we have reached home - otherwise we have failed in our mission.

We utilise basic probability calculations to compute our chances of succeeding. If we arrive at a landing that has M out-going edges, then the probability of choosing a particular edge is 1/M. Therefore, the probability we will survive is simply the sum of 1/M * f(edge i, K-1) for all out-going edges i. The base cases are obvious, probability of surviving is 1.0 if we reach an 'E' exit otherwise it's 0.0. The function below illustrates this concept:

double calcPr(int startNode, int K) {
   if (K == 0) {
      // check if it's the exit
      if (table[startNode] != 1) return 0.0;
      return 1.0;
   }
   double res = 0.0;
   // note: if there are no more valid out-going edges it defaults to a 
   //    0.0 probability
   int poss = adjlist[startNode].size();
   for (int i = 0; i < adjlist[startNode].size(); i++) {
      res += (1.0 / poss) * calcPr(adjlist[startNode][i], K-1);
   }
   return res;
}

However, we can optimise this by memoising it. Otherwise we may run into trouble with time limits for the larger cases. This is easily done by setting up a cache table and marking it with a sentinel value to determine whether the state has been computed or not. Since no probability can be negative, we can exploit this fact.

double dp[51][51];

double calcPr(int startNode, int K) {
   if (dp[startNode][K] >= 0.0) return dp[startNode][K];
   if (K == 0) {
      if (table[startNode] != 1) return dp[startNode][K] = 0;
      return dp[startNode][K] = 1.0;
   }
   double res = 0.0;
   int poss = adjlist[startNode].size();
   for (int i = 0; i < adjlist[startNode].size(); i++) {
      res += (1.0 / poss) * calcPr(adjlist[startNode][i], K-1);
   }
   return dp[startNode][K] = res;
}

Now, we have finished the first step of the problem. The second step involves calculating out of the many possible landings we have chosen and given that we did arrive home safely, what is the probability we chose startLanding? The easiest way to think of this is to think geometrically. If we let a circle C be the total sum of the probabilities that we arrived safely from all possible landings, we can scale this to 1.0 as we know for certainty that we arrived safely (given in the problem). Now we simply need to calculate how much probability was contributed by starting from startLanding. This is simply calcPr(startLanding, K) / totalSum. This can be visualised as:



double tot = 0.0;
for (int i = 0; i < 51; i++) for (int j = 0; j < 51; j++) dp[i][j] = -1.0;
for (int i = 0; i < landings.size(); i++) {
   tot += calcPr(i, K);
}
return calcPr(startLanding, K) / tot;

Combining all of the above steps, we yield the final code:

class ParkAmusement {
public:
   double getProbability(vector <string>, int, int);
};

vector<vector<int> > adjlist;
int table[51];
double dp[51][51];

double calcPr(int startNode, int K) {
   if (dp[startNode][K] >= 0.0) return dp[startNode][K];
   if (K == 0) {
      if (table[startNode] != 1) return dp[startNode][K] = 0;
      return dp[startNode][K] = 1.0;
   }
   double res = 0.0;
   int poss = adjlist[startNode].size();
   for (int i = 0; i < adjlist[startNode].size(); i++) {
      res += (1.0 / poss) * calcPr(adjlist[startNode][i], K-1);
   }
   return dp[startNode][K] = res;
}

double ParkAmusement::getProbability(vector <string> landings, int startLanding, 
      int K) {
   adjlist = vector<vector<int> >(landings.size(), vector<int>());
   for (int i = 0; i < landings.size(); i++) {
      if (landings[i][i] != '0') {
         // exit
         if (landings[i][i] == 'P') table[i] = 2;   // pond
         else table[i] = 1; // exit
         continue;
      }
      for (int j = 0; j < landings[i].size(); j++) {
         if (landings[i][j] == '1') adjlist[i].push_back(j);
      }
   }
   double tot = 0.0;
   for (int i = 0; i < 51; i++) for (int j = 0; j < 51; j++) dp[i][j] = -1.0;
   for (int i = 0; i < landings.size(); i++) {
      tot += calcPr(i, K);
   }
   return calcPr(startLanding, K) / tot;
}

Tuesday, December 8, 2009

TC SRM 454 D1-500/D2-1000 (NumbersAndMatches)

Category: Dynamic Programming
URL: http://www.topcoder.com/stat?c=problem_statement&pm=10709

We are given a number which is composed of a sequence of digits. Each digit is represented by a combination of matches represented in a certain position to represent the digit. We are given the option of using up to K different match moves and are asked to calculate the number of different integers we can construct without discarding or adding any new matches to our given number.

First, we decompose the problem into several parts. The first involves somehow representing each digit in our program. A simple way to do this is to use a 2D integer array in which index (i,j) is equal to 1 if the digit i has a match in position j. An alternative way is to use a bitmask with the bit representation of the j-th bit set to 1 if there is a match in position j.

Next, we need to somehow compute the transitional costs between two different digits. For example, we need to quantify how many moves we would use up if we convert say a 6 to a 3. Since we are only allowed a limited number of moves its obvious we need a way to keep track of how many we have used up. One possible way is to determine how many places the two numbers differ. Then the number of moves we need to make is equivalent to exactly half of this.

There is another issue we need to deal with, that is the situation where the number of matches used to construct two digits differ. For example, the digit 6 uses 6 matches whereas digit 3 only uses 5 matches. We need to keep track of how much surplus (or deficit) matches we get from transitioning from one digit to another. This is important as we cannot add or remove matches from the final number - i.e. the number of matches used in constructing a unique integer must be the same as the original integer. Now we need to keep track of two things when we consider transforming one digit to another: the number of moves it takes and the difference in matches we have used so far.

We can define the cost function as follows:



To calculate the surplus/deficit it's quite simple, just calculate the number of matches used for each digit and subtract one from the other depending on which direction the transformation is taken. The above cost function ensures that we don't double count the number of moves required from moving a match from a surplus digit to a deficit digit. A way to visualise this is we "precharge" the surplus with an additional move but we are also granted a free move since this surplus is already counted once so the actual number of moves decreases if we need to import more matches from another transformation.

Lastly, we need to somehow compute all the ways we can construct unique integers given that we keep the number of matches used the same (i.e. at the end of our decision algorithm the number of surplus matches is exactly 0) as well as using less than or equal to K moves (which we subtract along the way). This suggests a recursive algorithm from potentially changing the first digit to something else, depending on the transformation we are left with either more/less/equal number of matches and the number of available moves to be less or equal to what we started with. We must also handle the case where we "borrow" matches from the latter part of the integer, so we must also handle cases where our match balance is negative. This can be done by computing an offset to keep the total number positive (to use for our caching purposes).

To efficiently compute the number of ways we need to memoise or apply bottom up DP to our algorithm. Note that the number of unique states is rather small so computational time isn't an issue. An implementation of the above ideas is shown below (using memoisation):

class NumbersAndMatches {
public:
   long long differentNumbers(long long, int);
};

int mat[10][7] = {
  {1,1,1,0,1,1,1}, // 0
  {0,0,1,0,0,1,1}, // 1
  {1,0,1,1,1,0,1}, // 2
  {1,0,1,1,0,1,1}, // 3
  {0,1,1,1,0,1,0}, // 4
  {1,1,0,1,0,1,1}, // 5
  {1,1,0,1,1,1,1}, // 6
  {1,0,1,0,0,1,0}, // 7
  {1,1,1,1,1,1,1}, // 8
  {1,1,1,1,0,1,1}  // 9
};

// defines the total number of matches used for digit i
int tot[10] = {6, 3, 5, 5, 4, 5, 6, 3, 7, 6};   

long long dp[19][155][155];
pair<int,int> costTab[10][10];
string num;

#define SURPLUS_MOD 72

long long func(int idx, int moves, int surplus) {
   // offset due to the possibility of negative surplus
   int surplusMod = surplus + SURPLUS_MOD;
   if (idx >= num.size()) {
      // requires the surplus amount of matches to be 0 otherwise its not valid
      if (moves >= 0 && surplus == 0) return 1;
      return 0;   
   }
   if (dp[idx][moves][surplusMod] != -1) return dp[idx][moves][surplusMod];
   long long res = 0;
   // try each digit
   int curDig = num[idx] - '0';
   for (int dig = 0; dig <= 9; dig++) {
      int c = costTab[curDig][dig].first;
      int d = costTab[curDig][dig].second;
      if (moves - c >= 0) {
         res += func(idx+1, moves - c, surplus + d);
      }
   }
   return dp[idx][moves][surplusMod] = res;
}

long long NumbersAndMatches::differentNumbers(long long N, int K) {
   long long res = 0;
   memset(dp,-1,sizeof(dp));
   // compute the cost function for every digit -> digit transition
   for (int i = 0; i <= 9; i++) {
      for (int j = 0; j <= 9; j++) {
         int c = 0;
         for (int k = 0; k < 7; k++) {
            if (mat[i][k] != mat[j][k]) c++;
         }
         int s = tot[i] - tot[j];
         costTab[i][j] = make_pair((c+s)/2, s);
      }
   }
   // decompose the number into a string
   long long vN = N;
   while (vN > 0) {
      num += ((vN % 10) + '0');
      vN /= 10;
   }
   reverse(num.begin(),num.end());
   // calculate the number of ways
   return res = func(0, K, 0);
}

Thursday, November 26, 2009

TC SRM453.5 D1 450 (PlanarGraphShop)

Category: Dynamic Programming, Graph Theory
URL: http://www.topcoder.com/stat?c=problem_statement&pm=10412

This problem asks us to determine the minimum number of simple planar graphs we need to exact use of N gold coins. Each planar graph has an associated cost function based on the number of vertices and edges the graph has, more specifically: cost = (number of vertices)^3 + (number of edges)^2. We should immediately see a knapsack-like algorithm which strongly suggests using dynamic programming. The main problem lies with determining what number of vertices and edges constitutes a valid simple planar graph.

Luckily there is a simple theorem which we can use as an upper bound calculation to determine whether a pairing of (#V,#E) is a valid pair. For graphs that has at least 3 vertices: if e > 3v - 6 then this will imply nonplanarity. Hence, we can deduce that all e <= 3v - 6 has a simple planar graph representation. For graphs that has less than 3 vertices we can simply hand calculate the combinations. Once we know what our valid pairings are we can construct a set of feasible cost moves based on our cost function by iterating through all possible pairings of # vertices and # edges. We then apply a simple DP algorithm to calculate the minimum number of planar graphs needed.

An implementation is shown below:

class PlanarGraphShop {
public:
   int bestCount(int);
};

int dp[50001];
vector<int> moves;

#define INF 1 << 20

int func(int N) {
   if (N == 0) return 0;
   if (N < 0) return INF;
   if (dp[N] != -1) return dp[N];
   int res = INF;
   for (int i = 0; i < moves.size() && moves[i] <= N; i++) 
      res <?= 1 + func(N - moves[i]);
   return dp[N] = res;
}

int PlanarGraphShop::bestCount(int N) {
   set<int> costs;
   costs.insert(1); costs.insert(8); costs.insert(9);
   for (int v = 1; v < 37; v++) 
      for (int e = 0; e <= 3 * v - 6; e++)
         costs.insert(v*v*v + e*e);
   moves = vector<int>(costs.begin(),costs.end());
   memset(dp,-1,sizeof(dp));
   return func(N);
}

Tuesday, September 22, 2009

TC SRM424 D1 600 (StrengthOrIntellect)

Category: Dynamic Programming
URL: http://www.topcoder.com/stat?c=problem_statement&pm=10174

This problem asks us to determine the maximum number of missions we can complete given two character attributes: strength and intellect. We meet the pre-requirements to complete a mission if either one of the strength or intellect attributes are at least as high/equal to the mission's strength and intellect requirements respectively.

A greedy approach does not work here. You can informally argue by noting that when we are given two choices to increase strength or intellect, there is no clear decision on which option to choose. The local optimality does not hold for the global optimality, as if we choose one or the other we might miss out on earning a lot of points which will allow us to complete more missions. Another greedy heuristic might be to try both maximising strength and maximising intellect and choose the maximum number of missions from the two possible options. This also does not work because the optimal solution might not be on a extreme point of strength or intellect. Since the decisions we make need to be flexible enough to change between increasing strength and/or then increasing intellect or vice-versa, it's rather clear that no amount of heuristics will pass all the (corner) test cases.

The next obvious approach to consider is brute force. We can keep track of which missions we have accomplished and how much strength and intellect we have achieved so far, then for all possible mission transitions recurse along the mission and adjust the strength and intellect. In effect, this almost works if we memoise the solution. The problem is keeping track of the missions we have completed - the naive approach is using a bitmask which will definitely exceed the memory constraint (2^50 at worst case).

The trickest part of the problem is state compressing our DP algorithm. If we can somehow reduce the memory footprint whilst also keeping the correctness intact, then a viable solution is generated. Let's take out our bitmask and try to find a way to implicitly represent that information. A common reduction in this situation is to somehow find the relationship between our current DP state parameters (in this case, strength and intellect) and the DP state parameter(s) we want to get rid of (in this case, how we can implicitly represent the missions completed so far). In fact we can see that they have a direct relationship! If we are given a particular strength and intellect we can already uniquely determine which missions can be completed, as we take all viable missions as all the points are strictly positive (always something to gain, nothing to lose).

Now our DP state becomes simply:

D(S, I) = the maximum number of missions if we can get to strength S and intellect I.
D(S, I) = 0 if the state is unreachable (not enough points).

We can proceed to pre-compute the maximum number of missions for each state as there are only 1000 * 1000 = 1 million such states. We can also pre-compute whether or not a state is possible or not. A state is impossible if we don't have enough points to cover the increase from strength 1 to S, and intellect 1 to I. Of course, this doesn't fully dictate the feasibility of a state but it's a good start. To compute the number of excess (unused) points we simply sum up all the points gained from completing the missions and subtracting S-1 and I-1 from it (as we expended S-1 strength points and I-1 intellect points to get to the current state). This calculation is done as follows:

pair<int,int> preComp[1011][1011];
...
   for (int i = 0; i <= 1000; i++) {
      for (int j = 0; j <= 1000; j++) {
         int pointsGained = 0, missions = 0;
         for (int k = 0; k < strength.size(); k++) {
            if (strength[k] <= i || intellect[k] <= j) {
               pointsGained += points[k];
               missions++;
            }
         }
         int extraPoints = pointsGained - i - j + 2;
         // extraPoints < 0 --> inaccessible state set to 0 missions
         // the pair is (excess points, number of missions completed)
         preComp[i][j] = make_pair(extraPoints, extraPoints < 0 ? 0 : missions);
      }
   }
...

Next, we need to determine feasibility from our original DP state (1,1) - which is the strength and intellect we start with. We recursively call the next DP state if we have unused points by either deciding to upgrade our strength by 1 or our intellect by 1. Note the number of unique states is only 1 million so this should run in time. Combining this with a caching mechanism (i.e. memoisation) we come up with our final solution:

class StrengthOrIntellect {
public:
   int numberOfMissions(vector <int>, vector <int>, vector <int>);
};

pair<int,int> preComp[1011][1011];
int memo[1011][1011];

int func(int strength, int intellect) {
   // seen it before (caching mechanism)
   if (memo[strength][intellect] != -1) return memo[strength][intellect];
   memo[strength][intellect] = preComp[strength][intellect].second;
   int extra = preComp[strength][intellect].first;
   // if there are unused points, recurse on them otherwise the maximum
   // is simply just the number of missions we have obtained so far
   return memo[strength][intellect] >?= extra <= 0 ? memo[strength]
      [intellect] : max(func(strength+1, intellect), func(strength, intellect+1));
}

int StrengthOrIntellect::numberOfMissions(vector <int> strength, vector <int> 
      intellect, vector <int> points) {
   int res = 0;
   for (int i = 0; i <= 1000; i++) {
      for (int j = 0; j <= 1000; j++) {
         int pointsGained = 0, missions = 0;
         for (int k = 0; k < strength.size(); k++) {
            if (strength[k] <= i || intellect[k] <= j) {
               pointsGained += points[k];
               missions++;
            }
         }
         int extraPoints = pointsGained - i - j + 2;
         // extraPoints < 0 --> inaccessible state set to 0 missions
         preComp[i][j] = make_pair(extraPoints, extraPoints < 0 ? 0 : missions);
      }
   }
   memset(memo,-1,sizeof(memo));
   return res = func(1, 1);
}

Friday, September 11, 2009

TC TCO08 Round 2 D1 500 (CreaturesTraining)

Category: Dynamic Programming
URL: http://www.topcoder.com/stat?c=problem_statement&pm=8570

This problem asks us to upgrade a set of creatures to maximise the final score after D days. The score is defined as the sum of the power level of the creatures. This is a tricky problem to solve efficiently. The key observation to make is that you never downgrade creatures and you can also skip creatures. Therefore this suggests that we can go from the leftside to rightside in a linear fashion. We can do a rough sketch of the recursive state as:

D(idx, x) = the maximum number of points we can get when we are processing level idx of the creatures and we have x days left in credit to upgrade our creatures.

When we try to formulate our recurrence relation we run into a problem, the state does not fully describe the problem state as it does not keep track of which creatures we have upgraded so far. So we begin to re-formulate our DP state by incorporating the set of monsters we have upgraded so far, i.e. the modifications we have made to the count array. One obvious yet inefficient way is to keep track of the count array as part of the DP state. However, this is too large of a memory footprint and as a result causes our solution to become inefficient. However, it is a good starting point as we can re-use our observation before that we never look back at a previous level. With this in mind, we can simply just keep track of the current extra creatures we have upgraded so far and optimally decide how many of these to keep in each level. Therefore, our DP state becomes:

D(idx, x, carry) = the maximum number of points we can get when we are processing level idx of the creatures and we have x days left in credit to upgrade our creatures. Given that we are also carrying forth carry number of monsters from the previous level.

In each state we can choose between 0 to count[idx] + carry monsters to upgrade. Each of these are valid transitions/upgrades provided we have enough day credit to upgrade them. We need to check whether or not this state space fits within the memory constraints, as there are a maximum of 50 monster levels and 100 days and carries (as we cannot carry more than 100 as this will exceed the upper limit of day credits). So we are faced with a memory footprint of 50 * 100 * 100 = 500,000 64-bit integers which easily fits within the limit.

Lastly, we consider the base/terminating case, this is when we reach the last level of the creatures. We implicitly prune any transitions where we exceed the amount of remaining day credit, so we do not need to consider this in our terminating case.

Implementing the ideas above is rather simple and short:

class CreatureTraining {
public:
   long long maximumPower(vector <int>, vector <int>, int);
};

long long dp[51][111][111];
vector<int> _count;
vector<int> _power;

long long func(int idx, int day, int carry) {
   if (idx >= _count.size()) return dp[idx][day][carry] = 0;
   if (dp[idx][day][carry] != -1) return dp[idx][day][carry];
   long long res = 0;
   int mx = carry + _count[idx];
   for (int i = 0; i <= mx; i++) {
      // move a set of 0 -> mx monsters up to the next level
      // leave mx - i behind
      if (day - i >= 0) {
         res >?= func(idx+1, day - i, i) + (mx-i) * (long long)_power[idx];
      } else {
         break;
      }
   }
   return dp[idx][day][carry] = res;
}

long long CreatureTraining::maximumPower(vector <int> count, vector <int> power, 
      int D) {
   long long res = 0;
   memset(dp,-1,sizeof(dp));
   limitD = D; _count = count; _power = power;
   return res = func(0, D, 0);
}

TC TCO08 Qualification Round 1 D1 600 (PrimeSums)

Category: Dynamic Programming
URL: http://www.topcoder.com/stat?c=problem_statement&pm=8596

This problem continues with our DP-themed run of recent problems. It asks us to compute the number of subsets (including the empty subset) in which the total bag weight is a prime number. The constraints give us an upper bound of 500,000 (50 * 10,000) for the maximum prime. So the first task is to compute all the prime numbers under this number, this is an easy and efficient task if you know the Sieve of Eratosthenes. The next step is determining how many ways we can make a sub-bag of a weight W, we can use this in our last step by adding up all the weights W which are prime to yield our final answer. Before that though, we need to determine the number of sub-bags which totals a weight of W.

We notice that there are overlapping subproblems which screams to us: dynamic programming! To see this, imagine we have a bag of {1, 1, 2, 7}. To compute the number which weights are possible for all subsets: we start off with the first element. We can choose to use this element and place it into our bag, or we can reject this element and not put it into our bag. Either way, we have {1, 2, 7} as a sub-problem however if we don't do any caching or DP we will compute this sub-problem twice. This effect carries on from the first element to the n-th element and the number of sub-problems we are forced to compute grows exponentially. So we can simply memoise our solution, or can we?

Looking at the constraints we need to keep track of our current index inside the bag (so we know which element we are considering) and also the current weight of the bag (which has an upper bound of 500,000). Since the number of elements is 50 are the worst case, using a memoisation table of 50 * 500,000 64-bit integers is too large to fit into the memory limit. There must be another way to optimise our space complexity. As with all DP problems, we only need to keep sub-solutions to as much as our recurrence relation needs them. So let's define our recurrence relation in hopes of a good optimisation:

DP State: F(n,w) = the number of sub-bags we can construct using a set of 1..n bags (1-based index) that sums to exactly a weight of w.
DP Recurrence: F(n,w) = F(n-1,w - b[n]) + F(n-1,w)

The recurrence comes from the fact that we are faced with 2 decisions at a particular state, to take the item and place it into the bag (the first recursive call) and to skip the item and retain our weight (the second recursive call). As you can see from the depth of the recurrence relation is only at most 1 (i.e. we only need to keep track of the n-1'th pre-computed values), then we can reduce our memoisation table to 2 * 500,000 64-bit integers and use modulus arithmetic to implicitly swap the rows around. See the DP: State Compression article for more information on this technique.

As usual, we also need to define a base case. This is pretty simple, we are allowed to have an empty subset so naturally there is 1 way to have a weight of 0 without using any elements in the bag. Hence F(0,0) = 1. However, we encounter a major problem when we try and implement the solution. The fact there are non-unique elements in the bag can cause our dynamic programming algorithm to double-count the same sub-bag. For example, consider the bag {1, 1, 2, 7} if we skip the first bag and decide to use all the other bags we yield {1, 2, 7} with a weight of 10. However, we can also choose to use the first bag and skip the second bag to yield the same sub-bag {1, 2, 7} also with a weight of 10. Our current algorithm does not take this into account and hence will produce incorrect results for certain inputs.

A neat way to tackle this problem is to condense the same weights into one bag and keep track of how many elements we have of that weight. In the example used previously, we can make a mapping of { {1,2}, {2,1}, {7,1} } which has a pair in each element denoting the weight of the element and how many times it is inside our bag. Then we apply the same DP algorithm to this modified mapping and introduce additional decisions to include k = 1 to n of each element, where n is defined as the number of times it is in our bag. So we can choose to use 0, 1 or 2 elements of weight 1 but we are only allowed to choose 0 or 1 elements of either weight 2 and 7. Note that this does not cause us to double count as we only allow a transition of 2 elements of weight 1 once in our algorithm.

Lastly, we just need to add up our DP table for all the prime weights which is pretty self-explanatory.

class PrimeSums {
public:
   long long getCount(vector <int>);
};

int prime[500011];
long long dp[2][500011];

long long PrimeSums::getCount(vector <int> bag) {
   long long res = 0;
   sort(bag.begin(),bag.end());
   // sieve to generate the primes
   memset(prime,1,sizeof(prime));
   prime[0] = prime[1] = 0;
   for (int i = 2; i * i <= 500000; i++) {
      if (prime[i]) {
         for (int j = i * i; j <= 500000; j += i) {
            prime[j] = 0;
         }
      }
   }
   // merge identical bag elements together
   map<int,int> mappings;
   for (int i = 0; i < bag.size(); i++) mappings[bag[i]]++;
   
   int kth = 1;
   
   // base case
   dp[0][0] = 1;
   
   for (map<int,int>::iterator it = mappings.begin(); it != mappings.end(); it++, 
   kth++) {
      for (int i = 0; i <= 500000; i++) {
         dp[kth&1][i] = dp[!(kth&1)][i];
         // loop through possible weight arrangements based on the available set
         for (int j = 1; j <= it->second; j++) {
            if (i-j*it->first < 0) break;
            dp[kth&1][i] += dp[!(kth&1)][i-j*it->first];
         }
      }
   }
   // sum up all the number of sub-bags which have a prime weight
   for (int j = 0; j <= 500000; j++) {
      if (prime[j]) res += dp[!(kth&1)][j];
   }
   return res;
}

Thursday, September 10, 2009

TC TCO06 Round 1 D1 350 (SeparateConnections)

Category: Dynamic Programming, Matching in a General Graph
URL: http://www.topcoder.com/stat?c=problem_statement&pm=6095

This problem just asks us to find the maximum number of matches in a general graph (should be a straightforward reduction given the problem constraints). The most efficient way to solve this problem is to use a maximum matching algorithm (like Edmond's). However, this is not straightforward to code and certainly not required for a 350 pointer problem. Given the constraints, one easy way to implement a maximum matching algorithm is to use dynamic programming. We will be using a bitmask to efficiently keep track of the nodes that have already been matched so far.

For those new to bitmasks, here is a very short tutorial on it. As you may know, all decimal integer numbers (base 10) have a binary representation form (base 2). As a result, each bit inside a binary representation can represent up to 1-bit of information (duh). Besides representing numbers we can use each of these bits (0 or 1) as useful book-keeping information - in our case (and in many cases) we use them to mark whether or not we have visited or used a particular resource. So for example, the bit representation 01100 could mean we have used resource 2 and 3 (0-based index). Note that we go from left to right due to the fact that the left-most bit is also the most significant bit (MSB). This bit representation (01100) is simply known as 2^3 + 2^2 = 14 in decimal. Expanding on this, let's say in our problem we have used node 1, node 4 and node 5. This is represented as 2^5 + 2^4 + 2^1 = 32 + 16 + 2 = 50. There are several important bit operations we can use to manipulate and query the bits:

For starters, to get a '1' in the n-th position from the right (remember 0-based index) - you use the shift left operator '<<'. For example, 1 << 5 yields a bit representation equivalent to 0010 0000 (32).

To set a bit, we use the bitwise-OR operator denoted as '|'. The bitwise-OR operator works by applying the OR operation to each of the corresponding bits between two numbers. The OR operation is defined as:

0 OR 0 = 0
0 OR 1 = 1
1 OR 0 = 1
1 OR 1 = 1

For example to set the 5th bit (again 0-index) from the right, we first use the shift left operator to yield the correct number (i.e. 1 << 5). Let's say the bitmask before the operation is defined as 0100 0010 (68). Then applying 1 << 5 which has a representation of 0010 0000 (32) yields:

0100 0010 (68) OR 0010 0000 (32) Equals 0110 0010 (100)

To test whether a bit is set within a bitmask, we make use of the bitwise-AND operator denoted as '&'. The AND operation is defined as:

0 AND 0 = 0
0 AND 1 = 0
1 AND 0 = 0
1 AND 1 = 1

We use a similar method as we did for setting a bit in a particular position, let's say we wanted to test if the 5th bit was set. We use the shift left operator to yield the correct number (i.e. 1 << 5). Then we simply apply the bitwise-AND operator to the bitmask we wanted to test. If our mask was 0110 0010 then it would return 0010 0000. However if our mask was 0100 0010, then it would return 0000 0000 as it has no bits both turned on at the same position. Thus, we can use this to test whether or not a bitmask contains a bit in a particular position, if it the result is positive after the AND operation then it is inside the bitmask, if it's zero then it isn't inside the bitmask.

There are other bitmasking techniques which won't be discussed here as they aren't used for this problem.

Back to the problem, given that there are only 18 nodes at maximum. We can easily create a bitmask that keeps track of which nodes we have currently used up (there are 2^18 = 262144 different combinations). So for the actual DP state, we can iteratively introduce more nodes to the set and choose any adjacent nodes which haven't been used so far. Therefore, the state is D(n,m) where n is the index of which node we are processing and m is the bitmask of which nodes we have already used up. Next, we need to define the recurrence relation as well as any base cases. The recurrence relation can be found by selecting from a set of valid decisions at a particular node.

The two available options are to:
- Not used the node in the maximum matching (in hopes that it increases the available nodes later on)
- Use the node along with another node which is adjacent to it (increases the number of matchings by 1 and decreases the available set)

We define the recurrence as: D(n,m) = max { 1 + D(n+1, m | (1 << k), D(n+1, m) } where k is the set of all unused nodes which are adjacent to n. Take note, we do not consider matching node n, if it's already been matched by the time we reach index n. Then applying a simple memoization we yield a working solution:

class SeparateConnections {
public:
   int howMany(vector <string>);
};

int adjmat[19][19];
int memo[19][1<<19];
int sz;

int func(int idx, int mask) {
   if (idx >= sz) return 0;
   if (memo[idx][mask] != -1) return memo[idx][mask];
   int res = 0;
   if (!(mask & (1 << idx))) {
      for (int i = idx+1; i < sz; i++) {
         if (mask & (1 << i)) continue;
         if (!adjmat[idx][i]) continue;
         res >?= 1 + func(idx+1, mask | (1 << i));
      }
   }
   res >?= func(idx+1, mask);
   return memo[idx][mask] = res;
}

int SeparateConnections::howMany(vector <string> mat) {
   int res = 0;
   for (int i = 0; i < mat.size(); i++) {
      for (int j = 0; j < mat.size(); j++) {
         if (mat[i][j] == 'Y') adjmat[i][j] = adjmat[j][i] = 1;
      }
   }
   sz = mat.size();
   memset(memo,-1,sizeof(memo));
   return res = func(0,0) * 2;
}

However, one can also use a general matching algorithm like Edmond's Blossoming/Shrinking Algorithm to solve the problem. See below for such an implementation (based on Pape and Conradt variation of Edmond's which is incorrect for a certain set of inputs, but it still passed system tests for this problem). No further explanation would be provided for this algorithm as I'll leave this for a later blog entry in which we build up our matching algorithms.

#define MAX_VERTICES 21

int adjmat[MAX_VERTICES+1][MAX_VERTICES+1];
int table[MAX_VERTICES+1][MAX_VERTICES+1];
int mate[MAX_VERTICES+1];
int level[MAX_VERTICES+1];
int grandparent[MAX_VERTICES+1];

bool checkAncestor(int currentNode, int checkNode) {
       int v = grandparent[currentNode];
       if (v < 0) return false;   // not an ancestor
       if (v == checkNode) return true; // ancestor
       return checkAncestor(v, checkNode);
}

int augment(int curNode) {
       // BFS an augmenting path
       memset(level,1,sizeof(level));
       memset(grandparent,-1,sizeof(grandparent));
       queue<int> q;
       q.push(curNode);
       level[curNode] = 0;
       while (!q.empty()) {
               int x = q.front(); q.pop();
               for (int y = 1; y <= MAX_VERTICES; y++) {
                  // not part of an edge in G or seen the odd-vertex in T
                  if (!adjmat[x][y] || !level[y]) continue;   
                  if (mate[y] == 0) {
                     // found an augmenting path
                     int next = 0;
                     while (true) {
                        mate[y] = x;
                        next = mate[x];
                        mate[x] = y;
                        if (grandparent[x] < 0) break;
                        x = grandparent[x];
                        mate[next] = x;
                        y = next;
                     }
                     return 1;
                 } else if (mate[y] != 0 && !checkAncestor(x,y) && 
                         mate[y] != x) {
                     // push the node into the alternating tree T
                     q.push(mate[y]);
                     grandparent[mate[y]] = x;
                     level[y] = 0;
                 }
               }
       }
       return 0;
}

int maxMatching() {
   int res = 0;
   memset(mate,0,sizeof(mate));    // all currently unmatched
   memset(grandparent,-1,sizeof(grandparent));
   // start at the vertex i if it's single
   while (true) {
      int v = 0;
      for (int i = 1; i < MAX_VERTICES; i++) {
         if (!mate[i]) v += augment(i);
      }
      if (v == 0) break;
      res += v;
   }
   return res;
}

int SeparateConnections::howMany(vector <string> mat) {
   for (int i = 0; i < mat.size(); i++) {
      for (int j = 0; j < mat.size(); j++) {
         if (mat[i][j] == 'Y') adjmat[i+1][j+1] = adjmat[j+1][i+1] = 1;
      }
   }
   return maxMatching() * 2;
}

TC TCO05 Round 2 D1 500 (PartyGame)

Category: Dynamic Programming, Probability
URL: http://www.topcoder.com/stat?c=problem_statement&pm=4635

The problem asks us to compute the expected number of moves required to reach from left-most side of a row to the right-most side. The row consists of a number of R's and G's which control what actions we can take. We are also given a N-sided dice which allows us to randomly (uniformly) go 1 to N places to the right. If we land on an R we have to go back to the position we rolled from (i.e. no change in position, but a penalty of a move is applied). If we land on a G cell then all is good and we are allowed to keep our new position. If it's impossible to move from the left-most cell to the right-most cell then return -1. Also note that you can "over-hit" the target, anything which takes us beyond the last rightmost cell is considered a win.

The first thing to come into mind is dynamic programming. Why? Given a particular cell the expected number of moves does not change and is only dependent on the cells which come to the right of it. Another hint, is that there are overlapping subproblems as a result - since the number of possible paths is exponential it only makes sense to not re-calculate the same cells over and over. To apply DP to the problem we have already noted that it is sufficient to represent the expected number of cells by itself, this forms a 1-dimensional DP problem. The order in which we calculate should be obvious - we need to know all cells to the right of a particular cell, so we need to calculate from rightmost to leftmost cell.

Now defining the actual recurrence relation is more tricky, given a set of transitions from a cell some of them re-direct them back to our cell - how do we handle this case? Let's start off by defining f(x) as the expected number of moves for cell x. Now since the probability transitions are equal, i.e. 1 / N. If the next transition we land on for a given jump is 'G' then the expected number of moves grows by (1 / N) * (f(y) + 1) where y is the next cell. If the next transition we land on is 'R' then we have to go back to our current position, so the expected number of moves actually grows by (1 / N) * (f(x) + 1), where x is our current cell. Using an addition of all these expected moves defines our f(x). Hence we are left with something of the form:

f(x) = (1 / N) * (f(x) + 1) + (1 / N) * (f(x) + 1) + ... + (1 / N) * (f(x) + 1) + (1 / N) * (f(y) + 1) + (1 / N) * (f(z) + 1) + ...

Where the first part are all the cell transitions which lead them back to the current cell. The second part is the "good" transitions, in the example above these corresponds to the cells y and z. Now we can combine the first part together:

f(x) = (a / N) * (f(x) + 1) + (1 / N) * (f(y) + 1) + (1 / N) * (f(z) + 1) + ...

Where there are "a" number of cell transitions which lead it back to the current cell (i.e. 'R' cells). We need to formulate it in the form of f(x) = ... where the right hand side does not contain any f(x) terms. This requires a basic mathematical manipulation:

(1 - (a / N)) * f(x) = (1 / N) * (f(y) + 1) + (1 / N) * (f(z) + 1) + ...
f(x) = ((1 / N) * (f(y) + 1) + (1 / N) * (f(z) + 1) + ...) / (1 - (a / N))

This let's us calculate the expected number of moves including transitions which lead us back to our current cell. Hence combining all the concepts above, we come up with the following DP implementation:

class PartyGame {
public:
   double numberOfMoves(vector <string>, int);
};

double dp[2511];

double PartyGame::numberOfMoves(vector <string> field, int N) {
   string totalField;
   for (int i = 0; i < field.size(); i++) {
      totalField += field[i];
   }
   double res = 0.0;
   dp[totalField.size()-1] = totalField[totalField.size()-1] == 'R' ? 0.0 : 1.0;
   for (int i = totalField.size()-2; i >= 0; i--) {
      if (totalField[i] == 'R') continue;
      dp[i] = 0.0;
      int totBad = 0;
      double moveGood = 0.0;
      for (int j = 1; j <= N; j++) {
         if (i + j >= totalField.size()-1) {
            // instant win
            moveGood += (1.0 / N);
            continue;
         }
         if (totalField[i+j] == 'R') {
            // bad move (return to self)
            totBad++;
         } else {
            moveGood += (1.0 / N) * (dp[i+j] + 1);
         }
      }
      // solve equation
      dp[i] = (((double) totBad / N) + moveGood) / (1 - ((double)totBad / N));
   }
   
   return res = dp[0] > 1e50 ? -1.0 : dp[0];
}

TC TCO05 Round 1 D1 500 (FibonacciSum)

Category: Greedy, Dynamic Programming
URL: http://www.topcoder.com/stat?c=problem_statement&pm=4629

The problem asks us to compute the minimum number of Fibonacci numbers required to sum up to a given number n. If you are familiar with Zeckendorf's Theorem then you might be able to see a greedy strategy. Zeckendorf's Theorem states that every positive integer can be represented in a unique way as the sum of one or more distinct Fibonacci numbers such that it does not contain two consecutive Fibonacci numbers. A more obvious approach is to use dynamic programming, let D(n) denote the minimum number of Fibonacci terms to make up a sum of exactly n. Then we can define a recurrence relation such as: F(n) = min { 1 + F(n - T[k]) } where T[k] denotes the k-th Fibonacci term. Obviously F(n - T[k]) is only defined for n - T[k] >= 0. Additionally, the base case is defined as F(0) = 0, which acts as our stopping point as this means we have no remaining sum to make up for. Coding this is relatively simple:

class FibonacciSum {
public:
   int howMany(int);
};

int tab[101];
int memo[1000001];
int cnt;

#define INF (1 << 20)

int func(int n) {
   int res = INF;
   if (n == 0) return 0;
   if (memo[n] != -1) return memo[n];
   for (int i = cnt; i >= 0; i--) {
      if (n - tab[i] >= 0) res <?= 1 + func(n - tab[i]);
   }
   return memo[n] = res;
}

int FibonacciSum::howMany(int n) {
   tab[0] = tab[1] = 1;
   cnt = 2;
   memset(memo,-1,sizeof(memo));
   for (int i = 2; tab[i-1] + tab[i-2] < 1000000; i++, cnt++) {
      tab[i] = tab[i-1] + tab[i-2];
   }
   cnt--;
   return func(n);
}


We can also translate this into a bottom-up dynamic programming algorithm. To do this, we must consider the order of calculation - this is dictated by the recurrence relation we defined earlier. As this is a 1 dimensional DP problem, it's relatively simple: given an n value, we need to compute all values below n before considering n. This can be seen by the term inside the recursive call (n - T[k]), as T[k] is strictly positive for all k values we are guaranteed to only require all values below n to be calculated beforehand.

class FibonacciSum {
public:
   int howMany(int);
};

int dp[1000001];
int tab[101];

#define INF (1 << 20)

int FibonacciSum::howMany(int n) {
   tab[0] = tab[1] = 1;
   int cnt = 2;
   for (int i = 2; tab[i-1] + tab[i-2] < 1000000; i++, cnt++) {
      tab[i] = tab[i-1] + tab[i-2];
   }
   dp[0] = 0;
   for (int i = 1; i <= n; i++) {
      dp[i] = INF;
      for (int k = 1; k < cnt; k++) {
         if (i - tab[k] < 0) break;
         dp[i] = min(dp[i], 1 + dp[i - tab[k]]);
      }
   }
   return dp[n];
}

Friday, August 28, 2009

TC SRM447 D1 500 (PeopleYouMayKnow)

The problem statement can be accessed via:
http://www.topcoder.com/stat?c=problem_statement&pm=10580

The problem itself is pretty simple: it asks us to given a friends graph determine the number of friends to remove such that there are no paths from friend A to friend B which has a length less than 3.

The key observation is to remove all friends which have a shortest path between A and B which exceeds 3 as these will never conflict with our objective. To determine intermediate distances we can just use Floyd-Warshall's algorithm to calculate all-pairs shortest paths between friends. Then we can discard any nodes which violate dist[person1][i] + dist[i][person2] <= 3. Note that the problem fixes the constraint to not allow any paths with a length of less than 3 - this allows us to use a simple bipartite matching to determine the minimum cut of the network. For the general problem of having a length less than n, we can simply use a network flow algorithm (which we use for the implementation here).

We can divide the friends into two (one used as an in-node, the other used as an out-node) and run the maximum flow algorithm over the network. This concept is very similar to the one used in TC SRM360 D1 500 (Prince of Persia).

Implementation below:


class PeopleYouMayKnow {
public:
int maximalScore(vector <string>, int, int);
};

int augment(int cur, int flow, int ret);

int adjmat[201][201];
int used[201][201];
bool visited[201];
int t;

int maxflow() {
int flow = 0, thisFlow = 0;
memset(visited,false,sizeof(visited));
while ((thisFlow = augment(0, INT_MAX-1, 0)) != 0) {
memset(visited,false,sizeof(visited));
flow += thisFlow;
}
return flow;
}

int augment(int cur, int flow, int ret) {
if (cur == t-1) return flow;
visited[cur] = true;
for (int i = 0; i < t && ret == 0; i++) {
if (adjmat[cur][i] > used[cur][i] && !visited[i] &&
(ret = augment(i, min(flow,adjmat[cur][i]-used[cur][i]), 0)) > 0)
used[cur][i] += ret;
else if (used[i][cur] > 0 && !visited[i] &&
(ret = augment(i,min(flow, used[i][cur]), 0)) > 0)
used[i][cur] -= ret;
}
return ret;
}

int distTable[51][51];
int validNodes[51];
#define INF (1<<20)

int PeopleYouMayKnow::maximalScore(vector <string> friends, int A, int B) {
int res = 0;
int sz = friends.size();
for (int i = 0; i < sz; i++)
for (int j = 0; j < sz; j++)
if (i == j) distTable[i][j] = 0;
else distTable[i][j] = friends[i][j] == 'Y' ? 1 : INF;
for (int k = 0; k < sz; k++)
for (int i = 0; i < sz; i++)
for (int j = 0; j < sz; j++)
distTable[i][j] <?= distTable[i][k] + distTable[k][j];
// build the network - allowing infinite capacities between
// out -> in nodes and restricting capacities between nodes with 1 capacity
for (int i = 0; i < sz; i++)
for (int j = 0; j < sz; j++)
if (friends[i][j] == 'Y')
adjmat[i+1+sz][j+1] = INF;
for (int i = 0; i < sz; i++)
if (distTable[person1][i] + distTable[i][person2] > 3)
for (int j = 0; j < sz; j++) adjmat[j+1][i+1] = adjmat[i+1][j+1] = 0;
else
validNodes[i] = 1;
for (int i = 0; i < sz; i++)
if (validNodes[i]) adjmat[i+1][i+1+sz] = 1; else adjmat[i+1][i+1+sz] = 0;
adjmat[0][person1+1] = INF;
adjmat[person1+1][person1+1+sz] = INF;
adjmat[person2+1][person2+1+sz] = INF;
adjmat[person2+1+sz][2*sz+1] = INF;
t = 2*sz+2;
return res = maxflow();
}

Thursday, August 27, 2009

TC SRM 447 D2 900 (ImportsList)

The problem statement can be viewed:
http://www.topcoder.com/stat?c=problem_statement&pm=10579

This problem asks us to calculate the minimum number of scripts to be used within an import list. For a given script A, it needs access to a set of scripts S1, S2, ..., Sn. The relationships are transitive, so if script A has access to script B and script B has access to script C then script A also has access to script C. We must choose a combination such that a given script has access to all its required scripts whilst also choosing the minimum number of scripts to be imported.

The first observation to make is that the graph is directed, so access does not flow both ways. In fact, given the current constraints there are no cycles in the graph hence we have a directed acyclic graph (DAG). This gives us a huge hint to topologically sort (or linearise) the graph before doing further operations. A topological sort can be achieved with a simple DFS and recording the finished nodes in reverse order. Alright let's assume we now have a linearised DAG - how do we choose the minimum amount? Recall that for an DAG we can order the vertices such that all edges go from left to right. Therefore for a vertex V, we can greedily consider vertices from left to right and keeping track of which vertices can be accessed if we access vertex W.

This works for several reasons. Let's consider two vertices: W and X. Let's say that vertex X is to the right of vertex W. If vertex X is reachable from vertex W then vertex W can access at least as much as vertex X since we can simply access vertex X from W using an arbitrary set of traversals. Let's say that vertex X is not reachable from vertex W, then there could be some scripts which we can access between vertex W and X which we need to fulfill the requirements. If there isn't, then we can simply not use vertex W and we will consider vertex X eventually. Therefore, the optimality follows through the linearised order since we consider the largest set and choosing so does not cause conflicts later on.

A very simple way to represent what is accessible from a given vertex is to use a bitmask. Then the implementation is as follows:


class ImportsList {
public:
vector <int> importsCount(vector <string>);
};

#define N 52
int visited[N];
int order[N];
int cycleDetect;
int adjmat[N][N];
int dfsIdx;

void dfs(int node) {
if (visited[node]) { if (visited[node] == 1) cycleDetect = 1; return; }
visited[node] = 1;
for (int i = 0; i < N; i++) {
if (adjmat[node][i]) {
dfs(i);
}
}
visited[node] = 2;
order[dfsIdx] = node;
dfsIdx--;
}

long long mask[N];
long long procMask[N];

vector <int> ImportsList::importsCount(vector <string> requires) {
vector<int> res(requires.size(),0);
for (int i = 0; i < requires.size(); i++) {
long long m = 0;
for (int j = 0; j < requires.size(); j++) {
if (requires[i][j] == 'Y') { adjmat[i][j] = 1; m += (1LL << j); }
}
mask[i] = m;
}
dfsIdx = requires.size()-1;
for (int i = 0; i < requires.size(); i++) {
if (!visited[i]) {
dfs(i);
}
}
// greedy based on topological sorting
for (int i = requires.size()-1; i >= 0; i--) {
long long m = mask[order[i]];
long long p = m;
int num = 0;
for (int j = i+1; j < requires.size() && m != 0; j++) {
if (adjmat[order[i]][order[j]] && (m & (1LL << order[j]))) {
num++;
p |= procMask[order[j]];
m = m & ~(procMask[orderj]]);
}
}
procMask[order[i]] = p;
res[order[i]] = num;
}
return res;
}

Friday, August 21, 2009

TC MPR1 D1 550 (CantorDust)

The problem statement can be viewed:
http://www.topcoder.com/stat?c=problem_statement&pm=10418&rd=13935

This problem is pretty much exactly the same as the D2 1000. However, there is one additional case we need to handle - patterns which are empty. This is fairly simple to translate from our previous algorithm/implementation. We simply need to check if the pattern supplied is empty (i.e. no X's) and if so make a special case computation.

Now the question is, how do we extend our original algorithm to cover this case? Take note that if there were no X's in the fractal then the number of different occurrences of our blank pattern would be:

(size - row size of the pattern + 1) * (size - column size of the pattern + 1)

However there are X's which act as obstacles from reaching this upper bound. So we simply calculate the number of occurrences where our pattern does not fit into the fractal. This is rather simple, just calculate the number of occurrences in which our pattern does fit and subtract it from the total number of occurrences for the associated row/column.

The total number of occurrences for a row is: size - row size of the pattern + 1
The total number of occurrences for a column is: size - column size of the pattern + 1
The number of occurrences where the pattern doesn't fit (per row): (size - row size of the pattern + 1) - pattern does fit in the row
The number of occurrences where the pattern doesn't fit (per col): (size - column size of the pattern + 1) - pattern does fit in the column

Using the same multiplication principle, the total number of occurrences where the pattern doesn't fit in the 2D space is simply the product of the two. Lastly, we just need to exclude this many occurrences from our upper bound calculation to yield the answer:

Total number of occurrences of an empty pattern: ((size - row size of the pattern + 1) * (size - column size of the pattern + 1)) - ((size - row size of the pattern + 1) - pattern does fit in the row) * ((size - column size of the pattern + 1) - pattern does fit in the column)

The implementation is as follows:


class CantorDust {
public:
int occurrencesNumber(vector <string>, int);
};

string tab[10];
int visited[10];

string build(int num) {
if (num == 0) return "X";
if (visited[num]) return tab[num];
visited[num] = 1;
return tab[num] = build(num-1) + string(pow(3.0,(double)num-1),'.') + build(num-1);
}

int CantorDust::occurrencesNumber(vector <string> pattern, int t) {
int res = 0;
string str = build(t);
string emptyRow = string(pattern[0].size(),'.');
string emptyCol = string(pattern.size(),'.');
string rows = "";
int size = pow(3.0,(double)t);
for (int i = 0; i < pattern.size(); i++) {
if (pattern[i] == emptyRow) continue;
if (rows != "" && rows != pattern[i]) return 0;
rows = pattern[i];
}
string cols = "";
for (int i = 0; i < pattern[0].size(); i++) {
string col;
for (int j = 0; j < pattern.size(); j++) {
col += pattern[j][i];
}
if (col == emptyCol) continue;
if (cols != "" && cols != col) return 0;
cols = col;
}
bool isEmpty = false;
if (cols == "" && rows == "") {
cols = emptyCol;
rows = emptyRow;
isEmpty = true;
}
int numMatchRow = 0, numMatchCol = 0;
string::size_type pos = 0;
while ((pos = str.find(rows,pos)) != string::npos) {
numMatchRow++;
pos++;
}
pos = 0;
while ((pos = str.find(cols,pos)) != string::npos) {
numMatchCol++;
pos++;
}
if (isEmpty) {
int exclude = (size - emptyRow.size() + 1 -
numMatchRow) * (size - emptyCol.size() + 1 - numMatchCol);
return (size - emptyRow.size() + 1) *
(size - emptyCol.size() + 1) - exclude;
}
return res = numMatchRow * numMatchCol;
}

Thursday, August 20, 2009

TC Member Pilot Round 1 Division 2

This is the first member-driven "Single Round Match" for Topcoder (which also isn't rated). Division 2 had an interesting set of problems which proved to be slightly more difficult than the norm.

Division 2 250-pointer (TriFibonacci)
This problem can be viewed:
http://www.topcoder.com/stat?c=problem_statement&pm=10587&rd=13935

Basically this problem asks us to "fill in the blank" of a "Tri-Fibonacci" sequence which is defined as F(a) = F(a-1) + F(a-2) + F(a-3). The small catch here is that there can be sequences where this isn't well-defined (i.e. invalid). All the integers in the sequence must be strictly positive (i.e. zeroes are not allowed).

To solve this problem there are several options:
1. Brute force the blank and check whether or not it leads to a valid Tri-Fibonacci sequence. Due to the constraints this would take around 20 million operations which should run in time. If we can't find a valid sequence then we return -1.

2. Compute the blank based on the sequence given. Roughly speaking, we can divide this into two cases. One is where the number is within the first three terms and hence need to deduce the value by taking the sum of the other 2 and subtracting it from the next Tri-Fibonacci number.

For example, given {2, -1, 3, 7}. We compute 2+3=5, and then we can deduce the blank by noting 5+x=7 (by definition of the recurrence).

The other case is when the missing term is on the fourth term or beyond. Then it's easy to just use the given recurrence and compute the missing number from the last 3 numbers in the sequence before the blank. We then proceed to do a check on whether or not the unique number we obtained from the sequence is a valid Tri-Fibonacci. If it is, we accept and return the value. Otherwise we return -1.

The solution below uses approach 2:


class TriFibonacci {
public:
int complete(vector <int>);
};

bool checkFib(const vector<int>& A) {
for (int i = 3; i < A.size(); i++) {
if (A[i] != A[i-1] + A[i-2] + A[i-3]) return 0;
}
return 1;
}

int TriFibonacci::complete(vector <int> A) {
int pos = 0;
for (int i = 0; i < A.size(); i++) {
if (A[i] < 0) { pos = i; break; }
}
if (pos < 3) {
int s = 0;
for (int i = 0; i < 3; i++) if (A[i] >= 0) s += A[i];
if (A[3] - s <= 0) return -1;
A[pos] = A[3] - s;
return checkFib(A) ? A[pos] : -1;
}
A[pos] = A[pos-1] + A[pos-2] + A[pos-3];
return checkFib(A) ? A[pos] : -1;
}


Division 2 500-pointer (ErdosNumber)
This problem can be viewed:
http://www.topcoder.com/stat?c=problem_statement&pm=10377&rd=13935

This problem asks us to compute the Erdos Number based on a given set of publications. If you are interested in this number there are some facts on Wikipedia: http://en.wikipedia.org/wiki/Erd%C5%91s_number.

Again like most problems there are several ways to solve this. The bulk of the problem is really just parsing the string into a suitable data structure. This can be done using a map which maps a string (author names) and converts it into an integer index (for the adjacency matrix). We simply read each publication at a time and if it's a new author name we insert it into the map. Once we have read all the authors for a given publication we iterate in pairs and assign a weight of 1 to each edge between the authors. Then we simply need to run a shortest path algorithm from Erdos to all the other authors. There are several options available:

1. Most simplest is to use a Floyd-Warshall algorithm. This runs in O(n^3) time and since there are at most 100 vertices, it will easily run in time.
2. Use a breadth-first search from Erdos keeping track of the depth. This works as all the edges have a weight of 1. Reasonably simple to code and should be of choice if you don't know the Floyd-Warshall's All-Pairs shortest path algorithm. It's also more efficient.
3. Alternative methods include using Dijkstra's from Erdos to all vertices (Warshall's is more efficient and easier to code). You can also use Bellman-Ford's algorithm which will get the answer from the Erdos vertex to all other vertices.

A small implementation note is to keep track if a vertex is not reachable from the Erdos vertex.

See the below implementation which uses Floyd-Warshall's algorithm:


class ErdosNumber {
public:
vector <string> calculateNumbers(vector <string>);
};

#define INF (1<<20)
int adjmat[101][101];

vector <string> ErdosNumber::calculateNumbers(vector <string>
publications) {
vector<string> res;
map<string,int> authorMapping;
int authorIdx = 0;
for (int i = 0; i < 101; i++) for (int j = 0; j < 101; j++)
adjmat[i][j] = i == j ? 0 : INF;
for (int i = 0; i < publications.size(); i++) {
istringstream iss(publications[i]);
string author;
vector<string> assocs;
while (iss >> author) {
if (authorMapping.count(author) == 0) {
authorMapping.insert(make_pair(author,authorIdx++));
}
assocs.push_back(author);
}
for (int j = 0; j < assocs.size(); j++) {
for (int k = j+1; k < assocs.size(); k++) {
int first = authorMapping[assocs[j]];
int second = authorMapping[assocs[k]];
adjmat[first][second] = adjmat[second][first] = 1;
}
}
}
int erdosIdx = authorMapping["ERDOS"];
for (int k = 0; k < authorMapping.size(); k++) {
for (int i = 0; i < authorMapping.size(); i++) {
for (int j = 0; j < authorMapping.size(); j++) {
adjmat[i][j] <?= adjmat[i][k] + adjmat[k][j];
}
}
}
for (map<string,int>::iterator it = authorMapping.begin();
it != authorMapping.end(); it++) {
ostringstream oss;
oss << it->first;
int dist = adjmat[erdosIdx][it->second];
if (dist < INF) oss << " " << dist;
res.push_back(oss.str());
}
return res;
}


Division 2 1000-pointer (CantorDustEasy)
This problem can be viewed:
http://www.topcoder.com/stat?c=problem_statement&pm=10588&rd=13935

Finally the Division 2 Hard problem! This is a somewhat tricky problem and takes some hand-drawn diagrams to really see the picture. The problem is as given: given a recursive fractal algorithm determine the number of occurrences of a specified pattern. Given the low constraints of the pattern (50x50 at most) we can try and calculate the fractal up to 3^4 = 81 > 50 and attempt to recursively build on the answer. Why do we only need up to a time of 4? As each black square "sub-section" is divided into the four corners there is an empty space between them made up of 3^(t-1) blocks. This means that no pattern can possible "cross" over and given the constraint there must be at least 1 black square in the pattern - it makes this approach feasible. This is somewhat difficult to implement though so we attempt to make it easier.

On a side note, you can't naively enumerate all the patterns in a given fractal. The size at 9 is 3^18 which is nearly 400 million. We can't store this much let alone compute pattern matching on each cell which would certainly TLE.

So there is an easier approach but requires some effort to discover. If we draw a few diagrams we can see that all rows are either the same or they are empty (made up entirely of white cells). This is due to the symmetric property of the fractal, likewise - the columns are either the same or they are empty.

But how does this help us count the number of occurrences of a pattern? Firstly it helps eliminate infeasible patterns - if the non-empty rows/columns don't match then we can immediately discard it as a feasible pattern. The next observation we need to note is that our pattern now has either rows and columns that are the same or they are empty. Let's take a non-empty row R, since we know that all of the non-empty rows in the actual fractal F are the same. We can just count the number of R's we see in F. How does this help?

To illustrate the concept, take "X.X...X.X" as a sample row from F. Our pattern's row is: ".X.", using a string finding algorithm we find that it occurs twice in F. Hence if our original (non-collapsed) pattern only had 1 row of ".X." then we immediately retrieve the answer we have. Note this is a 1 dimensional problem - we need to convert this to our 2 dimensional problem. So we do the same thing with columns and determine how many occurrences. Let's say the count we have for the rows is "x" and the count we have for columns is "y". Then by the multiplication principle, the total number of occurrences in the 2D space is simply x * y (draw some diagrams to convince yourself - if we have non-empty rows that have x occurrences each then the total number is just how many times these x occurrences occur in a vertical fashion from the column pattern).

Combining this together we yield a simple implementation:


class CantorDustEasy {
public:
int occurrencesNumber(vector <string>, int);
};

string tab[10];
int visited[10];

string build(int num) {
if (num == 0) return "X";
if (visited[num]) return tab[num];
visited[num] = 1;
return tab[num] = build(num-1) + string(pow(3.0,(double)num-1),'.') + build(num-1);
}

int CantorDustEasy::occurrencesNumber(vector <string> pattern, int t) {
int res = 0;
string str = build(t);
string emptyRow = string(pattern[0].size(),'.');
string emptyCol = string(pattern.size(),'.');
string rows = "";
for (int i = 0; i < pattern.size(); i++) {
if (pattern[i] == emptyRow) continue;
if (rows != "" && rows != pattern[i]) return 0;
rows = pattern[i];
}
string cols = "";
for (int i = 0; i < pattern[0].size(); i++) {
string col;
for (int j = 0; j < pattern.size(); j++) {
col += pattern[j][i];
}
if (col == emptyCol) continue;
if (cols != "" && cols != col) return 0;
cols = col;
}
int numMatchRow = 0, numMatchCol = 0;
string::size_type pos = 0;
while ((pos = str.find(rows,pos)) != string::npos) {
numMatchRow++;
pos++;
}
pos = 0;
while ((pos = str.find(cols,pos)) != string::npos) {
numMatchCol++;
pos++;
}
return res = numMatchRow * numMatchCol;
}

Monday, August 17, 2009

TC SRM443 D1 600 (BinaryFlips)

This problem can be accessed via:
http://www.topcoder.com/stat?c=problem_statement&pm=10387

A fairly standard problem with somewhat tricky implementation due to the problem constraints. The problem gives us "A" number of zeroes and "B" and number of ones and using a swap of exactly "K" digits in each turn, determine the minimum number of moves to get it into a state of all ones (and returning -1 if it's impossible to do so). Looking at the constraints of around 100,000 means we can't naively try all combinations. A simple approach to take is to Breadth-First searching the moves and making sure we don't re-visit already processed nodes with lower number of moves (as they will never be optimal).

We need to make several observations to help simplify the problem. As we can choose any K digits to swap over in a given turn, this means we can view the problem as a set of 0's followed by a set of 1's for simplicity. This helps us determine that the state required for the BFS is simply the number of zeroes (or ones) and not a string of disordered 0's and 1's. A further optimisation is to note that the number of zeroes can be derived from the number of ones as the total size of the string does not change. Hence the maximum state space is simply 200,000 so we are fine for memory usage. The only problem left is the transitions between states inside the inner cycle of the BFS. As the K is as large as 100,000 we will run into time limit issues if we aren't efficient in which states we check/traverse to.

This is perhaps the trickest part of the problem - the other parts are fairly standard and simple. Now we ask ourselves, "is there a way to make sure we only consider non-visited nodes?". Using a set data structure which keeps all its elements unique allows us to only traverse the unvisited states. To make another optimisation we need to ensure that we only visit the reachable unvisited states - we calculate the state range in which we can reach and do a subset traversal of the set (instead of going through the whole set and determining whether it's a valid move). To determine the range of states we need to see that if we can traverse to x-number of 1's and we can also traverse to y-number of 1's then we can also traverse to all number of 1's between x and y with the same parity. To see this, arrange a state in a sequence of 0's and 1's and take an arbitrary K (K = 4 in this case):

Move to 9 1's: 0[0000]11111 -> 0[1111]11111
Move to 7 1's: 00[0001]1111 -> 00[1110]1111
Move to 5 1's: 000[0011]111 -> 000[1100]111
Move to 3 1's: 0000[0111]11 -> 0000[1000]11
Move to 1 1's: 00000[1111]1 -> 00000[0000]1

As you can see, we can only move a specific parity (odd/even) between the minimum and maximum range. This can be seen above where we use a sliding "window" mechanism for the inversion. The parity is determined by a combination of the window size (i.e. K) and the parity of the state. To see this connection, we build an informal proof:

Let the K-sized window be denoted by W.
Let P be the parity of the state S, where S is a sequence of 0's followed by 1's.
Then we define x to be the number of 0's in S and y to be the number of 1's in S.

By definition, P = y (mod 2). Furthermore we need to define x' and y' to be the number of 0's and 1's in the window W respectively. We take x'' and y'' to be the number of 0's and 1's in the state S excluding the digits inside window W. Then it follows that:

x'' = x - x' (mod 2)
y'' = y - y' (mod 2)
x' + y' = K

Parity of the inversed window W = x' (mod 2) since it's the number of 0's in the window (as they turn into 1's after the transform which becomes the parity). We now calculate P', the parity of S after being transformed by an arbitrary K-sized window:

P' = x' (mod 2) + y'' (mod 2) (parity of the inversed window plus the region of S - W)
P' = K - y' (mod 2) + y'' (mod 2)
P' = K - y' + y'' (mod 2)
P' = K - y' + 2y' + y'' (mod 2) (add 2y' - this can be done since it doesn't alter the odd/even-ness of the equation)
P' = K + y' + y'' (mod 2)
P' = K + y (mod 2)

Therefore the next state's parity is P' = K + y (mod 2). So we only need to compute this to determine whether the parity of the next state is odd/even. So it suffices to just calculate the minimum and maximum number of 1's and then traverse based on the parity. Using all the concepts above gives us an implementation which runs at around 100ms worst case:


class BinaryFlips {
public:
int minimalMoves(int, int, int);
};

int visited[200011];

int BinaryFlips::minimalMoves(int A, int B, int K) {
if (A == 0) return 0;
if (A + B < K) return -1;
queue<int> q;
int N = A + B;
memset(visited,-1,sizeof(visited));
set<int> evenVisit, oddVisit;
for (int i = 0; i <= N; i++) {
if (i & 1) oddVisit.insert(i);
else evenVisit.insert(i);
}
if (B & 1) oddVisit.erase(oddVisit.find(B));
else evenVisit.erase(evenVisit.erase(B));
q.push(B);
visited[B] = 0;
while (!q.empty()) {
int t = q.front(); q.pop();
if (t == N) return visited[t];
int numZeroes = N - t;
int numOnes = t;
int x = numOnes < K ? K - numOnes : numOnes - K;
int y = numZeroes < K ? numOnes + numZeroes * 2 - K : numOnes + K;
int s = min(x,y);
int st = max(x,y);
set<int>* visitSet = (t+K) & 1 ? &oddVisit : &evenVisit;
if (visitSet->size() == 0) continue;
set<int>::iterator it = visitSet->lower_bound(s);
while (it != visitSet->end() && *it <= st) {
if (visited[*it] < 0) {
q.push(*it);
visited[*it] = 1 + visited[t];
set<int>::iterator itTemp = it; itTemp++;
visitSet->erase(it);
it = itTemp;
} else {
it++;
}
}
}
return -1;
}