/* The following program calculates a minimum genarator and a maximum independent
   subset for a family of subpaths of a main Circuit. The notation "arc" in the program 
   and the comments is equivalent to that of "directed edge" in my thesis. */

 /* This program was implemented in 1997/98 by Joerg Foerster at the Eotvos University in 
   Budapest. */
   
/* e_mail: joergfoerster@hotmail.com  or  joergf@cs.elte.hu */



#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <unistd.h>


/* member of the dynamic field of directed paths */
typedef struct Path
    {   int                  left_node;  /* number of the left endnode */
        int                  right_node; /* number of the right endnode */
        int                  id;         /* -1 if freebie gen., >=0 otherwise */               
    } path;

typedef struct Double_Int
    {   int                  first;
        int                  last;
    } double_int;


typedef struct Chain_Tree
     {
       int                   id;
       int                   left_node;  /* number of the left endnode on the mainpath */
       int                   right_node; /* number of the right endnode -"- */
       struct Chain_Tree     *l_branch;
       struct Chain_Tree     *r_branch;
     } chain_tree;
      

/* member of the list E_R of essential pairs */
typedef struct E_pair
    {   int                  id;            /* index in the list of selected ess. pairs, -1 otherwise */
        int                  path_number;   /* path-index in the sorted p_field */
	int                  arc_number;    /* index of the ess. arc's left endnode */
        struct E_pair        *next;         /* pointer to the next member in the list */
	struct E_pair        *prev;         /* pointer to the previous - " - */
        struct E_pair        *next_on_path; /* pointer to the next essential member 
					       on the SAME path */
        struct E_pair        *prev_on_path; /* pointer to the previous - " - */
        struct Chain_Tree    *cover;
     } e_pair;


/* member of the list of forbidden ess. repr. subpaths */
typedef struct Forbidden
    {
        struct E_pair        *forb;   /* pointer to the corresponding E_pair struct */
        struct Forbidden     *next;   /* pointer to the next member in the list */
        struct Forbidden     *prev;   /* pointer to th previous - " - */
    } forbidden;

/* member of the bipartite graph's edge-list */
typedef struct Edge
    { 
      int          adjac;     /* index of the adjacent node */
      struct Edge  *next;     /* pointer to the next edge in the list */
    } edge;

/* structure for the nodes of the bip. graph */
typedef struct Node1
   {  int          match;   /* index of the matching-pair if matched, 0 othewise */ 
      int          root;    /* root of the alternating path */
      int          parent;  /* parent in the alt. path from the root */
      edge         *elist;  /* edge-list */
      edge         *last;   /* pointer to the last member of the edge-list */
    } node1;


/* structure for the nodes of the bip. graph */
typedef struct Node2
   {  int          match;   /* 1 if matched, 0 othewise */ 
      int          root;    /* root of the alternating path */
      int          parent;  /* parent in the alt. path */
   } node2;



/* the bipartite graph */
typedef struct Bip_Graph
    { 
      int          n_nodes;  /* number of nodes in the bip. graph */     
      int          n_edges;  /* number of edges in the bip. graph */
      node1        *nodes1;  /* dynamic field of nodes */
      node2        *nodes2;
    } bip_graph;                    

/* memeber of the dynamic field of chains constructed from the max. matching */
typedef struct Chain
     {
       int                   left_node;  /* number of the left endnode on the mainpath */
       int                   right_node; /* number of the right endnode -"- */
       int                   e_left;     /* index of the smallest member in E_R */
       int                   e_right;    /*      -"-     largest     -"-  */
     } chain;



/* helpful struct for the construction of the set of essential
   pairs (called E_R) in which we store the essential 'subintervals' of an interval */
typedef struct Interval
     {         
       int                   left_node;
       int                   right_node; 
       struct Interval       *next;
     } interval;


/* structure for the construction of E_R where we store the instantaneous
last essential arc and next poss arc for a certain path of the path-field */
typedef struct Constr_Essesntial
    {
       int                   nextposs; /* the next possible essential arc on a path */
       e_pair                *lastess; /* the last essential arc found on a path */
    } constr_essential;



/* bipartite matching algorithm in file bipmatch.c */
extern int *Bip_Match();
/* 'get resource usage' for runtime tests */
extern int getrusage ();


/***********************************************************************************/

/*int Sign(i)
  int i;

  {
  if(i>=0)
  return (0);  
  else
  return (1);
    
  }*/

/***********************************************************************************/

#define Modulo(a,b) ((a)%(b))
#define Sign(i) ((i)>=0 ? 0 : 1) 

/* int Modulo(int a, int b)
   {
   int c;
   
   c = a - (a / b) * b;
   
   return(c);
   } */
   
/***********************************************************************************/

/* gets two intersecting paths (resp. their endnodes) and returns 1 if path two is a 
   true subpath of path one and 0 otherwise. */
int Two_in_One(p_field, p1, p2)
path *p_field;
int p1;
int p2;


{
 int help, l1, r1, l2, r2;
 
 l1 = p_field[p1].left_node;
 r1 = p_field[p1].right_node;
 l2 = p_field[p2].left_node;
 r2 = p_field[p2].right_node;

 help = Sign(l2-l1) + Sign(r2-l2) + Sign(r1-r2) + Sign(l1-r1);
 
 if (help == 3)
   return(1);
 else
   return(0);
}

/***********************************************************************************/

/* Gets the numbers of three nodes a,b and c on the main circuit and decides 
   whether b is between a and c or equal to one of them according to the main 
   direction. The function returns 1 if the case is like the above mentioned
   and 0 otherwise. It also returns 1 if a = b = c.*/

int ABC ( a, b, c)
int a;
int b;
int c;

{
 int help;

 help = Sign(b-a) + Sign(c-b) + Sign(a-c);
 
 if ( ((help == 1) && (a!=c)) || ((a==b) && (b==c)) )
   return(1);
 else
   return(0);
}




/***********************************************************************************/


/* Sorts a field of subpaths first according to their left endnode from left to
   right, then sorts Paths with the same left endnode according to their right 
   endnode from right to left. */

void My_Qsort( p_field ,  left, right, circ_length ) 
path *p_field;
int left;
int  right;
int circ_length;

{
  int i, last, length_a, length_b;
  void My_Swap(path p_field[], int i, int j);

  if (left >= right)
    return;
  My_Swap( p_field, left, (left+right)/2);
  last = left;
  for (i=left+1; i<=right; i++)
    {
      if (p_field[i].left_node < p_field[left].left_node)
	My_Swap(p_field, ++last, i);

  /* If they have the same left endnode, sort them according to their right endnode
     such that the longer precedes the shorter one */
      if (p_field[i].left_node == p_field[left].left_node) 
	{
	  /* In the case of a main circuit the length of a subpath is a little more 
	     difficult to obtain.*/
	  if (p_field[i].left_node < p_field[i].right_node)
	    length_b = p_field[i].right_node - p_field[i].left_node;
	  else
	    length_b = p_field[i].right_node + (circ_length - p_field[i].left_node);

	  if (p_field[left].left_node < p_field[left].right_node)
	    length_a = p_field[left].right_node - p_field[left].left_node;
	  else
	    length_a = p_field[left].right_node + (circ_length - p_field[left].left_node);

	  if(length_a < length_b)   /* in this case we have to switch them */
	    My_Swap(p_field, ++last, i);
	}
    }
      My_Swap(p_field, left, last);
      My_Qsort(p_field, left, last-1, circ_length);
      My_Qsort(p_field, last+1, right, circ_length);
}

void My_Swap(path p_field[], int i, int j)
{
  path temp;

  temp = p_field[i];
  p_field[i] = p_field[j];
  p_field[j] = temp;
}

/***********************************************************************************/


/* Reads 1. the length of the main circuit
         2. cardinality of the family of subpaths
	 3. each subpath of the family (the number of endnodes)
         4. a freebie generator (if desired) 
   from a specified file */

path *Read_Input(circ_length, size, c)  /* returns a field of paths */
int *circ_length;                       /* a dynamic variable storing the length of the circ. */
int *size;                              /* size of the subfamily of subpaths */
char *c;                 
{  
  int i;
  int help, help2, subfam_size=0, gen_size=0;
  char file[50];
  FILE *fp1, *fp2=NULL;
  path *p_field;
  
  printf("Please type in the name of the inputfile: ");  /* reads the filename */
  scanf("%s",file);
  
  if (file[0] == '\0') fp1 = stdin;
  else fp1 = fopen(file,"r");
  
  if (fp1 == NULL) 
    {
      printf("\nRead_Input: file %s can't be opened\n",file);
      exit(0);
    }
  printf("Read_Input: Reading file %s.\n",file); 
    
  fscanf(fp1,"%d", circ_length);  /* length of the circuit */
  
  while(getc(fp1) != '\n') ;
  
  fscanf(fp1,"%d",size);          /* cardinality of the family of subpaths */
  subfam_size = *size;
  while(getc(fp1) != '\n') ;


  do
    {
      while(getchar() != '\n') ;
      printf("\nDo you wish to enter a freebie generator? (y/n):");
    }
  while ( ((*c=getchar()) != 'y') && (*c != 'n') );

  if (*c=='y')
    {
      printf("Please type in the name of the inputfile: ");  /* reads the filename */
      scanf("%s",file);
      
      if (file[0] == '\0') fp2 = stdin;
      else fp2 = fopen(file,"r");
      
      if (fp2 == NULL) 
	{
	  printf("\nRead_Input: file %s can't be opened\n",file);
	  exit(0);
	}
      printf("Read_Input: Reading file %s.\n",file); 
      
      fscanf(fp2,"%d",size);           /* cardinality of the freebie generator */
      gen_size = *size;
      while(getc(fp2) != '\n') ;
    }

  /*Allocating sufficient memory for a dynamic field of subpaths */
  p_field = (path*) malloc ( (subfam_size+gen_size) * sizeof(struct Path));
  

  for(i=0; i<subfam_size; i++)         /* reading the subpaths into the path-field */
    {
      if ( fscanf(fp1,"%d%d",&help,&help2) == EOF)
	{
	  printf("ERROR! Fewer distinct paths than expected!\n");
	  exit(0);
	}
      /*p_field[i].left_node = help;
      p_field[i].right_node = help2;*/

      p_field[i].left_node = Modulo(help, *circ_length);
      p_field[i].right_node = Modulo(help2, *circ_length);
    }
  fclose(fp1);

  i=0;
  My_Qsort(p_field, i, subfam_size-1, *circ_length);  /* sorting the path-family */  
  for(i=0; i<subfam_size; i++)     
    {
      p_field[i].id = i;                /* this way we mark them as the real paths and 
			                   save their index */
    }

  if( *c=='y')
    {
      for(i=subfam_size; i< (gen_size+subfam_size); i++) /* reading of the subpaths into the path-field */
	{
	  if ( fscanf(fp2,"%d%d",&help,&help2) == EOF)
	    {
	      printf("ERROR! Fewer distinct freebie subpaths than expected!\n");
	      exit(0);
	    }	  
	  p_field[i].left_node = help;
	  p_field[i].right_node = help2;
	  p_field[i].id = -1;       /* marking the members of the freebie generator, so we are 
                                       able to distinguish them from the real subpaths */
	}
      fclose(fp2);
    }

  *size = subfam_size + gen_size;   /* the altogether size of the path-family */
  return (p_field);
}

