/* An Implementation of the Enumeration of Regular Triangulations */
#include <stdio.h>
#include <fcntl.h>
#include <time.h>
#include "gnu_mp.h"


/*#define DEBUG*/

/*#define EXPERIMENT*/


/* for the whole program */ 
#define INDEX int
DTYPE BZERO;
int fd;
int fd_;
int fd__; /* for restart */
int stop_watch;
int sumof_time;

void bprintf(c)
MP_RAT *c;
{
  MP_INT num, den;

  mpz_init(&num);
  mpz_init(&den);

  mpq_get_num(&num, c);
  mpq_get_den(&den, c);

  mpz_out_str(stdout, 10, &num);
  printf("/");
  mpz_out_str(stdout, 10, &den);

  mpz_clear(&den);
  mpz_clear(&num);
}

void still()
{
/*  while (getchar() == EOF);*/
}

void watch_on()
{
  stop_watch = clock();
}

void watch_see()
{
  int t;
  char msg[BUFSIZ];
  write(fd, "[time]: ", 8);
  t = clock() - stop_watch;
  sprintf(msg, "%16d\n", (long) t);
  write(fd, msg, strlen(msg));
  free(msg);
  sumof_time += t;
}



/* for Gaussian elimination */
#define ELIM_ERROR    -1
#define ELIM_TERMINATE 0

DTYPE    *elim_matrix;
int    elim_numof_row;
int elim_numof_column;
INDEX      *row, *col;



/* for linear programming */
#define LP_CONTINUE   0
#define LP_UNBOUNDED  1
#define LP_TERMINATE -1

DTYPE     *tableau; /* 0th row represents the objective function. */
int sizeof_tableau;
int           rank;
INDEX       *basic;
INDEX    *nonbasic;
DTYPE    *equation;


/* for the main part */
DTYPE *matrix;
DTYPE *gale;

DTYPE *initial_mtrx;
DTYPE *initial_glmtrx;
DTYPE *volume_mtrx;
DTYPE *normal_mtrx;
DTYPE *normal_glmtrx;
DTYPE *signed_radon_mtrx;
DTYPE *signed_radon_glmtrx;
DTYPE *set1of_normals;
DTYPE *set2of_normals;
DTYPE *volvec1;
DTYPE *volvec2;
DTYPE *vector1;
DTYPE *vector2;
DTYPE *vector3;
DTYPE *vector4;

INDEX *global_Gamma;
INDEX *global_Gamma_;
INDEX *global_Delta_;

int    numof_rows;
int numof_columns;
int     dimension; 
int   codimension;
int       sizetri;
int      sizecand;

int     numof_reg;

#define Simplex INDEX*
#define Radon   INDEX*

#ifndef BUFSIZ
#define BUFSIZ 1024
#endif

#define SMPLXNUM 30

#define Tri    INDEX*
#define Cand   INDEX*
#define Normal DTYPE*




void print_matrix(mtrx, nr, nc)
DTYPE *mtrx;
int nr, nc;
{
  int i, j;
  for (i = 0; i < nr; i++) {
    for (j = 0; j < nc; j++) {
      bprintf(&mtrx[i * nc + j]);
      printf("\t");
    }
    printf("\n");
  }
}



int nonzero_in_column(c, d)
int c, d;
{
  int i;
  for (i = d; i < elim_numof_row; i++)
    if (bcmp(&elim_matrix[row[i] * elim_numof_column + col[c]], &BZERO) != 0)
      return i;
  return ELIM_ERROR;
}

void row_mul(factor, r)
DTYPE *factor;
int r;
{
  int i;
  DTYPE temp;
  bcreate(&temp);
  for (i = 0; i < elim_numof_column; i++) {
    bmul(&temp, &elim_matrix[row[r] * elim_numof_column + col[i]], factor); 
    bass(&elim_matrix[row[r] * elim_numof_column + col[i]], &temp);
  }
  bdestroy(&temp);
}

void row_sub(r1, r2)
int r1, r2;
{
  int i;
  DTYPE temp;
  bcreate(&temp);
  for (i = 0; i < elim_numof_column; i++) {
    bsub(&temp, &elim_matrix[row[r1] * elim_numof_column + col[i]], 
	 &elim_matrix[row[r2] * elim_numof_column + col[i]]);
    bass(&elim_matrix[row[r1] * elim_numof_column + col[i]], &temp);
  }
  bdestroy(&temp);
}

int search_nonzero_element(c)
int c;
{
  int r, c0;
  INDEX temp;
  c0 = c;
  r = nonzero_in_column(c, c0);
  while (c < elim_numof_column && (r = nonzero_in_column(c, c0)) == ELIM_ERROR)
    c++;
  if (c == elim_numof_column)
    return ELIM_ERROR;
  temp = row[c0]; row[c0] = row[r]; row[r] = temp;
  temp = col[c0]; col[c0] = col[c]; col[c] = temp;  
  return ELIM_TERMINATE;
}

int elim_column(c)
int c;
{
  int i;
  DTYPE temp;
  if (c >= elim_numof_column)
    return ELIM_ERROR;
  if (search_nonzero_element(c) != ELIM_ERROR) {
    bcreate(&temp);
    for (i = c + 1; i < elim_numof_row; i++) {
      if (bcmp(&elim_matrix[row[i] * elim_numof_column + col[c]], &BZERO) != 0) {
	 bdiv(&temp, &elim_matrix[row[c] * elim_numof_column + col[c]], 
	      &elim_matrix[row[i] * elim_numof_column + col[c]]);
	 row_mul(&temp, i);
	 row_sub(i, c);
      }
    }
    bdestroy(&temp);
    return ELIM_TERMINATE;
  } else 
    return ELIM_ERROR; /* matrix denotes a zero matrix. */
}

int elim_column_backward(c)
int c;
{
  int i;
  DTYPE temp;
  bcreate(&temp);
  binv(&temp, &elim_matrix[row[c] * elim_numof_column + col[c]]);
  row_mul(&temp, c);
  for (i = c - 1; i >= 0; i--) {
    if (bcmp(&elim_matrix[row[i] * elim_numof_column + col[c]], &BZERO) != 0) {
      binv(&temp, &elim_matrix[row[i] * elim_numof_column + col[c]]);
      row_mul(&temp, i);
      row_sub(i, c);
    }
  }
  bdestroy(&temp);
  return ELIM_TERMINATE;
}

int elim_column_for_volume(c, factor)
int c;
DTYPE *factor;
{
  int i;
  DTYPE temp, temp2;

  bint(factor, 1);
  if (c >= elim_numof_column)
    return ELIM_ERROR;
  if (search_nonzero_element(c) != ELIM_ERROR) {
    bcreate(&temp);
    bcreate(&temp2);
    for (i = c + 1; i < elim_numof_row; i++) {
      if (bcmp(&elim_matrix[row[i] * elim_numof_column + col[c]], &BZERO) != 0) {
	 bdiv(&temp, &elim_matrix[row[c] * elim_numof_column + col[c]], 
	      &elim_matrix[row[i] * elim_numof_column + col[c]]);
	 row_mul(&temp, i);
	 row_sub(i, c);
	 bmul(&temp2, factor, &temp);
	 bass(factor, &temp2);
      }
    }
    bdestroy(&temp);
    bdestroy(&temp2);
    return ELIM_TERMINATE;
  } else
    return ELIM_ERROR; /* matrix denotes a zero matrix. */
}

