polyparam.c

Go to the documentation of this file.

/***********************************************************************/
/*                Parametrized polyhedra V4.20                         */
/*                copyright 1995-2000 Vincent Loechner                 */
/*                copyright 1996-1997, Doran Wilde                     */
/***********************************************************************/
 
/********************* -----------USER #DEFS-------- ***********************/
/* These are mainly for debug purposes. You shouldn't need to change       */
/* anything for daily usage...                                             */
/***************************************************************************/
 
/* you may define each macro independently */
/* #define DEBUGPP      */
/* #define DEBUGPP3     */ /* initialization of domain, context, ... */
/* #define DEBUGPP31    */ /* even more init-domains */
/* #define DEBUGPP32    */ /* even even more... (Elim_Columns) */
/* #define DEBUGPP4     */ /* m-faces scan */
/* #define DEBUGPP41    */ /* inverse Di in scan */
/* #define DEBUGPP5     */ /* Compute_PDomains */
/********************* ---------END USER #DEFS------ ***********************/
 
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifdef DEBUGPP
#include <time.h>
#endif
 
#include <polylib/polylib.h>
 
static void traite_m_face(Polyhedron *, unsigned int *, unsigned int *);
static void scan_m_face(int, int, Polyhedron *, unsigned int *);
 
/*
 * Return the intersection of two polyhedral domains 'Pol1' and 'Pol2' such
 * that if the intersection is a polyhedron of lower dimension (a degenerate
 * polyhedron) than the operands, it is discarded from the resulting polyhedra
 * list.
 */
Polyhedron *PDomainIntersection(Polyhedron *Pol1, Polyhedron *Pol2,
                                unsigned NbMaxRays) {
 
  Polyhedron *p1, *p2, *p3, *d;
 
  if (!Pol1 || !Pol2)
    return (Polyhedron *)0;
  if ((Pol1->Dimension != Pol2->Dimension) || (Pol1->NbEq != Pol2->NbEq)) {
    fprintf(stderr,
            "? PDomainIntersection: operation on different dimensions\n");
    return (Polyhedron *)0;
  }
 
  POL_ENSURE_FACETS(Pol1);
  POL_ENSURE_VERTICES(Pol1);
  POL_ENSURE_FACETS(Pol2);
  POL_ENSURE_VERTICES(Pol2);
 
  d = (Polyhedron *)0;
  for (p1 = Pol1; p1; p1 = p1->next) {
    for (p2 = Pol2; p2; p2 = p2->next) {
      p3 = AddConstraints(p2->Constraint[0], p2->NbConstraints, p1, NbMaxRays);
      if (!p3)
        continue;
 
      /* If the new polyhedron 'p3' has lower dimension, discard it */
      if (p3->NbEq != Pol1->NbEq)
        Polyhedron_Free(p3);
 
      /* Otherwise add it to the new polyhderal domain 'd'. */
      else
        d = AddPolyToDomain(p3, d);
    }
  }
  return d;
} /* PDomainIntersection */
 
/*
 * Given polyhderal domains 'Pol1' and 'Pol2', return the difference of the
 * two domains with a modification that the resulting polyhedra in the new
 * domain don't have a 1 unit space around cut and the degenerate results
 * (of smaller dimension) are discarded.
 */
Polyhedron *PDomainDifference(Polyhedron *Pol1, Polyhedron *Pol2,
                              unsigned NbMaxRays) {
 
  Polyhedron *p1, *p2, *p3, *d;
  int i;
 
  if (!Pol1 || !Pol2)
    return (Polyhedron *)0;
  if ((Pol1->Dimension != Pol2->Dimension) || (Pol1->NbEq != Pol2->NbEq)) {
    fprintf(stderr, "? PDomainDifference: operation on different dimensions\n");
    return (Polyhedron *)0;
  }
 
  POL_ENSURE_FACETS(Pol1);
  POL_ENSURE_VERTICES(Pol1);
  POL_ENSURE_FACETS(Pol2);
  POL_ENSURE_VERTICES(Pol2);
 
  d = (Polyhedron *)0;
  for (p2 = Pol2; p2; p2 = p2->next) {
    for (p1 = Pol1; p1; p1 = p1->next) {
      for (i = 0; i < p2->NbConstraints; i++) {
 
        /* Add the constraint (-p2->Constraint[i]) >= 0 in 'p1' */
        p3 = SubConstraint(p2->Constraint[i], p1, NbMaxRays, 2);
        if (!p3)
          continue;
 
        /* If the new polyhedron 'p3' is empty or is a polyhedron of lower */
        /* dimension, discard it.                                          */
        if (emptyQ(p3) || p3->NbEq != Pol1->NbEq)
          Polyhedron_Free(p3);
 
        /* Otherwise add 'p3' to the new polyhderal domain 'd' */
        else
          d = AddPolyToDomain(p3, d);
      }
    }
    if (p2 != Pol2)
      Domain_Free(Pol1);
    Pol1 = d;
    d = (Polyhedron *)0;
  }
  return Pol1;
} /* PDomainDifference */
 
/*
 * Return 1 if matrix 'Mat' is full column ranked, otherwise return 0.
 */
static int TestRank(Matrix *Mat) {
 
  int i, j, k;
  Value m1, m2, m3, gcd, tmp;
 
  /* Initialize all the 'Value' variables */
  value_init(m1);
  value_init(m2);
  value_init(m3);
  value_init(gcd);
  value_init(tmp);
 
  for (k = 0; k < Mat->NbColumns; ++k) {
 
    /* If the digonal entry (k,k) is zero, search down the column(k) */
    /* starting from row(k) to find a non-zero entry                 */
    if (value_zero_p(Mat->p[k][k])) {
      for (j = k + 1; j < Mat->NbRows; ++j) {
 
        /* If a non-zero entry (j,k) is found */
        if (value_notzero_p(Mat->p[j][k])) {
 
          /* Exchange row(k) and row(j) */
          for (i = k; i < Mat->NbColumns; ++i) {
            value_assign(tmp, Mat->p[j][i]);
            value_assign(Mat->p[j][i], Mat->p[k][i]);
            value_assign(Mat->p[k][i], tmp);
          }
          break;
        }
      }
 
      /* If no non-zero entry is found then the matrix 'Mat' is not full */
      /* ranked. Return zero.                                            */
      if (j >= Mat->NbRows) {
 
        /* Clear all the 'Value' variables */
        value_clear(m1);
        value_clear(m2);
        value_clear(m3);
        value_clear(gcd);
        value_clear(tmp);
        return 0;
      }
    }
 
    /* Now Mat[k][k] is the pivot element */
    for (j = k + 1; j < Mat->NbRows; ++j) {
 
      /* Make every other entry (below row(k)) in column(k) zero */
      value_gcd(gcd, Mat->p[j][k], Mat->p[k][k]);
      for (i = k + 1; i < Mat->NbColumns; ++i) {
 
        /* pour tous les indices i > k */
        value_multiply(m1, Mat->p[j][i], Mat->p[k][k]);
        value_multiply(m2, Mat->p[j][k], Mat->p[k][i]);
        value_subtract(m3, m1, m2);
        value_division(Mat->p[j][i], m3, gcd);
      }
    }
  }
 
  /* Clear all the 'Value' variables */
  value_clear(m1);
  value_clear(m2);
  value_clear(m3);
  value_clear(gcd);
  value_clear(tmp);
 
  /* The matrix 'Mat' is full ranked, return 1 */
  return 1;
} /* TestRank */
 
/*
 * The Saturation matrix is defined to be an integer (int type) matrix. It is
 * a boolean matrix which has a row for every constraint and a column for
 * every line or ray. The bits in the binary format of each integer in the
 * saturation matrix stores the information whether the corresponding constr-
 * aint is saturated by ray(line) or not.
 */
typedef struct {
  unsigned int NbRows;
  unsigned int NbColumns;
  unsigned int **p;
  unsigned int *p_init;
} SatMatrix;
 
static SatMatrix *SMAlloc(int rows, int cols) {
 
  unsigned int **q, *p;
  int i;
 
  SatMatrix *result = (SatMatrix *)malloc(sizeof(SatMatrix));
  assert(result != NULL);
 
  result->NbRows = rows;
  result->NbColumns = cols;
  result->p = q = (unsigned int **)malloc(rows * sizeof(unsigned int *));
  assert(result->p != NULL);
  result->p_init = p =
      (unsigned int *)malloc(rows * cols * sizeof(unsigned int));
  assert(result->p_init != NULL);
 
  for (i = 0; i < rows; i++) {
    *q++ = p;
    p += cols;
  }
 
  return result;
} /* SMAlloc */
 
