[TECH] Simple comparision of real numbers

We know that the fixed precision representation of real numbers is done by mapping them into integer space I from real space R. However the underlying integer representation of the 32-bit floats and 64-bit doubles is a sign magnitude representation (that is the reason why we see -0.00 during output of some numerical algorithms). People seem to be adding a constant to convert the underlying sign magnitude form to 2's compliment. However following is a more logical conversion (since we know that 2's compliment is 1's compliment +1 , the motivation of 2's compliment is to avoid the -0 as we know). The following comparison routine compares two floating point numbers (add the tolerance in the terms of the number the real numbers to make it approximate comparison).

/*Less than or greater than*/
char CompareFloats(float a, float b){
    int aint = *(int *)&a;
    int bint = *(int *)&b;
    aint = (aint & (1UL<<31))? ~(aint^(1UL<<31))+1:aint;
    bint = (bint & (1UL<<31))? ~(bint^(1UL<<31))+1:bint;
    aint -= bint;
    printf("%f %s %f \n",a,(aint <0)?"<":">=",b);
    printf("There are %d real numbers \n",(aint<0)?-aint:aint);
    printf("between %f and %f (in 32 bit float representation)\n",a,b);
}

[TECH] Launchpad and Bazaar

I have now started to move all my development efforts to launchpad and bazaar. I really found launchpad a great place to maintain and release my work. I tried several alternatives before I moved to launchpad, I tried to use code.google but unfortunately its not as rich as launchpad and bazaar. I have released version 0.9 of libEditScript project whose aim to build a high performance and space efficient sequence alignment library. I have charted out the features for the next release, one thing I'm sure people would like to use is a JNI (Java Native Interface) so that people can use this space efficient code in their java code. I have seen that this problem of the need of edit script arises in several occasions, I have seen people using the space inefficient version O(n^2) of edit script computation. This is were we gain significantly in terms of the space saving. Also I have one more idea to make it more parallel in the sense I want to use some concurrent techniques to make this parallel, especially since I have circular queue I can easily use several threads to make this concurrent which I will add as a next step in the next release..

Checkout the libEditScript at https://launchpad.net/libeditscript. Next in my list would be to make a release of the weightedmatching library.