void volume(answer)
DTYPE *answer;
{
  int i, dim, flag;
  DTYPE factor, temp1, temp2;

  elim_matrix     = volume_mtrx;
  elim_numof_row    = dimension;
  elim_numof_column = dimension;

  row = (INDEX*) malloc(elim_numof_row    * sizeof(INDEX));
  col = (INDEX*) malloc(elim_numof_column * sizeof(INDEX));
  for (i = 0; i < elim_numof_row; i++)    row[i] = (INDEX) i;
  for (i = 0; i < elim_numof_column; i++) col[i] = (INDEX) i;

  dim = elim_numof_row;
  bcreate(&factor);
  bcreate(&temp1);
  bcreate(&temp2);
  bint(&factor, 1);
  for (i = 0; i < elim_numof_row; i++) {
    flag = elim_column_for_volume(i, &temp1);
    bmul(&temp2, &factor, &temp1);
    bass(&factor, &temp2);
    if (flag == ELIM_ERROR) {
      dim = i;
      i = elim_numof_row;
    }
  }
  binv(&temp1, &factor);
  bass(&factor, &temp1);
  if (dim != elim_numof_row) {
    bint(answer, 0);
  } else {
    for (i = 0; i < elim_numof_row; i++) {
      bmul(&temp1, &factor, &elim_matrix[row[i] * elim_numof_column + col[i]]);
      bass(&factor, &temp1);
    }
    bass(answer, &factor);
  }
  if (bcmp(answer, &BZERO) < 0) {
    bneg(&temp1, answer);
    bass(answer, &temp1);
  }
  bdestroy(&factor);
  bdestroy(&temp1);
  bdestroy(&temp2);
  free(row);
  free(col);
}



int gale_transform(mtrx, nor, noc, gale_mtrx)
DTYPE      *mtrx;
int     nor, noc;
DTYPE *gale_mtrx;
{
  int i, j, codim;

  elim_matrix      = mtrx;
  elim_numof_row    = nor;
  elim_numof_column = noc;

  row = (INDEX *) malloc(elim_numof_row    * sizeof(INDEX));
  col = (INDEX *) malloc(elim_numof_column * sizeof(INDEX));
  for (i = 0; i < elim_numof_row;    i++) 
    row[i] = (INDEX) i;
  for (i = 0; i < elim_numof_column; i++) 
    col[i] = (INDEX) i;

  codim = elim_numof_column - elim_numof_row;
  for (i = 0; i < elim_numof_row; i++) 
    if (elim_column(i) == ELIM_ERROR) {
      codim = elim_numof_column - i;
      i = elim_numof_row;
    }
  for (i = elim_numof_column - codim - 1; i >= 0; i--) 
    if (elim_column_backward(i) == ELIM_ERROR)
      return 0;

  for (i = 0; i < codim; i++)
    for (j = 0; j < noc - codim; j++) {
      bcreate(&gale_mtrx[i * noc + col[j]]);
      bass(&gale_mtrx[i * noc + col[j]], 
	   &mtrx[row[j] * elim_numof_column + col[i + noc - codim]]);
    }
  for (i = 0; i < codim; i++)
    for (j = 0; j < codim; j++) {
      bcreate(&gale_mtrx[i * elim_numof_column + col[j + noc - codim]]);
      bint(&gale_mtrx[i * elim_numof_column + col[j + noc - codim]], -(i == j));
    }
  free(row);
  free(col);
  return codim;
}



int check(matrix1, matrix2, nor1, nor2, noc)
DTYPE *matrix1, *matrix2;
int nor1, nor2, noc;
{
  int i, j, k, retval;
  DTYPE temp1, temp2, temp3;

  retval = 1;
  bcreate(&temp1);
  bcreate(&temp2);
  bcreate(&temp3);
  for (i = 0; i < nor1; i++)
    for (j = 0; j < nor2; j++) {
      bint(&temp1, 0);
      for (k = 0; k < noc; k++) {
	bmul(&temp2, &matrix1[i * noc + k], &matrix2[j * noc + k]);
	badd(&temp3, &temp1, &temp2);
	bass(&temp1, &temp3);
      }
      if (bcmp(&temp1, &BZERO) != 0)
	 retval = 0;
    }
  bdestroy(&temp3);
  bdestroy(&temp2);
  bdestroy(&temp1);
  return retval;
}



void load(file_name, matrix_p, nor_p, noc_p)
char *file_name;
DTYPE **matrix_p;
int *nor_p, *noc_p;
{
  int i, n, fd, temp;
  int numof_data, cursor;
  char buf[BUFSIZ];

  if ((fd = open(file_name, O_RDONLY)) == -1) {
    fprintf(stderr, "Can't open %s\n", file_name);
    exit(1);
  }

  cursor = 0;
  if (read(fd, buf, BUFSIZ) <= 0) {
    fprintf(stderr, "Read error.\n");
    exit(1);
  }
  if (sscanf(buf, "%d %d%n", nor_p, noc_p, &n) <=0) {
    fprintf(stderr, "Illegal file.\n");
    exit(1);
  }

  *matrix_p = (DTYPE *) malloc(*nor_p * *noc_p * sizeof(DTYPE));
  cursor += n; numof_data = 0; i = 0;
  while (numof_data < *nor_p) {
    if (sscanf(&buf[cursor], "%d%n", &temp, &n) <= 0) {
      if (read(fd, buf, BUFSIZ) <= 0) {
	 fprintf(stderr, "Few data in: %s.\n", file_name);
	 exit(1);
      }
      cursor = 0;
      if (sscanf(&buf[cursor], "%d%n", &temp, &n) <= 0) {
	 fprintf(stderr, "Illegal data.\n");
	 exit(1);
      }
    }
    bcreate(&(*matrix_p)[numof_data * *noc_p + i]);
    bint(&(*matrix_p)[numof_data * *noc_p + i], temp);
    i++; cursor += n;
    if (i == *noc_p) {
      numof_data++;
      i = 0;
    }
  }

  free(buf);
  close(fd);
  printf("Loading has finished.\n");
}



void random_matrix()
{
  int i;
  time_t t;

  srand((unsigned int) time(&t));
  matrix = (DTYPE *) malloc(numof_rows * numof_columns * sizeof(DTYPE));
  for (i = 0; i < numof_columns; i++) {
    bcreate(&matrix[i]);
    bint(&matrix[i], 1);
  }
  for (i = numof_columns; i < numof_rows * numof_columns; i++) {
    bcreate(&matrix[i]);
    bint(&matrix[i], (int) rand() % 1000 - 499);
  }
}



int exchange(ent, lea)
int ent, lea;
{
  int i, j;
  INDEX temp_index;
  DTYPE temp1, temp2;

  bcreate(&temp1);
  bcreate(&temp2);

  for (i = 0; i <= codimension; i++) 
    bint(&equation[i], 0);
  bass(&temp2, &tableau[lea * rank + ent]);
  for (i = 0; i < rank; i++) {
    bdiv(&temp1, &tableau[lea * rank + i], &temp2);
    bneg(&tableau[lea * rank + i], &temp1);
  }
  binv(&tableau[lea * rank + ent], &temp2);

  for (j = 0; j < sizeof_tableau; j++)
    if (j != lea) {
      for (i = 0; i < rank; i++)
	 if (i != ent) {
	   bmul(&temp1, &tableau[lea * rank + i], &tableau[j * rank + ent]);
	   badd(&equation[i], &temp1, &tableau[j * rank + i]);
	 } else
	   bmul(&equation[i], &tableau[lea * rank + i], &tableau[j * rank + ent]);
      for (i = 0; i < rank; i++)
	 bass(&tableau[j * rank + i], &equation[i]);
    }
  temp_index = basic[lea]; basic[lea] = nonbasic[ent]; nonbasic[ent] = temp_index;

  bdestroy(&temp2);
  bdestroy(&temp1);
  return LP_CONTINUE;
}