#if defined(DEBUGPP4)
static void SMPrint(SatMatrix *matrix) {
 
  unsigned int *p;
  int i, j;
  unsigned NbRows, NbColumns;
 
  fprintf(stderr, "%d %d\n", NbRows = matrix->NbRows,
          NbColumns = matrix->NbColumns);
  for (i = 0; i < NbRows; i++) {
    p = *(matrix->p + i);
    for (j = 0; j < NbColumns; j++)
      fprintf(stderr, " %10X ", *p++);
    fprintf(stderr, "\n");
  }
} /* SMPrint */
#endif
 
static void SMFree(SatMatrix *matrix) {
 
  free((char *)matrix->p_init);
  free((char *)matrix->p);
  free((char *)matrix);
  return;
} /* SMFree */
 
/* -------------------------------------------------------------------------
 * Shared Global Variables:
 * Used by procedures: Find_m_face, scan_m_face, Poly2Sat, traite_m_face,
 *                     count_sat
 * -------------------------------------------------------------------------
 */
static int m;     /* number of parameters */
static int m_dim; /* dimension of m-face */
static int n;     /* dimension (not including parameters) */
static int ws;    /* Working Space size */
static int nr;    /* (NbRays-1)/32 + 1 */
 
static Polyhedron *CEqualities; /* Equalities in the context */
static SatMatrix *Sat;          /* Saturation Matrix (row=constraint, col=ray)*/
static unsigned int *egalite;   /* Bool vector marking constraints in m-face  */
static Matrix *Xi, *Pi;         /* Xi and Pi */
static Matrix *PiTest;          /* Matrix used to test if Pi is full ranked? */
static Matrix *CTest;
static Matrix *PiInv;  /* Matrix inverse Pi, with the last col of   */
                       /* each line = denominator of the line       */
static Matrix *RaysDi; /* Constraint matrix for computing Di */
 
static int KD; /* Flag : keep the full domains in memory ? */
               /* 1 = yes; 0 = no, keep constraints only   */
 
static int nbPV;                  /* The number of parameterized vertices */
static Param_Vertices *PV_Result; /* List of parameterized vertices */
static Param_Domain *PDomains;    /* List of domains. */
 
#ifdef DEBUGPP
static int nbfaces;
#endif
 
/*
 * Add the constraints from the context polyhedron 'CEqualities' to the
 * constraints of polyhedra in the polyhedral domain 'D' and return the new
 * polyhedral domain. Polyhedral domain 'D' is subsequently deleted from memory
 */
static Polyhedron *Add_CEqualities(Polyhedron *D) {
 
  Polyhedron *d, *r, *tmp;
 
  if (!CEqualities)
    return D;
  else {
    if (!D || emptyQ(D)) {
      if (D)
        Domain_Free(D);
      return (Polyhedron_Copy(CEqualities));
    }
    r = AddConstraints(D->Constraint[0], D->NbConstraints, CEqualities, ws);
    tmp = r;
    for (d = D->next; d; d = d->next) {
      tmp->next =
          AddConstraints(d->Constraint[0], d->NbConstraints, CEqualities, ws);
      tmp = tmp->next;
    }
    Domain_Free(D);
    return (r);
  }
} /* Add_CEqualities */
 
#define INT_BITS (sizeof(unsigned) * 8)
 
unsigned int *int_array2bit_vector(unsigned int *array, int n) {
  int i, ix;
  unsigned bx;
  int words = (n + INT_BITS - 1) / INT_BITS;
  unsigned int *bv = (unsigned int *)calloc(words, sizeof(unsigned));
  assert(bv);
  for (i = 0, ix = 0, bx = MSB; i < n; ++i) {
    if (array[i])
      bv[ix] |= bx;
    NEXT(ix, bx);
  }
  return bv;
}
 
/*----------------------------------------------------------------------*/
/* traite_m_face                                                        */
/*       Given an m-face, compute the parameterized vertex              */
/* D       - The entire domain                                          */
/* mf      - Bit vector marking the lines/rays in the m-face            */
/* egalite - boolean vector marking saturated constraints in m-face     */
/*----------------------------------------------------------------------*/
static void traite_m_face(Polyhedron *D, unsigned int *mf,
                          unsigned int *egalite) {
  Matrix *Si;      /* Solution */
  Polyhedron *PDi; /* polyhedron Di */
  Param_Vertices *PV;
  int j, k, c, r;
  unsigned kx, bx;
 
#ifdef DEBUGPP
  ++nbfaces;
#endif
 
  /* Extract  Xi, Pi, and RaysDi from D */
  RaysDi->NbRows = 0;
  for (k = 0, c = 0, kx = 0, bx = MSB; k < D->NbRays; ++k) {
    if (mf[kx] & bx) { /* this ray is in the current m-face */
      if (c < m + 1) {
        int i;
 
        /* tester si cette nouvelle colonne est lin. indep. des autres */
        /* i.e. si gauss ne donne pas de '0' sur la colonne Pi */
        /* jusqu'a l'indice 'c'                                */
 
        /* construit PiTest */
        for (j = 0; j < m + 1; ++j) {
          for (i = 0; i < c; ++i)
 
            /* les c premieres colonnes */
            value_assign(PiTest->p[j][i], Pi->p[j][i]);
 
          /* la nouvelle */
          value_assign(PiTest->p[j][c], D->Ray[k][j + 1 + n]);
        }
        PiTest->NbColumns = c + 1;
        r = TestRank(PiTest);
        if (r /* TestRank(PiTest) */) {
 
          /* Ok, c'est lin. indep. */
          for (j = 0; j < n; j++)
            value_assign(Xi->p[j][c], D->Ray[k][j + 1]); /* Xi */
          for (j = 0; j < m; j++)
            value_assign(Pi->p[j][c], D->Ray[k][j + 1 + n]); /* Pi */
          value_assign(Xi->p[n][c], D->Ray[k][n + m + 1]);   /* const */
          value_assign(Pi->p[m][c], D->Ray[k][n + m + 1]);   /* const */
          c++;
        }
      }
 
      /* Status bit */
      value_assign(RaysDi->p[RaysDi->NbRows][0], D->Ray[k][0]);
      Vector_Copy(&D->Ray[k][n + 1], &RaysDi->p[RaysDi->NbRows][1], (m + 1));
      ++RaysDi->NbRows;
    }
    NEXT(kx, bx);
  }
 
#ifdef DEBUGPP41
  fprintf(stderr, "\nRaysDi=\n");
  Matrix_Print(stderr, P_VALUE_FMT, RaysDi);
  if (c < m + 1)
    fprintf(stderr, "Invalid ");
  fprintf(stderr, "Pi=\n");
  Matrix_Print(stderr, P_VALUE_FMT, Pi);
#endif
 
#ifdef DEBUGPP4
  if (c < m + 1)
    fprintf(stderr, "Eliminated because of no vertex\n");
#endif
 
  if (c < m + 1)
    return;
 
    /* RaysDi->numColumns = m+2; */ /* stays the same */
 
    /*  Xi->NbColumns = m+1;*/ /* VIN100: stays the same. was 'c'! */
    /*  Xi->NbRows = n+1; */   /* stays the same */
    /*  Pi->NbColumns = m+1;*/ /* VIN100: stays the same. was 'c'! */
    /*  Pi->NbRows = m+1; */   /* stays the same */
 
#ifdef DEBUGPP4
  fprintf(stderr, "Xi = ");
  Matrix_Print(stderr, P_VALUE_FMT, Xi);
  fprintf(stderr, "Pi = ");
  Matrix_Print(stderr, P_VALUE_FMT, Pi);
#endif
 
  /* (Right) invert Pi if POSSIBLE, if not then next m-face */
  /* Pi is destroyed                                        */
  if (!MatInverse(Pi, PiInv)) {
 
#ifdef DEBUGPP4
    fprintf(stderr, "Eliminated because of no inverse Pi\n");
#endif
 
    return;
  }
 
#ifdef DEBUGPP4
  fprintf(stderr, "FACE GENERATED!\n");
  fprintf(stderr, "PiInv = ");
  Matrix_Print(stderr, P_VALUE_FMT, PiInv);
#endif
 
  /* Compute  Si (now called Ti in the paper) */
  Si = Matrix_Alloc(Xi->NbRows, PiInv->NbColumns);
  rat_prodmat(Si, Xi, PiInv);
 
#ifdef DEBUGPP4
  fprintf(stderr, "Si = ");
  Matrix_Print(stderr, P_VALUE_FMT, Si);
#endif
 
  Si->NbRows--; /* throw out the last row = 0 ... 0 1 */
 
  /* Copy all of that into the PV structure */
  PV = (Param_Vertices *)malloc(sizeof(Param_Vertices));
  PV->next = PV_Result;
  PV->Vertex = Si;
  PV->Domain = NULL;
  PV->Facets = int_array2bit_vector(egalite, D->NbConstraints);
  PV_Result = PV;
  nbPV++; /* increment vertex count */
 
  /* Ok... now compute the parameter domain */
  PDi = Rays2Polyhedron(RaysDi, ws);
 
#ifdef DEBUGPP3
  fprintf(stderr, "RaysDi = ");
  Matrix_Print(stderr, P_VALUE_FMT, RaysDi);
  fprintf(stderr, "PDi = ");
  Polyhedron_Print(stderr, P_VALUE_FMT, PDi);
#endif
 
  if (KD == 0) {
 
    /* Add the equalities again to the domain */
    PDi = Add_CEqualities(PDi);
    PV->Domain = Polyhedron2Constraints(PDi);
    Polyhedron_Free(PDi);
  } else {
    Param_Domain *PD;
    PD = (Param_Domain *)malloc(sizeof(Param_Domain));
    PD->Domain = PDi;
    PD->F = NULL;
    PD->next = PDomains;
    PDomains = PD;
  }
  return;
} /* traite_m_face */
 