/***********************************************************************************/

/* Deletes all nodes (and adjacent arcs) which are no endpoint of any path and returns 
   a field in which the original values of each pathendpoint are stored for a later restruction. */ 

int *Reduce(p_field, subfam_size, circ_length)
path *p_field;     /* path-field */
int subfam_size;   /* size of the path-field */
int *circ_length;   /* pointer to the mainpath */

{
  extern void Screen_Subfam();   
  extern void My_Qsort();

  int *h_field, *red_unred, *unred_red;
  int i,hole_counter=1,  length;
  /* int red_counter = 0, l_old, l_new, max_right=0*/;

 
  length = *circ_length;
  printf("length = %d\n", length);

  h_field   = (int*) malloc ( length * sizeof(int) ); /* stores 1 if corresponding node is endpoint 
							 of some paths, 0 otherwise */
  unred_red = (int*) malloc ( length * sizeof(int) ); /* for the reduction: in block i of the field we
                                                         store the value of number i after the reduction */

  red_unred = (int*) malloc ( length * sizeof(int) ); /* to reverse the reduction: in block i of the field
                                                             - " -                   before the reduction */

  for (i=0; i<length; i++)
    { 
      /* initialization */
      red_unred[i] = i;
    }

  hole_counter = 0;          /* initialization */
  for (i=0; i<length; i++)
    h_field[i] = 0;          /* initialization */
  
  for(i=0; i< subfam_size; i++)
    {
      /* setting the values of h_field for the endpoints of the paths */
      h_field[p_field[i].left_node] = 1; 
      h_field[p_field[i].right_node] = 1;
    }
  for(i=0; i<length; i++)
    {
      if (h_field[i] == 0)  /* this is a node to delete (I call it a hole) */
	hole_counter++;
      else
	{
	  unred_red[i] = i - hole_counter;          /* setting the reduced values */   
	  red_unred[i-hole_counter] = red_unred[i]; /* storing the  original nonreduced values */
	}
    }
  /* and now we actually convert the unreduced endpoints into the reduced ones */
  for (i=0; i< subfam_size; i++)
    {
      /*l_old =  p_field[i].right_node -  p_field[i].left_node;*/
      p_field[i].left_node = unred_red[p_field[i].left_node];
      p_field[i].right_node = unred_red[p_field[i].right_node];
      /*if (p_field[i].right_node > max_right)
	max_right = p_field[i].right_node;   
      l_new = p_field[i].right_node -  p_field[i].left_node;*/
      /*red_counter += (l_old - l_new);*/
    }
  *circ_length = length - hole_counter;           /* modifying the main path's right node */
  printf("\nThe reduced length of the main circuit is: %d\n", *circ_length); 
  /*printf ("Number of deleted edges: %d\n", red_counter);*/
  free(unred_red);    /* we do not need this field anymore */
  return(red_unred);  /* returning the field to be able to reverse the reduction later */ 
}

/***********************************************************************************/

/* a simple ascii presentation of the subpaths on the screen */

void Screen_Subfam (circ_length, subfam_size, p_field)
int circ_length;
int subfam_size;
struct Path *p_field;

{  
  int i,j;
  
  printf("\nRepresentation of the family of subpaths on the screen:\n");
  printf("circ ( 0,%2d)|0",circ_length);
  for (i=1; i<=circ_length; i++)
    {
      printf("%4d",i);
    }
  printf("\n");
  for (i=0; i<subfam_size; i++)
    {
      printf("path (%2d,%2d)|",p_field[i].left_node,p_field[i].right_node);
      if ( p_field[i].left_node < p_field[i].right_node) 
	{
	  for (j=0; j < p_field[i].left_node; j++)
	    {
	      printf("    ");
	    }
	  
	  for (j=0; j< (p_field[i].right_node - p_field[i].left_node); j++)
	    {
	      printf("*-->");
	    }
	  printf("*\n");
	}
      else
	{
	  for (j=0; j < p_field[i].right_node; j++)
	    {
	      printf("*-->");
	    }
	  printf("*   ");
	   for (j=0; j< (p_field[i].left_node - p_field[i].right_node) -1; j++)
	    {
	      printf("    ");
	    }
	   for (j=p_field[i].left_node; j < circ_length; j++)
	     {
	       printf("*-->");
	     }
	   printf("*  (...)\n");
	}
    
    }
}


/***********************************************************************************/

/* Constructs a field of integer pairs with the length of the mainpath, storing for
each node i of the mainpath the index of the first and last path in p_field with
left endnode i. We actually construct the so called blocks defined in section 8.
Note that p_field is sorted such, that the paths in one block are in sequence. */ 
 
double_int *Reference_P_Field(p_field, subfam_size, circ_length)
path *p_field;
int subfam_size;
int circ_length;

{
  int i,j,k;


  register double_int *p_ref;

  /* allocating the sufficient memory */
  p_ref = (double_int*) malloc(circ_length * sizeof(struct Double_Int)); 

  /* initializing  the first empty blocks with -1s*/
  for(i=0; i < p_field[0].left_node; i++)
    {
      p_ref[i].first = -1;
      p_ref[i].last = -1;
    }

  /* setting the values for path number 0 */
  p_ref[p_field[0].left_node].first = 0;
  j = p_field[0].left_node;

  /* setting the values for all other paths */
  for(i=0; i < (subfam_size - 1); i++)
    {
      if (p_field[i+1].left_node > j) 
	{
	  p_ref[j].last = i;

	  /* initializing empty blocks with -1s */
	  for(k=j+1; k < p_field[i+1].left_node; k++)
	    {
	      p_ref[k].first = -1;
	      p_ref[k].last = -1;
	    }

	  j =  p_field[i+1].left_node;

	  p_ref[j].first = i+1;
	}
    }

  /* initializing remaining empty blocks with -1s */
  p_ref[j].last = subfam_size-1;
  for(k=j+1; k< circ_length; k++)
    {
      p_ref[k].first = -1;
      p_ref[k].last = -1;
    }

  return(p_ref);
}

/***********************************************************************************/

/* Constructs the structure E_R. E_R is a field of lists of essential pairs, whre the field
   has the length of the main circuit. A block of the field contains a list of essential 
   pairs, each of them consisting of a path and one of its arcs, called 'representing arc'. 
   In the list of block i of the field E_R, that is, in E_R[i] all essential pairs are contained,
   which are represented by arc number i.*/

e_pair **Construct_E_R(p_field, subfam_size, E_R_size, p_ref, circ_length)
path *p_field;
int subfam_size;
int *E_R_size;
double_int *p_ref;
int circ_length;