int pivot()
{
  int enter, leave, i;
  DTYPE temp1, temp2;

  bcreate(&temp1);
  bcreate(&temp2);

  enter = 1;
  while (enter < rank && bcmp(&tableau[enter], &BZERO) <= 0) 
    enter++;
  if (enter == rank)
    return LP_TERMINATE;
  bass(&temp1, &tableau[enter]);
  for (i = enter + 1; i < rank; i++) {
    if (bcmp(&tableau[i], &BZERO) > 0 && nonbasic[enter] > nonbasic[i]) {
      bass(&temp1, &tableau[i]);
      enter = i;
    }
  }
#ifdef DEBUG
  printf("Entering variable is x%d:\t pivot.\n", (int) nonbasic[enter]);
#endif
  leave = 1;
  while ((bcmp(&tableau[leave * rank], &BZERO) < 0 
	  || bcmp(&tableau[leave * rank + enter], &BZERO) >= 0)
	 && leave < sizeof_tableau)
    leave++;
  if (leave == sizeof_tableau)
    return LP_UNBOUNDED;

  bdiv(&temp2, &tableau[leave * rank], &tableau[leave * rank + enter]);
  for (i = 1; i < sizeof_tableau; i++) {
    if (bcmp(&tableau[i * rank], &BZERO) >= 0 
	&& bcmp(&tableau[i * rank + enter], &BZERO) < 0) {
      bdiv(&temp1, &tableau[i * rank], &tableau[i * rank + enter]);
      if (bcmp(&temp2, &temp1) < 0) {
	 bass(&temp2, &temp1);
	 leave = i;
      }
      else if (bcmp(&temp2, &temp1) == 0 && basic[leave] > basic[i])
	 leave = i;
    }
  }
#ifdef DEBUG
  printf("Leaving variable is x%d:\t pivot.\n", basic[leave]);
#endif
  bdestroy(&temp2);
  bdestroy(&temp1);
  return exchange(enter, leave);
}


void print_tableau()
{
  int i, j;

  for (i = 0; i < sizeof_tableau; i++) {
    if (i == 0)
      printf("w  = ");
    else
      printf("x%d = ", (int) basic[i]);
    bprintf(&tableau[i * rank + 0]);
    printf(" ");
    for (j = 1; j < rank; j++) {
      if (bcmp(&tableau[i * rank + j], &BZERO) >= 0)
	 printf("+");
      bprintf(&tableau[i * rank + j]);
      printf(" x%d ", (int) nonbasic[j]);
    }
    printf("\n");
  }
}


int linear_program(nor, noc, answer)
int  nor, noc;
DTYPE *answer;
{
  int i;
  int state, retval;
  int enter,  leave;
  char msg[BUFSIZ];

  sizeof_tableau = nor; rank = noc;
  basic = (INDEX *) malloc(sizeof_tableau * sizeof(INDEX));
  basic[0] = 0;
  for (i = 1; i < sizeof_tableau; i++) basic[i] = i + rank - 1;
  nonbasic = (INDEX *) malloc(rank * sizeof(INDEX));
  for (i = 1; i < rank; i++) nonbasic[i] = i;
#ifdef DEBUG
  print_tableau();
  printf(">>>\n");
#endif
  for (enter = 1; enter < rank; enter++) {
    leave = 1;
    while (leave < sizeof_tableau 
	   && (bcmp(&tableau[leave * rank + enter], &BZERO) == 0 || basic[leave] < rank))
      leave++;
    if (leave == sizeof_tableau) {
      fprintf(stderr, "Fatal error: initial_pivot.\n");
      sprintf(msg, "Fatal error: initial_pivot.\n");
      write(fd, msg, strlen(msg));
      exit(1);
    }
#ifdef DEBUG
    printf("%d <---> %d: initial_pivot.\n", nonbasic[enter], basic[leave]);
#endif
    exchange(enter, leave);
  }

  do {
#ifdef DEBUG
    print_tableau();
#endif
  } while ((state = pivot()) == LP_CONTINUE);

  for (i = 0; i < rank - 1; i++)
    bint(&answer[i], 0);
  for (i = 1; i < sizeof_tableau ; i++)
    if (basic[i] > 0 && basic[i] < rank)
      bass(&answer[basic[i] - 1], &tableau[i * rank]);

  switch (state) {
    case LP_TERMINATE :
      retval = 1;
      break;
    case LP_UNBOUNDED :
      retval = -1;
      break;
  }

  free(basic);
  free(nonbasic);
  return retval;
  free(msg);
}




/* main part */


void print_tuple(s, n)
INDEX *s;
int n;
{
  int i;
  printf("<");
  for (i = 0; i < n; i++) printf(" %2d", s[i] + 1);
  printf(">");
}

void print_simplex(s)
Simplex s;
{ printf("[S]"); print_tuple(s, dimension); }

void print_radon(r)
Radon r;
{ printf("[R]"); print_tuple(r, dimension + 1); }

void write_simplex(f, s)
int f;
Simplex s;
{ 
  int i;
  char msg[BUFSIZ];
  write(f, "[S]", 3); 
  write(f, "<", 1);
  for (i = 0; i < dimension; i++) {
    sprintf(msg, " %2d", s[i] + 1);
    write(f, msg, strlen(msg));
  }
  write(f, ">", 1);
  free(msg);
}

int read_simplex(f, s)
int f;
Simplex s;
{
  int i, j;
  char msg[BUFSIZ];
  read(f, msg, 4);
  if (msg[0] != '[') {
    free(msg);
    return 0;
  }
  for (i = 0; i < dimension; i++) {
    if (read(f, msg, 3) < 3) {
      free(msg);
      return 0;
    }
    sscanf(msg, " %2d", &j);
    s[i] = j - 1;
  }
  read(f, msg, 2);
  free(msg);
  return 1;
}



int member_simplex(smplx, elem)
Simplex smplx;
INDEX elem;
{
  int i;
  i = 0;
  while (i < dimension && smplx[i] != elem) i++;
  return (i != dimension);
}



int simplex_in_radon(smplx, rdn)
Simplex smplx;
Radon rdn;
{
  int i, j;
  if (smplx == NULL || rdn == NULL)
    return 0;
  i = 0; j = 0;
  do {
    if (smplx[i] == rdn[j]) {
      i++; j++;
    } else if (smplx[i] > rdn[j])
      j++;
    else
      return 0;
  } while ((i < dimension) && (j - i <= 1));
  return (i == dimension);
}



int lexico_simplex(smplx1, smplx2)
Simplex smplx1;
Simplex smplx2;
{
  int i;
  i = 0;
  while  (i < dimension && smplx1[i] == smplx2[i]) i++;
  return (i == dimension)? 0 : smplx1[i] - smplx2[i];
}

int lexico_radon(rdn1, rdn2)
Radon rdn1;
Radon rdn2;
{
  int i;
  i = 0;
  while  (i <= dimension && rdn1[i] == rdn2[i]) i++;
  return (i == dimension + 1)? 0 : rdn1[i] - rdn2[i];
}



int radon_from_simplex(rdn, smplx, elem)
Radon     rdn;
Simplex smplx;
INDEX    elem;
{
  int i;
  i = 0;
  while (i < dimension && smplx[i] < elem) {
    rdn[i] = smplx[i];
    i++;
  }
  if (i != dimension && smplx[i] == elem)
      return 0;
  rdn[i++] = elem;
  while (i < dimension + 1) {
    rdn[i] = smplx[i - 1];
    i++;
  }
  return 1;
}



int increment(tpl)
INDEX *tpl;
{
  int i, j;

  for (j = 0; j < codimension; j++)
    if (tpl[j] > dimension + j)
      return 0;
  if (tpl[0] == dimension)
    return 0;
  if (tpl[codimension - 1] < numof_columns - 1) {
    tpl[codimension - 1]++;
  } else {
    j = codimension - 2;
    while (tpl[j] == dimension + j) j--;
    tpl[j]++;
    for (i = j + 1; i < codimension; i++)
      tpl[i] = tpl[j] + i - j;
  }
  return 1;
}



int complement(tpl1, size, tpl2)
INDEX *tpl1;
int    size;
INDEX *tpl2;
{
  int i, c1, c2;
  c1 = 0; 
  c2 = 0;
#ifdef DEBUG
  printf(">>>: complement: ");
  print_tuple(tpl1, size);
#endif 
  for (i = 0; i < numof_columns; i++)
    if (c1 == size || i != tpl1[c1]) {
      tpl2[c2] = i;
      c2++;
    } else {
      c1++;
    }
  return (c2 == numof_columns - size);
}