/*----------------------------------------------------------------------*/
/* count_sat                                                            */
/*      count the number of saturated rays in the bit vector mf         */
/*      Uses nr from global area                                        */
/*----------------------------------------------------------------------*/
int cntbit[256] = {/* counts for 8 bits */
                   0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
                   1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
                   1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
                   2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
 
                   1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
                   2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
                   2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
                   3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
 
                   1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
                   2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
                   2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
                   3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
 
                   2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
                   3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
                   3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
                   4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8};
 
static int count_sat(unsigned int *mf) {
 
  register unsigned int i, tmp, cnt = 0;
 
  for (i = 0; i < nr; i++) {
    tmp = mf[i];
    cnt = cnt + cntbit[tmp & 0xff] + cntbit[(tmp >> 8) & 0xff] +
          cntbit[(tmp >> 16) & 0xff] + cntbit[(tmp >> 24) & 0xff];
  }
  return cnt;
} /* count_sat */
 
/* Returns true if all bits in part are also set in bv */
static int bit_vector_includes(unsigned int *bv, int len, unsigned int *part) {
  int j;
 
  for (j = 0; j < len; j++) {
#ifdef DEBUGPP4
    fprintf(stderr, "mf=%08X Sat=%08X &=%08X\n", part[j], bv[j],
            (part[j] & bv[j]));
#endif
    if ((part[j] & bv[j]) != part[j])
      return 0;
  }
  return 1;
}
 
/*----------------------------------------------------------------------*/
/* let D + E + L be the dimension of the polyhedron                     */
/* D = Dimension of polytope (ray space)                                */
/* L = Dimension of Lineality space (number of lines, bid)              */
/* E = Dimension of Affine hull (number of equations)                   */
/* n = number of data indices                                           */
/* m = number of parameters                                             */
/* full domain:                                                         */
/*     n + m = D + E + L                                                */
/* projected domains:                                                   */
/*     m = D_m + E_m + L_m                                              */
/*     n = D_n + E_n + L_n                                              */
/* What dimension M-face, when projected onto parameter space,          */
/* will give an L_m-face?                                               */
/* What are the conditions?                                             */
/*   - at least one vertex                                              */
/*   - number of rays >= D_m+1 after removal of redundants              */
/*                                                                      */
/* dim of face    nb saturated constraints   nb saturated lines,rays    */
/* -----------    ------------------------   -----------------------    */
/* (0+L)-face     all E eqns + >=D ineq      all L lines + 1 ray        */
/* (M+L)-face     all E eqns + >=(D-M) ineq  all L lines + >=(M+1) rays */
/* (D+L)-face     all E eqns + 0 ineq        all L lines + >=(D+1) rays */
/*----------------------------------------------------------------------*/
/*----------------------------------------------------------------------*/
/* scan_m_face                                                          */
/*      pos : the next candidate constraint position                    */
/*    nb_un : number of saturated constraints needed to finish a face   */
/*        D : the source polyhedron (context included )                 */
/*       mf : bit-array marking rays which are saturated so far         */
/* From Global area:                                                    */
/* ----------------                                                     */
/*        n : number of data indices                                    */
/*        m : number of parameters                                      */
/*  egalite : boolean vector marking saturated constraints in m-face    */
/*      Sat : Saturation Matrix (row=constraints, col=rays)             */
/*       ws : working space size                                        */
/*       nr : (NbRays-1)/32 + 1                                         */
/*                                                                      */
/* Recursive function to find the rays and vertices of each m-face      */
/*----------------------------------------------------------------------*/
static void scan_m_face(int pos, int nb_un, Polyhedron *D, unsigned int *mf) {
  /* pos   - the next candidate constraint position */
  /* nb_un - the number of constraints needed to finish a face */
  /* D     - the source polyhedron */
  /* mf    - (bit vector) marks rays that are saturated so far */
 
  unsigned int *new_mf;
 
#ifdef DEBUGPP4
  fprintf(stderr, "Start scan_m_face(pos=%d, nb_un=%d, n=%d, m=%d\n", pos,
          nb_un, n, m);
  fprintf(stderr, "mf = ");
  {
    int i;
    for (i = 0; i < nr; i++)
      fprintf(stderr, "%08X", mf[i]);
    fprintf(stderr, "\nequality = [");
    for (i = 0; i < D->NbConstraints; i++)
      fprintf(stderr, " %1d", egalite[i]);
    fprintf(stderr, "]\n");
  }
#endif
 
  if (nb_un == 0) { /* Base case */
    int i;
 
    /*********** ELIMINATION OF REDUNDANT FACES ***********/
    /* if all these vertices also verify a previous constraint */
    /* which is NOT already selected, we eliminate this face */
    /* This keeps the lexicographically greatest selection */
    for (i = 0; i < pos - 1; i++) {
      if (egalite[i])
        continue; /* already selected */
 
        /* if Sat[i] & mf == mf then it's redundant */
#ifdef DEBUGPP4
      fprintf(stderr, "Sat[%d]\n", i);
#endif
      if (bit_vector_includes(Sat->p[i], nr, mf)) {
#ifdef DEBUGPP4
        fprintf(stderr, "Redundant with constraint %d\n", i);
#endif
        return; /* it is redundant */
      }
    }
    /********* END OF ELIMINATION OF DEGENERATE FACES *********/
    /* Now check for other constraints that are verified */
    for (i = pos; i < D->NbConstraints; ++i) {
      if (bit_vector_includes(Sat->p[i], nr, mf))
        egalite[i] = 1;
    }
    /* if we haven't found a constraint verified by all */
    /* the rays, its OK, it's a new face. */
    traite_m_face(D, mf, egalite);
    for (i = pos; i < D->NbConstraints; ++i)
      egalite[i] = 0;
    return;
  }
 
  /* See if there are enough constraints left to finish */
  if ((pos + nb_un) > D->NbConstraints)
    return;
 
  /* Recurring part of the procedure */
  /* Add the pos'th constraint, compute new saturation vector */
  {
    int k;
    new_mf = (unsigned int *)malloc(nr * sizeof(unsigned int));
    for (k = 0; k < nr; k++)
      new_mf[k] = mf[k] & Sat->p[pos][k];
  }
#ifdef DEBUGPP4
  fprintf(stderr, "new_mf = ");
  {
    int i;
    for (i = 0; i < nr; i++) {
      fprintf(stderr, "%08X", new_mf[i]);
    }
    fprintf(stderr, "\ncount(new_mf) = %d\n", count_sat(new_mf));
  }
#endif
 
  {
    int c;
    c = count_sat(new_mf);
    /* optimization : at least m_dim+1 rays must be saturated to add this
     * constraint */
    if (c > m_dim) {
      int redundant = 0;
 
      egalite[pos] = 1; /* Try it with the pos-th constraint */
 
      /* If this constraint does not change anything,
       * it is redundant with respect to the selected
       * equalities and the remaining inequalities.
       * Check whether it is redundant with respect
       * to just the selected equalities.
       */
      if (c == count_sat(mf)) {
        int i, c, j;
 
        for (i = 0, c = 0; i < D->NbConstraints; ++i) {
          if (egalite[i] == 0 || egalite[i] == -1)
            continue;
          for (j = 0; j < D->Dimension + 1; ++j)
            value_assign(CTest->p[j][c], D->Constraint[i][j + 1]);
          ++c;
        }
        CTest->NbColumns = c;
#ifdef DEBUGPP41
        Matrix_Print(stderr, P_VALUE_FMT, CTest);
#endif
        redundant = !TestRank(CTest);
      }
 
      /* Do not decrement nb_un if equality is redundant. */
      if (redundant) {
        egalite[pos] = -1; /* Don't use in further redundance test
                            */
        scan_m_face(pos + 1, nb_un, D, new_mf);
      } else {
        scan_m_face(pos + 1, nb_un - 1, D, new_mf);
      }
    }
  }
  free(new_mf);
  egalite[pos] = 0; /* Try it without the pos-th constraint */
  if ((pos + nb_un) >= D->NbConstraints)
    return;
  scan_m_face(pos + 1, nb_un, D, mf);
  return;
} /* scan_m_face */
 
