source/ehrhart/eval_ehrhart.c

   1 /************************************************/
   2 /*  eval_ehrhart.c                              */
   3 /* functions to evaluate an Ehrhart polynomial. */
   4 /* written by Emmanuel Jeannot (c) 1997.        */
   5 /*  Emmanuel.Jeannot@ens-lyon.fr                */
   6 /*  http://www.ens-lyon.fr/~ejeannot            */
   7 /*                                              */
   8 /* modified 1998, 2000, Vincent Loechner        */
   9 /* (ArithmetiqueLib, Param_Names)               */
  10 /************************************************/
  11
  12 #include <stdio.h>
  13 #include <math.h>
  14 #include <assert.h>
  15 #include <stdlib.h>
  16
  17 #include <polylib/polylib.h>
  18
  19 /* #define EVAL_EHRHART_DEBUG */
  20
  21 /********************************************************/
  22 /* function in domain                                   */
  23 /*    check if the parameters in list_args              */
  24 /*    verifies the constraints of Polyhedron P          */
  25 /********************************************************/
  26 int in_domain(Polyhedron *P, Value *list_args) {
  27
  28   int col,row;
  29   Value v; /* value of the constraint of a row when
  30                parameters are instanciated*/
  31   Value tmp;
  32
  33   value_init(v);
  34   value_init(tmp);
  35
  36   /*P->Constraint constraint matrice of polyhedron P*/
  37   for(row=0;row<P->NbConstraints;row++) {
  38     value_assign(v,P->Constraint[row][P->Dimension+1]); /*constant part*/
  39     for(col=1;col<P->Dimension+1;col++) {
  40       value_multiply(tmp,P->Constraint[row][col],list_args[col-1]);
  41       value_addto(v,v,tmp);
  42     }
  43     if (value_notzero_p(P->Constraint[row][0])) {
  44
  45       /*if v is not >=0 then this constraint is not respected */
  46       if (value_neg_p(v)) {
  47         value_clear(v);
  48         value_clear(tmp);
  49         return 0;
  50       }
  51     }
  52     else {
  53
  54       /*if v is not = 0 then this constraint is not respected */
  55       if (value_notzero_p(v)) {
  56         value_clear(v);
  57         value_clear(tmp);
  58         return 0;
  59       }
  60     }
  61   }
  62
  63   /*if not return before this point => all
  64     the constraints are respected */
  65   value_clear(v);
  66   value_clear(tmp);
  67   return 1;
  68 } /* in_domain */
  69
  70 /****************************************************/
  71 /* function compute enode                           */
  72 /*     compute the value of enode p with parameters */
  73 /*     list "list_args                              */
  74 /*     compute the polynomial or the periodic       */
  75 /****************************************************/
  76
  77 static double compute_enode(enode *p, Value *list_args) {
  78
  79   int i;
  80   Value m, param;
  81   double res=0.0;
  82
  83   if (!p)
  84     return(0.);
  85
  86   value_init(m);
  87   value_init(param);
  88
  89   if (p->type == polynomial) {
  90     if (p->size > 1)
  91                 value_assign(param,list_args[p->pos-1]);
  92
  93     /* Compute the polynomial using Horner's rule */
  94     for (i=p->size-1;i>0;i--) {
  95       res +=compute_evalue(&p->arr[i],list_args);
  96       res *=VALUE_TO_DOUBLE(param);
  97     }
  98     res +=compute_evalue(&p->arr[0],list_args);
  99   }
 100   else if (p->type == periodic) {
 101     value_assign(param,list_args[p->pos-1]);
 102
 103     /* Choose the right element of the periodic */
 104     value_absolute(m,param);
 105     value_set_si(param,p->size);
 106     value_modulus(m,m,param);
 107     res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
 108   }
 109   value_clear(m);
 110   value_clear(param);
 111   return res;
 112 } /* compute_enode */
 113
 114 /*************************************************/
 115 /* return the value of Ehrhart Polynomial        */
 116 /* It returns a double, because since it is      */
 117 /* a recursive function, some intermediate value */
 118 /* might not be integral                         */
 119 /*************************************************/
 120
 121 double compute_evalue(evalue *e,Value *list_args) {
 122
 123   double res;
 124
 125   if (value_notzero_p(e->d)) {
 126     if (value_notone_p(e->d))
 127       res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
 128     else
 129       res = VALUE_TO_DOUBLE(e->x.n);
 130   }
 131   else
 132     res = compute_enode(e->x.p,list_args);
 133   return res;
 134 } /* compute_evalue */
 135
 136
 137 /****************************************************/
 138 /* function compute_poly :                          */
 139 /* Check for the good validity domain               */
 140 /* return the number of point in the Polyhedron     */
 141 /* in allocated memory                              */
 142 /* Using the Ehrhart pseudo-polynomial              */
 143 /****************************************************/
 144 Value *compute_poly(Enumeration *en,Value *list_args) {
 145
 146   Value *tmp;
 147   /*    double d; int i; */
 148
 149   tmp = (Value *) malloc (sizeof(Value));
 150   assert(tmp != NULL);
 151   value_init(*tmp);
 152   value_set_si(*tmp,0);
 153
 154   if(!en)
 155     return(tmp);        /* no ehrhart polynomial */
 156   if(en->ValidityDomain) {
 157     if(!en->ValidityDomain->Dimension) { /* no parameters */
 158       value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
 159       return(tmp);
 160     }
 161   }
 162   else
 163     return(tmp);  /* no Validity Domain */
 164   while(en) {
 165     if(in_domain(en->ValidityDomain,list_args)) {
 166
 167 #ifdef EVAL_EHRHART_DEBUG
 168       Print_Domain(stdout,en->ValidityDomain);
 169       print_evalue(stdout,&en->EP);
 170 #endif
 171
 172       /*                        d = compute_evalue(&en->EP,list_args);
 173                                 i = d;
 174                                 printf("(double)%lf = %d\n", d, i ); */
 175       value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
 176       return(tmp);
 177     }
 178     else
 179       en=en->next;
 180   }
 181   value_set_si(*tmp,0);
 182   return(tmp); /* no compatible domain with the arguments */
 183 } /* compute_poly */
 184
 185
 186
 187