{ int i, k, m, size, first_orig, last_orig;
  e_pair *runptr, *insptr;
  e_pair **E_R, **Ins_field;
  short DONE[circ_length];
 
  register int j;
  register int l;
  register constr_essential *h_field;
 

  
  E_R   = (e_pair**) malloc ( circ_length * sizeof(e_pair*) );

  Ins_field = (e_pair**) malloc ( circ_length * sizeof(e_pair*) );
  /* Ins_field[i] holds a pointer to the essential pair in list E_R[i] behind which
     we have to insert the next member of the list. At first this will be the last member
     of the current list, but after a 'certain time' we will have to start inserting the new 
     members 'at the top' of the current list again. This is the case when we find the first 
     member, which belongs to the top of the current collumn, since it is 'smaller' the the
     other ones. */

  /*for(i=0; i<circ_length; i++)
  {
    printf("Block %d goes from path %d to path %d.\n",i,p_ref[i].first, p_ref[i].last);
  }*/


  for(i=0; i<circ_length; i++)
    {
      /* construction and init. of a dummy-element as the first member 
	 of each column's list */     
      E_R[i] = (e_pair*) malloc ( sizeof (struct E_pair) );
      E_R[i]->next = NULL;
      E_R[i]->prev = NULL;
      Ins_field[i] = E_R[i];  /* inserting the next member behind the dummy element */
      DONE[i] = 0;
   }
     
  /* allocating field for the storage of nextposs info and the list of 
     essential intervals for each path */ 
  h_field = (constr_essential*) malloc(subfam_size * sizeof(constr_essential));

  for (j=0; j<subfam_size; j++)
    h_field[j].lastess = NULL;        /* initializing the lastess pointers, pointing to the latest
					 essential pair found on the current j'th path */

   for (i=0; i<circ_length; i++)
    {
      /*printf("\ni=%d\n",i);
	printf("p_ref[i].first= %d, p_ref[i].last= %d\n",p_ref[i].first, p_ref[i].last);*/
     first_orig = p_ref[i].first;
     last_orig = p_ref[i].last;


      if(p_ref[i].first != -1)        /* considering nonempty blocks only */
	{
	  size = p_ref[i].last - p_ref[i].first + 1;   /* the size (that is, the number of paths) of the block */

	  /* initializing the nextposs info for the last member of the current block */
	  if(p_field[p_ref[i].last].id != -1)  /* it is not a member of the freebie gen. */
	    h_field[p_ref[i].last].nextposs = p_field[p_ref[i].last].left_node;
	  else                                 /* it is a member of the freebie gen. */
	    h_field[p_ref[i].last].nextposs = p_field[p_ref[i].last].right_node;

          /* initializing the nextposs info for all other members of the current block */
	  for(j=p_ref[i].last; j>p_ref[i].first; j--)
	    {
	      if(p_field[j-1].id != -1)
		h_field[j-1].nextposs = p_field[j].right_node;    /* this path is not in the freebie generator */
	      else
		h_field[j-1].nextposs = p_field[j-1].right_node;  /* this path is in the freebie generator */
	    }
	  k=Modulo(i+1, circ_length);
	  while ( (p_ref[k].first == -1) && 
		  ( ABC(p_field[p_ref[i].first].left_node, k, p_field[p_ref[i].first].right_node-1) ) )
	    k = Modulo(k+1, circ_length);             /* skipping empty blocks */
      
	  do
	    {
	      /*printf("   k=%d\n",k);*/

	      while ( ( p_ref[i].first <= p_ref[i].last) && 
		      ( ABC(p_field[p_ref[i].last].left_node, p_field[p_ref[i].last].right_node, k) ) )
		{

                     /* we may delete all paths from block i with a right node less or equal to k, which is the 
		     left endnode of all members of the currently considered 'k'th block of subpaths */   
		  j=p_ref[i].last;
		  if ( ABC(p_field[j].left_node, h_field[j].nextposs, p_field[j].right_node-1) )
		    {
		      /* before deleting them, we have to put the remaining essential 
			 pairs of these paths into E_R */
		      for (m = h_field[j].nextposs; 
			   ABC(p_field[j].left_node, m, p_field[j].right_node-1);
			   m = Modulo(m+1, circ_length) )
			{
			  /*printf("j=%d\n",j);
			  printf("h_field[10].nextposs=%d\n",h_field[10].nextposs);*/

			  /* Now we have to insert the remaining essential arcs found along the 
			     j'th path into E_R at a position according to the  predefined order */
					
			  /* inserting an element into the structure of E_R */
			  insptr = (e_pair*) malloc( sizeof(struct E_pair) );  /* allocating memory */
			  (*E_R_size)++;                        /* resetting the cardinality of E_R */
			  insptr->path_number = p_field[j].id;  /* setting the path number */
			  insptr->arc_number =  m;              /* setting the arc number  */ 
			  insptr->next_on_path = NULL;          /* initialization */
			  insptr->prev_on_path = NULL;          /* initialization */		  
			  insptr->cover = NULL;                   /* initialization */
			  if (h_field[j].lastess != NULL )  /* setting the next_on_path pointer */
			    {
			      h_field[j].lastess->next_on_path = insptr;
			      insptr->prev_on_path = h_field[j].lastess;
			    }

			  h_field[j].lastess = insptr;			  

			  if ( (m <= i-1) && (DONE[m] != 1) )
			       {
				 Ins_field[m] = E_R[m];  /* here we have to set the "insert-position" for the column
							    m back to the beginning of the "column-list" */ 
				 DONE[m] = 1;
			       }

                          /* double-linking into the list */
			  insptr->next = Ins_field[m]->next;
			  Ins_field[m]->next = insptr;       
			  insptr->prev = Ins_field[m];  
			  if(insptr->next != NULL)
			    insptr->next->prev = insptr;
			  Ins_field[m] = Ins_field[m]->next; /* setting the "ins.-position" one ahead */
		
			  
			} /* for (m = h_field[j].nextposs; m < p_field[l].left_node ; m++) */ 
		      
		    } /* if ( h_field[j].lastess < p_field[j].right_node) */
		  p_ref[i].last--; /* this is the next path in block i we have to consider */

		}  /*while ( ( p_ref[i].first <= p_ref[i].last) && 
		     ( ABC(p_field[p_ref[i].last].left_node, p_field[p_ref[i].last].right_node, k) ) )*/


	      j = p_ref[i].last;   /* help-indices j and l for the merging of the to blocks */
	
	      l = p_ref[k].last;
	      /*printf("j=%d, p_ref[i].first= %d, l=%d p_ref[k].first=%d\n",j, p_ref[i].first,l, p_ref[k].first);*/
	      /* Merging the blocks i and k formally. i is the currently considered block. */
	      while ( (j >= p_ref[i].first) && (l >= p_ref[k].first) )
		{

		  if( ABC(p_field[j].left_node,p_field[j].right_node, p_field[l].right_node-1) )
		    {
		      j--; /* seting j back */		     
		    }
		  else
		    {
		      /*printf("1. in else; i=%d, j=%d, k=%d, l=%d\n",i,j,k,l);*/
		      if ( (ABC(p_field[j].left_node, h_field[j].nextposs, p_field[l].left_node)) 
			   && (h_field[j].nextposs != p_field[l].left_node) )
			{
			  /* an ess. interval on path j was found */

			  for (m = h_field[j].nextposs; 
			       ABC(p_field[j].left_node, m, p_field[l].left_node-1);
			       m = Modulo(m+1, circ_length) )       		 
			    {
			      /* Now we have to insert the essential arcs found along the i'th path
				 into E_R at a position according to the predefined order */
			      
			      /*printf("jl=%d, m=%d, ll=%d\n",p_field[j].left_node, m, p_field[l].left_node-1);*/		      
			      insptr = (e_pair*) malloc( sizeof(struct E_pair) ); /* allocating memory */
			      (*E_R_size)++;                       /* resetting the cardinality of E_R */
			      insptr->path_number = p_field[j].id; /* setting the path number */
			      insptr->arc_number =  m;             /* setting the arc number  */ 
			      insptr->next_on_path = NULL;         /* initialization */
			      insptr->prev_on_path = NULL;         /* initialization */
			      insptr->cover = NULL;               /* initialization */

			      if (h_field[j].lastess != NULL )     /* setting the next_on_path pointer */
				{
				  h_field[j].lastess->next_on_path = insptr;
				  insptr->prev_on_path = h_field[j].lastess;
				}

			      h_field[j].lastess = insptr;
			      
			      /* looking for the proper position to insert the element */
    			     
			      if ( (m <= i-1) && (DONE[m] != 1) )
				   {
				     Ins_field[m] = E_R[m];  /* here we have to set the "insert-position" for the column
								m back to the beginning of the "column-list" */ 
				     DONE[m] = 1;
				   }

			      /* double-linking into the list */
			      insptr->next = Ins_field[m]->next;
			      Ins_field[m]->next = insptr;       
			      insptr->prev = Ins_field[m];  
			      if(insptr->next != NULL)
				insptr->next->prev = insptr;
			      Ins_field[m] = Ins_field[m]->next; /* setting the "ins.-position" one ahead */

			    } /*  for ... ( (ABC(p_field[j].left_node, m, p_field[l].left_node))
				  &&  (h_field[j].nextposs != p_field[l].left_node) ) */
			   h_field[j].nextposs=m;  /* resetting the nextposs value */
      		       
			} /* if (p_field[l].left_node > h_field[j].nextposs) */
	    
                     if( ABC(p_field[j].left_node, h_field[j].nextposs, p_field[l].right_node-1) )
			h_field[j].nextposs = p_field[l].right_node; /* setting nextposs info for path j */

		      l--;   /* setting l back */
		    }
		}   /*while ( (j >= p_ref[i].first) && (l >= p_ref[k].first) )*/

	      k = Modulo(k+1, circ_length);   /* looking for the next block k to compare with block i */ 
	      while ( (ABC(p_field[p_ref[i].first].left_node, k, p_field[p_ref[i].first].right_node-1)) && 
		      (p_ref[k].first == -1) )
 		k = Modulo(k+1, circ_length); /* skipping empty blocks */
	    }
	  while (p_ref[i].first <= p_ref[i].last);  /* as long as block i is not empty */   	       
      
	}/* if(p_ref[i].first != -1) */
      p_ref[i].first = first_orig;
      p_ref[i].last = last_orig;

      } /*for (i=0; i < circ_length; i++)*/

   for (i=0; i < circ_length; i++)
     { 
       runptr = E_R[i];   /* skipping the dummy element */ 
       E_R[i] = E_R[i]->next;
       if (E_R[i]!= NULL)
	 E_R[i]->prev = NULL;
       free(runptr);
     }

   /*for(i=0; i<circ_length; i++)
    {
      runptr = E_R[i];
      printf("i=%d\n",i);
      while (runptr != NULL)
	{
	  printf(" (p:%d, a:%d) ",runptr->path_number, runptr->arc_number);
	  runptr = runptr->next;
	}
      printf("\n");
    }*/

   return(E_R);    /* returning a pointer to the first member of the list */

} /*Constuct_E_R*/

/***********************************************************************************/

/* This function removes the members of the freebie generator from p_field, since
   we do not need them anymore */
void Remove_Free_Gen(p_field, subfam_size)
path *p_field;


int *subfam_size;

{
  int i,counter=0;

  for (i=0; i<*subfam_size; i++)
    {
      if(p_field[i].id == -1)
	counter++;
      else
	p_field[i-counter] = p_field[i];
    }

  *subfam_size = *subfam_size - counter;   /* resetting the size of the subfamily of subpaths */
}
      

/***********************************************************************************/
/* An ascii representation of the structure E_R of essential pairs on the screen.*/

void Screen_E_R2(circ_length, p_field, subfam_size, E_R)
int circ_length;
path *p_field;
int subfam_size;
e_pair **E_R;