/*
 * Create a saturation matrix with rows correspond to the constraints and
 * columns correspond to the rays of the polyhedron 'Pol'. Global variable
 * 'nr' is set in the function.
 */
static SatMatrix *Poly2Sat(Polyhedron *Pol, unsigned int **L) {
 
  SatMatrix *Sat;
  int i, j, k, kx;
  unsigned int *Temp;
  Value *p1, *p2, p3, tmp;
  unsigned Dimension, NbRay, NbCon, bx;
 
  /* Initialize all the 'Value' variables */
  value_init(p3);
  value_init(tmp);
 
  NbRay = Pol->NbRays;
  NbCon = Pol->NbConstraints;
  Dimension = Pol->Dimension + 1; /* Homogeneous Dimension */
 
  /* Build the Sat matrix */
  nr = (NbRay - 1) / (sizeof(int) * 8) + 1; /* Set globally */
  Sat = SMAlloc(NbCon, nr);
  Temp = (unsigned int *)malloc(nr * sizeof(unsigned int));
  memset(Sat->p_init, 0, nr * NbCon * sizeof(int));
  memset(Temp, 0, nr * sizeof(unsigned int));
  kx = 0;
  bx = MSB;
  for (k = 0; k < NbRay; k++) {
    for (i = 0; i < NbCon; i++) {
      p1 = &Pol->Constraint[i][1];
      p2 = &Pol->Ray[k][1];
      value_set_si(p3, 0);
      for (j = 0; j < Dimension; j++) {
        value_multiply(tmp, *p1, *p2);
        value_addto(p3, p3, tmp);
        p1++;
        p2++;
      }
      if (value_zero_p(p3))
        Sat->p[i][kx] |= bx;
    }
    Temp[kx] |= bx;
    NEXT(kx, bx);
  }
 
  /* Set 'L' to an array containing ones in every bit position of its */
  /* elements.                                                        */
  *L = Temp;
 
  /* Clear all the 'Value' variables */
  value_clear(p3);
  value_clear(tmp);
 
  return Sat;
} /* Poly2Sat */
 
/*
 * Create a parametrized polyhedron with zero parameters. This function was
 * first written by Xavier Redon, and was later modified by others.
 */
Param_Polyhedron *GenParamPolyhedron(Polyhedron *Pol, Matrix *Rays) {
  Param_Polyhedron *result;
  int nbRows, nbColumns;
  int i, size, rays;
 
  nbRows = Pol->NbRays;
  nbColumns = Pol->Dimension + 2;
 
  /* Count the number of rays */
  for (i = 0, rays = 0; i < nbRows; i++)
    if (value_notone_p(Pol->Ray[i][0]) ||
        value_zero_p(Pol->Ray[i][nbColumns - 1]))
      ++rays;
 
  /* Initialize the result */
  result = (Param_Polyhedron *)malloc(sizeof(Param_Polyhedron));
  result->nbV = nbRows - rays;
  result->V = NULL;
  result->Constraints = Polyhedron2Constraints(Pol);
  result->Rays = Rays;
 
  /* Build the parametric vertices */
  for (i = 0; i < nbRows; i++) {
    Matrix *vertex;
    Param_Vertices *paramVertex;
    int j;
 
    if (value_notone_p(Pol->Ray[i][0]) ||
        value_zero_p(Pol->Ray[i][nbColumns - 1]))
      continue;
 
    vertex = Matrix_Alloc(nbColumns - 2, 2);
    for (j = 1; j < nbColumns - 1; j++) {
      value_assign(vertex->p[j - 1][0], Pol->Ray[i][j]);
      value_assign(vertex->p[j - 1][1], Pol->Ray[i][nbColumns - 1]);
    }
    paramVertex = (Param_Vertices *)malloc(sizeof(Param_Vertices));
    paramVertex->Vertex = vertex;
 
    /* There is one validity domain : universe of dimension 0 */
    paramVertex->Domain = Matrix_Alloc(1, 2);
    value_set_si(paramVertex->Domain->p[0][0], 1);
    value_set_si(paramVertex->Domain->p[0][1], 1);
    paramVertex->Facets = NULL;
    paramVertex->next = result->V;
    result->V = paramVertex;
  }
 
  /* Build the parametric domains (only one here) */
  if (nbRows > 1)
    size = (nbRows - 1) / (8 * sizeof(int)) + 1;
  else
    size = 1;
  result->D = (Param_Domain *)malloc(sizeof(Param_Domain));
  result->D->next = NULL;
  result->D->Domain = Universe_Polyhedron(0);
  result->D->F = (unsigned int *)malloc(size * sizeof(int));
  memset(&result->D->F[0], 0xFF, size * sizeof(int));
 
  return result;
} /* GenParamPolyhedron */
 
/*----------------------------------------------------------------------*/
/* PreElim_Columns                                                      */
/* function being called before Elim_Columns                            */
/* Equalities in E are analysed to initialize ref and p.                */
/* These two vectors are used to construct the new constraint matrix    */
/* PreElim_Columns returns the transformation matrix to re-convert the  */
/* resulting domains in the same format (by adding empty columns)       */
/* in the parameter space                                               */
/*----------------------------------------------------------------------*/
Matrix *PreElim_Columns(Polyhedron *E, int *p, int *ref, int m) {
 
  int i, j, l;
  Matrix *T;
 
  /* find which columns to eliminate */
  /* p contains, for each line in E, the column to eliminate */
  /* (i.e. the corresponding parameter index to eliminate) */
  /* 0 <= i <= E->NbEq, and  1 <= p[i] <= E->Dimension */
  memset(p, 0, sizeof(int) * (E->NbEq));
 
#ifdef DEBUGPP32
  fprintf(stderr, "\n\nPreElim_Columns starting\n");
  fprintf(stderr, "E =\n");
  Polyhedron_Print(stderr, P_VALUE_FMT, E);
#endif
 
  for (l = 0; l < E->NbEq; ++l) {
    if (value_notzero_p(E->Constraint[l][0])) {
      fprintf(stderr, "Internal error: Elim_Columns (polyparam.c), expecting "
                      "equalities first in E.\n");
      free(p);
      return (NULL);
    }
    for (i = 1; value_zero_p(E->Constraint[l][i]); ++i) {
      if (i == E->Dimension + 1) {
        fprintf(stderr, "Internal error: Elim_Columns (polyparam.c), expecting "
                        "non-empty constraint in E.\n");
        free(p);
        return (NULL);
      }
    }
    p[l] = i;
 
#ifdef DEBUGPP32
    fprintf(stderr, "p[%d] = %d,", l, p[l]);
#endif
  }
 
  /* Reference vector: column ref[i] in A corresponds to column i in M */
  for (i = 0; i < E->Dimension + 2 - E->NbEq; ++i) {
    ref[i] = i;
    for (j = 0; j < E->NbEq; ++j)
      if (p[j] <= ref[i])
        ref[i]++;
 
#ifdef DEBUGPP32
    fprintf(stderr, "ref[%d] = %d,", i, ref[i]);
#endif
  }
 
  /* Size of T : phdim-nbEq * phdim */
  T = Matrix_Alloc(m + 1 - E->NbEq, m + 1);
  for (i = 0; i < T->NbColumns; i++)
    for (j = 0; j < T->NbRows; j++)
      if (ref[E->Dimension - m + j + 1] == E->Dimension - m + i + 1)
        value_set_si(T->p[j][i], 1);
      else
        value_set_si(T->p[j][i], 0);
 
#ifdef DEBUGPP32
  fprintf(stderr, "\nTransMatrix =\n");
  Matrix_Print(stderr, P_VALUE_FMT, T);
#endif
 
  return (T);
 
} /* PreElim_Columns */
 
