source/ehrhart/homogenization.c

   1 /*
   2     This file is part of PolyLib.
   3
   4     PolyLib is free software: you can redistribute it and/or modify
   5     it under the terms of the GNU General Public License as published by
   6     the Free Software Foundation, either version 3 of the License, or
   7     (at your option) any later version.
   8
   9     PolyLib is distributed in the hope that it will be useful,
  10     but WITHOUT ANY WARRANTY; without even the implied warranty of
  11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12     GNU General Public License for more details.
  13
  14     You should have received a copy of the GNU General Public License
  15     along with PolyLib.  If not, see <http://www.gnu.org/licenses/>.
  16 */
  17
  18 /** homogenization.c
  19     copyright 2004-2005 Bavo Nootaert
  20 **/
  21
  22 #include <assert.h>
  23 #include <stdlib.h>
  24 #include <string.h>
  25 #include <polylib/polylib.h>
  26 #include <polylib/homogenization.h>
  27
  28 static evalue *dehomogenize_periodic(enode *en);
  29 static evalue *dehomogenize_polynomial(enode *en);
  30
  31 Polyhedron *homogenize(Polyhedron *P, unsigned MAXRAYS)
  32 {
  33     Matrix M, *M2;
  34     /* Pretend P is a Matrix for a second */
  35     M.NbRows = P->NbConstraints;
  36     M.NbColumns = P->Dimension+2;
  37     M.p_Init = P->p_Init;
  38     M.p = P->Constraint;
  39     M2 = AddANullColumn(&M);
  40     P = Constraints2Polyhedron(M2, MAXRAYS);
  41     Matrix_Free(M2);
  42     return P;
  43 }
  44
  45 /** dehomogenize an evalue. The last parameter (nb_param) is replaced by 1.
  46     This function is mutually recursive with dehomogenize_enode.
  47 **/
  48 void dehomogenize_evalue(evalue *ep, int nb_param){
  49   evalue *w;
  50
  51   /** cannot dehomogenize rationals **/
  52   if (value_zero_p(ep->d)){
  53
  54     /** we need to replace the last parameter **/
  55     if (ep->x.p->pos == nb_param){
  56       if (ep->x.p->type == periodic && ep->x.p->size > 1){
  57         w = dehomogenize_periodic(ep->x.p);
  58       }
  59       else{
  60         w = dehomogenize_polynomial(ep->x.p);
  61       }
  62       free_evalue_refs(ep);
  63       memcpy(ep, w, sizeof(evalue));
  64       free(w);
  65     }
  66     else{
  67       /** Not the last parameter. Recurse **/
  68       dehomogenize_enode(ep->x.p, nb_param);
  69     }
  70
  71   }
  72 }
  73
  74 /** dehomogenize all evalues in an enode.
  75     This function is mutually recursive with dehomogenize_evalue.
  76 **/
  77 void dehomogenize_enode(enode *p, int nb_param){
  78   evalue *temp;
  79   int i;
  80   for (i = 0; i < p->size; i++){
  81     dehomogenize_evalue(&p->arr[i], nb_param);
  82   }
  83 }
  84
  85
  86 /** return the 1st element of an enode representing a periodic **/
  87 static evalue *dehomogenize_periodic(enode *en){
  88   evalue *w;
  89   assert(en->type == periodic);
  90   assert(en->size > 1);
  91   assert(value_notzero_p(en->arr[1].d));
  92   w = (evalue*)malloc(sizeof(evalue));
  93   value_init(w->d); value_init(w->x.n);
  94   value_assign(w->d, en->arr[1].d); value_assign(w->x.n, en->arr[1].x.n);
  95   return w;
  96 }
  97
  98 /** dehomogenize a polynomial. Assume the enode contains a polynomial in
  99     one variable, the homogenous parameter.
 100     Returns an new evalue, representing a rational.
 101  **/
 102 static evalue *dehomogenize_polynomial(enode *en){
 103   evalue *enn;
 104   evalue *ev;
 105   int i;
 106   double som;
 107   Value num, den, gcd, f1, f2;
 108   assert(en->type == polynomial);
 109   som = 0;
 110   value_init(num); value_init(den); value_init(gcd);
 111   value_init(f1); value_init(f2);
 112   value_set_si(den, 1);
 113
 114   /** enumerate over all coefficients (which are either periodic or rational,
 115       but not polynomial) **/
 116   for (i = 0; i < en->size; i++){
 117     if (value_zero_p(en->arr[i].d)){
 118       if (en->arr[i].x.p->size > 1)
 119         ev = &en->arr[i].x.p->arr[1];
 120       else
 121         ev = &en->arr[i].x.p->arr[0];
 122     }
 123     else{
 124       ev = &en->arr[i];
 125     }
 126     /** add ev (fraction) to num/den **/
 127     value_multiply(f1, den, ev->x.n);
 128     value_multiply(f2, num, ev->d);
 129     value_addto(num, f1, f2);
 130     value_multiply(den, den, ev->d);
 131   }
 132
 133   /** simplify num/den **/
 134   value_gcd(gcd, num, den);
 135   value_divexact(num, num, gcd);
 136   value_divexact(den, den, gcd);
 137
 138   /** create new evalue representing num/den**/
 139   enn = (evalue*)malloc(sizeof(evalue));
 140   value_init(enn->d); value_init(enn->x.n);
 141   value_assign(enn->d, den);
 142   value_assign(enn->x.n, num);
 143
 144   /** cleanup **/
 145   value_clear(gcd);
 146   value_clear(f1); value_clear(f2);
 147   value_clear(num); value_clear(den);
 148
 149   return enn;
 150 }
 151
 152 /** dehomogenize a polyhedron. Assume the polyhedron p is homogenous.
 153     Returns a new polyhedron.
 154 **/
 155 Polyhedron *dehomogenize_polyhedron(Polyhedron *p, int maxRays){
 156   Matrix *constr, *constrh;
 157   Polyhedron *ph;
 158   int i;
 159   constr = Polyhedron2Constraints(p);
 160   constrh = Matrix_Alloc(constr->NbRows, constr->NbColumns - 1);
 161   for (i = 0; i < constr->NbRows; i++){
 162     Vector_Copy(constr->p[i], constrh->p[i], constr->NbColumns - 1);
 163   }
 164   ph = Constraints2Polyhedron(constrh, maxRays);
 165   Matrix_Free(constr); Matrix_Free(constrh);
 166   return ph;
 167 }
 168
 169 /** dehomogenize an enumeration. Replaces each validity domain and
 170     Ehrhart polynomial in the Enumeration en with the dehomogenized form.
 171  **/
 172 void dehomogenize_enumeration(Enumeration* en, int nb_params, int maxRays){
 173   Enumeration *en2;
 174   Polyhedron *vd;
 175   for (en2 = en; en2; en2 = en2->next) {
 176     vd = dehomogenize_polyhedron(en2->ValidityDomain, maxRays);
 177     Polyhedron_Free(en2->ValidityDomain);
 178     en2->ValidityDomain = vd;
 179     dehomogenize_evalue(&en2->EP, nb_params);
 180   }
 181 }