int radon_to_simplex(rdn, elem, smplx)
Radon rdn;
INDEX elem;
Simplex smplx;
{
  int i, j;
  i = 0;
  while (i <= dimension && rdn[i] != elem) i++;
  if (i == dimension + 1)
    return 0;
  for (j = 0; j < i; j++)
    smplx[j] = rdn[j];
  for (j = i; j < dimension; j++)
    smplx[j] = rdn[j + 1];
  return 1;
}



void print_tri(rt)
Tri rt;
{
  int i, c;
  i = 0; 
  c = 1;
  while (rt[i] != -1) {
    print_simplex(&rt[i]);
    i += dimension;
    if (c == 3) {
      printf("\n");
      c = 0;
    }
    c++;
  }
  printf("\n");
}

void print_cand(cd, nm)
Cand cd;
Normal nm;
{
  int i, ii;
  i = 0;
  ii = 0;
  while (cd[i] != -1) {
    print_radon(&cd[i]);
    print_matrix(&nm[ii], 1, codimension);
    i += (dimension + 1);
    ii += codimension;
  }
}

void write_tri(f, rt)
int f;
Tri rt;
{
  int i, c;
  i = 0;
  c = 1;
  while (rt[i] != -1) {
#ifndef EXPERIMENT
    write_simplex(f, &rt[i]);
    if (c == 3) {
      write(f, "\n", 1);
      c = 0;
    }
    c++;
#endif
    i += dimension;
  }
  write(f, "\n", 1);
}

void read_tri(f, rt)
int f;
Tri rt;
{
  int i, j;
  i = 0;
  while (read_simplex(f, &rt[i]))
    i += dimension;
  for (j = 0; j < dimension; j++)
    rt[i + j] = -1;
}



int find_tri(rt, smplx)
Tri rt;
Simplex smplx;
{
  int i, j;
  i = 0;
  while (rt[i] != -1) {
    j = 0;
    while  (j < dimension && rt[i + j] == smplx[j]) j++;
    if (j == dimension) 
      return 1;
    i += dimension;
  }
  return 0;
}



void assignment_tri(rt1, rt2)
Tri rt1;
Tri rt2;
{
  int i;
  i = 0;
  while (rt1[i] != -1 || rt2[i] != -1) {
    rt1[i] = rt2[i];
    i++;
  }
}

void assignment_cand(cd1, nm1, cd2, nm2)
Cand   cd1;
Normal nm1;
Cand   cd2;
Normal nm2;
{
  int i, j, jj;
  j = 0;
  jj = 0;
  while (cd1[j] != -1 || cd2[j] != -1) {
    for (i = 0; i <= dimension; i++)
      cd1[j + i] = cd2[j + i];
    for (i = 0; i < codimension; i++)
      bass(&nm1[jj + i], &nm2[jj + i]);
    j += (dimension + 1);
    jj += codimension;
  }
}

void delete_for_flip(rt, r)
Tri  rt;
Radon r;
{
  int i, j, k;
  i = 0;
  while (rt[i] != -1) {
    if (simplex_in_radon(&rt[i], r)) {
      j = i + dimension;
      while (rt[j] != -1) {
	rt[j - dimension] = rt[j];
	j++;
      }
      for (k = 0; k < dimension; k++)
	rt[j + k - dimension] = -1;
    } else
      i += dimension;
  }
}



int equal_tri(rt1, rt2)
Tri rt1;
Tri rt2;
{
  int i;
  i = 0;
  while (rt1[i] != -1) {
    if (lexico_simplex(&rt1[i], &rt2[i]) != 0) 
      return 0;
    i += dimension;
  }
  return 1;
}



void add_tri(rt, s) 
Tri rt;
Simplex s;
{
  int i, j, k, flag;
  i = 0;
  while (rt[i] != -1) {
    if ((flag = lexico_simplex(&rt[i], s)) > 0) {
      k = i;
      while (rt[k] != -1) 
	k += dimension;
      for (j = k - 1; j >= i; j--) 
	rt[j + dimension] = rt[j];
      for (j = 0; j < dimension; j++) 
	rt[i + j] = s[j];
      return;
    } else if (flag == 0)
      return;
    i += dimension;
  }
  for (j = 0; j < dimension; j++)
    rt[i + j] = s[j];
  return;
}

void add_cand(cd, nm, r, v)
Cand cd;
Normal nm;
Radon r;
DTYPE *v;
{
  int i, ii, j, k, kk, flag;
  i = 0;
  ii = 0;
  while (cd[i] != -1) {
    if ((flag = lexico_radon(&cd[i], r)) > 0) {
      k = i; 
      kk = ii;
      while (cd[k] != -1) {
	k += (dimension + 1);
	kk += codimension;
      }
      for (j = k - 1; j >= i; j--) 
	cd[j + dimension + 1] = cd[j];
      for (j = 0; j <= dimension; j++) 
	cd[i + j] = r[j];
      for (j = kk - 1; j >= ii; j--)
	bass(&nm[j + codimension], &nm[j]);
      for (j = 0; j < codimension; j++)
	bass(&nm[ii + j], &v[j]);
      return;
    } else if (flag == 0)
      return;
    i += (dimension + 1);
    ii += codimension;
  }
  for (j = 0; j <= dimension; j++)
    cd[i + j] = r[j];
  for (j = 0; j < codimension; j++)
    bass(&nm[ii + j], &v[j]);
  return;
}



void random_vector(vector, size)
DTYPE *vector;
int size;
{
  int i;
  time_t t;
  srand((unsigned int) time(&t));
  for(i = 0; i < size; i++)
    bint(&vector[i], (int) rand() % 1000 - 499);
}



void initial_triangulation(Delta)
Tri Delta;
{
  int i, j;
  int flag, dim, codim;
  int sign1, sign2;
  INDEX   *cone;
  Simplex smplx;
  DTYPE *vector;

  vector = (DTYPE *) malloc(codimension * sizeof(DTYPE));
  for (i = 0; i < codimension; i++)
    bcreate(&vector[i]);
  cone  = (INDEX *) malloc(codimension * sizeof(INDEX)); 
  smplx = (INDEX *) malloc(dimension   * sizeof(INDEX));

  for (i = 0; i < codimension; i++)
    cone[i] = i;
  do {
    flag = 1;
    printf(">>>: initial_triangulation...\n");
    random_vector(vector, codimension);
    print_matrix(vector, 1, codimension);
    do {
      flag = 1;
      for (i = 0; i < codimension; i++)
	 for (j = 0; j < codimension; j++)
	   bass(&initial_mtrx[i + j * (codimension + 1)], &gale[cone[i] + j * numof_columns]);
      for (i = 0; i < codimension; i++)
	 bass(&initial_mtrx[codimension + i * (codimension + 1)], &vector[i]);
      codim = gale_transform(initial_mtrx, codimension, codimension + 1, initial_glmtrx);
      dim = codimension + 1 - codim;
      if (dim != codimension) {
	 /* Caution!!!: The chosen vector may be contained
	    in the interiors of some other cones. */
	 /* A degenerate cone is detected. */
	 flag = 0;
      }
      if (flag == 1 && bcmp(&initial_glmtrx[codimension], &BZERO) != 0) {
	sign2 = bcmp(&initial_glmtrx[codimension], &BZERO);
	for (i = 0; i < codimension; i++) {
	  sign1 = bcmp(&initial_glmtrx[i], &BZERO);
	  if (sign1 == 0) {
	    /* We must choose another vector. */
	    flag = -1;
	    i = codimension + 1;
	  } else if ((sign1 > 0 && sign2 > 0) 
		     || (sign1 < 0 && sign2 < 0))
	    i = codimension + 1;
	}
	if (i == codimension) {
	  complement(cone, codimension, smplx);
	  add_tri(Delta, smplx);
	}
      }
    } while (flag != -1 && increment(cone));
  } while (flag == -1);

  for (i = 0; i < codimension; i++)
    bdestroy(&vector[i]);
  free(vector);
  free(cone);
  free(smplx);
  printf("...<<<: initial_triangulation.\n");
}