{
  int i, j, flag=0;
  int block[subfam_size+1][circ_length+1];  /* this block helps the representation */ 
  e_pair *runptr;

  printf("subfam_size = %d, circ_lenght = %d\n",subfam_size, circ_length);

  for(i=0; i<subfam_size; i++)
    for(j=0;j<circ_length; j++)
      block[i][j] = 0;    /* initializing */
  
      /* initializing the position of the subpaths */
  for(i=0; i<subfam_size; i++)
    {
      if (p_field[i].left_node < p_field[i].right_node)
	{
	  for(j=p_field[i].left_node; j<p_field[i].right_node; j++)
	    block[i][j] = 1;
	}
      else
	{
	  for(j=0; j<p_field[i].right_node; j++)
	    block[i][j] = 1;
	  for(j=p_field[i].left_node; j<circ_length; j++)
	    block[i][j] = 1;
	}
   } 

  /*j=p_field[i].left_node;
  while( ABC(p_field[i].left_node, j, p_field[i].right_node-1) )
    { 
      block[i][j] = 1;
      j = Modulo(j+1, circ_length);
    }*/
   

  for (i=0; i<circ_length; i++)
    {
      runptr = E_R[i];
      /*printf("i=%d , ",i);*/
      while (runptr != NULL)
	{
	  /*printf("runptr: (%d, %d)",runptr->path_number, runptr->arc_number);*/
	  block[runptr->path_number][runptr->arc_number] = 2;
	  runptr = runptr->next;
        }
      /*printf("\n");*/
    }
  
  printf("\nRepresentation of the essential pairs on the screen:\n");
  printf("circ ( 0,%2d)|0",circ_length);
  for (i=1; i<=circ_length; i++)
    {
      printf("%4d",i);
    }
  printf("\n");


  for(i=0; i<subfam_size; i++)
    {
      printf("path (%2d,%2d)|",p_field[i].left_node,p_field[i].right_node);  
      if(p_field[i].right_node == 0)
	printf("*");

      for(j=0;j<circ_length; j++)
	{
	  if ( block[i][j] == 0)
	    {
	      if ( (p_field[i].right_node < p_field[i].left_node) && 
		   (j < circ_length-1) && (block[i][j+1] != 0) )
		{
		  printf("   ");
		  flag = 0;
		}
	      else
		printf("    ");
	    }
	  else
	    {
	      if (flag == 0)
		{
		  printf("*");
		  flag = 1;
		}
	      if (block[i][j] == 1)
		printf("-->*");
	      else
		printf("==>*");
	    }
	}
      flag = 0;
      printf("\n");
    }
}



/***********************************************************************************/

/* A simple ascii presentation of E_R on the screen  */
void Screen_E_R(circ_length, p_field, subfam_size, E_R)
int circ_length;
path p_field[];
int subfam_size;
e_pair **E_R;

{
  int i, j, k, FOUND;
  e_pair *runptr=NULL;

  printf("\nRepresentation of the essential pairs on the screen:\n");
  printf("circ ( 0,%2d)|0",circ_length);
  for (i=1; i<=circ_length; i++)
    {
      printf("%4d",i);
    }
  printf("\n");
  for (i=0; i<subfam_size; i++)
    {
      printf("path (%2d,%2d)|",p_field[i].left_node,p_field[i].right_node);
      if ( p_field[i].left_node < p_field[i].right_node) 
	{
	  for (j=0; j < p_field[i].left_node; j++)
	    {
	      printf("    ");
	    }
	  
	  /* looking for the first essential arc on path i  */
	  k=p_field[i].left_node;
	  FOUND = 0;
	  while ( ( k < p_field[i].right_node ) && (FOUND == 0) ) 
	    {
	      runptr = E_R[k];

	      while( (runptr != NULL) && (runptr->path_number != i) )
		{
		  runptr = runptr->next;
		}

	      if (runptr != NULL)
		FOUND = 1;

	      k++;
	    }

	  if( FOUND == 1 )   /* at least one was found */
	    {
	      for (j=0; j < (runptr->arc_number - p_field[i].left_node); j++)
		{
		  printf("*-->");
		  
		}
	      /*printf("*==>");*/
	      printf("*%2d>", runptr->id);
	      
	      while( runptr->next_on_path != NULL )
		{
		  for(j=0; j < ( runptr->next_on_path->arc_number 
				 - runptr->arc_number - 1) ; j++)
		    {
		      printf("*-->");
		    }
		  /*printf("*==>");*/
		  printf("*%2d>", runptr->next_on_path->id);
		  
		  runptr = runptr->next_on_path;
		}
	      for (j = runptr->arc_number + 1; j<p_field[i].right_node; j++)
		{
		  printf("*-->");
		}
	    }
	  else  /* to if( runptr->path_number == i ) */
	    {
	      for (j=p_field[i].left_node; j<p_field[i].right_node; j++)
		{
		  printf("*-->");
		}
	    }
	  printf("*\n");

	}
      else   /* to if ( p_field[i].left_node < p_field[i].right_node) */ 
	{
	  /* now in this case we have to consider both parts of the path (0...right endnode)
	     and (left endnode...circ_length) separately */


	  
	  /* looking for the first essential arc on the first considered part (0...right endnode)
	   of the path */
	  k=0;
	  FOUND = 0;
	  while ( ( k < p_field[i].right_node ) && (FOUND == 0) ) 
	    {
	      runptr = E_R[k];
	      while( (runptr != NULL) && (runptr->path_number != i) )
		{
		  runptr = runptr->next;
		}

	      if (runptr != NULL)
		FOUND = 1;

	      k++;
	    }

	  if( FOUND == 1 )   /* at least one was found */
	    {
	      for (j=0; j < runptr->arc_number; j++)
		{
		  printf("*-->");
		}
	      /*printf("*==>");*/
	      printf("*%2d>", runptr->id);
	      
	      while( runptr->next_on_path != NULL )
		{
		  for(j=0; j < ( runptr->next_on_path->arc_number 
				 - runptr->arc_number - 1) ; j++)
		    {
		      printf("*-->");
		    }
		  /*printf("*==>");*/
		  printf("*%2d>", runptr->next_on_path->id);
		  
		  runptr = runptr->next_on_path;
		}
	      for (j = runptr->arc_number + 1; j<p_field[i].right_node; j++)
		{
		  printf("*-->");
		}
	      printf("*   ");
	    }
	  else
	    {
	      /* that is, no essential arc was found in the first considered part of the path */
	      for (j=0; j < p_field[i].right_node; j++)
		{
		  printf("*-->");
		}
	      printf("*   ");
	    }

	  /* filling the space between the to parts of the path */
	  for (j=0; j< (p_field[i].left_node - p_field[i].right_node) -1; j++)
	    {
	      printf("    ");
	    }
	  

	  /* looking for essential arcs in the second considered part (left endnode...circ_length)
	     of the path */
	  k=p_field[i].left_node;
	  FOUND = 0;
	  while ( ( k < circ_length ) && (FOUND == 0) ) 
	    {
	      runptr = E_R[k];

	      while( (runptr != NULL) && (runptr->path_number != i) )
		{
		  runptr = runptr->next;
		}

	      if (runptr != NULL)
		FOUND = 1;

	      k++;
	    }

	  if( FOUND == 1 )   /* at least one was found */
	    {
	      for (j=0; j < (runptr->arc_number - p_field[i].left_node) ; j++)
		{
		  printf("*-->");
		}
	      /*printf("*==>");*/
	      printf("*%2d>", runptr->id);
	      
	      while( (runptr->next_on_path != NULL) && 
		     (runptr->next_on_path->arc_number >= p_field[i].left_node) )
		{
		  for(j=0; j < ( runptr->next_on_path->arc_number 
				 - runptr->arc_number - 1) ; j++)
		    {
		      printf("*-->");
		    }
		  /*printf("*==>");*/
		  printf("*%2d>", runptr->next_on_path->id);
		  
		  runptr = runptr->next_on_path;
		}
	      for (j = runptr->arc_number + 1; j<circ_length; j++)
		{
		  printf("*-->");
		}
	      printf("* (...)\n");
	    }
	  else
	    {
	      for (j=p_field[i].left_node; j < circ_length; j++)
		{
		  printf("*-->");
		}
	      printf("* (...)\n");
	    }
 	}
    }
}

/***********************************************************************************/
/* NOT WORKING!*/
/* a siple ascii presentation of the data in E_R */
void Print_E_R(E_R)
e_pair *E_R;

{
  e_pair *runptr;

  runptr=E_R;

  while( runptr->next != NULL)
    {


      printf(" (path, arc, next_on_path->path, next_on_path->arc: ");
      printf("(%d, %d, " ,runptr->path_number, runptr->arc_number);


      if (runptr->next_on_path != NULL)
	{
	  printf("%d, %d)\n",runptr->next_on_path->path_number, 
		 runptr->next_on_path->arc_number);
	}
      else
	{
	  printf("NULL, NULL)\n");
	}

      printf(" ( id: %d;  prev: %d;  next: %d )\n ", runptr->id, runptr->prev->id, runptr->next->id);  

      runptr = runptr->next;
    }
}

/***********************************************************************************/

/* Executes phase 1 of the algorithm, that is, a crossfree subset K of selected 
ess. repr. paths and a subset T of forbidden ess. repr. paths is being constructed */

forbidden *Crossfree(p_field, E_R, E_R_size, circ_length)
path *p_field;
e_pair **E_R;
int *E_R_size;
int circ_length;