/*----------------------------------------------------------------------*/
/* Elim_Columns                                                         */
/* Eliminate columns from A, using the equalities in E.                 */
/* ref and p are precalculated by PreElim_Columns, using E;             */
/* these two vectors are used to construct the new constraint matrix    */
/*----------------------------------------------------------------------*/
Polyhedron *Elim_Columns(Polyhedron *A, Polyhedron *E, int *p, int *ref) {
 
  int i, l, c;
  Matrix *M, *C;
  Polyhedron *R;
  Value tmp1, tmp2;
 
  /* Initialize all the 'Value' variables */
  value_init(tmp1);
  value_init(tmp2);
 
#ifdef DEBUGPP32
  fprintf(stderr, "\nElim_Columns starting\n");
  fprintf(stderr, "A =\n");
  Polyhedron_Print(stderr, P_VALUE_FMT, A);
#endif
 
  /* Builds M = constraint matrix of A, useless columns zeroed */
  M = Polyhedron2Constraints(A);
  for (l = 0; l < E->NbEq; ++l) {
    for (c = 0; c < M->NbRows; ++c) {
      if (value_notzero_p(M->p[c][p[l]])) {
 
        /* A parameter to eliminate here ! */
        for (i = 1; i < M->NbColumns; ++i) {
          value_multiply(tmp1, M->p[c][i], E->Constraint[l][p[l]]);
          value_multiply(tmp2, E->Constraint[l][i], M->p[c][p[l]]);
          value_subtract(M->p[c][i], tmp1, tmp2);
        }
      }
    }
  }
 
#ifdef DEBUGPP32
  fprintf(stderr, "\nElim_Columns after zeroing columns of A.\n");
  fprintf(stderr, "M =\n");
  Matrix_Print(stderr, P_VALUE_FMT, M);
#endif
 
  /* Builds C = constraint matrix, useless columns eliminated */
  C = Matrix_Alloc(M->NbRows, M->NbColumns - E->NbEq);
  for (l = 0; l < C->NbRows; ++l)
    for (c = 0; c < C->NbColumns; ++c) {
      value_assign(C->p[l][c], M->p[l][ref[c]]);
    }
 
#ifdef DEBUGPP32
  fprintf(stderr, "\nElim_Columns after eliminating columns of A.\n");
  fprintf(stderr, "C =\n");
  Matrix_Print(stderr, P_VALUE_FMT, C);
#endif
 
  R = Constraints2Polyhedron(C, ws);
  Matrix_Free(C);
  Matrix_Free(M);
  value_clear(tmp1);
  value_clear(tmp2);
  return (R);
} /* Elim_Columns */
 
static Polyhedron *Recession_Cone(Polyhedron *P, unsigned nvar,
                                  unsigned MaxRays) {
  int i;
  Matrix *M = Matrix_Alloc(P->NbConstraints, 1 + nvar + 1);
  Polyhedron *R;
  for (i = 0; i < P->NbConstraints; ++i)
    Vector_Copy(P->Constraint[i], M->p[i], 1 + nvar);
  R = Constraints2Polyhedron(M, MaxRays);
  Matrix_Free(M);
  return R;
}
 
/* Compute lines/unidirectional rays of the (non parametric) polyhedron */
/* Input :
 *   D1 : combined polyhedron,
 * Output :
 *   Rays : non parametric ray matrix
 *   return value : number of lines
 */
static int ComputeNPLinesRays(int n, Polyhedron *D1, Matrix **Rays) {
  int i, j, nbr, dimfaces;
  Polyhedron *RC; /* Recession Cone */
 
  RC = Recession_Cone(D1, n, ws);
 
  /* get the rays/lines from RC */
  nbr = 0;
  for (i = 0; i < RC->NbRays; i++)
    if (value_zero_p(RC->Ray[i][n + 1]))
      nbr++;
  *Rays = Matrix_Alloc(nbr, n + 2);
  for (i = 0, j = 0; j < nbr; i++)
    if (value_zero_p(RC->Ray[i][n + 1]))
      Vector_Copy(RC->Ray[i], (*Rays)->p[j++], n + 2);
 
  dimfaces = RC->NbBid;
  Polyhedron_Free(RC);
 
#ifdef DEBUGPP31
  fprintf(stderr, "Rays = ");
  Matrix_Print(stderr, P_VALUE_FMT, *Rays);
  fprintf(stderr, "dimfaces = %d\n", dimfaces);
#endif
 
  return dimfaces;
}
 
/*
 * Given a polyhedron 'Di' in combined data and parameter space and a context
 * polyhedron 'C' representing the constraints on the parameter space, create
 * a list of parameterized vertices and assign values to global variables:
 * m,n,ws,Sat,egalite,mf,Xi,Pi,PiInv,RaysDi,CEqualities.
 */
