/*
  The following source is an approximation of the code
  contained in Appendix A of

    H. Le Verge.
    A note on Chernikova's algorithm.
    Publication interne 635,
    IRISA, Campus de Beaulieu,
    Rennes, France, 1992.

  This source has been obtained by

    1) converting Le Verge's listings to electronic form with a scanner;
    2) converting the images into an approximation of the text with
       an Optical Character Recognition (OCR) program;
    3) performing some manual corrections.

  All these steps are error-prone and no guarantee of accuracy
  can thus be provided.  USE AT YOUR OWN RISK!!!

  Here is an approximation of the constraint system that is used
  in Le Verge's paper to illustrate the use of the program:

    5 5
    1  -2   2  -1   0
    1   4   2  -1   0
    1   2   6  -3 -20
    1  -6  10  -5 -24
    1   0   0   0   1
*/

#include <malloc.h>
#include <stdio.h>
#include <stdlib.h>

#define EXCH(mat, l1, l2, aux) \
aux = mat.p[l1]; \
mat.p[l1] = mat.p[l2]; \
mat.p[l2] = aux;

/* #define NULL 0 */
#define TOP 2147483647
#define true 1
#define false 0

/* data structures */
typedef struct vector {
  int size;
  int *p;
} vector;

typedef struct matrix {
  int nbrows;
  int nbcolumns;
  int **p;
  int *p_init;
} matrix;

char *My_Alloc(int nbelement, int element_size)
{
  char *p;
  if ((p = malloc((unsigned)(nbelement*element_size))) == NULL) {
    printf("out of memory error\n");
    exit(1);
  }
  return p;
} /* My_Alloc */

void Alloc_Vector(vector *vec, int size)
{
  vec->size = size;
  vec->p = (int*) My_Alloc(size, sizeof(int));
} /* Alloc_Vector */

void Alloc_Matrix(matrix *mat, int nbrows, int nbcolumns)
{
  int i, *p, **q;
  mat->nbrows = nbrows;
  mat->nbcolumns = nbcolumns;
  q = (int**) My_Alloc(nbrows, sizeof(int *));
  mat->p = q;
  p = (int*) My_Alloc(nbrows*nbcolumns, sizeof(int));
  mat->p_init = p;
  for (i = 0; i<nbrows; i++) {
    mat->p[i] = p;
    p = p+nbcolumns;
  }
} /* Alloc_Matrix */

void Free_Vector(vector *vec)
{
  free((char *)vec->p);
} /* Free_Vector */

void Free_Matrix(matrix *mat)
{
  free((char *)mat->p);
  free((char *)mat->p_init);
} /* Free-Matrix */

void Read_Matrix(matrix *mat)
{
  int nbrows, nbcolumns, i, j;

  (void) scanf("%d %d", &nbrows, &nbcolumns);
  Alloc_Matrix(mat, nbrows, nbcolumns);
  for (i = 0; i < nbrows; i++)
    for (j = 0; j < nbcolumns; j++)
      (void) scanf("%d", &(mat->p[i][j]));
} /* Read-Matrix */

void Write_Matrix(matrix *mat, int dim)
{
  int nbrows, i, j;

  (void) printf ("%d %d\n", nbrows = mat->nbrows, dim);
  for (i = 0; i < nbrows; i++) {
    for (j = 0; j < dim; j++)
      (void) printf("%4d ", mat->p[i][j]);
    (void) printf("\n");
  }
} /* Write-Matrix */

void Gcd(int a, int b, int *g)
{
  int aux;
  if ((a==b) && (b==0))
    *g=0;
  else {
    a = abs(a);
    b = abs(b);
    while (a) {
      aux = b % a;
      b = a;
      a = aux;
    }
    *g=b;
  }
} /* Gcd */

/* finds the minimum (absolute value) of vector vec */
void Vector_Min_Not_0(vector *vec, int *min, int *index)
{
  int m, size, i, ind, aux;

  size = vec->size;
  m = TOP;
  ind = -1;
  for (i = 0; i < size; i++)
    if (vec->p[i] != 0)
      if (m > (aux = abs(vec->p[i]))) {
	ind = i;
	m = aux;
      }
  if (ind == -1)
    *min = 1;
  else
    *min = m;
  *index = ind;
} /* Vector_Min_Not_0 */

