/*
 *  plainGARSP.c
 *  
 *  A clean implementation of the GARSP algorithm.
 *
 *  Created by Michael Feiri on Wed Jan 02 2002.
 *  Copyright (c) 2002 Michael Feiri. All rights reserved.
 *
 */

#include <stdio.h>

#define NUMWORDS 8	/* 8 words (8*32=256) -> we can create rulers up to a length of 255 */
                        /* but thanks to "principle #1" we can safely get rulers up to (256*2)-1=511 */
#define MAXMARKS 19	/* depth of stacks -> we can place up to 18 marks (plus the "0" mark means OGR-19) */

unsigned int plainGARSP (int maxmarks, int maxlength); /* prototype */

unsigned int plainGARSP (int maxmarks, int maxlength) {

    int depth;                    /* current mark we are placing (corresponds to the D in the PDF) */
    unsigned int count[MAXMARKS]; /* current ruler length (e.g. LENGTH[])*/
    unsigned int list[MAXMARKS][NUMWORDS];   /* list bitmap */
    unsigned int dist[MAXMARKS][NUMWORDS];   /* dist bitmap */
    unsigned int comp[MAXMARKS][NUMWORDS];   /* dist bitmap */
    unsigned int nodes=0;

    unsigned long bit[200];                     /* which bit of LIST to update */
    unsigned char first[65536];   /* first blank in 16 bit bitmap, range: 1..16 */
    int i,j,k,m,s;                     /* counters */
   
    /* Start at the very beginning (one mark set at position zero) */
    depth = 1;
    count[0] = 0;
    count[1] = 0;

    /* Lookup table for the injection of a new marker bit */
    for( i=1; i<200; i++) {
        bit[i] = 0x80000000 >> ((i-1)%32);
    }

    /* Lookup table to find first zero bit in 16 bits */
    k = 0; m = 0x8000;
    for( i=1; i <= 16; i++) {
        for (j = k; j < k+m; j++) first[j] = i;
        k += m;
        m >>= 1;
    }
    first[0xffff] = 17;     /* just in case we use it */

    list[depth][0]=0;
    list[depth][1]=0;
    list[depth][2]=0;
    list[depth][3]=0;
    list[depth][4]=0;
    list[depth][5]=0;
    list[depth][6]=0;
    list[depth][7]=0;

    dist[depth][0]=0;
    dist[depth][1]=0;
    dist[depth][2]=0;
    dist[depth][3]=0;
    dist[depth][4]=0;
    dist[depth][5]=0;
    dist[depth][6]=0;
    dist[depth][7]=0;

    comp[depth][0]=0;
    comp[depth][1]=0;
    comp[depth][2]=0;
    comp[depth][3]=0;
    comp[depth][4]=0;
    comp[depth][5]=0;
    comp[depth][6]=0;
    comp[depth][7]=0;

    /* search through the search space **** */
    while (1) {

    nodes++;

   /* Find the next available mark location for this level */

    if( comp[depth][0] < 0xffff0000 ) {
        s = first[comp[depth][0]>>16];
    } else if( comp[depth][0] < 0xffffffff ) {
        s = 16 + first[comp[depth][0] & 0x0000ffff];
    
    } else if( comp[depth][1] < 0xffff0000 ) {
        s = 32 + first[comp[depth][1]>>16];
    } else if( comp[depth][1] < 0xffffffff ) {
        s = 48 + first[comp[depth][1] & 0x0000ffff];

    } else if( comp[depth][2] < 0xffff0000 ) {
        s = 64 + first[comp[depth][2]>>16];
    } else if( comp[depth][2] < 0xffffffff ) {
        s = 80 + first[comp[depth][2] & 0x0000ffff];

    } else if( comp[depth][3] < 0xffff0000 ) {
        s = 96 + first[comp[depth][3]>>16];
    } else if( comp[depth][3] < 0xffffffff ) {
        s = 112 + first[comp[depth][3] & 0x0000ffff];

    } else if( comp[depth][4] < 0xffff0000 ) {
        s = 128 + first[comp[depth][4]>>16];
    } else if( comp[depth][4] < 0xffffffff ) {
        s = 144 + first[comp[depth][4] & 0x0000ffff];

    } else if( comp[depth][5] < 0xffff0000 ) {
        s = 160 + first[comp[depth][5]>>16];
    } else if( comp[depth][5] < 0xffffffff ) {
        s = 176 + first[comp[depth][5] & 0x0000ffff];

    } else if( comp[depth][6] < 0xffff0000 ) {
        s = 192 + first[comp[depth][6]>>16];
    } else if( comp[depth][6] < 0xffffffff ) {
        s = 208 + first[comp[depth][6] & 0x0000ffff];

    } else if( comp[depth][7] < 0xffff0000 ) {
        s = 224 + first[comp[depth][7]>>16];
    } else /*if( comp[depth][7] < 0xffffffff )*/ {
        s = 240 + first[comp[depth][7] & 0x0000ffff];
    }

        /*if the length of our ruler does not exceeds the maximum length we can proceed */
        if ((count[depth] += s) <= maxlength) {
        
            /* shift by s (LIST to the right, COMP to the left) */
            while (s>=32) {
                comp[depth][0] = comp[depth][1];
                comp[depth][1] = comp[depth][2];
                comp[depth][2] = comp[depth][3];
                comp[depth][3] = comp[depth][4];
                comp[depth][4] = comp[depth][5];
                comp[depth][5] = comp[depth][6];
                comp[depth][6] = comp[depth][7];
                comp[depth][7] = 0;  
                list[depth][7] = list[depth][6];   
                list[depth][6] = list[depth][5];   
                list[depth][5] = list[depth][4];   
                list[depth][4] = list[depth][3];   
                list[depth][3] = list[depth][2];   
                list[depth][2] = list[depth][1];   
                list[depth][1] = list[depth][0];   
                list[depth][0] = 0;
                s -= 32;
            }
            if (s!=0) {
                int ss=32-s;
                comp[depth][0] = (comp[depth][0] << s) | (comp[depth][1] >> ss);   
                comp[depth][1] = (comp[depth][1] << s) | (comp[depth][2] >> ss);   
                comp[depth][2] = (comp[depth][2] << s) | (comp[depth][3] >> ss);   
                comp[depth][3] = (comp[depth][3] << s) | (comp[depth][4] >> ss);  
                comp[depth][4] = (comp[depth][4] << s) | (comp[depth][5] >> ss);   
                comp[depth][5] = (comp[depth][5] << s) | (comp[depth][6] >> ss);   
                comp[depth][6] = (comp[depth][6] << s) | (comp[depth][7] >> ss);
                comp[depth][7] <<= s;
                list[depth][7] = (list[depth][7] >> s) | (list[depth][6] << ss);
                list[depth][6] = (list[depth][6] >> s) | (list[depth][5] << ss);
                list[depth][5] = (list[depth][5] >> s) | (list[depth][4] << ss);
                list[depth][4] = (list[depth][4] >> s) | (list[depth][3] << ss);
                list[depth][3] = (list[depth][3] >> s) | (list[depth][2] << ss);
                list[depth][2] = (list[depth][2] >> s) | (list[depth][1] << ss);
                list[depth][1] = (list[depth][1] >> s) | (list[depth][0] << ss);
                list[depth][0] >>= s;
            }
            
            /* did we find a ruler with the desired number of marks golomb? */
            if (depth == (maxmarks)) {
                int i;
                printf("FOUND A RULER!!!\n");
                for(i=0; i<=maxmarks; i++ ) printf("%4d",count[i]);
                printf("\n");
            } else { /* if we didnt distribute all marks we need to go deeper... */
                
                /* prepare LIST[D+1]=LIST+(new_mark) */
                list[depth+1][0] = list[depth][0];
                list[depth+1][1] = list[depth][1];
                list[depth+1][2] = list[depth][2];
                list[depth+1][3] = list[depth][3];
                list[depth+1][4] = list[depth][4];
                list[depth+1][5] = list[depth][5];
                list[depth+1][6] = list[depth][6];
                list[depth+1][7] = list[depth][7];
                {
                int d = count[depth]-count[depth-1];
                if( d <= 32 ) { 
                    list[depth+1][0] = list[depth][0] + bit[ d ];
                } else if( d <= 64 ) { 
                    list[depth+1][1] = list[depth][1] + bit[ d ];
                } else if( d <= 96 ) { 
                    list[depth+1][2] = list[depth][2] + bit[ d ];
                } else if( d <= 128 ) { 
                    list[depth+1][3] = list[depth][3] + bit[ d ];
                } else if( d <= 160 ) { 
                    list[depth+1][4] = list[depth][4] + bit[ d ];
                } else if( d <= 192 ) { 
                    list[depth+1][5] = list[depth][5] + bit[ d ];
                } else if( d <= 224 ) { 
                    list[depth+1][6] = list[depth][6] + bit[ d ];
                } else /*if( d <= 256 )*/ { 
                    list[depth+1][7] = list[depth][7] + bit[ d ];
                }
                }
                
                /* prepare the new DIST[D+1]=DIST|LIST[D+1] */
                dist[depth+1][0] = list[depth+1][0] | dist[depth][0];
                dist[depth+1][1] = list[depth+1][1] | dist[depth][1];
                dist[depth+1][2] = list[depth+1][2] | dist[depth][2];
                dist[depth+1][3] = list[depth+1][3] | dist[depth][3];
                dist[depth+1][4] = list[depth+1][4] | dist[depth][4];
                dist[depth+1][5] = list[depth+1][5] | dist[depth][5];
                dist[depth+1][6] = list[depth+1][6] | dist[depth][6];
                dist[depth+1][7] = list[depth+1][7] | dist[depth][7];
            
                /* prepare COMP[D+1]=COMP|DIST[D+1] */
                comp[depth+1][0] = dist[depth+1][0] | comp[depth][0];
                comp[depth+1][1] = dist[depth+1][1] | comp[depth][1];
                comp[depth+1][2] = dist[depth+1][2] | comp[depth][2];
                comp[depth+1][3] = dist[depth+1][3] | comp[depth][3];
                comp[depth+1][4] = dist[depth+1][4] | comp[depth][4];
                comp[depth+1][5] = dist[depth+1][5] | comp[depth][5];
                comp[depth+1][6] = dist[depth+1][6] | comp[depth][6];
                comp[depth+1][7] = dist[depth+1][7] | comp[depth][7];
                
                count[depth+1] = count[depth]; /* pass on our current length */
                depth++;		       /* one level deeper */

                continue; /* go back to the start of our while(1) loop */
            }
        }
        /*if the length of our ruler exceeds the maximum length we need to go up one depth */
        else {
            if (--depth == 0) return nodes;/*break; // EXIT if we reached end of search space */
        }
    }

    return nodes;
}

#ifndef AS_MODULE

#include <sys/time.h>
#include <sys/resource.h>

double dtime()
{
struct rusage rusage;
double q;

 getrusage(RUSAGE_SELF,&rusage);

 q = (double)(rusage.ru_utime.tv_sec);
 q = q + (double)(rusage.ru_utime.tv_usec) * 1.0e-06;
	
 return q;
}

int main (int argc, const char * argv[]) {

    double starttime, benchtime;
    unsigned int nodes;
    
    printf("\nplainGARSP started...\n");
    
    starttime = dtime();
    nodes = plainGARSP(10,72);
    benchtime = dtime() - starttime;
    
    printf("Nodes per second = %.3f ( %u nodes/ %.3f seconds)\n",nodes/benchtime,nodes,benchtime);    
 
    return 0;
}

#endif