int signed_radon(rdn, plus, minus, numof_plus, numof_minus)
Radon        rdn;
INDEX      *plus;
INDEX     *minus;
int  *numof_plus;
int *numof_minus;
{
  int i, j, p, m, sgn;
  int dim, codim;

  for (i = 0; i <= dimension; i++)
    for (j = 0; j < dimension; j++)
      bass(&signed_radon_mtrx[j * (dimension + 1) + i], &matrix[j * numof_columns + rdn[i]]);
  codim = gale_transform(signed_radon_mtrx, dimension, dimension + 1, signed_radon_glmtrx);
  dim = dimension + 1 - codim;
  if (dim != dimension) 
    return 0;

  p = 0; m = 0;
  for (i = 0; i <= dimension; i++) {
    sgn = bcmp(&signed_radon_glmtrx[i], &BZERO);
    if (sgn > 0)
      plus[p++]  = rdn[i];
    else if (sgn < 0)
      minus[m++] = rdn[i];
  }

  *numof_plus  = p;
  *numof_minus = m;
  return 1;
}



int flippable(rdn, Delta)
Radon rdn;
Tri Delta;
{
  int ii, jj, k, p, m;
  INDEX *pls, *mns;
  Simplex smplx;

  smplx = (INDEX *) malloc(dimension * sizeof(INDEX));
  pls   = (INDEX *) malloc(dimension * sizeof(INDEX));
  mns   = (INDEX *) malloc(dimension * sizeof(INDEX));
  if (!signed_radon(rdn, pls, mns, &p, &m)) 
    return 0;
  ii = 0;
  for (k = 0; k < p; k++)
    if (radon_to_simplex(rdn, pls[k], smplx))
      if (find_tri(Delta, smplx))
	 ii++;
  if (p == ii) {
    free(pls);
    free(mns);
    free(smplx);
    return 1;
  }
  jj = 0;
  for (k = 0; k < m; k++)
    if (radon_to_simplex(rdn, mns[k], smplx))
      if (find_tri(Delta, smplx))
	 jj++;

  free(mns);
  free(pls);
  free(smplx);

  return (m == jj);
}



void adjust_normal(vector, Delta, rdn)
DTYPE *vector;
Tri Delta;
Radon rdn;
{
  int i, cursor;
  char msg[BUFSIZ];
  DTYPE temp1, temp2, temp3;

  bcreate(&temp1);
  bcreate(&temp2);
  bcreate(&temp3);

  i = 0;
  while (Delta[i] != -1 && !simplex_in_radon(&Delta[i], rdn))
    i += dimension;
  if (Delta[i] == -1) {
    fprintf(stderr, "Can't find any appropriate simplex: adjust_normal.\n");
    sprintf(msg, "Can't find any appropriate simplex: adjust_normal.\n");
    write(fd, msg, strlen(msg));
    exit(1);
  }
  cursor = 0;
  while (cursor < dimension + 1 && member_simplex(&Delta[i], rdn[cursor]))
    cursor++;
  bint(&temp2, 0);
  for (i = 0; i < codimension; i++) {
    bmul(&temp1, &vector[i], &gale[i * numof_columns + rdn[cursor]]);
    badd(&temp3, &temp2, &temp1);
    bass(&temp2, &temp3);
  }
  if (bcmp(&temp2, &BZERO) < 0)
    for (i = 0; i < codimension; i++) {
      bneg(&temp1, &vector[i]);
      bass(&vector[i], &temp1);
    }

  bdestroy(&temp1);
  bdestroy(&temp2);
  bdestroy(&temp3);
  free(msg);
  return;
}



int radon_to_normal(rdn, Delta, vec)
Radon rdn;
Tri Delta;
DTYPE *vec;
{
  int i, j, cd;
  INDEX *nu;

#ifdef DEBUG
  printf(">>>: radon_to_normal: ");
  print_radon(rdn);
  printf("in\n");
  print_tri(Delta);
#endif
  nu   = (INDEX *) malloc((codimension - 1) * sizeof(INDEX));
  complement(rdn, dimension + 1, nu);
  for (i = 0; i < codimension - 1; i++)
    for (j = 0; j < codimension; j++)
      bass(&normal_mtrx[i * codimension + j], &gale[j * numof_columns + nu[i]]);
#ifdef DEBUG
  print_tuple(nu, codimension - 1);
  print_matrix(normal_mtrx, codimension - 1, codimension);
#endif
  cd = gale_transform(normal_mtrx, codimension - 1, codimension, normal_glmtrx);
  if (cd == 1) {
    for (i = 0; i < codimension; i++) 
      bass(&vec[i], &normal_glmtrx[i]);
  } else {
#ifdef DEBUG
    print_matrix(normal_glmtrx, cd, codimension);
#endif
    free(nu);
    return 0;
  }
  adjust_normal(vec, Delta, rdn);
  free(nu);
#ifdef DEBUG
  printf("... <<<: radon_to_normal.\n");
#endif
  return 1;
}



void make_candidate(Delta, Gamma, Nu)
Tri  Delta;
Cand Gamma;
Normal  Nu;
{
  int i, j;
  Radon rho;
  char msg[BUFSIZ];

  rho   = (INDEX *) malloc((dimension + 1) * sizeof(INDEX));
  for (i = 0; i < codimension; i++)
    bint(&vector1[i], 0);

#ifdef DEBUG
  printf(">>>: make_candidate...\n");
#endif
  j = 0;
  while (Delta[j] != -1) {
    for (i = 0; i < numof_columns; i++) 
      if (radon_from_simplex(rho, &Delta[j], i)) 
	if (flippable(rho, Delta)) {
#ifdef DEBUG
	  print_radon(rho);
	  printf("is flippable in\n");
	  print_tri(Delta);
#endif
	  if (radon_to_normal(rho, Delta, vector1)) {
	    add_cand(Gamma, Nu, rho, vector1);
	  } else {
	    fprintf(stderr, "Fatal error: radon_to_normal.\n");
	    sprintf(msg, "Fatal error: radon_to_normal.\n");
	    write(fd, msg, strlen(msg));
	    exit(1);
	  }
	}
    j += dimension;
  }
#ifdef DEBUG
  print_cand(Gamma, Nu);
  printf("...<<<: make_candidate.\n");
#endif
  free(msg);
  free(rho);
}



int representative(Gamma, Nu, rdn, nrml)
Cand  Gamma;
Normal   Nu;
Radon   rdn;
DTYPE *nrml;
{
  int i, j, jj;
  char msg[BUFSIZ];

  j = 0; 
  jj = 0;
  while (Gamma[j] != -1 && lexico_radon(rdn, &Gamma[j]) != 0) {
    j += (dimension + 1);
    jj += codimension;
  }
  for (i = 0; i < codimension; i++)
    bass(&nrml[i], &Nu[jj + i]);
  j = 0;
  jj = 0;
  while (Gamma[j] != -1) {
    i = 0;
    while (i < codimension && bcmp(&Nu[jj + i], &nrml[i]) == 0)
      i++;
    if (i == codimension) {
      free(msg);
      return (lexico_radon(rdn, &Gamma[j]) == 0);
    }
    j += (dimension + 1);
    jj += codimension;
  }
  fprintf(stderr, "Fatal error: representative.\n");
  sprintf(msg, "Fatal error: radon_to_normal.\n");
  write(fd, msg);
  exit(1);
}



