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 Domain 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
  32   if( !P )
  33           return( 0 );
  34
  35   POL_ENSURE_INEQUALITIES(P);
  36
  37   value_init(v);
  38
  39   /* P->Constraint constraint matrice of polyhedron P */
  40   for(row=0;row<P->NbConstraints;row++) {
  41     value_assign(v,P->Constraint[row][P->Dimension+1]); /*constant part*/
  42     for(col=1;col<P->Dimension+1;col++) {
  43       value_addmul(v, P->Constraint[row][col], list_args[col-1]);
  44     }
  45     if (value_notzero_p(P->Constraint[row][0])) {
  46
  47       /*if v is not >=0 then this constraint is not respected */
  48       if (value_neg_p(v)) {
  49         value_clear(v);
  50         return( in_domain(P->next, list_args) );
  51       }
  52     }
  53     else {
  54
  55       /*if v is not = 0 then this constraint is not respected */
  56       if (value_notzero_p(v)) {
  57         value_clear(v);
  58         return( in_domain(P->next, list_args) );
  59       }
  60     }
  61   }
  62
  63   /* if not return before this point => all the constraints are respected */
  64   value_clear(v);
  65   return 1;
  66 } /* in_domain */
  67
  68 /****************************************************/
  69 /* function compute enode                           */
  70 /*     compute the value of enode p with parameters */
  71 /*     list "list_args                              */
  72 /*     compute the polynomial or the periodic       */
  73 /****************************************************/
  74
  75 static double compute_enode(enode *p, Value *list_args) {
  76
  77   int i;
  78   Value m, param;
  79   double res=0.0;
  80
  81   if (!p)
  82     return(0.);
  83
  84   value_init(m);
  85   value_init(param);
  86
  87   if (p->type == polynomial) {
  88     if (p->size > 1)
  89                  value_assign(param,list_args[p->pos-1]);
  90
  91     /* Compute the polynomial using Horner's rule */
  92     for (i=p->size-1;i>0;i--) {
  93       res +=compute_evalue(&p->arr[i],list_args);
  94       res *=VALUE_TO_DOUBLE(param);
  95     }
  96     res +=compute_evalue(&p->arr[0],list_args);
  97   }
  98   else if (p->type == periodic) {
  99     value_assign(m,list_args[p->pos-1]);
 100
 101     /* Choose the right element of the periodic */
 102     value_set_si(param,p->size);
 103     value_pmodulus(m,m,param);
 104     res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
 105   }
 106   value_clear(m);
 107   value_clear(param);
 108   return res;
 109 } /* compute_enode */
 110
 111 /*************************************************/
 112 /* return the value of Ehrhart Polynomial        */
 113 /* It returns a double, because since it is      */
 114 /* a recursive function, some intermediate value */
 115 /* might not be integral                         */
 116 /*************************************************/
 117
 118 double compute_evalue(evalue *e,Value *list_args) {
 119
 120   double res;
 121
 122   if (value_notzero_p(e->d)) {
 123     if (value_notone_p(e->d))
 124       res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
 125     else
 126       res = VALUE_TO_DOUBLE(e->x.n);
 127   }
 128   else
 129     res = compute_enode(e->x.p,list_args);
 130   return res;
 131 } /* compute_evalue */
 132
 133
 134 /****************************************************/
 135 /* function compute_poly :                          */
 136 /* Check for the good validity domain               */
 137 /* return the number of point in the Polyhedron     */
 138 /* in allocated memory                              */
 139 /* Using the Ehrhart pseudo-polynomial              */
 140 /****************************************************/
 141 Value *compute_poly(Enumeration *en,Value *list_args) {
 142
 143   Value *tmp;
 144   /*        double d; int i; */
 145
 146   tmp = (Value *) malloc (sizeof(Value));
 147   assert(tmp != NULL);
 148   value_init(*tmp);
 149   value_set_si(*tmp,0);
 150
 151   if(!en)
 152     return(tmp);        /* no ehrhart polynomial */
 153   if(en->ValidityDomain) {
 154     if(!en->ValidityDomain->Dimension) { /* no parameters */
 155       value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
 156       return(tmp);
 157     }
 158   }
 159   else
 160     return(tmp);  /* no Validity Domain */
 161   while(en) {
 162     if(in_domain(en->ValidityDomain,list_args)) {
 163
 164 #ifdef EVAL_EHRHART_DEBUG
 165       Print_Domain(stdout,en->ValidityDomain,NULL);
 166       print_evalue(stdout,&en->EP,NULL);
 167 #endif
 168
 169       /*                        d = compute_evalue(&en->EP,list_args);
 170                                 i = d;
 171                                 printf("(double)%lf = %d\n", d, i ); */
 172       value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
 173       return(tmp);
 174     }
 175     else
 176       en=en->next;
 177   }
 178   value_set_si(*tmp,0);
 179   return(tmp); /* no compatible domain with the arguments */
 180 } /* compute_poly */
 181
 182
 183
 184