/* computes the gcd of all components of vec */
void Vector_Gcd(int *vec, int size, int *g)
{
  vector vec_A;
  int min, index, *p, not_all_zero, i;

  Alloc_Vector(&vec_A, size);
  p = vec;
  for (i = 0; i < size; i++)
    vec_A.p[i] = *p++;
  do {
    Vector_Min_Not_0(&vec_A, &min, &index);
    if (min == 1)
      break;
    not_all_zero = false;
    for (i = 0; i < size; i++)
      if (i != index)
	not_all_zero |= (vec_A.p[i] %= min);
  } while (not_all_zero);
  *g = min;
  Free_Vector(&vec_A);
} /* Vector_Gcd */

/* normalize a vector in such a way that gcd(vec) = 1 */
void Vector_Normalize(int *vec, int size, int dimension)
{
  int gcd, i, *p;
  Vector_Gcd(vec, dimension, &gcd);
  if (gcd >= 2) {
    p = vec;
    for (i = 0; i < size; i++)
      *p++ /= gcd;
  }
} /* Vector_Normalize */

/* computes r3 as combination of r1 and r2 in such a way that r3[k] = 0 */
void combine(int *r1, int *r2, int *r3, int k, int nbcolumns, int dim)
{
  int aux, a1, a2, j;

  Gcd(r1[k], r2[k], &aux);
  a1 = *(r1+k) / aux;
  a2 = *(r2+k) / aux;
  for (j = 0; j < nbcolumns; j++)
    r3[j] = a2*r1[j] - a1*r2[j];
  Vector_Normalize(r3, nbcolumns, dim);
} /* combine */