{
  e_pair *run1, *run2;
  e_pair *run3;
  forbidden *T, *last;
  int T_counter=0, i=0;


  /* constructing a double-linked list of forbidden ess. pairs */
  T = (forbidden*) malloc( sizeof(struct Forbidden) );
  T->forb = NULL;                           /* a dummy element */
  T->next = NULL;
  T->prev = NULL;
  last = T;

  /* going along the field E_R */
  for (i=0;i<circ_length; i++)
    {
      /* going along the list in block E_R[i] */
      run1 = E_R[i];
      while ( run1 != NULL ) 
	{
	  /*printf("Forb. by (%d, %d): ",run1->path_number, run1->arc_number);*/

	  /* Forbidding to the left: */
          if (run1->prev != NULL)
	    {
	      run2 = run1->prev;
	      run3 = run2->prev_on_path;
	      while ( ( run3 != NULL ) &&
		      ( ABC(p_field[run1->path_number].left_node, run3->arc_number, run1->arc_number) ) )
		{
		 
		
		  /* 'Oneway' unlinking from the list of essential pairs. This means, that the pointers
		     from the previous resp. following member to run3 in the list are relinked as if run3 
		     would have been deleted, but the pointers from the run3 element to the others of the 
		     list are kept linked. */

                  if (run3->prev != NULL)
		    run3->prev->next = run3->next;      /* resetting the pointer of the previous pair if existent */
		  else
		    E_R[run3->arc_number] = run3->next; /* resetting the corresponding E_R pointer if run3 is the 
							   topmost member of the column */
		  
   		  if(run3->next != NULL)                /* are still set */
		    run3->next->prev = run3->prev; 

		  run3->prev = run3->next;              /* now the run3->prev pointer points to the union of 
							   the two pairs run3 and run1 */
		  run3->next = run2;                    /* now the run3->next pointer points to the intersection of 
							   the two pairs run3 and run1 */
		  
		  /*printf("(%d, %d) with intersect.= (%d,%d) union = (%d,%d) , ",run3->path_number, run3->arc_number, 
		    run3->next->path_number, run3->next->arc_number, run3->prev->path_number, run3->prev->arc_number);*/
		  

		  T->prev = (forbidden*) malloc ( sizeof(struct Forbidden) );
		  T->prev->next = T;                /* inserting a new element at the end of T */  
		  T = T->prev;
		  T->forb = run3;                   /* putting into the double-linked list of T */
		  T_counter++;                      /* counting the cardinality of T */
		 
		  run3 = run3->prev_on_path;        /* setting run3 on the path one ahead */ 
		}

	      run2->prev_on_path = run3;            /* resetting the prev_on_path resp. next_on_path */
	      if(run3 != NULL)                      /* pointers avoiding the deleted elements */
		run3->next_on_path = run2;
	    } /* if (run1->prev != NULL) */


	  /* Forbidding to the right: */
	  run2 = run1->next;
	  if (run2 != NULL) 
	    {
	      /* Go along the second ess. pair's (run2's) path and forbid all essential
		 pairs 'on it', which have their repr. arc in the intersection of 
		 the two paths (the paths of run1 and run2) */
	      run3 = run2->next_on_path;
	      while ( (run3 != NULL) && 
		      (ABC(run2->arc_number, run3->arc_number, p_field[run1->path_number].right_node-1)) )
		{
		  /* setting the 'forbidden_by' pointer of the pair */	  
		  /*if (run3->prev->forb_by != NULL) 
		    run3->forb_by = run3->prev->forb_by;
		    else
		    run3->forb_by = run1;*/
		
		  /* 'oneway' unlinking from the list of essential pairs */
		  run3->prev->next = run3->next;    /* the pointers of the forbidden pair */
		  if(run3->next != NULL)            /* are still set */
		    run3->next->prev = run3->prev;   

		  
		  run3->next = run2;                /* this points now to the inersection of the two pairs run1 and run3,
						     while the run3->prev pointer points to their union */

		  /*printf("(%d, %d) with intersect.= (%d,%d) union = (%d,%d) , ",run3->path_number, run3->arc_number, 
			 run3->next->path_number,run3->next->arc_number, run3->prev->path_number, run3->prev->arc_number);*/
		 
		  T->prev  = (forbidden*) malloc ( sizeof(struct Forbidden) );
		  T->prev->next = T;  
		  T = T->prev;
		  T->forb = run3;                   /* putting into the double-linked list of T */
		  T_counter++;                      /* counting the cardinality of T */
		 
	
		  run3 = run3->next_on_path;
		}
	      run2->next_on_path = run3;
	      if (run3 != NULL)
		run3->prev_on_path = run2;
	    }
	  
	  run1 = run2;
	  /*printf("\n");*/
	} /* while (run1 != NULL) */


    } /* for ... */





  *E_R_size = *E_R_size - T_counter;            /* resetting the cardinality of E_R */
  printf("The number of forbidden pairs is %d\n",T_counter);
  T->prev = NULL;
  last = last->prev;                                  /* skipping the dummy-element */
  if (last != NULL)
    last->next = NULL;
  else
    T = NULL;

  return (T);                                   /* returning the double-linked list T */

}
/***********************************************************************************/

/* Here we are setting the pointers in the blocks of the field PER, which will then
   point to the first non-outcrossed essential member of the according path or to NULL */
void Construct_PER( E_R, circ_length, PER )
e_pair **E_R;
int circ_length;
e_pair **PER;
     
{
  int i;
  e_pair *run;
  
  for(i=0;i<circ_length; i++)
    {
      run = E_R[i];
      while (run != NULL)
	{
	  if( run->prev_on_path == NULL )
	    PER[run->path_number] = run;
	    
	  run = run->next;
	}
    }
}


/***********************************************************************************/

/* Constructs a dynamic field of pointers to the elements of E_R. Denote that 
   the pointers to members of E_R have the same idex in the field as the member's id. 
   We will use this coincidence later, when we have to transform back the results of the 
   Dilwoth-algorithm into the original system of subpaths. */

e_pair **Reference_E_R(E_R, circ_length, E_R_size) 
e_pair **E_R;
int circ_length;
int E_R_size;

{
  int i=1, j;
  e_pair *run; 
  e_pair **E_R_field;

  E_R_field = (e_pair**) malloc ( (E_R_size +1) * sizeof(e_pair*) );
 
  for (j=0; j<circ_length; j++)
    {
      run = E_R[j];
      while (run != NULL)
	{
	  
	  E_R_field[i] = run;
	  run->id = i;
	  i++;
	  run = run->next;
	}
    }

  return(E_R_field);  /* returns a field of pointers to E_R */
}

/***********************************************************************************/

/* initializes the bipartite graph */
bip_graph *Init_Bip_Graph(E_R_size)
int E_R_size;
{
  bip_graph *G;
  int i;

  G = (bip_graph*) malloc( sizeof( struct Bip_Graph ) );
  G->n_nodes = E_R_size;          /* number of nodes (in each partition) */
  G->n_edges = 0;                 /* number of edges still 0 */
  G->nodes1  = (node1*) malloc ( (E_R_size + 1) * sizeof(struct Node1));  /* field of nodes in partition 1 */
  G->nodes2  = (node2*) malloc ( (E_R_size + 1) * sizeof(struct Node2));  /* field of nodes in partition 2 */

  for ( i=0; i<= G->n_nodes; i++)
    {
      /* It is sufficient for our our purposes to store an edge-list for
	 the nodes in partition one only. */
      G->nodes1[i].match = 0;                                   /* no matching edges yet */
      G->nodes1[i].elist = (edge*) malloc(sizeof(struct Edge)); /*dummy element for the edge-list */
      G->nodes1[i].elist->next = NULL;
      G->nodes1[i].elist->adjac = -1;
      G->nodes1[i].last = G->nodes1[i].elist;
                                    
      G->nodes2[i].match = 0;
      
    }
  return(G);                                                    /* returning the graph */
}

/***********************************************************************************/
/* Initializes the dynamic field bottom, which is a helpful structure to
   locate all neighboring edges of the bip. graph, see next routine. */

void Init_Bottom(E_R, p_field, bottom, circ_length)
e_pair **E_R;
path *p_field;
e_pair **bottom;
int circ_length;

{
  int i, max;
  e_pair *run;

  for(i=0; i<p_field[0].left_node; i++)
    bottom[i] = NULL;

  for(i=p_field[0].left_node; i<circ_length; i++)
    {
      if (E_R[i] != NULL)
	{
	  run = E_R[i];
	  max = run->path_number;
	  while ( (run->next!= NULL) && (run->next->path_number >= max) )
	    {
	      run = run->next;
	      max = run->path_number;
	    }  

	  if( i < p_field[run->path_number].left_node)
	    {
	      bottom[i] = run;
	    } 
	  else
	    bottom[i] = NULL;
	}
      else
	bottom[i] = NULL;
    }
}


/***********************************************************************************/

/* Finds all neighboring edges of the bipartite graph, that is, finds all comparable pairs
   of K following one another. (see step 5 in section 8) */     
bip_graph *Compare ( p_field, E_R, E_R_size, PER, subfam_size, circ_length )
path *p_field;
e_pair **E_R;
int E_R_size;
e_pair **PER; /* PER is such a field of pointers, that PER[i] points to the first 
		 member of K on the 'i'th subpath */ 
int subfam_size;
int circ_length;