void flipping(Delta, Gamma, Nu, nrml, Delta_, Gamma_, Nu_)
Tri   Delta;
Cand  Gamma;
Normal   Nu;
DTYPE *nrml;
Tri  Delta_;
Cand Gamma_;
Normal  Nu_;
{
  int i, j, jj, p, m, flag;
  INDEX *pls, *mns;
  DTYPE temp;
  Simplex smplx;
  char msg[BUFSIZ];

  smplx = (INDEX *) malloc(dimension * sizeof(INDEX));
  pls   = (INDEX *) malloc(dimension * sizeof(INDEX));
  mns   = (INDEX *) malloc(dimension * sizeof(INDEX));
  bcreate(&temp);

#ifdef DEBUG
  printf(">>>: flipping...\n ");
#endif
  j = 0;
  jj = 0;
  assignment_tri(Delta_, Delta);
  while (Gamma[j] != -1) {
    i = 0;
    while (i < codimension && bcmp(&nrml[i], &Nu[jj + i]) == 0)
      i++;
    if (i == codimension) {
      signed_radon(&Gamma[j], pls, mns, &p, &m);
      flag = 1;
      radon_to_simplex(&Gamma[j], pls[0], smplx);
      if (!find_tri(Delta, smplx)) {
	radon_to_simplex(&Gamma[j], mns[0], smplx);
	if (!find_tri(Delta, smplx)) {
	  fprintf(stderr, "Fatal error: flipping.\n");
	  sprintf(msg, "Fatal error: flipping.\n");
	  write(fd, msg, strlen(msg));
	  exit(1);
	}
	flag = -1;
      }
      delete_for_flip(Delta_, &Gamma[j]);
      if (flag == -1) {
	for (i = 0; i < p; i++) {
	  radon_to_simplex(&Gamma[j], pls[i], smplx);
	  add_tri(Delta_, smplx);
	}
      } else {
	for (i = 0; i < m; i++) {
	  radon_to_simplex(&Gamma[j], mns[i], smplx);
	  add_tri(Delta_, smplx);
	}
      }
    }
    j += (dimension + 1);
    jj += codimension;
  }
  assignment_cand(Gamma_, Nu_, Gamma, Nu);
  j = 0;
  jj = 0;
  while (Gamma_[j] != -1) {
    i = 0;
    while (i < codimension && bcmp(&Nu_[jj + i], &nrml[i]) == 0)
      i++;
    if (i == codimension)
      for (i = 0; i < codimension; i++) {
	bneg(&temp, &Nu_[jj + i]);
	bass(&Nu_[jj + i], &temp);
      }
    j += (dimension + 1);
    jj += codimension;
  }
  
#ifdef DEBUG
  print_tri(Delta);
  printf(" to...\n");
  print_tri(Delta_);
  printf("...<<<: flipping.\n");
#endif
  free(msg);
  free(mns);
  free(pls);
  free(smplx);
  bdestroy(&temp);
}



void volume_vector(Delta, volvec)
Tri Delta;
DTYPE *volvec;
{
  int i, j, k;
  DTYPE vlm, temp;

  bcreate(&vlm);
  bcreate(&temp);
  j = 0;
  while (Delta[j] != -1) {
    for (i = 0; i < dimension; i++)
      for (k = 0; k < dimension; k++)
	bass(&volume_mtrx[k * dimension + i], &matrix[k * numof_columns + Delta[j + i]]);
    volume(&vlm);
    for (i = 0; i < dimension; i++) {
      badd(&temp, &volvec[Delta[j + i]], &vlm);
      bass(&volvec[Delta[j + i]], &temp);
    }
    j += dimension;
  }
  bdestroy(&vlm);
  bdestroy(&temp);
}



int comparison_tri(Dlt1, Dlt2)
Tri Dlt1;
Tri Dlt2;
{
  int i, flag;

  for (i = 0; i < numof_columns; i++) { 
    bint(&volvec1[i], 0);
    bint(&volvec2[i], 0);
  }
  volume_vector(Dlt1, volvec1);
  volume_vector(Dlt2, volvec2);
#ifdef DEBUG
  print_matrix(volvec1, 1, numof_columns);
  print_matrix(volvec2, 1, numof_columns);
#endif
  i = 0;
  flag = bcmp(&volvec1[0], &volvec2[0]);
  while (i < numof_columns && flag == 0) {
#ifdef DEBUG
    printf(" %d , ", flag);
#endif 
    flag = bcmp(&volvec1[i], &volvec2[i]);
    i++;
  }
#ifdef DEBUG
  printf("(%d): comparison_tri.\n", flag);
#endif
  return flag;
}



void make_tableau(Nu, size, obj)
Normal  Nu;
int   size;
DTYPE *obj;
{
  int i, j, k;
  INDEX *nu;
  DTYPE temp1, temp2;

  nu = (INDEX *) malloc((codimension - 1) * sizeof(INDEX));
  bcreate(&temp1);
  bcreate(&temp2);
  k = 0;
  for (i = 0; i < size; i++) {
    bint(&temp1, 0);
    for (j = 0; j < codimension; j++) {
      bass(&tableau[(i + 1) * (codimension + 1) + j + 1], &Nu[k + j]);
      badd(&temp2, &Nu[k + j], &temp1);
      bass(&temp1, &temp2);
    }
    bneg(&tableau[(i + 1) * (codimension + 1)], &temp1);
    k += codimension;
  }

  /* Set the objective function. */
  for (i = 0; i < codimension; i++) {
    bneg(&tableau[i + 1], &obj[i]);
    bass(&obj[i], &tableau[i + 1]);
  }
  bint(&tableau[0], 0);

  /* Set an artificial constraint. */
  for (i = 0; i < codimension; i++)
    bint(&tableau[(size + 1) * (codimension + 1) + i + 1], -1);
  bint(&tableau[(size + 1) * (codimension + 1)], codimension * 2);

  free(nu);
  bdestroy(&temp2);
  bdestroy(&temp1);
}



int regularity(Gamma, Nu, nrml)
Cand  Gamma;
Normal   Nu;
DTYPE *nrml;
{
  int i, j, count, retval;
  DTYPE temp1, temp2, temp3;

  j = 0;
  count = 0;
  while (Gamma[j] != -1) {
    j += (dimension + 1);
    count++;
  }
  for (i = 0; i < codimension; i++)
    bint(&vector2[i], 0);
  bcreate(&temp1);
  bcreate(&temp2);
  bcreate(&temp3);

#ifdef DEBUG
  printf(">>>: regularity...\n");
#endif
  make_tableau(Nu, count, nrml);
  linear_program(count + 2, codimension + 1, vector2);

  /* Check if the optimal solution lies on the facet
     corresponding to the flipped Radon partition. */
  retval = 1;
  bint(&temp1, 1);
  for (i = 0; i < codimension; i++) {
    bsub(&temp2, &vector2[i], &temp1);
    bass(&vector2[i], &temp2);
  }
  bint(&temp2, 0);
  for (i = 0; i < codimension; i++) {
    bmul(&temp1, &nrml[i], &vector2[i]);
    badd(&temp3, &temp2, &temp1);
    bass(&temp2, &temp3);
  }
  if (bcmp(&temp2, &BZERO) == 0)
    retval = 0;

#ifdef DEBUG
  printf("...<<<: regularity.\n");
#endif
  bdestroy(&temp1);
  bdestroy(&temp2);
  bdestroy(&temp3);
  return retval;
}