Param_Polyhedron *Find_m_faces(Polyhedron **Di, Polyhedron *C, int keep_dom,
                               int working_space, Polyhedron **CEq,
                               Matrix **CT) {
 
  unsigned int *mf;
  int i, j, dimfaces;
  Polyhedron *D = *Di, *C1, /* true context */
      *D1;                  /* the combined polyhedron, including context C */
  Matrix *Rays,             /* lines/rays (non parametric) */
      *M;
  Param_Polyhedron *res;
  int *p, *ref;
 
  CEqualities = NULL;
 
  if (CT) {
    *CEq = NULL;
    *CT = NULL;
  }
 
  if (!D || !C)
    return (Param_Polyhedron *)0;
 
  ws = working_space;
  m = C->Dimension;
  n = D->Dimension - m;
  if (n < 0) {
    fprintf(stderr, "Find_m_faces: ?%d parameters of a %d-polyhedron !\n", m,
            n);
    return (Param_Polyhedron *)0;
  }
  if (m == 0)
    return GenParamPolyhedron(D, Matrix_Alloc(0, 2));
 
  /* Add constraints from Context to D -> result in D1 */
  C1 = align_context(C, D->Dimension, ws);
 
#ifdef DEBUGPP31
  fprintf(stderr, "m = %d\n", m);
  fprintf(stderr, "D = ");
  Polyhedron_Print(stderr, P_VALUE_FMT, D);
  fprintf(stderr, "C1 = ");
  Polyhedron_Print(stderr, P_VALUE_FMT, C1);
#endif
 
  D1 = DomainIntersection(D, C1, ws);
 
#ifdef DEBUGPP31
  fprintf(stderr, "D1 = ");
  Polyhedron_Print(stderr, P_VALUE_FMT, D1);
#endif
 
  Domain_Free(C1);
 
  if (!D1)
    return (NULL);
  if (emptyQ(D1)) {
    Polyhedron_Free(D1);
    return (NULL);
  }
 
  /* Compute the true context C1 */
  /* M : lines in the direction of the first n indices (index space) */
  M = Matrix_Alloc(n, D1->Dimension + 2);
  for (i = 0; i < n; i++)
    value_set_si(M->p[i][i + 1], 1);
  C1 = DomainAddRays(D1, M, ws);
  Matrix_Free(M);
 
#ifdef DEBUGPP31
  fprintf(stderr, "True context C1 = ");
  Polyhedron_Print(stderr, P_VALUE_FMT, C1);
#endif
 
  dimfaces = ComputeNPLinesRays(n, D1, &Rays);
 
  /* CEqualities contains the constraints (to be added again later) */
  /* *CT is the transformation matrix to add the removed parameters */
  if (!CT) {
    if (C1->NbEq == 0) {
      Polyhedron_Free(C1);
      C1 = NULL;
    } else {
      Polyhedron *CEq1, /* CEqualities, in homogeneous dim */
          *D2;          /* D1, (temporary) simplified */
 
      /* Remove equalities from true context C1 and from D1             */
      /* Compute CEqualities = matrix of equalities in C1, projected in */
      /* the parameter space                                            */
      M = Matrix_Alloc(C1->NbEq, m + 2);
      for (j = 0, i = 0; i < C1->NbEq; ++i, ++j) {
        while (value_notzero_p(C1->Constraint[j][0]))
          ++j;
        value_assign(M->p[i][0], C1->Constraint[j][0]);
        Vector_Copy(&C1->Constraint[j][D->Dimension - m + 1], &M->p[i][1],
                    (m + 1));
      }
      CEqualities = Constraints2Polyhedron(M, ws);
      Matrix_Free(M);
      CEq1 = align_context(CEqualities, D->Dimension, ws);
 
      /* Simplify D1 and C1 (remove the equalities) */
      D2 = DomainSimplify(D1, CEq1, ws);
      Polyhedron_Free(D1);
      Polyhedron_Free(C1);
      Polyhedron_Free(CEq1);
      D1 = D2;
      C1 = NULL;
    }
  } else {            /* if( CT  ) */
    Polyhedron *CEq1, /* CEqualities */
        *D2;          /* D1, (temporary) simplified */
 
    /* Suppress all useless constraints in parameter domain */
    /* when CT is not NULL (ehrhart) */
    /* Vin100, march 01 */
    CEq1 = C1;
    M = Matrix_Alloc(C1->NbConstraints, m + 2);
    for (i = 0; i < C1->NbConstraints; ++i) {
      value_assign(M->p[i][0], C1->Constraint[i][0]);
      Vector_Copy(&C1->Constraint[i][D->Dimension - m + 1], &M->p[i][1],
                  (m + 1));
    }
    CEqualities = Constraints2Polyhedron(M, ws);
    Matrix_Free(M);
 
    D2 = DomainSimplify(D1, CEq1, ws);
    Polyhedron_Free(D1);
    D1 = D2;
    C1 = Universe_Polyhedron(D2->Dimension);
 
    /* if CT is not NULL, the constraints are eliminated                */
    /* *CT will contain the transformation matrix to come back to the   */
    /* original dimension (for a polyhedron, in the parameter space)    */
    if (CEq1->NbEq) {
      m -= CEq1->NbEq;
      p = (int *)malloc(sizeof(int) * (CEq1->NbEq));
    } else
      p = NULL;
    ref = (int *)malloc(sizeof(int) * (CEq1->Dimension + 2 - CEq1->NbEq));
    *CT = PreElim_Columns(CEq1, p, ref, CEqualities->Dimension);
    D2 = Elim_Columns(D1, CEq1, p, ref);
    if (p)
      free(p);
    free(ref);
 
#ifdef DEBUGPP3
    fprintf(stderr, "D2\t Dim = %3d\tNbEq = %3d\tLines = %3d\n", D2->Dimension,
            D2->NbEq, D2->NbBid);
    C2 = Elim_Columns(C1, CEq1, p, ref);
    fprintf(stderr, "C2\t Dim = %3d\tNbEq = %3d\tLines = %3d\n", C2->Dimension,
            C2->NbEq, C2->NbBid);
    Polyhedron_Free(C2);
#endif
 
    Polyhedron_Free(D1);
    Polyhedron_Free(C1);
    D1 = D2;
    C1 = NULL;
    *CEq = CEqualities;
 
#ifdef DEBUGPP3
    fprintf(stderr, "Polyhedron CEq = ");
    Polyhedron_Print(stderr, P_VALUE_FMT, *CEq);
    fprintf(stderr, "Matrix CT = ");
    Matrix_Print(stderr, P_VALUE_FMT, *CT);
#endif
 
    Polyhedron_Free(CEq1);
    CEqualities = NULL; /* don't simplify ! */
 
    /* m changed !!! */
    if (m == 0) {
      /* return the new D1 too */
      *Di = D1;
 
      return GenParamPolyhedron(D1, Rays);
    }
  }
 
#ifdef DEBUGPP3
  fprintf(stderr, "Polyhedron D1 (D AND C) = ");
  Polyhedron_Print(stderr, P_VALUE_FMT, D1);
  fprintf(stderr, "Polyhedron CEqualities = ");
  if (CEqualities)
    Polyhedron_Print(stderr, P_VALUE_FMT, CEqualities);
  else
    fprintf(stderr, "NULL\n");
#endif
 
  KD = keep_dom;
  PDomains = NULL;
  PV_Result = NULL;
  nbPV = 0;
 
  if (emptyQ(D1)) {
    Polyhedron_Free(D1);
    Matrix_Free(Rays);
    return NULL;
  }
 
  /* mf : a bit array indicating which rays are part of the m-face */
  /* Poly2Sat initializes mf to all ones */
  /* set global variable nr to size (number of words) of mf */
  Sat = Poly2Sat(D1, &mf);
 
#ifdef DEBUGPP4
  fprintf(stderr, "Sat = ");
  SMPrint(Sat);
 
  fprintf(stderr, "mf = ");
  for (i = 0; i < nr; i++)
    fprintf(stderr, "%08X", mf[i]);
  fprintf(stderr, "\n");
#endif
 
  /* A boolean array saying which constraints are part of the m-face */
  egalite = (unsigned int *)malloc(sizeof(int) * D1->NbConstraints);
  memset(egalite, 0, sizeof(int) * D1->NbConstraints);
 
  for (i = 0; i < D1->NbEq; i++)
    egalite[i] = 1;
 
  Xi = Matrix_Alloc(n + 1, m + 1);
  Pi = Matrix_Alloc(m + 1, m + 1);
  PiTest = Matrix_Alloc(m + 1, m + 1);
  CTest = Matrix_Alloc(D->Dimension + 1, D->NbConstraints);
  PiInv = Matrix_Alloc(m + 1, m + 2);
  RaysDi = Matrix_Alloc(D1->NbRays, m + 2);
  m_dim = m;
 
  /* m_dim has to be increased by the dimension of the smallest faces
   * of the (non parametric) polyhedron
   */
  m_dim += dimfaces;
 
  /* if the smallest face is of smaller dimension than m_dim,
   * then increase m_dim (I think this should never happen --Vincent)
   */
#ifdef DEBUGPP3
  if (m_dim < D1->NbBid)
    fprintf(stderr, "m_dim (%d) < D1->NbBid (%d)\n", m_dim, D1->NbBid);
#endif
  if (m_dim < D1->NbBid)
    m_dim = D1->NbBid;
 
#ifdef DEBUGPP
  nbfaces = 0;
#endif
#ifdef DEBUGPP3
  fprintf(stderr, "m_dim = %d\n", m_dim);
  fprintf(stderr,
          "Target: find faces that saturate %d constraints and %d rays/lines\n",
          D1->Dimension - m_dim, m_dim + 1);
#endif
 
  /* D1->NbEq constraints already saturated ! */
  scan_m_face(D1->NbEq, (D1->Dimension - m_dim - D1->NbEq), D1, mf);
 
  /* pos, number of constraints needed     */
 
#ifdef DEBUGPP
  fprintf(stderr, "Number of m-faces: %d\n", nbfaces);
#endif
 
  Matrix_Free(RaysDi);
  Matrix_Free(PiInv);
  Matrix_Free(PiTest);
  Matrix_Free(CTest);
  Matrix_Free(Pi);
  Matrix_Free(Xi);
  free(egalite);
  free(mf);
  SMFree(Sat);
 
  /*    if(CEqualities && keep_dom==0) {
           Domain_Free(CEqualities);
        }
  */
 
  res = (Param_Polyhedron *)malloc(sizeof(Param_Polyhedron));
  res->nbV = nbPV;
  res->V = PV_Result;
  res->D = PDomains;
  res->Constraints = Polyhedron2Constraints(D1);
  res->Rays = Rays;
 
  if (CT) /* return the new D1 too ! */
    *Di = D1;
  else
    Domain_Free(D1);
 
  return (res);
} /* Find_m_faces */
 
/*
 * Given parametric domain 'PD' and number of parametric vertices 'nb_domains',
 * find the vertices that belong to distinct sub-domains.
 */
void Compute_PDomains(Param_Domain *PD, int nb_domains, int working_space) {
 
  unsigned bx;
  int i, ix, nv;
  Polyhedron *dx, *d1, *d2;
  Param_Domain *p1, *p2, *p2prev, *PDNew;
 
  if (nb_domains == 0) {
 
#ifdef DEBUGPP5
    fprintf(stderr, "No domains\n");
#endif
 
    return;
  }
 
  /* Already filled out by GenParamPolyhedron */
  if (!PD->next && PD->F)
    return;
 
  /* Initialization */
  nv = (nb_domains - 1) / (8 * sizeof(int)) + 1;
 
#ifdef DEBUGPP5
  fprintf(stderr, "nv = %d\n", nv);
#endif
 
  for (p1 = PD, i = 0, ix = 0, bx = MSB; p1; p1 = p1->next, i++) {
 
    /* Assign a bit array 'p1->F' of suitable size to include the vertices */
    p1->F = (unsigned *)malloc(nv * sizeof(unsigned));
 
    /* Set the bit array to zeros */
    memset(p1->F, 0, nv * sizeof(unsigned));
    p1->F[ix] |= bx; /* Set i'th bit to one */
    NEXT(ix, bx);
  }
 
#ifdef DEBUGPP5
  fprintf(stderr, "nb of vertices=%d\n", i);
#endif
 
  /* Walk the PD list with two pointers */
  ix = 0;
  bx = MSB;
  for (p1 = PD; p1; p1 = p1->next) {
    for (p2prev = p1, p2 = p1->next; p2; p2prev = p2, p2 = p2->next) {
 
      /* Find intersection */
      dx = PDomainIntersection(p1->Domain, p2->Domain, working_space);
 
      if (!dx || emptyQ(dx)) {
 
#ifdef DEBUGPP5
        fprintf(stderr, "Empty dx (p1 inter p2). Continuing\n");
#endif
        if (dx)
          Domain_Free(dx);
        continue;
      }
 
#ifdef DEBUGPP5
      fprintf(stderr, "Begin PDomainDifference\n");
      fprintf(stderr, "p1=");
      Polyhedron_Print(stderr, P_VALUE_FMT, p1->Domain);
      fprintf(stderr, "p2=");
      Polyhedron_Print(stderr, P_VALUE_FMT, p2->Domain);
#endif
      d1 = PDomainDifference(p1->Domain, p2->Domain, working_space);
      d2 = PDomainDifference(p2->Domain, p1->Domain, working_space);
 
#ifdef DEBUGPP5
      fprintf(stderr, "p1\\p2=");
      Polyhedron_Print(stderr, P_VALUE_FMT, d1);
      fprintf(stderr, "p2\\p1=");
      Polyhedron_Print(stderr, P_VALUE_FMT, d2);
      fprintf(stderr, "END PDomainDifference\n\n");
#endif
      if (!d1 || emptyQ(d1) || d1->NbEq != 0) {
 
#ifdef DEBUGPP5
        fprintf(stderr, "Empty d1\n");
#endif
        if (d1)
          Domain_Free(d1);
        Domain_Free(dx);
 
        if (!d2 || emptyQ(d2) || d2->NbEq != 0) {
 
#ifdef DEBUGPP5
          fprintf(stderr, "Empty d2 (deleting)\n");
#endif
          /* dx = p1->Domain = p2->Domain */
          if (d2)
            Domain_Free(d2);
 
          /* Update p1 */
          for (i = 0; i < nv; i++)
            p1->F[i] |= p2->F[i];
 
          /* Delete p2 */
          p2prev->next = p2->next;
          Domain_Free(p2->Domain);
          free(p2->F);
          free(p2);
          p2 = p2prev;
        } else { /* d2 is not empty --> dx==p1->domain */
 
#ifdef DEBUGPP5
          fprintf(stderr, "p2 replaced by d2\n");
#endif
          /* Update p1 */
          for (i = 0; i < nv; i++)
            p1->F[i] |= p2->F[i];
 
          /* Replace p2 with d2 */
          Domain_Free(p2->Domain);
          p2->Domain = d2;
        }
      } else { /* d1 is not empty */
        if (!d2 || emptyQ(d2) || d2->NbEq != 0) {
 
#ifdef DEBUGPP5
          fprintf(stderr, "p1 replaced by d1\n");
#endif
          if (d2)
            Domain_Free(d2);
 
          /* dx = p2->domain */
          Domain_Free(dx);
 
          /* Update p2 */
          for (i = 0; i < nv; i++)
            p2->F[i] |= p1->F[i];
 
          /* Replace p1 with d1 */
          Domain_Free(p1->Domain);
          p1->Domain = d1;
        } else { /*d2 is not empty-->d1,d2,dx are distinct */
 
#ifdef DEBUGPP5
          fprintf(stderr, "Non-empty d1 and d2\nNew node created\n");
#endif
          /* Create a new node for dx */
          PDNew = (Param_Domain *)malloc(sizeof(Param_Domain));
          PDNew->F = (unsigned int *)malloc(nv * sizeof(int));
          memset(PDNew->F, 0, nv * sizeof(int));
          PDNew->Domain = dx;
 
          for (i = 0; i < nv; i++)
            PDNew->F[i] = p1->F[i] | p2->F[i];
 
          /* Replace p1 with d1 */
          Domain_Free(p1->Domain);
          p1->Domain = d1;
 
          /* Replace p2 with d2 */
          Domain_Free(p2->Domain);
          p2->Domain = d2;
 
          /* Insert new node after p1 */
          PDNew->next = p1->next;
          p1->next = PDNew;
        }
      }
    } /* end of p2 scan */
    if (p1->Domain->next) {
      Polyhedron *C = DomainConvex(p1->Domain, working_space);
      Domain_Free(p1->Domain);
      p1->Domain = C;
    }
  } /* end of p1 scan */
} /* Compute_PDomains */
 