{  e_pair *run1, *run2;
   edge *help_ptr;
   int i, k, edgecounter=0, blocker, blocker_old;
   e_pair **bottom;

   register bip_graph *G;


   /*for(i=0; i<subfam_size; i++)
     {
       if(PER[i] != NULL)
	 printf("PER[%d]=(%d, %d)\n",i, PER[i]->path_number, PER[i]->arc_number);
     }*/
   G = Init_Bip_Graph(E_R_size);                         /* initializing graph G */ 


   /* Bottom is a special field of pointers to the members of K, where bottom[e] points to the
      actual 'bottom-most member' of K with arc e, for this is a candidate for a following
      element to the currently considered member of K */ 
   bottom = (e_pair**) malloc (circ_length * sizeof(e_pair*)+1); 

   /*if (bottom == NULL)
     printf("bottom-field could not be allocated!\n");*/
   
   /*memset(bottom,0,length * sizeof(e_pair*));*/
   Init_Bottom(E_R, p_field, bottom, circ_length);

   /*printf("Bottom initialized to:\n");
   for(i=0; i<circ_length; i++)
     {
       if(bottom[i] != NULL)
	 printf("bottom[%d] = %d, ", i, bottom[i]->id);
     }
     printf("\n");*/

   /* We take the subpaths one by one according to their enumeration and for each subpath 
      we go along its members of K from left to right and take then the next subpath */
   for( k=0; k<subfam_size; k++ )
     {  
       /*printf("k=%d\n",k);*/
       run1 = PER[k];    /* we go along the field PER, where PER[k] points to the 
			    leftmost member of K on subpath number k */

       /* updating the bottom field with the members of K on the currently considered path */
       run2 = run1;
       while(run2 != NULL)
	 {
	   bottom[run2->arc_number] = run2;
	   run2 = run2->next_on_path;
	 }
       
       /* looking for edges of the bip. graph with member run1 on one end */
       while ( run1 != NULL )
	 { 

	   /* looking for a directed edge 'above' run1's edge */
	   blocker = Modulo(run1->arc_number+1, circ_length);  /* initializing the blocker variable */
	   run2 = run1->prev;
	   if (run2 != NULL) 
	     {
	       /* following element to run1 was found */


	       blocker = p_field[run2->path_number].right_node; /* blocking members of K according
								   to the rule in step 5 section 8 */
	       /*printf("run1->id = %d\n",run1->id);
	       printf("run2->id=%d\n", run2->id);*/	     
	       /* inserting the corresponding neighboring edge into G */


	       G->nodes1[run1->id].last->next = (edge*) malloc( sizeof(struct Edge) );

    	       if ( G->nodes1[run1->id].last->next == NULL)
		 printf("Ohoooo!\n");

	      G->nodes1[run1->id].last = G->nodes1[run1->id].last->next;
	      
	       if ( G->nodes1[run1->id].last == NULL)
		 printf("Ohoooo2!\n");


	       /*printf("run1->id = %d\n",run1->id);
		 printf("run2->id=%d\n", run2->id);*/

	       G->nodes1[run1->id].last->adjac = run2->id;	   
	       edgecounter++;                                   /* counting the edges */
	       G->nodes1[run1->id].last->next = NULL;
	     }
	     
	   i=blocker; /* blocker was either set already or it has 
			 the value of the column next to run1 */

	   if(run1->id == 5)
	     {
	       /*printf("blocker3 = %d\n",blocker);
		 if (bottom[9] != NULL)
		 printf("bottom[9] = %d\n", bottom[9]->id);*/
	     }
  
	   /*printf("1. blocker = %d\n",blocker);*/

	   blocker = run1->arc_number;
	   blocker_old = blocker;
	   while ( ABC(blocker_old, i, p_field[run1->path_number].right_node-1) )
	     {

	       /*printf("run1->arc_number: %d, path: (%d,%d), i = %d\n", run1->arc_number, p_field[run1->path_number].left_node,
		 p_field[run1->path_number].right_node, i);*/
	       if (bottom[i] != NULL)
		 {
		   /* a following element to run1 was found */
		   blocker_old = blocker;
		   blocker = p_field[bottom[i]->path_number].right_node; /* blocking members of K according 
									    to the rule of step 5 in section 8 */  

		   /* inserting the corresponding neighboring edge into G */
		   G->nodes1[run1->id].last->next = (edge*) malloc( sizeof( struct Edge ) );

		   if ( G->nodes1[run1->id].last->next == NULL)
		     printf("Ohoooo3!\n");


		   G->nodes1[run1->id].last =  G->nodes1[run1->id].last->next;
		   G->nodes1[run1->id].last->adjac = bottom[i]->id;	   

		  
		   edgecounter++;                                /* counting the edges */
		   G->nodes1[run1->id].last->next = NULL;
		   i = blocker;                                  /* considering next unblocked member of K */
		 }
	       else
		 {
		   i = Modulo(i+1, circ_length);                 /* considering next unblocked member of K */
		 }
	     }

	   run1 = run1->next_on_path;                            /* go along the members of K on actual path k */
	 } /* while ( run1 != NULL ) */

       if (k+1 < subfam_size)
	 {
	   i = p_field[k].left_node; 
	   while ( (ABC(p_field[k].left_node, i, p_field[k+1].left_node)) 
		   && (i !=  p_field[k+1].left_node) )
	     {
	       bottom[i] = NULL;
	       i = Modulo(i+1, circ_length);
	     }
	 }
       
     }  /* for (k=0; k<subfam_size; k++)  */

   G->n_edges = edgecounter;
   for (i=1; i<=G->n_nodes; i++)                 /* removing dummy-element from the beginning of the edge-lists */
     {
       help_ptr = G->nodes1[i].elist;
       G->nodes1[i].elist = help_ptr->next;
       free(help_ptr); 
     }
   return(G);                                    /* returning the graph */            
} /*Compare*/


/***********************************************************************************/

/* writes the edges of the bip. graph in the file Graph.in in Dimacs format if wanted */
void Write_Bip_Graph(G)
bip_graph *G;

{
  char file[]="Graph.in";
  edge *run;
  int i;
  FILE *fp3 = NULL;
 

  fp3 = fopen(file,"w");
  if (fp3 == NULL) 
    {
      printf("\nWrite_Bip_Graph(G): file %s can't be opened\n",file);
      exit(0);
    }  
  /*else
    {
      printf("\nfile %s opened successfully!\nfp3->_cnt = %d\n",file, fp3->_cnt);
      if (fp3->_ptr != NULL)
	printf("fp3->_ptr = %c\n",fp3->_ptr);
      else
	printf("fp3->_ptr = NULL\n");
      if (fp3->_base != NULL)
	printf("fp3->_base = %c\n",fp3->_base);
      else
	printf("fp3->_base = NULL\n");
      
      printf("fp3->_flag = %c\n",fp3->_flag);
      printf("fp3->_file = %c\n",fp3->_file);
    }*/


  fprintf(fp3,"%c",'p');
  fprintf(fp3,"%s"," edge ");
  fprintf(fp3,"%d %d\n", G->n_nodes, G->n_edges);

  for (i=1; i<=G->n_nodes; i++)
	{
	  run = G->nodes1[i].elist;
	  while ( run!=NULL)
	    {
	      fprintf(fp3,"%c ",'e');
	      fprintf(fp3,"%d %d\n",i, run->adjac );
	      run = run->next;
	    }
	}

  fclose(fp3);
}




/***********************************************************************************/

/* derives a maximum independent subset and writes  
   it to a file named Independent and returns its cardinality */ 
int  Max_Indep_Subs (p_field, E_R_field, G, red_unred, R)
path *p_field;
e_pair **E_R_field;
bip_graph *G;
int *red_unred;
int *R;    /* help field from the extern function Bip_Match() */

{
  int i, number=0;
  char file[] = "Independent";
  FILE *fp;


  fp = fopen (file,"w");
  if (fp == NULL) 
    {
      printf("\nMax_Indep_Subs: file %s can't be opened\n",file);
      exit(0);
    }

  fprintf(fp,"The following subfamily of represented paths\n");
  fprintf(fp,"is independent with maximal cardinality:\n");  

  /* using the info of the latest alternating forest found in Bip_Match() */
  for(i=1; i<=G->n_nodes; i++)
    {
      if ( (R[i] == 1) &&
	   ( (G->nodes1[i].match == 0) ||
	     ( (G->nodes1[i].match != 0) && (R[G->nodes1[i].match] == 3) ) ) )
	{
	 
	  /* writing the unreduced values using the red_unred field constructed in function Reduce() */
	  fprintf(fp,"Path(%d, %d) ", red_unred[p_field[E_R_field[i]->path_number].left_node],
		  red_unred[p_field[E_R_field[i]->path_number].right_node]);
	  
	  fprintf(fp,"Arc:%d, Ess.id = %d \n", red_unred[E_R_field[i]->arc_number], i); 
	
	  number++;
     	}
    }
printf("\n");
  fclose(fp);
  return (number); /* returns sigma(F) (see Gyori's theorem for subpaths) */
}

/***********************************************************************************/

/* constructs a minimum Covering of K from the max. matching on G */
chain_tree *Min_Chain_Decomp(p_field, circ_length, E_R_field, E_R_size, G, Gyori_size)
path *p_field;
int circ_length;
e_pair **E_R_field;
int E_R_size;
bip_graph *G;
int Gyori_size;  /* cardinality of the independent subsyst. */

{
  int i=1, j, k=0;
  chain_tree *chain_field; 

  /* the user may uncomment the folowing commands and obtain 
     an output of the minimum chaindecomposition of K */

  /*printf("\nA minimum chaindecomposition of K is:\n");*/
  chain_field = (chain_tree*) malloc ( (Gyori_size) * sizeof(chain_tree) );

  for( i=1; i<=G->n_nodes; i++)
    {
      if ( (G->nodes1[i].match != 0) && (G->nodes2[i].match == 0) )
	{
	  /* that is, beginning of a chain was found in partition1 */
	  chain_field[k].left_node = Modulo(E_R_field[i]->arc_number, circ_length);

	  E_R_field[i]->cover = chain_field + k;
	  
	  j = G->nodes1[i].match;
	  E_R_field[j]->cover = chain_field + k;
	  /*printf ("%d. chain:{%d, %d",k,i,j);*/
	  while ( G->nodes1[j].match != 0)
	    {
	      /* go along the chain */
	      /*printf(", %d, %d",j,G->nodes1[j].match);*/ 
	      j = G->nodes1[j].match;
	      E_R_field[j]->cover = chain_field + k;

	    }

	  chain_field[k].right_node = Modulo(E_R_field[j]->arc_number+1, circ_length);
	  chain_field[k].id = k;
	  chain_field[k].l_branch = NULL;
	  chain_field[k].r_branch = NULL;
     
	  /*printf ("}, Endnodes: (%d, %d)\n", chain_field[k].left_node, chain_field[k].right_node);*/
	  k++;
	}
      if ( (G->nodes1[i].match == 0) && (G->nodes2[i].match == 0) )
	{
	  /* Chain of length one */   
	  chain_field[k].left_node = Modulo(E_R_field[i]->arc_number, circ_length);
	  chain_field[k].right_node = Modulo(E_R_field[i]->arc_number+1, circ_length);
	  chain_field[k].id = k;
	  chain_field[k].l_branch = NULL;
	  chain_field[k].r_branch = NULL;

	  E_R_field[i]->cover = chain_field + k;
	  
	  /*printf ("%d. chain {%d} with Endnodes(%d, %d)\n",k, i, 
		  chain_field[k].left_node, chain_field[k].right_node);*/
	  k++;
	}
    }

  /*for(i=1; i<=E_R_size; i++)
    printf(" %d covered by the %d. chain ",i ,E_R_field[i]->cover);
  printf("\n");*/

  return(chain_field); /* returning the covering of K */
}


/***********************************************************************************/


int Covered(pair, chain_part, help_chain, Gyori_size, p_field)
e_pair *pair;
chain *chain_part;
chain_tree *help_chain;
int Gyori_size;
path *p_field;