int adjacency(Delta, ind, Adj)
Tri Delta;
int ind;
Tri Adj;
{
  int i, retval;
  Radon rho;

#ifdef DEBUG
  printf("\n>>>: adjacency...\n");
#endif
  for (i = 0; i < codimension; i++)
    bint(&vector3[i], 0);
  rho = (INDEX *) malloc((dimension + 1) * sizeof(INDEX));
  for (i = 0; i < sizecand * (dimension + 1); i++) {
    global_Gamma[i] = -1;
    global_Gamma_[i] = -1;
  }

  make_candidate(Delta, global_Gamma, set1of_normals);
  if (global_Gamma[ind * (dimension + 1)] == -1) {
    free(rho);
    return 0;
  }
  for (i = 0; i <= dimension; i++)
    rho[i] = global_Gamma[ind * (dimension + 1) + i];

  retval = -1;
  if (representative(global_Gamma, set1of_normals, rho, vector3)) {
    flipping(Delta, global_Gamma, set1of_normals, vector3, Adj, global_Gamma_, set2of_normals);
    if (regularity(global_Gamma_, set2of_normals, vector3)) 
      retval = 1;
    else
      retval = -1;
  } else {
    retval = 2;
#ifdef DEBUG
    printf("Not the representative: adjacency.\n");
#endif
  }
#ifdef DEBUG
  printf("...<<<: adjacency.\n");
#endif
  free(rho);
  return retval;
}



int local_search(Delta, Best)
Tri Delta;
Tri  Best;
{
  int i, j;
  char msg[BUFSIZ];

  for (i = 0; i < codimension; i++)
    bint(&vector4[i], 1);

  assignment_tri(Best, Delta);
  for (i = 0; i < sizetri * dimension; i++)
    global_Delta_[i] = -1;
  for (i = 0; i < sizecand * (dimension + 1); i++) {
    global_Gamma[i] = -1;
    global_Gamma_[i] = -1;
  }

#ifdef DEBUG
  printf("\n>>>: local_search...\n");
#endif

  make_candidate(Delta, global_Gamma, set1of_normals);
  j = 0;
  while (global_Gamma[j] != -1) {
    if (representative(global_Gamma, set1of_normals, &global_Gamma[j], vector4)) {
      flipping(Delta, global_Gamma, set1of_normals, vector4, 
	       global_Delta_, global_Gamma_, set2of_normals);
      if (regularity(global_Gamma_, set2of_normals, vector4)) 
	if (comparison_tri(Best, global_Delta_) < 0) {
	  assignment_tri(Best, global_Delta_);
#ifdef DEBUG
	  printf("Exchange!!!\n");
#endif
	}
    } else {
#ifdef DEBUG
      printf("Not the representative: local_search.\n");
#endif      
    }
    j += (dimension + 1);
  }
#ifdef DEBUG
  print_tri(Delta);
  printf(" to...\n");
  print_tri(Best);
  printf("[%d]...<<<: local_search.\n", !equal_tri(Best, Delta));
#endif
  free(msg);
  return (!equal_tri(Best, Delta));
}