/*
 * Given a polyhedron 'Din' in combined data and parametre space, a context
 * polyhedron 'Cin' representing the constraints on the parameter space and
 * a working space size 'working_space', return a parametric polyhedron with
 * a list of parametric vertices and their defining domains.
 */
Param_Polyhedron *Polyhedron2Param_Vertices(Polyhedron *Din, Polyhedron *Cin,
                                            int working_space) {
 
  Param_Polyhedron *result;
 
  POL_ENSURE_FACETS(Din);
  POL_ENSURE_VERTICES(Din);
  POL_ENSURE_FACETS(Cin);
  POL_ENSURE_VERTICES(Cin);
 
#ifdef DEBUGPP
  fprintf(stderr, "Polyhedron2Param_Vertices algorithm starting at : %.2fs\n",
          (float)clock() / CLOCKS_PER_SEC);
#endif
 
  /***************** Scan the m-faces ****************/
  result = Find_m_faces(&Din, Cin, 0, working_space, NULL, NULL);
 
#ifdef DEBUGPP
  fprintf(stderr, "nb of points : %d\n", result->nbV);
#endif
 
#ifdef DEBUGPP
  fprintf(stderr, "end main loop : %.2fs\n", (float)clock() / CLOCKS_PER_SEC);
#endif
 
  return (result);
} /* Polyhedron2Param_Vertices */
 
/*
 * Free the memory allocated to a list of parametrized vertices
 */
void Param_Vertices_Free(Param_Vertices *PV) {
 
  Param_Vertices *next_pv;
 
  while (PV) {
    next_pv = PV->next;
    if (PV->Vertex)
      Matrix_Free(PV->Vertex);
    if (PV->Domain)
      Matrix_Free(PV->Domain);
    if (PV->Facets)
      free(PV->Facets);
    free(PV);
    PV = next_pv;
  }
} /* Param_Vertices_Free */
 
/*
 * Print a list of parametrized vertices *
 */
void Print_Vertex(FILE *DST, Matrix *V, char **param_names) {
  int l, v;
  int first;
  Value gcd, tmp;
 
  value_init(gcd);
  value_init(tmp);
 
  fprintf(DST, "[");
  for (l = 0; l < V->NbRows; ++l) {
 
    /* Variables */
    first = 1;
    fprintf(DST, " ");
    for (v = 0; v < V->NbColumns - 2; ++v) {
      if (value_notzero_p(V->p[l][v])) {
        value_gcd(gcd, V->p[l][v], V->p[l][V->NbColumns - 1]);
        value_divexact(tmp, V->p[l][v], gcd);
        if (value_posz_p(tmp)) {
          if (!first)
            fprintf(DST, "+");
          if (value_notone_p(tmp)) {
            value_print(DST, VALUE_FMT, tmp);
          }
        } else { /* V->p[l][v]/gcd<0 */
          if (value_mone_p(tmp))
            fprintf(DST, "-");
          else {
            value_print(DST, VALUE_FMT, tmp);
          }
        }
        value_divexact(tmp, V->p[l][V->NbColumns - 1], gcd);
        if (value_notone_p(tmp)) {
          fprintf(DST, "%s/", param_names[v]);
          value_print(DST, VALUE_FMT, tmp);
        } else
          fprintf(DST, "%s", param_names[v]);
        first = 0;
      }
    }
 
    /* Constant */
    if (value_notzero_p(V->p[l][v]) || first) {
      if (value_posz_p(V->p[l][v]) && !first)
        fprintf(DST, "+");
      value_gcd(gcd, V->p[l][v], V->p[l][V->NbColumns - 1]);
      value_divexact(tmp, V->p[l][v], gcd);
      value_print(DST, VALUE_FMT, tmp);
      value_divexact(tmp, V->p[l][V->NbColumns - 1], gcd);
      if (value_notone_p(tmp)) {
        fprintf(DST, "/");
        value_print(DST, VALUE_FMT, tmp);
        fprintf(DST, " ");
      }
    }
    if (l < V->NbRows - 1)
      fprintf(DST, ", ");
  }
  fprintf(DST, " ]");
  value_clear(gcd);
  value_clear(tmp);
  return;
} /* Print_Vertex */
 
/*----------------------------------------------------------------------*/
/* VertexCT                                                             */
/* convert a paramvertex from reduced space to normal m-space           */
/*----------------------------------------------------------------------*/
Matrix *VertexCT(Matrix *V, Matrix *CT) {
 
  Matrix *Vt;
  int i, j, k;
 
  if (CT) {
 
    /* Have to transform the vertices to original dimension */
    Vt = Matrix_Alloc(V->NbRows, CT->NbColumns + 1);
    for (i = 0; i < V->NbRows; ++i) {
      value_assign(Vt->p[i][CT->NbColumns], V->p[i][V->NbColumns - 1]);
      for (j = 0; j < CT->NbColumns; j++) {
        for (k = 0; k < CT->NbRows; k++)
          if (value_notzero_p(CT->p[k][j]))
            break;
        if (k < CT->NbRows)
          value_assign(Vt->p[i][j], V->p[i][k]);
        else
          value_set_si(Vt->p[i][j], 0);
      }
    }
    return (Vt);
  } else
    return (NULL);
} /* VertexCT */
 
/*
 * Print the validity Domain 'D' of a parametric polyhedron
 */
void Print_Domain(FILE *DST, Polyhedron *D, char **pname) {
  int l, v;
  int first;
 
  POL_ENSURE_FACETS(D);
  POL_ENSURE_VERTICES(D);
 
  for (l = 0; l < D->NbConstraints; ++l) {
    fprintf(DST, "         ");
    first = 1;
    for (v = 1; v <= D->Dimension; ++v) {
      if (value_notzero_p(D->Constraint[l][v])) {
        if (value_one_p(D->Constraint[l][v])) {
          if (first)
            fprintf(DST, "%s ", pname[v - 1]);
          else
            fprintf(DST, "+ %s ", pname[v - 1]);
        } else if (value_mone_p(D->Constraint[l][v]))
          fprintf(DST, "- %s ", pname[v - 1]);
        else {
          if (value_pos_p(D->Constraint[l][v]) && !first)
            fprintf(DST, "+ ");
          value_print(DST, VALUE_FMT, D->Constraint[l][v]);
          fprintf(DST, "%s ", pname[v - 1]);
        }
        first = 0;
      }
    }
    if (value_notzero_p(D->Constraint[l][v])) {
      if (value_pos_p(D->Constraint[l][v]) && !first)
        fprintf(DST, "+");
      fprintf(DST, " ");
      value_print(DST, VALUE_FMT, D->Constraint[l][v]);
    }
    fprintf(DST, (value_notzero_p(D->Constraint[l][0])) ? " >= 0" : " = 0");
    fprintf(DST, "\n");
  }
  fprintf(DST, "\n");
 
  if (D->next) {
    fprintf(DST, "UNION\n");
    Print_Domain(DST, D->next, pname);
  }
  return;
} /* Print_Domain */
 