int chernikova(matrix *bid, matrix *uni,
	       /* matrices of bidirectional and unidirectional rays */
	       int nbbidray, int nbray,
	       /* number of bidirectional and unidirectional rays */
	       vector *ineq) /* vector of flags inequality/equation */
{
  vector commonzero, inequality;
  matrix ray, bidray;

  int nbcolumns, index_non_zero, nbcommonconstraints, max_rays;
  int i, j, k, l, m, *p, bound, inf_bound, equal_bound, sup_bound, redundant;
  int dimension; /* dimension of the space */

  bidray = *bid;
  ray = *uni;
  inequality = *ineq;
  nbcolumns = bidray.nbcolumns; /* width of both matrices uni and bid */
  max_rays = ray.nbrows;	/* size of matrix uni */
  dimension = nbbidray;
  Alloc_Vector(&commonzero, nbcolumns);
  for (k = dimension; k < nbcolumns; k++) {
    /* finds if exists, a bidirectional ray which does not saturate the
       current constraint */
    index_non_zero = -1;
    for (i = 0; i < nbbidray; i++)
      if (bidray.p[i][k] != 0) {
	index_non_zero = i;
	break;
      }
    if (index_non_zero != -1) {
      /* discards index_non_zero bidirectional ray */
      nbbidray--;
      if (nbbidray != index_non_zero) {
	p = bidray.p[nbbidray];
	bidray.p[nbbidray] = bidray.p[index_non_zero];
	bidray.p[index_non_zero] = p;
      }
      /* Computes the new lineality space */
      for (i = 0; i < nbbidray;i++)
	if (bidray.p[i][k] != 0)
	  combine (bidray.p[i], bidray.p[nbbidray], bidray.p[i], k,
		   nbcolumns, dimension);
      /* add the positive part of index-non-zero bidirectional ray to
	 the set of unidirectional rays */
      if (nbray == max_rays)
	return 1;
      if (bidray.p[nbbidray][k] < 0)
	for (j = 0; j < nbcolumns; j++)
	  ray.p[nbray][j] = -bidray.p[nbbidray][j];
      else
	for (j = 0; j < nbcolumns; j++)
	  ray.p[nbray][j] = bidray.p[nbbidray][j];
      /* Computes the new pointed cone */
      for (i = 0; i < nbray; i++)
	if (ray.p[i][k] != 0)
	  combine(ray.p[i], ray.p[nbray], ray.p[i], k, nbcolumns, dimension);
      if (inequality.p[k])
	nbray++;
    }
    else {
      /* sort rays : 0 <= i < equal_bound : saturates the constraint
                   : equal-bound <= i < sup_bound : verifies the constraint
                   : sup_bound <= i < bound : does not verify */
      equal_bound = 0;
      sup_bound = 0;
      inf_bound = nbray;
      while (inf_bound > sup_bound) {
	if (ray.p[sup_bound][k] == 0) {
	  EXCH(ray, equal_bound, sup_bound, p);
	  equal_bound++;
	  sup_bound++;
	}
	else if (ray.p[sup_bound][k] < 0) {
	  inf_bound--;
	  EXCH(ray, inf_bound, sup_bound, p);
	}
	else
	  sup_bound++;
      }
      /* Computes only the new pointed cone */
      bound = nbray;
      for (i = equal_bound; i < sup_bound; i++)
	for(j = sup_bound; j < bound; j++) {
	  /* computes the set of common saturated constraints */
	  nbcommonconstraints = 0;
	  for (l = dimension; l < k; l++)
	    if ((ray.p[i][l] == 0) && (ray.p[j][l]== 0)) {
	      commonzero.p[nbcommonconstraints] = l;
	      nbcommonconstraints++;
	    }
	  if (nbcommonconstraints+nbbidray >= dimension-2) {
	    /* check whether a ray m saturates the same set of constraints */
	    redundant = false;
	    for (m = 0; m < bound; m++)
	      if ((m != i) && (m != j)) {
		for (l = 0; l < nbcommonconstraints; l++)
		  if (ray.p[m][commonzero.p[l]] != 0)
		    break;
		if (l == nbcommonconstraints) {
		  /* the combination of ray i and j will generate a
		     non extremal ray so ... */
		  redundant = true;
		  break;
		}
	      }
	    if (!redundant) {
	      if (nbray == max_rays)
		return 1;
	      /* computes the new ray */
	      combine(ray.p[j], ray.p[i], ray.p[nbray], k, nbcolumns,
		      dimension);
	      nbray++;
	    }
	  }
	}
      /* Eliminates all non extremal rays */
      if (inequality.p[k])
	j = sup_bound;
      else
	j = equal_bound;
      i = nbray;
      while ((j < bound) && (i > bound)) {
	i--;
	EXCH(ray, i, j, p);
	j++;
      }
      if (j == bound)
	nbray = i;
      else
	nbray = j;
    }
  }
  Free_Vector(&commonzero);

  bid->nbrows = nbbidray;
  uni->nbrows = nbray;
  return 0;
} /* chernikova */

int main(int argc, char **argv)
{
  matrix mat, ray, bidray;
  vector inequality; /* vector of flags equation/inequality */
  int nbconstraints, dimension, nbcolumns, nbcolumns2, i, j;

  Read_Matrix(&mat);

  nbconstraints = mat.nbrows;
  nbcolumns = mat.nbcolumns;
  dimension = nbcolumns - 1;
  nbcolumns2 = dimension + nbconstraints;

  Alloc_Matrix(&ray, atoi(argv[1]), nbcolumns2);
  Alloc_Matrix(&bidray, dimension, nbcolumns2);

  /* set the identity matrix as bidirectional rays */
  for (i=0;i<dimension;i++) {
    for (j = 0; j < dimension; j++)
      bidray.p[i][j]=0;
    bidray.p[i][i]=1;
  }

  for (i = 0; i < nbconstraints; i++)
    for (j = 1; j < nbcolumns; j++)
      bidray.p[j-1][i+dimension] = mat.p[i][j];

  Alloc_Vector(&inequality, nbcolumns2);
  for (i = 0; i < nbconstraints; i++)
    inequality.p[i+dimension] = mat.p[i][0];

  if (chernikova(&bidray, &ray, dimension, 0, &inequality) != 0) {
    printf("not enough rays available\n");
    exit(2);
  }

  Free_Vector(&inequality);

  printf ("bidirectional rays\n");
  Write_Matrix(&bidray, mat.nbcolumns-1);
  printf("unidirectional rays\n");
  Write_Matrix(&ray, mat.nbcolumns-1);

  Free_Matrix(&bidray);
  Free_Matrix(&ray);

  return 0;
} /* main */