void regular_triangulations(restart)
int restart;
{
  int i, flag;
  char msg[BUFSIZ];
  Tri Opt;
  Tri Delta;
  Tri Delta_;
  Tri Delta0;

  watch_on();

  printf(">>>: >>>: Original matrix.\n");
  print_matrix(matrix, numof_rows, numof_columns);

  gale = (DTYPE *) malloc((numof_columns - 1) * numof_columns * sizeof(DTYPE));
  for (i = 0; i < (numof_columns - 1) * numof_columns; i++)
    bcreate(&gale[i]);
  codimension = gale_transform(matrix, numof_rows, numof_columns, gale);
  dimension = numof_columns - codimension;
  if (!check(matrix, gale, dimension, codimension, numof_columns)) {
    fprintf(stderr, "Can't compute Gale transform.\n");
    exit(1);
  }
  printf(">>>: >>>: Gale transform.\n");
  print_matrix(gale, codimension, numof_columns);


  sizetri    = SMPLXNUM + 1;
  sizecand   = SMPLXNUM * dimension + 1; 
  tableau      = (DTYPE *) malloc((sizecand + 2) * (codimension + 1) * sizeof(DTYPE));
  initial_mtrx    = (DTYPE *) malloc(codimension * (codimension + 1) * sizeof(DTYPE));
  initial_glmtrx  = (DTYPE *) malloc(codimension * (codimension + 1) * sizeof(DTYPE));
  volume_mtrx     = (DTYPE *) malloc(dimension * dimension * sizeof(DTYPE));
  normal_mtrx     = (DTYPE *) malloc((codimension - 1) * codimension * sizeof(DTYPE));
  normal_glmtrx   = (DTYPE *) malloc((codimension - 1) * codimension * sizeof(DTYPE));
  signed_radon_mtrx  =  (DTYPE *) malloc(dimension * (dimension + 1) * sizeof(DTYPE));
  signed_radon_glmtrx = (DTYPE *) malloc(dimension * (dimension + 1) * sizeof(DTYPE));
  set1of_normals      = (DTYPE *) malloc(sizecand * codimension * sizeof(DTYPE));
  set2of_normals      = (DTYPE *) malloc(sizecand * codimension * sizeof(DTYPE));
  volvec1  = (DTYPE *) malloc(numof_columns * sizeof(DTYPE));
  volvec2  = (DTYPE *) malloc(numof_columns * sizeof(DTYPE));
  equation = (DTYPE *) malloc((codimension + 1) * sizeof(DTYPE));
  vector1  = (DTYPE *) malloc(codimension * sizeof(DTYPE));
  vector2  = (DTYPE *) malloc(codimension * sizeof(DTYPE));
  vector3  = (DTYPE *) malloc(codimension * sizeof(DTYPE));  
  vector4  = (DTYPE *) malloc(codimension * sizeof(DTYPE));
  global_Gamma  = (INDEX *) malloc(sizecand * (dimension + 1) * sizeof(INDEX));
  global_Gamma_ = (INDEX *) malloc(sizecand * (dimension + 1) * sizeof(INDEX));
  global_Delta_ = (INDEX *) malloc(sizetri  * dimension * sizeof(INDEX));

  for (i = 0; i < rank; i++)
    bcreate(&equation[i]);
  for (i = 0; i < (sizecand + 2) * (codimension + 1); i++)
    bcreate(&tableau[i]);
  for (i = 0; i < codimension * (codimension + 1); i++) {
    bcreate(&initial_mtrx[i]);
    bcreate(&initial_glmtrx[i]);
  }
  for (i = 0; i < dimension * dimension; i++)
    bcreate(&volume_mtrx[i]);
  for (i = 0; i < (codimension - 1) * codimension; i++) {
    bcreate(&normal_mtrx[i]);
    bcreate(&normal_glmtrx[i]);
  }
  for (i = 0; i < dimension * (dimension + 1); i++) {
    bcreate(&signed_radon_mtrx[i]);
    bcreate(&signed_radon_glmtrx[i]);
  }
  for (i = 0; i < sizecand * codimension; i++) {
    bcreate(&set1of_normals[i]);
    bcreate(&set2of_normals[i]);
  }
  for (i = 0; i < numof_columns; i++) {
    bcreate(&volvec1[i]);
    bcreate(&volvec2[i]);
  }
  for (i = 0; i <= codimension; i++) 
    bcreate(&equation[i]);
  for (i = 0; i < codimension; i++) {
    bcreate(&vector1[i]);
    bcreate(&vector2[i]);
    bcreate(&vector3[i]);
    bcreate(&vector4[i]);
  }

  Opt         = (INDEX *) malloc(sizetri * dimension * sizeof(INDEX));
  Delta       = (INDEX *) malloc(sizetri * dimension * sizeof(INDEX));
  Delta_      = (INDEX *) malloc(sizetri * dimension * sizeof(INDEX));
  Delta0      = (INDEX *) malloc(sizetri * dimension * sizeof(INDEX));
  for (i = 0; i < sizetri * dimension; i++) {
    Opt[i] = -1;
    Delta[i] = -1;
    Delta_[i] = -1;
    Delta0[i] = -1;
  }
  for (i = 0; i < sizecand * (dimension + 1); i++) {
    global_Gamma[i] = -1;
    global_Gamma_[i] = -1;
  }
  for (i = 0; i < sizetri * dimension; i++)
    global_Delta_[i] = -1;

#ifndef EXPERIMENT
  if (!restart) {
#endif
    initial_triangulation(Opt);
#ifndef EXPERIMENT
    sprintf(msg, ">>>: >>>: Initial triangulation.\n");
    write(fd, msg, strlen(msg));
    write_tri(fd, Opt);
#endif
    while (local_search(Opt, Delta)) {
      assignment_tri(Opt, Delta);
#ifndef EXPERIMENT
      sprintf(msg, "***--- To reach the optimum...\n");
      write(fd, msg, strlen(msg));
      write_tri(fd, Opt);
#endif
    } 
#ifndef EXPERIMENT
  } else 
    read_tri(fd__, Opt);
#endif

#ifndef EXPERIMENT
  sprintf(msg, "****** [1] Optimal triangulation ******\n");
  write(fd, msg, strlen(msg));
  write_tri(fd, Opt);
#endif
#ifdef DEBUG
  printf("\n****** [%d] Optimal triangulation ******\n", numof_reg);
  print_tri(Opt);
#endif

#ifndef EXPERIMENT
  if (!restart) 
    assignment_tri(Delta, Opt);
  else {
    read_tri(fd__, Delta);
    sprintf(msg, "\n****** [%d] Regular Triangulation ******\n", numof_reg);
    write(fd, msg, strlen(msg));
    write_tri(fd, Delta);
  }
#endif
    
  i = 0;
  do {
    while ((flag = adjacency(Delta, i, Delta_)) != 0) {
#ifndef EXPERIMENT
      sprintf(msg, "***------ [%d]th adjacent.", i+1);
      write(fd, msg, strlen(msg));
      if ((i + 1) % 3 == 0) 
	write(fd, "\n", 1);
#endif
#ifdef DEBUG
      printf("***------ [%d]th adjacent.", i+1);
      if ((i + 1) % 3 == 0)
	printf("\n");
#endif
      i++;
      if (flag == 1) {
	local_search(Delta_, Delta0);
	if (equal_tri(Delta0, Delta)) {
	  assignment_tri(Delta, Delta_);
	  numof_reg++;
#ifndef EXPERIMENT
	  sprintf(msg, "\n****** [%d] Regular Triangulation ******\n", numof_reg);
	  write(fd, msg, strlen(msg));
	  write_tri(fd, Delta);
#endif
#ifdef DEBUG
	  printf("\n****** [%d] Regular Triangulation ******\n", numof_reg);
	  print_tri(Delta);
#endif
	  i = 0;
	}
      } else if (flag == -1) {
#ifndef EXPERIMENT
/*	sprintf(msg, "****** NON-REGULAR OR NOT A TRIANGULATION ******\n");
	write(fd_, msg, strlen(msg));
	write_tri(fd_, Delta_);
	write(fd,  msg, strlen(msg));*/
#endif
#ifdef DEBUG
	printf("****** NON-REGULAR OR NOT A TRIANGULATION ******\n");
/*	print_tri(Delta_);*/
#endif
      } else if (flag == 2) {
      }
#ifdef DEBUG
      still();
#endif
    }
    if (!equal_tri(Opt, Delta)) {
      assignment_tri(Delta_, Delta);
      local_search(Delta_, Delta);
      i = 0;
      do {
	adjacency(Delta, i, Delta0);
	i++;
      } while (!equal_tri(Delta0, Delta_));
#ifndef EXPERIMENT
      sprintf(msg, "*** ascend the search tree ***\n");
      write(fd, msg, strlen(msg));
#endif
#ifdef DEBUG
      printf("*** ascend the search tree ***\n");
#endif
    } else 
      Delta_[0] = -1;
#ifdef DEBUG
    still();
#endif
  } while (Delta_[0] != -1);
#ifdef EXPERIMENT
  sprintf(msg, "%d regular trinagulations.\n", numof_reg);
  write(fd, msg, strlen(msg));
#endif

  for (i = 0; i < numof_rows * numof_columns; i++)
    bdestroy(&matrix[i]);
  for (i = 0; i < codimension * numof_columns; i++)
    bdestroy(&gale[i]);
  for (i = 0; i < (sizecand + 2) * (codimension + 1); i++)
    bdestroy(&tableau[i]);
  for (i = 0; i < codimension * (codimension + 1); i++) {
    bdestroy(&initial_mtrx[i]);
    bdestroy(&initial_glmtrx[i]);
  }
  for (i = 0; i < dimension * dimension; i++)
    bdestroy(&volume_mtrx[i]);
  for (i = 0; i < (codimension - 1) * codimension; i++) {
    bdestroy(&normal_mtrx[i]);
    bdestroy(&normal_glmtrx[i]);
  }
  for (i = 0; i < dimension * (dimension + 1); i++) {
    bdestroy(&signed_radon_mtrx[i]);
    bdestroy(&signed_radon_glmtrx[i]);
  }
  for (i = 0; i < sizecand * codimension; i++) {
    bdestroy(&set1of_normals[i]);
    bdestroy(&set2of_normals[i]);
  }
  for (i = 0; i < numof_columns; i++) {
    bdestroy(&volvec1[i]);
    bdestroy(&volvec2[i]);
  }
  for (i = 0; i <= codimension; i++) 
    bdestroy(&equation[i]);
  for (i = 0; i < codimension; i++) {
    bdestroy(&vector1[i]);
    bdestroy(&vector2[i]);
    bdestroy(&vector3[i]);
    bdestroy(&vector4[i]);
  }

  free(matrix);
  free(gale);
  free(tableau);
  free(initial_mtrx);
  free(initial_glmtrx);
  free(volume_mtrx);
  free(normal_mtrx);
  free(signed_radon_mtrx);
  free(signed_radon_glmtrx);
  free(set1of_normals);
  free(set2of_normals);
  free(volvec1);
  free(volvec2);
  free(equation);
  free(vector1);
  free(vector2);
  free(vector3);
  free(vector4);

  free(global_Gamma);
  free(global_Gamma_);
  free(global_Delta_);

  free(Delta);
  free(Delta_);
  free(Delta0);
  free(Opt);

  free(msg);

  watch_see();
  return;
}



int main(argc, argv)
int argc;
char *argv[];
{
  int i;
  char msg[BUFSIZ];

  bcreate(&BZERO);
  bint(&BZERO, 0);
  sumof_time  = 0;

#ifndef EXPERIMENT
  numof_reg = 1;
  if (argc == 5) {
    if ((fd__ = open(argv[3], O_RDONLY)) == -1) {
      fprintf(stderr, "Can't open %s.\n", argv[4]);
      exit(1);
    }
    sscanf(argv[4], "%d", &numof_reg);
  }

  if (argc < 3) {
    fprintf(stderr, "Few arguments.\n");
    exit(1);
  }
  load(argv[1], &matrix, &numof_rows, &numof_columns);

  if ((fd = open(argv[2], O_WRONLY | O_CREAT, 0644)) == -1) {
    fprintf(stderr, "Can't open %s.\n", argv[2]);
    exit(1);
  }
#endif
#ifdef EXPERIMENT
  numof_reg = 1;
  if (argc < 4) {
    fprintf(stderr, "Few arguments.\n");
    exit(1);
  }

  if ((fd = open(argv[3], O_WRONLY | O_CREAT, 0644)) == -1) {
    fprintf(stderr, "Can't open %s.\n", argv[2]);
    exit(1);
  }
  sscanf(argv[1], "%d", &numof_rows);
  sscanf(argv[2], "%d", &numof_columns);
#endif

#ifdef EXPERIMENT
  for (i = 0; i < 20; i++) {
    random_matrix();
#endif
    if (argc == 5) 
      regular_triangulations(1);
    else
      regular_triangulations(0);
#ifdef EXPERIMENT
    numof_reg = 1;
  }
#endif

#ifdef EXPERIMENT
  sprintf(msg, "Total: %d msc.\n%d msc. per reg. tri.\n", sumof_time, (long) sumof_time / 20);
  write(fd, msg, strlen(msg));
#endif  

  free(msg);
  bdestroy(&BZERO);
  close(fd);
  exit(0);
  return 0;
}