/*
 * Given a list of parametrized vertices and an array of parameter names, Print
 * a list of parametrized vertices in a comprehensible format.
 */
void Param_Vertices_Print(FILE *DST, Param_Vertices *PV,
                          char **param_names) {
  Polyhedron *poly;
 
  while (PV) {
    fprintf(DST, "Vertex :\n");
    Print_Vertex(DST, PV->Vertex, param_names);
 
    /* Pour le domaine : */
    fprintf(DST, "   If :\n");
    poly = Constraints2Polyhedron(PV->Domain, 200);
    Print_Domain(DST, poly, param_names);
    Domain_Free(poly);
    PV = PV->next;
  }
  return;
} /* Param_Vertices_Print */
 
/*
 * Given a polyhedron 'Din' in combined data and parametre space, a context
 * polyhedron 'Cin' representing the constraints on the parameter space and
 * a working space size 'working_space', return a parametric polyhedron with
 * a list of distinct validity domains and a complete list of valid vertices
 * associated to each validity domain.
 */
Param_Polyhedron *Polyhedron2Param_Domain(Polyhedron *Din, Polyhedron *Cin,
                                          int working_space) {
 
  Param_Polyhedron *result;
  Param_Domain *D;
 
  POL_ENSURE_FACETS(Din);
  POL_ENSURE_VERTICES(Din);
  POL_ENSURE_FACETS(Cin);
  POL_ENSURE_VERTICES(Cin);
 
  if (emptyQ(Din) || emptyQ(Cin))
    return NULL;
 
#ifdef DEBUGPP
  fprintf(stderr, "Polyhedron2Param_Polyhedron algorithm starting at : %.2fs\n",
          (float)clock() / CLOCKS_PER_SEC);
#endif
 
  /* Find the m-faces, keeping the corresponding domains */
  /* in the linked list PDomains */
  result = Find_m_faces(&Din, Cin, 1, working_space, NULL, NULL);
 
#ifdef DEBUGPP
  if (result)
    fprintf(stderr, "Number of vertices : %d\n", result->nbV);
  fprintf(stderr, "Vertices found at : %.2fs\n",
          (float)clock() / CLOCKS_PER_SEC);
#endif
 
  /* Processing of PVResult and PDomains */
  if (result && Cin->Dimension > 0) /* at least 1 parameter */
    Compute_PDomains(result->D, result->nbV, working_space);
  if (result && CEqualities)
    for (D = result->D; D; D = D->next)
      D->Domain = Add_CEqualities(D->Domain);
  Polyhedron_Free(CEqualities);
 
#ifdef DEBUGPP
  fprintf(stderr, "domains found at : %.2fs\n",
          (float)clock() / CLOCKS_PER_SEC);
#endif
 
  return (result);
} /* Polyhedon2Param_Domain */
 
/*
 *
 */
Param_Polyhedron *Polyhedron2Param_SimplifiedDomain(Polyhedron **Din,
                                                    Polyhedron *Cin,
                                                    int working_space,
                                                    Polyhedron **CEq,
                                                    Matrix **CT) {
 
  Param_Polyhedron *result;
 
  assert(CEq != NULL);
  assert(CT != NULL);
 
  POL_ENSURE_FACETS(*Din);
  POL_ENSURE_VERTICES(*Din);
  POL_ENSURE_FACETS(Cin);
  POL_ENSURE_VERTICES(Cin);
 
#ifdef DEBUGPP
  fprintf(stderr, "Polyhedron2Param_Polyhedron algorithm starting at : %.2fs\n",
          (float)clock() / CLOCKS_PER_SEC);
#endif
 
  /* Find the m-faces, keeping the corresponding domains */
  /* in the linked list PDomains */
  result = Find_m_faces(Din, Cin, 1, working_space, CEq, CT);
 
#ifdef DEBUGPP
  if (result)
    fprintf(stderr, "Number of vertices : %d\n", result->nbV);
  fprintf(stderr, "Vertices found at : %.2fs\n",
          (float)clock() / CLOCKS_PER_SEC);
#endif
 
  /* Processing of PVResult and PDomains */
  if (result && Cin->Dimension > 0) /* at least 1 parameter */
    Compute_PDomains(result->D, result->nbV, working_space);
 
    /* Removed this, Vin100, March 01 */
    /*  if(result && CEqualities )
          for(D=result->D;D;D=D->next)
          D->Domain = Add_CEqualities(D->Domain);
    */
 
#ifdef DEBUGPP
  fprintf(stderr, "domains found at : %.2fs\n",
          (float)clock() / CLOCKS_PER_SEC);
#endif
 
  return (result);
} /* Polyhedron2Param_SimplifiedDomain */
 
/*
 * Free the memory allocated to a list of validity domain of a parametrized
 * polyhedron.
 */
void Param_Domain_Free(Param_Domain *PD) {
 
  Param_Domain *next_pd;
 
  while (PD) {
    free(PD->F);
    Domain_Free(PD->Domain);
    next_pd = PD->next;
    free(PD);
    PD = next_pd;
  }
  return;
} /* Param_Domain_Free */
 
/*
 * Free the memory allocated to a parametric polyhedron 'P'
 */
void Param_Polyhedron_Free(Param_Polyhedron *P) {
 
  if (!P)
    return;
  Param_Vertices_Free(P->V);
  Param_Domain_Free(P->D);
  if (P->Constraints)
    Matrix_Free(P->Constraints);
  if (P->Rays)
    Matrix_Free(P->Rays);
  free(P);
  return;
} /* Param_Polyhedron_Free */
 
/*
 * Scales the parametric polyhedron such that all vertices are integer.
 */
void Param_Polyhedron_Scale_Integer(Param_Polyhedron *PP, Polyhedron **P,
                                    Value *det, unsigned MaxRays) {
  int i;
  int nb_param, nb_vars;
  Vector *denoms;
  Param_Vertices *V;
  Value global_var_lcm;
  Matrix *expansion;
 
  value_set_si(*det, 1);
  if (!PP->nbV)
    return;
 
  nb_param = PP->D->Domain->Dimension;
  nb_vars = PP->V->Vertex->NbRows;
 
  /* Scan the vertices and make an orthogonal expansion of the variable
     space */
  /* a- prepare the array of common denominators */
  denoms = Vector_Alloc(nb_vars);
  value_init(global_var_lcm);
 
  /* b- scan the vertices and compute the variables' global lcms */
  for (V = PP->V; V; V = V->next)
    for (i = 0; i < nb_vars; i++)
      value_lcm(denoms->p[i], denoms->p[i], V->Vertex->p[i][nb_param + 1]);
 
  value_set_si(global_var_lcm, 1);
  for (i = 0; i < nb_vars; i++) {
    value_multiply(*det, *det, denoms->p[i]);
    value_lcm(global_var_lcm, global_var_lcm, denoms->p[i]);
  }
 
  /* scale vertices */
  for (V = PP->V; V; V = V->next)
    for (i = 0; i < nb_vars; i++) {
      Vector_Scale(V->Vertex->p[i], V->Vertex->p[i], denoms->p[i],
                   nb_param + 1);
      Vector_Normalize(V->Vertex->p[i], nb_param + 2);
    }
 
  /* the expansion can be actually writen as global_var_lcm.L^{-1} */
  /* this is equivalent to multiply the rows of P by denoms_det */
  for (i = 0; i < nb_vars; i++)
    value_division(denoms->p[i], global_var_lcm, denoms->p[i]);
 
  /* OPT : we could use a vector instead of a diagonal matrix here (c- and
   * d-).*/
  /* c- make the quick expansion matrix */
  expansion = Matrix_Alloc(nb_vars + nb_param + 1, nb_vars + nb_param + 1);
  for (i = 0; i < nb_vars; i++)
    value_assign(expansion->p[i][i], denoms->p[i]);
  for (i = nb_vars; i < nb_vars + nb_param + 1; i++)
    value_assign(expansion->p[i][i], global_var_lcm);
 
  /* d- apply the variable expansion to the polyhedron */
  if (P)
    *P = Polyhedron_Preimage(*P, expansion, MaxRays);
 
  Matrix_Free(expansion);
  value_clear(global_var_lcm);
  Vector_Free(denoms);
}