{

  e_pair *prev, *next;
  chain_tree *run;


  prev = pair->prev;
  next = pair->next;
 
  run = prev->cover;  /* here run cannot be NULL, for proof see the paper */

  /* as long as there is a left branch at all, that is, the current cover of the pair prev is
     no longer existent, since its right endnode has been exchanged by some other */
  while(run->l_branch != NULL)
   {
     if ( (ABC(p_field[prev->path_number].left_node, run->l_branch->left_node, prev->arc_number)) &&
	  (ABC(prev->arc_number+1, run->l_branch->right_node, p_field[prev->path_number].right_node)) )
       {
	 /* in this case we take the left branch of the tree of coveredges */
	 run = run->l_branch;
       }
     else
       {
	 /* otherwise we take the right branch of the tree of coveredges */
	 run = run->r_branch; 
       }
   }
  prev->cover = run;

  run = next->cover;
  /* as long as there is a left branch at all, that is, the current cover of the pair next is
     no longer existent, since its right endnode has been exchanged by some other */
  while(run->l_branch != NULL)
    {


/*if ( (run->r_branch != NULL) && (run->l_branch != NULL) )
  {
    printf("run = (%d,%d),\n",run->left_node, run->right_node);
      printf("l_brach->id = %d, r_branch->id = %d\n",run->l_branch->id, run->r_branch->id);
      printf("l_branch = (%d, %d)\n",run->l_branch->left_node, run->l_branch->right_node);
      printf("r_branch = (%d, %d)\n",run->r_branch->left_node, run->r_branch->right_node);
  }
else
  {
    printf("Hier ist der Fehler: run = (%d,%d)\n",run->left_node, run->right_node);
    exit(0);
  }*/

      if ( (ABC(p_field[next->path_number].left_node, run->l_branch->left_node, next->arc_number)) &&
	   (ABC(next->arc_number+1, run->l_branch->right_node, p_field[next->path_number].right_node)) )
	{
	  /* in this case we take the left branch of the tree of coveredges */
	 run = run->l_branch;
	}
      else
       {
	 /* otherwise we take the right branch of the tree of coveredges */
	 run = run->r_branch; 
       }
    }

  next->cover = run;
  
  /* if the cover of the intersection covers the considered pair, then set the cover id */
  if ( (ABC(p_field[pair->path_number].left_node, prev->cover->left_node, pair->arc_number)) &&
       (ABC(pair->arc_number+1, prev->cover->right_node, p_field[pair->path_number].right_node)) )
    {
      pair->cover = prev->cover;
    }
  else
    /* if the cover of the union covers the considered pair, then set the cover id */
    if ( (ABC(p_field[pair->path_number].left_node, next->cover->left_node, pair->arc_number)) &&
	 (ABC(pair->arc_number+1, next->cover->right_node, p_field[pair->path_number].right_node)) )
      {
	pair->cover = next->cover;
      }

  /*if (pair->cover == NULL)
  {
    printf("Exchanging the right endnodes of (%d,%d) and (%d,%d)\n", prev->cover->left_node,
	 prev->cover->right_node, next->cover->left_node, next->cover->right_node);

    if ( (pair->path_number == 0) && (pair->arc_number == 5) )
      printf("prev->id = %d , next->id = %d\n",prev->id, next->id);

  }
else
  printf("the pair (%d, %d)'s cover = (%d, %d)\n",pair->path_number, pair->arc_number, 
	 pair->cover->left_node, pair->cover->right_node);*/



  if (pair->cover == NULL)
    return(0);
  else
    return(1);
}

/***********************************************************************************/

/* Executes phase 3 of the algorithm, that is, if a forbidden ess. pair is not 
   covered, then this will be corrected. */ 
void Phase_3 ( p_field, E_R_field, chain_part, Gyori_size, N)
path *p_field;
e_pair **E_R_field;
chain_tree *chain_part;
int Gyori_size;
forbidden *N;

{
  int corr_count= 0, help, i;
  path *help_chain;
  chain_tree *prev_cov, *next_cov;


  help_chain = (path*) malloc(Gyori_size * sizeof(struct Path));
  /* help_chain is an up to date storage of the currently existing coveredges */

  for (i=0; i<Gyori_size; i++)
    {
      help_chain[i].left_node = chain_part[i].left_node;
      help_chain[i].right_node = chain_part[i].right_node;

      /*printf("chain[%d] = (%d, %d)\n",chain_part[i].id, chain_part[i].left_node, chain_part[i].right_node);*/

    }
 

  while ( N != NULL )
    { 
     
      /* printf("\nNo covering chain assigned yet to pair: (%d, %d)!!!\n", N->forb->path_number, N->forb->arc_number); */



      /*if (N->forb->arc_number == 1764)
	printf("N = (%d, %d)!!!\n", N->forb->path_number, N->forb->arc_number);*/

      if ( (N->forb->cover == NULL) && (!Covered(N->forb, chain_part, help_chain, Gyori_size, p_field)) )
	{
	  /*printf("\nThis pair is really uncovered: (%d, %d)!!!\n", N->forb->path_number, N->forb->arc_number);*/
	  /* doing the exchange step */
	  corr_count++;  /* counts the number of correction-loops */
	  prev_cov = N->forb->prev->cover;
	  next_cov = N->forb->next->cover;


	  /*printf("prev_cov->id = %d, next_cov->id = %d\n",prev_cov->id, next_cov->id);*/



          /* exchanging the right endnodes of the corresponding coveredges */
	  help = help_chain[prev_cov->id].right_node; 
	  help_chain[prev_cov->id].right_node = help_chain[next_cov->id].right_node;
	  help_chain[next_cov->id].right_node = help;



	  /* constructing the tree of exchanged coveredges */
	  prev_cov->l_branch = (chain_tree*) malloc (sizeof(struct Chain_Tree));
	  prev_cov->l_branch->l_branch = NULL;
	  prev_cov->l_branch->r_branch = NULL;
	  prev_cov->l_branch->id = prev_cov->id;
	  prev_cov->l_branch->left_node = help_chain[prev_cov->id].left_node;
	

	  /*for (i=0; i<Gyori_size; i++)
	    {	    
	      printf("help_chain[%d] = (%d, %d)\n",i, help_chain[i].left_node, help_chain[i].right_node);
	    }
	    printf("\n");*/

	  /* the folowing right endnode of the coveredge has been changed already */
	  prev_cov->l_branch->right_node = help_chain[prev_cov->id].right_node;
	  next_cov->l_branch = prev_cov->l_branch;
	


	  prev_cov->r_branch = (chain_tree*) malloc (sizeof(struct Chain_Tree));
	  prev_cov->r_branch->l_branch = NULL;
	  prev_cov->r_branch->r_branch = NULL;
	  prev_cov->r_branch->id = next_cov->id;
	  prev_cov->r_branch->left_node = help_chain[next_cov->id].left_node;
	  /* the folowing right endnode of the coveredge has been changed already */
	  prev_cov->r_branch->right_node = help_chain[next_cov->id].right_node;
	  next_cov->r_branch = prev_cov->r_branch;



	  /* both branches cover the currently considered forbidden pair */
	  /* we always take the coveredge in the right branch */
	  N->forb->cover = next_cov->r_branch;

	 
	}
	
      N = N->next;
    }

  /* Now we restore the chain part by copiing the actual set of coveredges to it, 
     which is stored in help_chain */
  for (i=0; i<Gyori_size; i++)
    {
      chain_part[i].left_node = help_chain[i].left_node;
      chain_part[i].right_node = help_chain[i].right_node;
    }



  printf("Number of Corrections: %d\n",corr_count);
  
}

/***********************************************************************************/

/* Writes the gen. system into a file called Generator */ 
void Write_Generating_System(chain_part, Gyori_size, red_unred)
chain_tree *chain_part;
int Gyori_size;
int *red_unred;

{
  int i;
  FILE *fp;
  char file[]="Generator";

  fp = fopen(file,"w");
  if (fp == NULL) 
    {
      printf("\nWrite_Generating_System: file %s can't be opened\n",file);
      exit(0);
    }

  fprintf(fp,"Size: %d\n", Gyori_size);
 
  /* writing the unreduced values using the red_unred field constructed in function Reduce() */
  for(i=0; i<Gyori_size; i++)
    {
      fprintf(fp, "(%d %d)\n", red_unred[chain_part[i].left_node], red_unred[chain_part[i].right_node]);
      /*fprintf(fp, "(%d %d)\n", chain_part[i].left_node, chain_part[i].right_node);*/

    }
  fclose(fp);
}



/***********************************************************************************/

/* tests the generating system, only for debugging reasons */ 
int Test_If_Generates(p_field, circ_length, subfam_size, chain_part, Gyori_size)
path *p_field;
int circ_length;
int subfam_size;
chain_tree *chain_part;
int Gyori_size;

{
  int i, j, help_right,  counter=0;
 

  printf("Testing the Generator...\n");
  for(i=0; i<subfam_size; i++)
    {
 
      j=0;
      while (chain_part[j].left_node < p_field[i].left_node ) 
	j++;

      help_right = p_field[i].left_node;

      counter = 0;
      while ( (ABC(p_field[i].left_node, chain_part[j].left_node, p_field[i].right_node-1)) &&
	       (counter < Gyori_size) )
	{
	  /*printf("chain[%d]: (%d,%d), p_field[%d] = (%d,%d)\n",j, chain_part[j].left_node, 
	    chain_part[j].right_node, i, p_field[i].left_node, p_field[i].right_node); */

	  if ( ( ABC(chain_part[j].left_node, help_right, chain_part[j].right_node ) ) &&
	       ( ABC(help_right, chain_part[j].right_node, p_field[i].right_node) ) )
	    {
	      help_right =  chain_part[j].right_node;
	    }
	  else
	    {
	      if ( ABC(p_field[i].left_node, help_right, chain_part[j].left_node-1) )
		{
		  printf("\nArc %d in path (%d, %d) Nr. %d not generated!!!\n", 
			 help_right, p_field[i].left_node, p_field[i].right_node, i);



		  printf("chain[%d]: (%d,%d), p_field[%d] = (%d,%d)\n",j, chain_part[j].left_node, 
			chain_part[j].right_node, i, p_field[i].left_node, p_field[i].right_node);
		  printf("chain[%d]: (%d,%d)\n",j-1, chain_part[j-1].left_node, 
			 chain_part[j-1].right_node);
		 
		  exit(0);
		}
	    }
	  j = Modulo(j+1, Gyori_size);
	  counter++;
	}

    }

  printf("\nGenerator-Test successful!\n\n");
  return(1);
      
}