[TECH] A non-recursive algorithm to compute the Edit Script in O(min(n1,n2)+log(min(n1,n2)) space

Recently I was working on an area efficient and synthesizable design to compute the edit script (minimum cost sequence of INSERT, DELETE and CHANGE) operations to transform string S₁ to S₂. After a little bit of thought I have the following idea to build a non-recursive version of the Hirschberg's algorithm, since the algorithm is non-recursive we can build an efficient digital circuit with this idea. We use a simple circular queue and apply DFS (Depth First Search) and we can prove that the capacity of this queue at any stage of the algorithm is θ(log(min(n₁,n₂)). The proof is simple we choose the geometrically decreasing string which has smaller length from the given strings(min(n₁,n₂)). Since we do a depth first search and the depth of subproblem tree is θ(log(min(n₁,n₂)), so we will have atmost θ(log(min(n₁,n₂)) subproblems in the circular queue at any stage of the algorithm.

The core non-recursive algorithm starts off with the initial problem and finds the minimum alignment split ('q') between (S₁[1.....n₁/2] and S₂ creates two sub problems and appends in the circular queue (currently it uses BFS). Download the linear space implementation from here Following is the outline of the core non recursive implementation.

     91     CQueue bfs_list;
     92     ESSubProblem es_sub;
     93 
     94     InitializeCQ(&bfs_list,sizeof(ESSubProblem));
     95     /*Put the toplevel subproblem into the CQueue*/
     96     es_sub.sub_s1 = s1; es_sub.sub_n1 = n1;
     97     es_sub.sub_s2 = s2; es_sub.sub_n2 = n2;
     98 
     99     AppendCQ(&bfs_list,&es_sub);
    100     while(DequeueCQ(&bfs_list,&es_sub)){
    101         sub_s1 = es_sub.sub_s1; sub_s2 = es_sub.sub_s2;
    102         sub_n1 = es_sub.sub_n1; sub_n2 = es_sub.sub_n2;
    103 
    104         if(sub_n1 <=1){
    105             /*If sub_n1 is 1 or 0 we cannot split it any more*/
    106             EditDistanceLS(sub_s1,sub_n1,sub_s2,sub_n2);
    107             i = sub_n1; j = sub_n2;
    108             /*Figure out the edit operations*/
    109             while(i || j){
    110                 op = GetOp(i,j,(i?sub_s1[i-1]:'\0'),
    111                     (j?sub_s2[j-1]:'\0'));
    112                 switch(op){
    113                     case 'I':
    114                             EditScriptS2[(&sub_s2[j-1]) - s2] = 'I';
    115                             j--;
    116                             break;
    117                     case 'D':
    118                             EditScriptS1[(&sub_s1[i-1]) - s1] = 'D';
    119                             i--;
    120                             break;
    121                     case 'C':
    122                             EditScriptS1[(&sub_s1[i-1]) - s1] =
    123                                 (sub_s1[i-1] == sub_s2[j-1])?':':'C';
    124                                 i--; j--;
    125                 }
    126             }
    127         }else {
    128             /*Add the two subproblems*/
    129             q = FindMinSplit(sub_s1,sub_n1,sub_s2,sub_n2);
    130             /*SUB PROBLEM 1*/
    131             es_sub.sub_n1 = CELING(sub_n1,2);
    132             es_sub.sub_n2 = q;
    133             AppendCQ(&bfs_list,&es_sub);
    134 
    135             /*SUB PROBLEM 2*/
    136             es_sub.sub_s1 = &(sub_s1[CELING(sub_n1,2)]);
    137             es_sub.sub_s2 = &(sub_s2[q]);
    138             es_sub.sub_n1 = sub_n1/2;
    139             es_sub.sub_n2 = sub_n2-q;
    140             AppendCQ(&bfs_list,&es_sub);
    141         }
    142     }

[TECH] A Graph with θ(log(n)) approximation ratio for greedy vertex cover

We know that the greedy algorithm (which picks the vertex with next highest degree) to solve the vertex cover problem is not optimal. However we have 2-Approximation algorithm for the vertex cover problem (based on a integer linear programming formulation). Actually it was not very straight forward to for me prove that the greedy algorithm has an approximation ratio greater than 2. However after some thinking the following construction indeed proves that the greedy algorithms approximation ratio is > 2. In fact we can easily see (from the figure) that the greedy algorithm will first pick the node with highest degree i.e n (the magenta). Then it deletes all the edges adjacent on the magenta node and continues and clearly the greedy algorithm will have all the nodes in the blue bounding box. However its clear that if we choose the nodes in the oval region we get a minimum vertex cover of size n. The size of greedy vertex cover is n/2 + n/3 + n/4 .... n/n ≈ θ(nlog(n)) (for large n) and hence the approximation ratio for this graph is θ(log(n)).

[TECH] On The Uniqueness Of a Perfect Match In Undirected Acyclic graphs.

I was studying about several results related to Tiling Theory. Tiling Theory adds formalism to the question " Given a simple rectangular region (R) can we find a valid tiling by using tiles from a set (T) ". If our tile set (T) is a singleton set with a unit square then we can tile any region. However what if the the tile set contains only a 2x1 rectangle (domino)? what if the tile set contains a set of polymino's (tetris game)?. Although the question may seem very simple there is a rich combinatorial work based on this question. In fact the study about the tiling problem originated from Wang's Conjecture about periodic tiling. Actually the algorithmic research about tiling theory was dormant for quite some time, even the proof for 2x1 dominoes is fairly recent see this . I was reading about the proof and was really impressed by its elegance. Often people consider simple things useless or don't have enough respect for simple facts, however I always find that good results are based on several very simple facts. In fact I would like to highlight simple ideas in every algorithm I study about.

In this proof to prove that the problem is NP-complete for 2x1 domino's a simple fact that " There is exactly only one perfect matching (if at all one exists) in a Tree " is used. This can be proved intuitively as follows.

Let M be the perfect match in the Tree T, if M is not a unique perfect match then there should be some other perfect match lets call it M^'. If we observe the leaf nodes in T which have only one edge adjacent on them so both M and M^' should agree on the matched edges which have leaf nodes of T as the vertices . We can now delete all the leaf nodes and their adjacent edges and continue the argument until we don't have any leaf nodes. Thus from this argument its clear that both M and M^' are the same and hence there is only a unique perfect match if at all if one exists in a tree T.

[TECH] Converting Market Matrix Sparse Representation to Row Compressed Sparse Representation

I wanted to use some good benchmarks to illustrate the performance gains of the Fast Dual Update Algorithm for pre-ordering sparse matrices. University of Florida Sparse Library has a rich collection of sparse matrices from various domains. Unfortunately all the matrices are in Market Matrix format but not in Row Compressed format which I currently use.

The following is a handy script which converts the Market Matrix format to Row Compressed Format mm2rc.pl. The script uses a simple logic based on external sorting.

[TECH] Computing Eigen Values and Eigen Vectors using LAPACK

I have used the LAPACK to compute the eigen values, till recently I needed the eigen vectors both left and right. I thought this could be really handy if you need to compute eigen values and eigen vectors. Download both library and this file from here

/*Computing both the eigen values and eigen vectors using LAPACK.
 *vamsi.krishnak (at) gmail (dot) com
 **/
#include "f2c.h"
#include "clapack.h"
#include<stdio.h>
#include<math.h>
#include<malloc.h>
#include<assert.h>

int main(){
    integer dim=0;
    unsigned int i;
    doublereal *A = NULL;
    real *eig_values = NULL;
    real *work;
    char JOBVL='V'; /*Left eigen vectors*/
    char JOBVR='V'; /*No right eigen vectors*/
    integer iwork[1];
    integer lwork,liwork=1,info,lda;
    doublereal *WR=NULL;
    doublereal *WI=NULL;
    doublereal *VL=NULL;
    doublereal *VR=NULL;
    integer LDVR=1;
    integer LDVL=1;
    integer LDA;
    doublereal *WORK=NULL;
    integer LWORK;


    printf("Please Enter the dimension of the array:\n");
    scanf("%d",&dim); LDA=dim; lda=dim; LDVL=dim;LDVR=dim;
    printf("Please enter the Matrix A\n");
    eig_values = (real *) malloc(sizeof(real)*dim);
    lwork = 10*dim; 
    work = (real *) malloc(sizeof(real)*lwork);
    A = (doublereal *) malloc(sizeof(doublereal)*(dim*dim));
    printf("Please Enter the elements of the matrix\n");
    i = 0;
    
    /*Double real and imaginary parts*/
    WR = (doublereal *)malloc(sizeof(doublereal)*dim);
    WI = (doublereal *)malloc(sizeof(doublereal)*dim);
    VL = (doublereal *)malloc(sizeof(doublereal)*(dim*dim));
    VR = (doublereal *)malloc(sizeof(doublereal)*(dim*dim));
    LWORK = 6*dim;
    WORK = (doublereal *)malloc(sizeof(doublereal)*LWORK);  
    assert(WI && WR && VL && VR);
    /**************/
    do{
        scanf("%lf",&A[i]);
        i++;
    }while(i<(dim*dim));
    /*Compute the Eigen values*/
    //ssyev_("N","U",&dim,A,&dim,eig_values,work,&lwork/*,iwork,&liworki*/,&info);
    dgeev_(&JOBVL,&JOBVR,&dim,A,&LDA,WR,WI,VL,&LDVL,VR,&LDVR,WORK,&LWORK,&info);
    if(!info){
        for(i=0;i<dim;i++){
            printf("%lf+i%lf\n",WR[i],WI[i]);
        }
        printf("===Right Eigen Vectors===\n");
        int k; int conj;
        for(i=0;i<dim;i++){
            /*Print i^th right eigen vector*/
            if(fabs(WI[i] > 0)){
                /*eigen vector i and eigen vector i+1*/
                for(conj=0;conj<2;conj++){
                    printf("[ ");
                    for(k=0;k<dim;k++){
                        if(!conj){
                            printf("(%lf)+i(%lf) ",VL[i*dim+k],VL[(i+1)*dim+k]);

                        }else{
                            printf("(%lf)-i(%lf) ",VL[i*dim+k],VL[(i+1)*dim+k]);
                        }
                    }
                    printf("  ]\n");
                }
                i++;
            }else{
                /*real eigen vector*/
                printf("[  ");
                for(k=0;k<dim;k++){
                    printf("%lf ",VL[i*dim+k]);
                }
                printf("] \n");
            }
            printf("\n");
        }
    }
}

[TECH] Cauchy's Inequality [continued...]

Just adding the solution to one more problem continuing from last post.

[TECH] Cauchy's inequality a great tool for approximation algorithms

Its been really long since I made a post, all these day's I'm doing a lot of theory and trying to come up with approximation algorithms to the Border Length Minimization Problem. After a long time I did prove the problem is NP-hard by reducing the T.S.P problem to the B.L.M.P problem. I'm now trying to find some approximation algorithms for the problem. I thought if I could strengthen my skills on inequalities it would be really great and I found a great problem book for inequalities by Michael Steele which has a good collections of problems on inequalities. In this post I provide my solutions to some of the problems which seemed interesting.

[TECH] The shopping problem is as hard as TSP.

The Shopping Problem : Let I={i₁,i₂...i_k} be the items on your shopping list. Let S = {s₁,s₂,....,s_m} be the list of stores each store s_i sells only a subset of the items on your shopping list, and also the cost the items sold by the shops may vary i.e item i_j may be sold with different prices in different shops, so let C[i,j] be the cost of item i in store j. Since the shops are located in different places let D[p,q] be the cost of reaching store q from store p. So our aim is to start at home (H) and buy all the items we want and return to home with minimum possible cost. As I mentioned in the previous post this problem is was mentioned in Google Code Jam practice problems, however I did not see many submissions to this problem.

Firstly the problem is NP-hard because if items sold by each of the stores are disjoint and they are singleton sets then we should visit all the stores in the right order to minimize the total distance traveled, this is nothing by TSP. The following is the dynamic programming formulation for this problem. Let S(V,j) be a optimal solution to buy all the items in set V ending at store j. Then the following is the DP formulation.
S(V,j) = min { S(V-{k},p) + D[p,j] + C[k,j] } ; k \in V, 1<=p<=m Final Answer = min { S(I,j) + D[j,Home] } Initialization of the DP is simple you may want to look at the code below.

Further complexity to the problem can be added by saying that the items are classified as PerishableP and Non-PerishableN. So when ever you buy a perishable item you need to go back home to put in the refrigerator. The above DP formulation can be extended to this easily as follows. Let S(V,j,P) be the optimal sub problem of shopping all the items in the set V ending at the store j and having a Perishable item in the shopping cart, on the same lines S(V,j,N) optimal sub problem of shopping all the items in the set V ending at store j with NO perishable item in the shopping cart. Then the following formulation will give the optimal solution.

S(V,j,P) = min {S(V-{k},p,P) + D[p,Home] + D[Home,j] + C[k,j] } or
              min {S(V-{k},p,N) + D[p,j] + C[k,j] } 
If k is a perishable item

We can write similar formulations for Non-Perishable items also, get the code here implementing the shopping problem.

  1 /*Vamsi Kundeti (vamsi(dot)krishnak(at)gmail*/
  2 #include<stdio.h>
  3 #include<stdlib.h>
  4 #include<string.h>
  5 #include<assert.h>
  6 #include<math.h>
  7 #include "List.h"
  8 unsigned int **GRAPH;
  9 //#define VERBOSE_MODE 1
 10 unsigned int item_bits;
 11 char **item_names;
 12 unsigned int item_subset;
 13 /*has a unsigned int indicating
 14 *if this shop sells the item or not
 15 */
 16 unsigned int *shop_sells;
 17 /*dist[i][j] indicates the distance 
 18  *between i and j*/
 19 double **dist;
 20 /*cost[i][j] indicates the cost of
 21 *item i in shop j
 22 */
 23 double **cost;
 24 /*The number of the items*/
 25 unsigned int icount;
 26 /*The number of test cases*/
 27 unsigned int tcount;
 28 unsigned int caseno;
 29 /*The number of shops**/
 30 unsigned int scount;
 31 typedef struct _pair_ {
 32     int x;
 33     int y;
 34 }Coor;
 35 Coor *shop_coor;
 36 void EnumSubSetsBFS(unsigned int n);
 37 unsigned int LookUpItem(const char *key){
 38     unsigned int i=0;
 39     for(i=0;i<icount;i++){
 40         if(!strcmp(key,item_names[i])){
 41             return i+1;
 42         }
 43     }
 44     /*Ignore this item*/
 45     return 0;
 46 }
 47 unsigned int gas_price;
 48 unsigned int perish;
 49 inline char CheckPerish(unsigned int i){
 50     return (perish & (1<<(i)))?1:0;
 51 }
 52 void ReadInput(void){
 53     unsigned int i,j;
 54     unsigned int len;
 55     unsigned int icost;
 56     char item_buffer[64];
 57     char delim; 
 58     scanf("%u",&tcount);
 59     while(tcount--){
 60         perish=0;
 61         scanf("%u %u %u",&icount,&scount,&gas_price);
 62         assert(icount && scount);
 63         /*Read the items to be bought*/
 64         item_names = (char **) malloc(sizeof(char *)*icount);
 65         for(i=0;i<icount;i++){
 66             scanf("%s",&item_buffer);
 67             len = strlen(item_buffer);
 68             len = (len > 64)?64:len;
 69             item_names[i] = (char *) malloc(sizeof(char)*(len+1));
 70             assert(memcpy(item_names[i],item_buffer,(len+1))==item_names[i]);
 71             if(len && (item_names[i])[len-1]=='!'){
 72                 perish = perish | (1 << (i));
 73                 (item_names[i])[len-1]='\0';
 74             }
 75         }
 76         shop_coor = (Coor *) malloc(sizeof(Coor)*(scount+1));
 77         cost = (double **) malloc(sizeof(double *)*(scount+1));
 78         shop_sells = (unsigned int *)malloc(sizeof(unsigned int)*(scount+1));
 79         shop_coor[0].x = 0.0;
 80         shop_coor[0].y = 0.0;
 81         for(i=1;i<scount+1;i++){
 82             cost[i] = (double *) malloc(sizeof(double)*icount);
 83             scanf("%d %d",&(shop_coor[i].x),&(shop_coor[i].y));
 84             shop_sells[i] = 0;
 85             for(j=0;j<icount;j++){
 86                 cost[i][j] = 0;
 87             }
 88             while(1){
 89                 scanf(" %[^:]:%u",item_buffer,&icost);
 90                 if(j=LookUpItem(item_buffer)){
 91                     cost[i][j-1] = icost;
 92                     /*Also set the sells bit_set*/
 93                     shop_sells[i] = shop_sells[i] | (1<<(j-1));
 94                 }
 95                 scanf("%c",&delim);
 96                 if(delim == '\n'){
 97                     break;
 98                 }
 99             }
100         }
101         /*Now Compute the distance matrix*/
102         dist = (double **)malloc(sizeof(double *)*(scount+1));
103         for(i=0;i<scount+1;i++){
104             dist[i] = (double *)malloc(sizeof(double)*(scount+1));
105             for(j=0;j<scount+1;j++){
106                 dist[i][j] = sqrtf((shop_coor[i].x - shop_coor[j].x)*
107                     (shop_coor[i].x - shop_coor[j].x) +(shop_coor[i].y - shop_coor[j].y)*
108                         (shop_coor[i].y - shop_coor[j].y));
109             }
110         }
111 #if 0
112         printf("Items..\n");
113         for(i=0;i<icount;i++){
114             printf("%s ",item_names[i]);
115         }
116         printf("\n");
117         printf("Cost Matrix\n");
118         for(i=0;i<scount+1;i++){
119             printf("(%d,%d)",shop_coor[i].x,shop_coor[i].y);
120             for(j=0;j<icount;j++){
121                 if(cost[i][j]){
122                     printf(" %s[%c]:%lf",item_names[j],((CheckPerish(j))?'P':'U'),cost[i][j]);
123                 }
124             }
125             printf("\n");
126         }
127         printf("\n");
128 #endif
129         /*Make the call to the dynamic programming routine*/
130         EnumSubSetsBFS(icount);
131 
132         /*Freeup the stuff*/
133         for(i=0;i<icount;i++){
134             free(item_names[i]); 
135         }
136         free(item_names);
137         free(shop_coor);
138         free(shop_sells);
139         for(i=1;i<scount+1;i++){
140             free(cost[i]);
141         }
142         free(cost);
143         for(i=0;i<scount+1;i++){
144             free(dist[i]);
145         }
146         free(dist);
147     }
148 }
149 
150 void CreateGraph(unsigned int *n){
151     unsigned int i,j;
152     printf("Enter the Graph Size\n");
153     scanf("%u",n);
154     printf("Enter the weighted graph\n");
155     GRAPH = (unsigned int **) malloc(sizeof(unsigned int *)*(*n));
156     for(i=0;i<*n;i++){
157         GRAPH[i] = (unsigned int *) malloc(sizeof(unsigned int)*(*n));
158         for(j=0;j<*n;j++){
159             scanf("%u",&(GRAPH[i][j]));
160         }
161     }
162 }
163 
164 
165 unsigned int subset_count=0;
166 int main(){
167     /*Create the Graph*/
168     unsigned int n;
169     caseno=0;
170     ReadInput();
171     //EnumSubSetsBFS(n);
172 
173 }
174 /*Returns the new tail of the list*/
175 inline List* ListAppend(List **tail,void *_data){
176     List *new_node = (List *)malloc(sizeof(new_node)*1);
177     assert(new_node);
178     if(!(*tail)){
179         (*tail) = new_node;
180         (*tail)->_data = _data;
181         (*tail)->next = NULL;
182         return *tail;
183     }
184     new_node->_data = _data;
185     new_node->next = NULL;
186     (*tail)->next = new_node;
187     return new_node;
188 }
189 typedef struct _subset_{
190     /*if i^th bit is set then subset has i^th
191      *element*/
192     unsigned int ebits;
193     unsigned int max; /*Maximum element in the subset*/
194     unsigned int size; /*Size of subset*/
195     /*Pointer to the sub-problem S({},i,0) S({},i,1)*/
196     void *sub_problem;
197 }SubSet;
198 
199 /*Enumerate the Subset using BFS, n is 
200  *the set size [1,2...n] */
201 void EnumSubSetsBFS(unsigned int n){
202     unsigned int check_sub_set,k,i,j;
203     unsigned int current_level = 0;
204     List *bfs_list = NULL;
205     List *tail = NULL; List *free_node,*ptr;
206     SubSet *sub_set,*new_sub_set;
207     SubSet *one_sets = NULL;
208 
209     one_sets = (SubSet *)malloc(sizeof(SubSet)*icount); 
210     assert(one_sets);
211     /*Add the 1-subsets i=items*/
212     for(i=0;i<icount;i++){
213         one_sets[i].ebits = 0; one_sets[i].ebits |= (1<<i);
214         one_sets[i].max = i; one_sets[i].size = 1;
215         /****SUBPROBLEM SPECIFIC*****/
216         one_sets[i].sub_problem = (double **)malloc(sizeof(double *)*2);
217         /* Perishable S({i},X,0) */
218         ((double **)one_sets[i].sub_problem)[0] = (double *) malloc(sizeof(double)*(scount+1));
219         /* Unperishable S({i},X,1) */
220         ((double **)one_sets[i].sub_problem)[1] = (double *) malloc(sizeof(double)*(scount+1));
221         double **sub_problem_i = (double **)one_sets[i].sub_problem;
222 
223         /*j is the shop*/
224         for(j=1;j<scount+1;j++){
225             /*if item i can be bought at shop j
226              *depending on perishable or not initialize sub_problem[0] 
227              *sub_problem[1]
228              */
229              sub_problem_i[0][j] = HUGE_VAL; sub_problem_i[1][j] = HUGE_VAL;
230              if(shop_sells[j]&(1<<i)){
231                  if(CheckPerish(i)){
232                      sub_problem_i[0][j] = dist[0][j]*(gas_price) + cost[j][i];
233                  }else{
234                      sub_problem_i[1][j] = dist[0][j]*(gas_price) + cost[j][i];
235                  }
236              }
237 #ifdef VERBOSE_MODE
238             printf("S({%u},%u,0) = %lf\n",i+1,j+1,sub_problem_i[0][j]);
239             printf("S({%u},%u,1) = %lf\n",i+1,j+1,sub_problem_i[1][j]);
240 #endif
241         }
242         /*****************************/
243         if(!tail){
244             ListAppend(&bfs_list,(void *)&(one_sets[i]));
245             tail = bfs_list;
246             continue;
247         }
248         tail = ListAppend(&tail,(void *)&(one_sets[i]));
249     }
250 
251 
252 
253     /*Now do the BFS*/
254     while(bfs_list){
255         sub_set = (SubSet *)bfs_list->_data; 
256         free_node = bfs_list;
257         bfs_list=bfs_list->next;
258         subset_count++;
259         assert(sub_set->max < icount);
260         for(i=((sub_set->max)+1);i<icount;i++){
261             new_sub_set = (SubSet *) malloc(sizeof(SubSet)*1);
262             new_sub_set->size = (sub_set->size)+1;
263             new_sub_set->max = i;
264             new_sub_set->ebits = (sub_set->ebits | (1<<i));
265             /*This subset enumeration is critical to Subset based D.P*/
266 
267             /*****<begin>**SUBPROBLEM SPECIFIC******/
268             /* S(X,j) = min_k {S(X-k,k) + d[k][j]}*/
269             ptr = free_node;
270             double current_min = HUGE_VAL;
271             double **new_sub_problem;
272             double **sub_problem; double local_min;
273             double local_min_perish,local_min_unperish;
274             unsigned int k1,perish,shop1;
275             List *ptr_prev;
276             
277             new_sub_set->sub_problem = (double **) malloc(sizeof(double *)*2);
278             ((double **)new_sub_set->sub_problem)[0] = (double *) malloc(sizeof(double)*(scount+1));
279             ((double **)new_sub_set->sub_problem)[1] = (double *) malloc(sizeof(double)*(scount+1));
280             new_sub_problem = (double **)new_sub_set->sub_problem;
281 
282             for(perish=0;perish<=1;perish++){
283                 for(j=1;j<scount+1;j++){
284                     /*Objective here is to fill up S({new_sub_set},j,0) 
285                      *and S({new_sub_set},j,1)*/
286                     ptr = free_node;
287                     current_min = HUGE_VAL;
288                     for(k=1;k<=new_sub_set->size;k++){
289                         do{
290                             check_sub_set = ((SubSet *)ptr->_data)->ebits; 
291                             check_sub_set ^= new_sub_set->ebits; 
292                             k1 = (int)log2(check_sub_set);
293                             ptr_prev = ptr;
294                             ptr=ptr->next;
295                         }while(check_sub_set & (check_sub_set-1));
296 
297                         /*j is the shop and k1 is the item which you want
298                          *to buy at this shop, depending on
299                          *perishablity create the value of current min.
300                         */
301                         if(!(shop_sells[j] & check_sub_set)){
302                              continue; /*The shop don't sell this*/
303                         }
304                         sub_problem = (double **)
305                             ((SubSet *)ptr_prev->_data)->sub_problem;
306                         /*item k1 is sold by shop j*/
307                         local_min = HUGE_VAL;
308                         if((perish && !CheckPerish(k1)) 
309                                 || (!perish && CheckPerish(k1))){
310                             for(shop1=1;shop1<scount+1;shop1++){
311                                 /*S{{new_sub_set}-{k1},shop1} + 
312                                 d(shop1,h) + d(h,j) + c[k1][j]*/
313                                 if(j != shop1 || perish){
314                                     local_min_perish = 
315                                         sub_problem[0][shop1] + 
316                                         (dist[shop1][0] + dist[0][j])*
317                                             (gas_price) + cost[j][k1];
318                                 }else if(!perish){
319                                     local_min_perish = sub_problem[0][shop1] + 
320                                         cost[j][k1];
321                                 }
322                                 local_min_unperish = sub_problem[1][shop1] + 
323                                     (dist[shop1][j])*(gas_price) + cost[j][k1];
324                                 local_min_perish = 
325                                  (local_min_perish < local_min_unperish)? 
326                                     local_min_perish:local_min_unperish;
327 
328                                 if(local_min_perish < local_min){
329                                     local_min = local_min_perish;
330                                 }
331                             }
332                         }
333                         if(local_min < current_min){
334                             current_min = local_min; 
335                         }
336                     }
337                     new_sub_problem[perish][j] = current_min;
338                 }
339             }
340 
341 #ifdef VERBOSE_MODE
342             printf("Printing S(V,i) size=%u \n",sub_set->size);
343             /*Print SubSet*/
344             unsigned int i1;
345             for(i1=1;i1<scount+1;i1++){
346                 printf("S({");
347                 unsigned int j1;
348                 /*Get the elements back*/
349                 check_sub_set = new_sub_set->ebits;
350                 for(j1=0;j1<icount;j1++){
351                     if(check_sub_set & (1<<j1)){
352                         printf(" %u ",j1+1);
353                     }
354                 }
355                 printf("},%u,0)=%lf\n",i1+1,((double **)(new_sub_set->sub_problem))[0][i1]);
356                 printf("S({ABOVE},%u,1)=%lf\n",i1+1,((double **)(new_sub_set->sub_problem))[1][i1]);
357             }
358 #endif
359 
360             /**********<END>****************/
361             tail = ListAppend(&tail,(void *)new_sub_set);
362         }
363 
364         if(sub_set->size == icount){
365             double final_answer=HUGE_VAL;
366             double local_min; double **sub_problem;
367             sub_problem = (double **) sub_set->sub_problem;
368             for(j=1;j<scount+1;j++){
369                 local_min = (sub_problem[0][j]<sub_problem[1][j])?
370                     sub_problem[0][j]:sub_problem[1][j];
371                 local_min += (dist[j][0])*(gas_price);
372                 if(local_min < final_answer){
373                     final_answer = local_min;
374                 }
375             }
376             printf("Case #%u: %0.7lf\n",++caseno,final_answer);
377         }
378 
379         if(!(sub_set >= one_sets && sub_set <= &(one_sets[icount-1]))){
380             free(((double **)sub_set->sub_problem)[0]);
381             free(((double **)sub_set->sub_problem)[1]);
382             free(sub_set->sub_problem);
383             free(sub_set); 
384         }
385         free(free_node);
386     }
387     for(i=0;i<icount;i++){
388         free(((double **)one_sets[i].sub_problem)[0]);
389         free(((double **)one_sets[i].sub_problem)[1]);
390         free(one_sets[i].sub_problem);
391     }
392     free(one_sets);
393 //    printf("SubsetCount=%u\n",subset_count);
394 }
395 

[TECH] A O(log(n)) time and O(1) space algorithm to find out Nth Fibonacci number.

We often come across Fibonacci numbers in several contexts, if you see Knuth's "Concrete Mathematics" in the chapter on "Generating functions" one interesting example illustrates the possible brick patterns to build a wall whose width is 'n' using bricks of dimension 2x1. Also see the UVA problem here

We know a simple dynamic programming based algorithm can lead us to a O(n) time algorithm with O(1). But we can do much better than that in fact we can do it in O(log(n)) time and O(1) space. We can rewrite the dynamic programming as matrix formulation as follows.

Fn Fn-1

=

1 1 1 0

Fn-1 Fn-2

This is really great we can get the closed form for F_n by repeated substitution. In fact the following is the closed form for F_n is as follows.

Fn Fn-1

=

1 1 1 0

n-1

1 0

So the problem of finding the n^th Fibonacci number is reduces to evaluating X^n-1 where X is a 2x2 matrix as shown above. Now the question is how we can evaluate X^n-1 as efficiently as possible. In fact to compute any a^b we need exactly |_{_}log(n)_{_}| multiplications by making use of the bit representation of b = (2^i₁ + 2^i₂ ...) where i_j are the positions of bits which have 1 in the binary representation of b. I used this trick to get a O(log(n)) time and O(1) time algorithm to compute the n^th Fibonacci number. Please see the code below get this here

  1 /*Find the n^th fibonacci number in O(log(n)) 
  2  *time, using the following idea of matrix formulation 
  3  * 
  4  * [   F_n  ]    [ 1      1  ] [F_{n-1}]
  5  * [        ] = [         ] [       ]
  6  * [ F_{n-1}]    [ 1      0  ] [F_{n-2}]
  7  *
  8  * vamsik(at)engr(dot)uconn July 11, 2008
  9  */
 10 #include<stdio.h>
 11 #include<stdlib.h>
 12 /*Fib returns the n^th Fibonacci number*/
 13 unsigned int Fib(unsigned int n){
 14     unsigned int deg=1;
 15     unsigned int fib[2][2] = {{1,1},{1,0}};
 16     unsigned int fibAns[2][2];
 17     unsigned char initAns=0;
 18     unsigned int mul[4];
 19     if(n<=2){ return (n<2)?n:2;}
 20     while(deg < n-1){
 21          mul[0] = fib[0][0]*fib[0][0] + fib[0][1]*fib[1][0];
 22          mul[1] = fib[0][0]*fib[0][1] + fib[0][1]*fib[1][1];
 23          mul[2] = fib[1][0]*fib[0][0] + fib[1][1]*fib[1][0];
 24          mul[3] = fib[1][0]*fib[0][1] + fib[1][1]*fib[1][1];
 25          fib[0][0] = mul[0]; fib[0][1] = mul[1]; 
 26          fib[1][0] = mul[2]; fib[1][1] = mul[3];
 27          deg <<=1;
 28          if(deg & n-1){ /*i^th bit set*/
 29              if(!initAns){
 30                  fibAns[0][0] = fib[0][0]; fibAns[0][1] = fib[0][1]; 
 31                  fibAns[1][0] = fib[1][0]; fibAns[1][1] = fib[1][1];
 32                 initAns=1;
 33              }else{
 34                  mul[0] = fibAns[0][0]*fib[0][0] + fibAns[0][1]*fib[1][0];
 35                  mul[1] = fibAns[0][0]*fib[0][1] + fibAns[0][1]*fib[1][1];
 36                  mul[2] = fibAns[1][0]*fib[0][0] + fibAns[1][1]*fib[1][0];
 37                  mul[3] = fibAns[1][0]*fib[0][1] + fibAns[1][1]*fib[1][1];
 38                  fibAns[0][0] = mul[0]; fibAns[0][1] = mul[1]; 
 39                  fibAns[1][0] = mul[2]; fibAns[1][1] = mul[3];
 40              }
 41          }
 42     }
 43     if(n-1 & 1){
 44         return fibAns[0][0] + fibAns[0][1];
 45     }
 46     return fibAns[0][0];
 47 }
 48 int main(int argc,char **argv){
 49     unsigned int n;
 50     while(1){
 51         scanf("%u",&n);
 52         if(!n){
 53             exit(0);
 54         }
 55         printf("%u\n",n,Fib(n));
 56     }
 57 }
 58 

[TECH] Illustration of how to use the BFS based subset enumeration algorithmic framework to solve TSP.

In the previous post I posted on the idea of using BFS to enumerate the subsets incrementally based on the size (see the previous post). The power of this core framework is so immense let me show how we can implement the TSP based on the framework. The extra code is between the lines 82-90 (TSP Initialization) and 113-175. The same core can be easily adapted to solve several sub-set based Dynamic programming optimizations like the job-scheduling, vehicle routing problem (VRP). If you have not looked into one of the Google Code Jam 2008 (practice) problems the problem#4 which is a shopping problem can be easily solved by this core framework. Check this code TSP.C (use -std=c99 to compile).

  1 /*******************************
  2  * Subset Enumeration Algorithmic 
  3  * Framework and its illustration in 
  4  * finding TSP. 
  5  ******************************/
  6 #include<stdio.h>
  7 #include<stdlib.h>
  8 #include<string.h>
  9 #include<assert.h>
 10 #include<math.h>
 11 #include "List.h"
 12 unsigned int **GRAPH;
 13 unsigned int subset_count=0;
 14 //#define VERBOSE_MODE 1
 15 void EnumSubSetsBFS(unsigned int n);
 16 
 17 void CreateGraph(unsigned int *n){
 18     unsigned int i,j;
 19     printf("Enter the Graph Size\n");
 20     scanf("%u",n);
 21     if(*n>32){
 22         fprintf(stderr,"TODO: The Bitset currently supports only 32 bits\n");
 23         exit(1);
 24     }
 25     printf("Enter the weighted graph\n");
 26     GRAPH = (unsigned int **) malloc(sizeof(unsigned int *)*(*n));
 27     for(i=0;i<*n;i++){
 28         GRAPH[i] = (unsigned int *) malloc(sizeof(unsigned int)*(*n));
 29         for(j=0;j<*n;j++){
 30             scanf("%u",&(GRAPH[i][j]));
 31         }
 32     }
 33 }
 34 
 35 int main(){
 36     /*Create the Graph*/
 37     unsigned int n;
 38     CreateGraph(&n);
 39     EnumSubSetsBFS(n);
 40 
 41 }
 42 /*Returns the new tail of the list*/
 43 inline List* ListAppend(List **tail,void *_data){
 44     List *new_node = (List *)malloc(sizeof(new_node)*1);
 45     assert(new_node);
 46     if(!(*tail)){
 47         (*tail) = new_node;
 48         (*tail)->_data = _data;
 49         (*tail)->next = NULL;
 50         return *tail;
 51     }
 52     new_node->_data = _data;
 53     new_node->next = NULL;
 54     (*tail)->next = new_node;
 55     return new_node;
 56 }
 57 typedef struct _subset_{
 58     /*if i^th bit is set then subset has i^th
 59      *element*/
 60     unsigned int ebits;
 61     unsigned int max; /*Maximum element in the subset*/
 62     unsigned int size; /*Size of subset*/
 63     /*Pointer to the sub-problem */
 64     void *sub_problem;
 65 }SubSet;
 66 /*Enumerate the Subset using BFS, n is 
 67  *the set size [1,2...n] */
 68 void EnumSubSetsBFS(unsigned int n){
 69     unsigned int check_sub_set,k,i,j;
 70     unsigned int current_level = 0;
 71     List *bfs_list = NULL;
 72     List *tail = NULL; List *free_node,*ptr;
 73     SubSet *sub_set,*new_sub_set;
 74     SubSet *one_sets = NULL;
 75 
 76     one_sets = (SubSet *)malloc(sizeof(SubSet)*n); 
 77     assert(one_sets);
 78     /*Add the 1-subsets*/
 79     for(i=1;i<n;i++){
 80         one_sets[i].ebits = 0; one_sets[i].ebits |= (1<<i);
 81         one_sets[i].max = i; one_sets[i].size = 1;
 82         /****SUBPROBLEM SPECIFIC*****/
 83         one_sets[i].sub_problem = (float *)malloc(sizeof(float)*n);
 84         for(j=1;j<n;j++){
 85             if(i==j){ continue; }
 86             /* TSP({i},j) */
 87             ((float *)one_sets[i].sub_problem)[j] = GRAPH[0][i]+GRAPH[i][j];
 88             printf("S({%u},%u) = %f\n",i+1,j+1,((float *)one_sets[i].sub_problem)[j]);
 89         }
 90         /*****************************/
 91         if(!tail){
 92             ListAppend(&bfs_list,(void *)&(one_sets[i]));
 93             tail = bfs_list;
 94             continue;
 95         }
 96         tail = ListAppend(&tail,(void *)&(one_sets[i]));
 97     }
 98     /*Now do the BFS*/
 99     while(bfs_list){
100         sub_set = (SubSet *)bfs_list->_data; 
101         free_node = bfs_list;
102         bfs_list=bfs_list->next;
103         subset_count++;
104         assert(sub_set->max < n);
105         for(i=((sub_set->max)+1);i<n;i++){
106             new_sub_set = (SubSet *) malloc(sizeof(SubSet)*1);
107             new_sub_set->size = (sub_set->size)+1;
108             new_sub_set->max = i;
109             new_sub_set->ebits = (sub_set->ebits | (1<<i));
110             
111             
112             
113             /*****<begin>**SUBPROBLEM SPECIFIC******/
114             /* S(X,j) = min_k {S(X-k,k) + d[k][j]}*/
115             ptr = free_node;
116             float current_min = HUGE_VAL;
117             float *new_sub_problem;
118             float *sub_problem; float local_min;
119             unsigned int k1,jstart,jend;
120             List *ptr_prev;
121             
122             new_sub_set->sub_problem = (float *) malloc(sizeof(float)*n);
123             new_sub_problem = (float *)new_sub_set->sub_problem;
124 
125             if(new_sub_set->size == n-1){
126                 jstart = 0; jend = 1;
127             }else{
128                 jstart = 1; jend = n;
129             }
130 
131             for(j=jstart;j<jend;j++){
132                 ptr = free_node;
133                 current_min = HUGE_VAL;
134                 if(new_sub_set->ebits & (1<<j)){
135                     /*Avoid computation of redundant subproblems like
136                      *S({....,x_i,x_j,x_k,....},x_i) */
137                     new_sub_problem[j]=HUGE_VAL;
138                     continue;
139                 }
140 
141                 for(k=1;k<=new_sub_set->size;k++){
142                     do{
143                         check_sub_set = ((SubSet *)ptr->_data)->ebits; 
144                         check_sub_set ^= new_sub_set->ebits; 
145                         k1 = (int)log2(check_sub_set);
146                         ptr_prev = ptr;
147                         ptr=ptr->next;
148                     }while(check_sub_set & (check_sub_set-1));
149                     sub_problem = (float *)((SubSet *)ptr_prev->_data)->sub_problem;
150                     local_min = sub_problem[k1] + GRAPH[k1][j];
151 
152                     if(local_min < current_min){
153                         current_min = local_min; 
154                         new_sub_problem[j] = current_min;
155                     }
156 
157                 }
158             }
159             /*Print Solution to the Subproblem */
160             unsigned int i1,istart,iend;
161             unsigned int j1;
162             
163             for(i1=jstart;i1<jend;i1++){
164                 check_sub_set = new_sub_set->ebits;
165                 if(check_sub_set & (1<<i1)){ continue; }
166                 printf("S(%u,{ ",i1+1);
167                 /*Get the elements back*/
168                 for(j1=0;j1<n;j1++){
169                     if(check_sub_set & (1<<j1)){
170                         printf(" %u ",j1+1);
171                     }
172                 }
173                 printf("})=%f\n",((float *)(new_sub_set->sub_problem))[i1]);
174             }
175             /**********<END>***SUBPROBLEM SPECIFIC*************/
176 
177 
178             tail = ListAppend(&tail,(void *)new_sub_set);
179         }
180         if(!(sub_set >= one_sets && sub_set <= &(one_sets[n-1]))){
181             free(sub_set->sub_problem);
182             free(sub_set); 
183         }
184         free(free_node);
185     }
186     for(i=0;i<n;i++){
187         free(one_sets[i].sub_problem);
188     }
189     free(one_sets);
190     printf("SubsetCount=%u\n",subset_count);
191 }
192 

Sample Input.
---------test5.txt----------
5
0 3 1 5 4 
1 0 5 4 3 
5 4 0 2 1 
3 1 3 0 3 
5 2 4 1 0
---------------------
[vamsik@abadon ~]$ ./a.out < test5.txt 
Enter the Graph Size
Enter the weighted graph
S({2},3) = 8.000000
S({2},4) = 7.000000
S({2},5) = 6.000000
S({3},2) = 5.000000
S({3},4) = 3.000000
S({3},5) = 2.000000
S({4},2) = 6.000000
S({4},3) = 8.000000
S({4},5) = 8.000000
S({5},2) = 6.000000
S({5},3) = 8.000000
S({5},4) = 5.000000
S(4,{  2  3 })=9.000000
S(5,{  2  3 })=8.000000
S(3,{  2  4 })=10.000000
S(5,{  2  4 })=9.000000
S(3,{  2  5 })=10.000000
S(4,{  2  5 })=7.000000
S(2,{  3  4 })=4.000000
S(5,{  3  4 })=6.000000
S(2,{  3  5 })=4.000000
S(4,{  3  5 })=3.000000
S(2,{  4  5 })=6.000000
S(3,{  4  5 })=8.000000
S(5,{  2  3  4 })=7.000000
S(4,{  2  3  5 })=8.000000
S(3,{  2  4  5 })=10.000000
S(2,{  3  4  5 })=4.000000
S(1,{  2  3  4  5 })=5.000000
SubsetCount=15