/***********************************************************************************/
/***********************************************************************************/

void main()
{ 
  struct Path *p_field;
  int circ_length;
  int i,subfam_size, E_R_size=0, Gyori_size;
  e_pair **E_R;
  bip_graph *G;
  e_pair **E_R_field, **PER;
  chain_tree *chain_part;
  forbidden *N;
  int *red_unred;
  int *R;
  char c;
  double_int *p_ref;
  struct rusage *usage0, *usage1, *usage2;
  double utime1, stime1, utime2, stime2, utime3, stime3, utime4, stime4 ;

  /* to mesure the run-time */
  usage0 = (struct rusage*) malloc(sizeof(struct rusage));
  usage1 = (struct rusage*) malloc(sizeof(struct rusage));
  usage2 = (struct rusage*) malloc(sizeof(struct rusage));

  
  getrusage(RUSAGE_SELF, usage0);

  p_field = Read_Input( &circ_length, &subfam_size, &c); 
  printf("\nNo. of distinct subpaths = %d.  ", subfam_size);
  printf("length of the circuit: %d\n",circ_length);

  if ( (subfam_size <= 30) && (circ_length <= 30) )
    {
      Screen_Subfam(circ_length, subfam_size, &p_field[0]); /* simple ascii output */
    }
  


  i=0;
  /*Sorting the family of subpaths */
  My_Qsort(p_field, i, subfam_size-1, circ_length);
  if (c == 'y')
    for(i=0; i<subfam_size-1; i++)     
      {
	/* making sure that the mebers of the freebie generator get in
	   the right position to mark all covered elements as NOT f-essential */
	if ( (p_field[i].id == -1) && 
	     (p_field[i].left_node == p_field[i+1].left_node) &&
	     (p_field[i].right_node == p_field[i+1].right_node) ) 
	  My_Swap(p_field,i,i+1);               
      }

  /*if ( (subfam_size <= 30) && (circ_length <= 30) )
    {
      Screen_Subfam(circ_length, subfam_size, &p_field[0]);  simple ascii output 
    }*/
  
  printf("\nDeleting redundant edges...\n"); 
  red_unred = Reduce(p_field, subfam_size, &circ_length);
  printf("Reduction ready.\n");

  if ( (subfam_size <= 30) && (circ_length <= 30) )
    {
      Screen_Subfam(circ_length, subfam_size, &p_field[0]); 
    }
  
  p_ref = Reference_P_Field(p_field, subfam_size, circ_length);
  printf("\nConstructing the set of essential pairs...\n");
  getrusage(RUSAGE_SELF, usage1);
  E_R = Construct_E_R(p_field, subfam_size, &E_R_size, p_ref, circ_length);

 
  getrusage(RUSAGE_SELF, usage2);

  utime1 =  (double) usage2->ru_utime.tv_sec +  0.000001 * usage2->ru_utime.tv_usec ;
  utime1 -= (double) usage1->ru_utime.tv_sec +  0.000001 * usage1->ru_utime.tv_usec ;

  stime1 =  (double) usage2->ru_stime.tv_sec +  0.000001 * usage2->ru_stime.tv_usec ;
  stime1 -= (double) usage1->ru_stime.tv_sec +  0.000001 * usage1->ru_stime.tv_usec ;
  /*printf("Usertime: %lf sec , Systemtime: %lf sec\n",utime1, stime1);*/
  printf("Construction ready.\n");

  /*if (c == 'y')
    Remove_Free_Gen(p_field, &subfam_size);*/

  if ( (subfam_size <= 30) && (circ_length <= 30) )
    {
      printf("\nMarking the essential ( ==> ) arcs:\n");
      Screen_E_R2(circ_length, p_field, subfam_size, E_R);
    }
  printf("%d essential pairs were found.\n", E_R_size);


  printf("\nOutcrossing...\n");
                                          /* initializing */
  N = Crossfree(p_field, E_R, &E_R_size, circ_length);

  printf("Outcrossing ready.\n");

  /*PER is a field of pointers, where PER[i] points to the 
    first member K on the 'i'th subpath*/ 
  PER = (e_pair**) malloc (subfam_size * sizeof(e_pair*) );   /* allocating memory */
  for (i=0; i<subfam_size; i++)
    PER[i] = NULL;  
  Construct_PER( E_R, circ_length, PER );

  printf("\nConstructing the reference-field for E_R\n");
  E_R_field = Reference_E_R(E_R, circ_length, E_R_size);

  if ( (subfam_size <= 30) && (circ_length <= 30) )
    {
      printf("\nRemaining essential ( ==> ) arcs:\n");
      Screen_E_R(circ_length, p_field, subfam_size, E_R);
    }


  printf("\nConstructing bipartite graph...\n");
  getrusage(RUSAGE_SELF, usage1);
  G = Compare( p_field, E_R, E_R_size, PER, subfam_size, circ_length);

  getrusage(RUSAGE_SELF, usage2);

  utime2 =  (double) usage2->ru_utime.tv_sec +  0.000001 * usage2->ru_utime.tv_usec ;
  utime2 -= (double) usage1->ru_utime.tv_sec +  0.000001 * usage1->ru_utime.tv_usec ;

  stime2 =  (double) usage2->ru_stime.tv_sec +  0.000001 * usage2->ru_stime.tv_usec ;
  stime2 -= (double) usage1->ru_stime.tv_sec +  0.000001 * usage1->ru_stime.tv_usec ;
  /*printf("Usertime: %lf sec , Systemtime: %lf sec\n",utime2, stime2);*/
  printf("Bipartite graph ready.\n"); 

  printf("\nNumber of nodes in the bip. graph = %d.\n",G->n_nodes);
  printf("Number of edges in the bip. graph = %d.\n",G->n_edges);

  /*if ( ( (main_path.right_node - main_path.left_node) <= 320) &&
       (subfam_size <= 250) ) 
    {
      do
	{
	  while(getchar() != '\n') ;
	  printf("\nShow Interval system graphically? (y/n):");
	}
      while ( ((c1=getchar()) != 'y') && (c1 != 'n') );
      if (c1=='y')
	{  
	  Show_E_R(p_field, subfam_size, main_path, PER, N, G, E_R_field);
	}
    }*/
 
  Write_Bip_Graph(G); /*if wanted*/ 
  printf("\nMatching...\n");
  getrusage(RUSAGE_SELF, usage1);
  R = Bip_Match(G);
  getrusage(RUSAGE_SELF, usage2);

  utime3 =  (double) usage2->ru_utime.tv_sec +  0.000001 * usage2->ru_utime.tv_usec ;
  utime3 -= (double) usage1->ru_utime.tv_sec +  0.000001 * usage1->ru_utime.tv_usec ;

  stime3 =  (double) usage2->ru_stime.tv_sec +  0.000001 * usage2->ru_stime.tv_usec ;
  stime3 -= (double) usage1->ru_stime.tv_sec +  0.000001 * usage1->ru_stime.tv_usec ;
  /*printf("Usertime: %lf sec , Systemtime: %lf sec\n",utime3, stime3);*/
  printf("Matching completed.\n");

  printf("\nConstructing the maximum independent subsystem...\n");  
  Gyori_size = Max_Indep_Subs (p_field, E_R_field, G, red_unred, R);
  printf("Ready.\n");

  chain_part = Min_Chain_Decomp(p_field, circ_length, E_R_field, E_R_size, G, Gyori_size);
  printf ("\nsigma(F) = gamma(F) = %d\n", Gyori_size);




  /*run = E_R[1764];
  while (run != NULL)
    {
      printf("Pair: (%d, %d), id = %d, cover = (%d, %d)\n", run->path_number, run->arc_number,
	     run->id, chain_part[run->cover].left_node, chain_part[run->cover].right_node);
      run = run->next;
    }*/




  printf ("\nCorrecting the chain partition...\n");
  Phase_3 ( p_field, E_R_field, chain_part, Gyori_size, N);
  printf ("Ready.\n");

  printf("Writing the gen. system.\n"); 
  Write_Generating_System(chain_part, Gyori_size, red_unred);
  printf("Written.\n");

  getrusage(RUSAGE_SELF, usage2);

  utime4 =  (double) usage2->ru_utime.tv_sec +  0.000001 * usage2->ru_utime.tv_usec ;
  utime4 -= (double) usage0->ru_utime.tv_sec +  0.000001 * usage0->ru_utime.tv_usec ;

  stime4 =  (double) usage2->ru_stime.tv_sec +  0.000001 * usage2->ru_stime.tv_usec ;
  stime4 -= (double) usage0->ru_stime.tv_sec +  0.000001 * usage0->ru_stime.tv_usec ;

  printf("\nSWAPS: %ld\n", usage2->ru_nswap);
  printf("Entire              Usertime: %8.4f sec , Systemtime: %8.4f sec\n",utime4, stime4);

  printf("\nPERCENT USAGE:\n");
  printf("Constr_E_R        : Usertime:  %7.4f %%   , Systemtime:  %7.4f %%\n", 100 * (double) utime1 / utime4,
	 0.01 * (double) stime1 / stime4); 
  printf("Constr.Bip.Graph  : Usertime:  %7.4f %%   , Systemtime:  %7.4f %%\n", 100 * (double) utime2 / utime4,
	 0.01 * (double) stime2 / stime4); 
  printf("Bip.Matching      : Usertime:  %7.4f %%   , Systemtime:  %7.4f %%\n\n", 100 * (double) utime3 / utime4,
	 0.01 * (double) stime3 / stime4);  

  Test_If_Generates(p_field, circ_length,  subfam_size, chain_part, Gyori_size);  /* for debugging reasons only */

}