?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
compute modulo if argument is known
[barvinok.git] / ev_operations.c
blobc86fd1757fed95c340d85235e67c937737eb367e
1 #include <assert.h>
2 #include <math.h>
3
4 #include "ev_operations.h"
5
6 void evalue_set_si(evalue *ev, int n, int d) {
7 value_set_si(ev->d, d);
8 value_init(ev->x.n);
9 value_set_si(ev->x.n, n);
10 }
11
12 void evalue_set(evalue *ev, Value n, Value d) {
13 value_assign(ev->d, d);
14 value_init(ev->x.n);
15 value_assign(ev->x.n, n);
16 }
17
18 void aep_evalue(evalue *e, int *ref) {
19
20 enode *p;
21 int i;
22
23 if (value_notzero_p(e->d))
24 return; /* a rational number, its already reduced */
25 if(!(p = e->x.p))
26 return; /* hum... an overflow probably occured */
27
28 /* First check the components of p */
29 for (i=0;i<p->size;i++)
30 aep_evalue(&p->arr[i],ref);
31
32 /* Then p itself */
33 if (p->type != modulo)
34 p->pos = ref[p->pos-1]+1;
35 return;
36 } /* aep_evalue */
37
38 /** Comments */
39 void addeliminatedparams_evalue(evalue *e,Matrix *CT) {
40
41 enode *p;
42 int i, j;
43 int *ref;
44
45 if (value_notzero_p(e->d))
46 return; /* a rational number, its already reduced */
47 if(!(p = e->x.p))
48 return; /* hum... an overflow probably occured */
49
50 /* Compute ref */
51 ref = (int *)malloc(sizeof(int)*(CT->NbRows-1));
52 for(i=0;i<CT->NbRows-1;i++)
53 for(j=0;j<CT->NbColumns;j++)
54 if(value_notzero_p(CT->p[i][j])) {
55 ref[i] = j;
56 break;
57 }
58
59 /* Transform the references in e, using ref */
60 aep_evalue(e,ref);
61 free( ref );
62 return;
63 } /* addeliminatedparams_evalue */
64
65 struct fixed_param {
66 int pos;
67 Value v;
68 };
69
70 void _reduce_evalue (evalue *e, int n, struct fixed_param *fixed) {
71
72 enode *p;
73 int i, j, k;
74
75 if (value_notzero_p(e->d))
76 return; /* a rational number, its already reduced */
77 if(!(p = e->x.p))
78 return; /* hum... an overflow probably occured */
79
80 /* First reduce the components of p */
81 for (i=0; i<p->size; i++)
82 _reduce_evalue(&p->arr[i], n, fixed);
83
84 if (p->type==periodic) {
85
86 /* Try to reduce the period */
87 for (i=1; i<=(p->size)/2; i++) {
88 if ((p->size % i)==0) {
89
90 /* Can we reduce the size to i ? */
91 for (j=0; j<i; j++)
92 for (k=j+i; k<e->x.p->size; k+=i)
93 if (!eequal(&p->arr[j], &p->arr[k])) goto you_lose;
94
95 /* OK, lets do it */
96 for (j=i; j<p->size; j++) free_evalue_refs(&p->arr[j]);
97 p->size = i;
98 break;
99
100 you_lose: /* OK, lets not do it */
101 continue;
102 }
103 }
104
105 /* Try to reduce its strength */
106 if (p->size == 1) {
107 value_clear(e->d);
108 memcpy(e,&p->arr[0],sizeof(evalue));
109 free(p);
110 }
111 }
112 else if (p->type==polynomial) {
113 for (k = 0; k < n; ++k) {
114 if (fixed[k].pos == p->pos) {
115 evalue v;
116 value_init(v.d);
117 value_set_si(v.d, 1);
118 value_init(v.x.n);
119 value_assign(v.x.n, fixed[k].v);
120 for (i=p->size-1;i>=1;i--) {
121 emul(&v, &p->arr[i]);
122 eadd(&p->arr[i], &p->arr[i-1]);
123 free_evalue_refs(&(p->arr[i]));
124 }
125 p->size = 1;
126 free_evalue_refs(&v);
127 }
128 }
129
130 /* Try to reduce the degree */
131 for (i=p->size-1;i>=1;i--) {
132 if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
133 break;
134 /* Zero coefficient */
135 free_evalue_refs(&(p->arr[i]));
136 }
137 if (i+1<p->size)
138 p->size = i+1;
139
140 /* Try to reduce its strength */
141 if (p->size == 1) {
142 value_clear(e->d);
143 memcpy(e,&p->arr[0],sizeof(evalue));
144 free(p);
145 }
146 }
147 else if (p->type==modulo) {
148 if (value_notzero_p(p->arr[0].d)) {
149 evalue v;
150 value_init(v.d);
151 value_set_si(v.d, 1);
152 value_init(v.x.n);
153 value_set_si(p->arr[0].d, p->pos);
154 mpz_fdiv_r(v.x.n, p->arr[0].x.n, p->arr[0].d);
155
156 for (i=p->size-1;i>=2;i--) {
157 emul(&v, &p->arr[i]);
158 eadd(&p->arr[i], &p->arr[i-1]);
159 free_evalue_refs(&(p->arr[i]));
160 }
161 p->size = 2;
162 free_evalue_refs(&v);
163 }
164
165 /* Try to reduce the degree */
166 for (i=p->size-1;i>=2;i--) {
167 if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
168 break;
169 /* Zero coefficient */
170 free_evalue_refs(&(p->arr[i]));
171 }
172 if (i+1<p->size)
173 p->size = i+1;
174
175 /* Try to reduce its strength */
176 if (p->size == 2) {
177 value_clear(e->d);
178 memcpy(e,&p->arr[1],sizeof(evalue));
179 free_evalue_refs(&(p->arr[0]));
180 free(p);
181 }
182 }
183 else if (p->type == relation) {
184 if (p->size == 3 && EVALUE_IS_ZERO(p->arr[2])) {
185 free_evalue_refs(&(p->arr[2]));
186 p->size = 2;
187 }
188 if (p->size == 2 && EVALUE_IS_ZERO(p->arr[1])) {
189 free_evalue_refs(&(p->arr[1]));
190 free_evalue_refs(&(p->arr[0]));
191 evalue_set_si(e, 0, 1);
192 }
193 }
194 } /* reduce_evalue */
195
196 void reduce_evalue (evalue *e) {
197 if (value_notzero_p(e->d))
198 return; /* a rational number, its already reduced */
199
200 if (e->x.p->type == partition) {
201 int i;
202 int n;
203 struct fixed_param *fixed = 0;
204 unsigned dim = -1;
205 for (i = 0; i < e->x.p->size/2; ++i) {
206 n = 0;
207 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
208 dim = D->Dimension;
209 if (!D->next && D->NbEq) {
210 int j, k;
211 if (!fixed) {
212 fixed = (struct fixed_param*)
213 malloc(D->Dimension * sizeof(*fixed));
214 for (j = 0; j < dim; ++j)
215 value_init(fixed[j].v);
216 }
217 for (j = 0; j < D->NbEq; ++j) {
218 for (k = 0; k < D->Dimension; ++k)
219 if (value_notzero_p(D->Constraint[j][k+1])) {
220 fixed[n].pos = k+1;
221 if (value_one_p(D->Constraint[j][k+1]))
222 value_oppose(fixed[n].v, D->Constraint[j][dim+1]);
223 else if (value_mone_p(D->Constraint[j][k+1]))
224 value_assign(fixed[n].v, D->Constraint[j][dim+1]);
225 else
226 assert(0);
227 ++n;
228 }
229 }
230 }
231 _reduce_evalue(&e->x.p->arr[2*i+1], n, fixed);
232 }
233 if (fixed) {
234 int j;
235 for (j = 0; j < dim; ++j)
236 value_clear(fixed[j].v);
237 free(fixed);
238 }
239 } else
240 _reduce_evalue(e, 0, 0);
241 }
242
243 void print_evalue(FILE *DST,evalue *e,char **pname) {
244
245 if(value_notzero_p(e->d)) {
246 if(value_notone_p(e->d)) {
247 value_print(DST,VALUE_FMT,e->x.n);
248 fprintf(DST,"/");
249 value_print(DST,VALUE_FMT,e->d);
250 }
251 else {
252 value_print(DST,VALUE_FMT,e->x.n);
253 }
254 }
255 else
256 print_enode(DST,e->x.p,pname);
257 return;
258 } /* print_evalue */
259
260 void print_enode(FILE *DST,enode *p,char **pname) {
261
262 int i;
263
264 if (!p) {
265 fprintf(DST, "NULL");
266 return;
267 }
268 switch (p->type) {
269 case evector:
270 fprintf(DST, "{ ");
271 for (i=0; i<p->size; i++) {
272 print_evalue(DST, &p->arr[i], pname);
273 if (i!=(p->size-1))
274 fprintf(DST, ", ");
275 }
276 fprintf(DST, " }\n");
277 break;
278 case polynomial:
279 fprintf(DST, "( ");
280 for (i=p->size-1; i>=0; i--) {
281 print_evalue(DST, &p->arr[i], pname);
282 if (i==1) fprintf(DST, " * %s + ", pname[p->pos-1]);
283 else if (i>1)
284 fprintf(DST, " * %s^%d + ", pname[p->pos-1], i);
285 }
286 fprintf(DST, " )\n");
287 break;
288 case periodic:
289 fprintf(DST, "[ ");
290 for (i=0; i<p->size; i++) {
291 print_evalue(DST, &p->arr[i], pname);
292 if (i!=(p->size-1)) fprintf(DST, ", ");
293 }
294 fprintf(DST," ]_%s", pname[p->pos-1]);
295 break;
296 case modulo:
297 fprintf(DST, "( ");
298 for (i=p->size-1; i>=1; i--) {
299 print_evalue(DST, &p->arr[i], pname);
300 if (i >= 2) {
301 fprintf(DST, " * ");
302 if (i > 2)
303 fprintf(DST, "(");
304 fprintf(DST, "(");
305 print_evalue(DST, &p->arr[0], pname);
306 fprintf(DST, ") mod %d", p->pos);
307 if (i>2)
308 fprintf(DST, ")^%d + ", i-1);
309 else
310 fprintf(DST, " + ", i-1);
311 }
312 }
313 fprintf(DST, " )\n");
314 break;
315 case relation:
316 fprintf(DST, "[ ");
317 print_evalue(DST, &p->arr[0], pname);
318 fprintf(DST, "= 0 ] * \n");
319 print_evalue(DST, &p->arr[1], pname);
320 if (p->size > 2) {
321 fprintf(DST, " +\n [ ");
322 print_evalue(DST, &p->arr[0], pname);
323 fprintf(DST, "!= 0 ] * \n");
324 print_evalue(DST, &p->arr[2], pname);
325 }
326 break;
327 case partition:
328 for (i=0; i<p->size/2; i++) {
329 Print_Domain(DST, EVALUE_DOMAIN(p->arr[2*i]), pname);
330 print_evalue(DST, &p->arr[2*i+1], pname);
331 }
332 break;
333 default:
334 assert(0);
335 }
336 return;
337 } /* print_enode */
338
339 static int mod_term_smaller(evalue *e1, evalue *e2)
340 {
341 if (value_notzero_p(e1->d)) {
342 if (value_zero_p(e2->d))
343 return 1;
344 return value_lt(e1->x.n, e2->x.n);
345 }
346 if (value_notzero_p(e2->d))
347 return 0;
348 if (e1->x.p->pos < e2->x.p->pos)
349 return 1;
350 else if (e1->x.p->pos > e2->x.p->pos)
351 return 0;
352 else
353 return mod_term_smaller(&e1->x.p->arr[0], &e2->x.p->arr[0]);
354 }
355
356 static void eadd_rev(evalue *e1, evalue *res)
357 {
358 evalue ev;
359 value_init(ev.d);
360 evalue_copy(&ev, e1);
361 eadd(res, &ev);
362 free_evalue_refs(res);
363 *res = ev;
364 }
365
366 static void eadd_rev_cst (evalue *e1, evalue *res)
367 {
368 evalue ev;
369 value_init(ev.d);
370 evalue_copy(&ev, e1);
371 eadd(res, &ev.x.p->arr[ev.x.p->type==modulo]);
372 free_evalue_refs(res);
373 *res = ev;
374 }
375
376 struct section { Polyhedron * D; evalue E; };
377
378 void eadd_partitions (evalue *e1,evalue *res)
379 {
380 int n, i, j;
381 Polyhedron *d, *fd;
382 struct section *s;
383 s = (struct section *)
384 malloc((e1->x.p->size/2+1) * (res->x.p->size/2+1) *
385 sizeof(struct section));
386 assert(s);
387
388 n = 0;
389 for (j = 0; j < e1->x.p->size/2; ++j) {
390 assert(res->x.p->size >= 2);
391 fd = DomainDifference(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
392 EVALUE_DOMAIN(res->x.p->arr[0]), 0);
393 if (!emptyQ(fd))
394 for (i = 1; i < res->x.p->size/2; ++i) {
395 Polyhedron *t = fd;
396 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
397 Domain_Free(t);
398 if (emptyQ(fd))
399 break;
400 }
401 if (emptyQ(fd)) {
402 Domain_Free(fd);
403 continue;
404 }
405 value_init(s[n].E.d);
406 evalue_copy(&s[n].E, &e1->x.p->arr[2*j+1]);
407 s[n].D = fd;
408 ++n;
409 }
410 for (i = 0; i < res->x.p->size/2; ++i) {
411 fd = EVALUE_DOMAIN(res->x.p->arr[2*i]);
412 for (j = 0; j < e1->x.p->size/2; ++j) {
413 Polyhedron *t;
414 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
415 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
416 if (emptyQ(d)) {
417 Domain_Free(d);
418 continue;
419 }
420 t = fd;
421 fd = DomainDifference(fd, EVALUE_DOMAIN(e1->x.p->arr[2*j]), 0);
422 if (t != EVALUE_DOMAIN(res->x.p->arr[2*i]))
423 Domain_Free(t);
424 value_init(s[n].E.d);
425 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
426 eadd(&e1->x.p->arr[2*j+1], &s[n].E);
427 s[n].D = d;
428 ++n;
429 }
430 if (!emptyQ(fd)) {
431 s[n].E = res->x.p->arr[2*i+1];
432 s[n].D = fd;
433 ++n;
434 } else
435 free_evalue_refs(&res->x.p->arr[2*i+1]);
436 if (fd != EVALUE_DOMAIN(res->x.p->arr[2*i]))
437 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
438 value_clear(res->x.p->arr[2*i].d);
439 }
440
441 free(res->x.p);
442 res->x.p = new_enode(partition, 2*n, -1);
443 for (j = 0; j < n; ++j) {
444 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
445 value_clear(res->x.p->arr[2*j+1].d);
446 res->x.p->arr[2*j+1] = s[j].E;
447 }
448
449 free(s);
450 }
451
452 static void explicit_complement(evalue *res)
453 {
454 enode *rel = new_enode(relation, 3, 0);
455 assert(rel);
456 value_clear(rel->arr[0].d);
457 rel->arr[0] = res->x.p->arr[0];
458 value_clear(rel->arr[1].d);
459 rel->arr[1] = res->x.p->arr[1];
460 value_set_si(rel->arr[2].d, 1);
461 value_init(rel->arr[2].x.n);
462 value_set_si(rel->arr[2].x.n, 0);
463 free(res->x.p);
464 res->x.p = rel;
465 }
466
467 void eadd(evalue *e1,evalue *res) {
468
469 int i;
470 if (value_notzero_p(e1->d) && value_notzero_p(res->d)) {
471 /* Add two rational numbers */
472 Value g,m1,m2;
473 value_init(g);
474 value_init(m1);
475 value_init(m2);
476
477 value_multiply(m1,e1->x.n,res->d);
478 value_multiply(m2,res->x.n,e1->d);
479 value_addto(res->x.n,m1,m2);
480 value_multiply(res->d,e1->d,res->d);
481 Gcd(res->x.n,res->d,&g);
482 if (value_notone_p(g)) {
483 value_division(res->d,res->d,g);
484 value_division(res->x.n,res->x.n,g);
485 }
486 value_clear(g); value_clear(m1); value_clear(m2);
487 return ;
488 }
489 else if (value_notzero_p(e1->d) && value_zero_p(res->d)) {
490 switch (res->x.p->type) {
491 case polynomial:
492 /* Add the constant to the constant term of a polynomial*/
493 eadd(e1, &res->x.p->arr[0]);
494 return ;
495 case periodic:
496 /* Add the constant to all elements of a periodic number */
497 for (i=0; i<res->x.p->size; i++) {
498 eadd(e1, &res->x.p->arr[i]);
499 }
500 return ;
501 case evector:
502 fprintf(stderr, "eadd: cannot add const with vector\n");
503 return;
504 case modulo:
505 eadd(e1, &res->x.p->arr[1]);
506 return ;
507 case partition:
508 assert(EVALUE_IS_ZERO(*e1));
509 break; /* Do nothing */
510 case relation:
511 /* Create (zero) complement if needed */
512 if (res->x.p->size < 3 && !EVALUE_IS_ZERO(*e1))
513 explicit_complement(res);
514 for (i = 1; i < res->x.p->size; ++i)
515 eadd(e1, &res->x.p->arr[i]);
516 break;
517 default:
518 assert(0);
519 }
520 }
521 /* add polynomial or periodic to constant
522 * you have to exchange e1 and res, before doing addition */
523
524 else if (value_zero_p(e1->d) && value_notzero_p(res->d)) {
525 eadd_rev(e1, res);
526 return;
527 }
528 else { // ((e1->d==0) && (res->d==0))
529 assert(!((e1->x.p->type == partition) ^
530 (res->x.p->type == partition)));
531 if (e1->x.p->type == partition) {
532 eadd_partitions(e1, res);
533 return;
534 }
535 if (e1->x.p->type == relation &&
536 (res->x.p->type != relation ||
537 mod_term_smaller(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
538 eadd_rev(e1, res);
539 return;
540 }
541 if (res->x.p->type == relation) {
542 if (e1->x.p->type == relation &&
543 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
544 if (res->x.p->size < 3 && e1->x.p->size == 3)
545 explicit_complement(res);
546 if (e1->x.p->size < 3 && res->x.p->size == 3)
547 explicit_complement(e1);
548 for (i = 1; i < res->x.p->size; ++i)
549 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
550 return;
551 }
552 if (res->x.p->size < 3)
553 explicit_complement(res);
554 for (i = 1; i < res->x.p->size; ++i)
555 eadd(e1, &res->x.p->arr[i]);
556 return;
557 }
558 if ((e1->x.p->type != res->x.p->type) ) {
559 /* adding to evalues of different type. two cases are possible
560 * res is periodic and e1 is polynomial, you have to exchange
561 * e1 and res then to add e1 to the constant term of res */
562 if (e1->x.p->type == polynomial) {
563 eadd_rev_cst(e1, res);
564 }
565 else if (res->x.p->type == polynomial) {
566 /* res is polynomial and e1 is periodic,
567 add e1 to the constant term of res */
568
569 eadd(e1,&res->x.p->arr[0]);
570 } else
571 assert(0);
572
573 return;
574 }
575 else if (e1->x.p->pos != res->x.p->pos ||
576 (res->x.p->type == modulo &&
577 !eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
578 /* adding evalues of different position (i.e function of different unknowns
579 * to case are possible */
580
581 switch (res->x.p->type) {
582 case modulo:
583 if(mod_term_smaller(res, e1))
584 eadd(e1,&res->x.p->arr[1]);
585 else
586 eadd_rev_cst(e1, res);
587 return;
588 case polynomial: // res and e1 are polynomials
589 // add e1 to the constant term of res
590
591 if(res->x.p->pos < e1->x.p->pos)
592 eadd(e1,&res->x.p->arr[0]);
593 else
594 eadd_rev_cst(e1, res);
595 // value_clear(g); value_clear(m1); value_clear(m2);
596 return;
597 case periodic: // res and e1 are pointers to periodic numbers
598 //add e1 to all elements of res
599
600 if(res->x.p->pos < e1->x.p->pos)
601 for (i=0;i<res->x.p->size;i++) {
602 eadd(e1,&res->x.p->arr[i]);
603 }
604 else
605 eadd_rev(e1, res);
606 return;
607 }
608 }
609
610
611 //same type , same pos and same size
612 if (e1->x.p->size == res->x.p->size) {
613 // add any element in e1 to the corresponding element in res
614 if (res->x.p->type == modulo)
615 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
616 i = res->x.p->type == modulo ? 1 : 0;
617 for (; i<res->x.p->size; i++) {
618 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
619 }
620 return ;
621 }
622
623 /* Sizes are different */
624 switch(res->x.p->type) {
625 case polynomial:
626 case modulo:
627 /* VIN100: if e1-size > res-size you have to copy e1 in a */
628 /* new enode and add res to that new node. If you do not do */
629 /* that, you lose the the upper weight part of e1 ! */
630
631 if(e1->x.p->size > res->x.p->size)
632 eadd_rev(e1, res);
633 else {
634
635 if (res->x.p->type == modulo)
636 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
637 i = res->x.p->type == modulo ? 1 : 0;
638 for (; i<e1->x.p->size ; i++) {
639 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
640 }
641
642 return ;
643 }
644 break;
645
646 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
647 case periodic:
648 {
649 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
650 of the sizes of e1 and res, then to copy res periodicaly in ne, after
651 to add periodicaly elements of e1 to elements of ne, and finaly to
652 return ne. */
653 int x,y,p;
654 Value ex, ey ,ep;
655 evalue *ne;
656 value_init(ex); value_init(ey);value_init(ep);
657 x=e1->x.p->size;
658 y= res->x.p->size;
659 value_set_si(ex,e1->x.p->size);
660 value_set_si(ey,res->x.p->size);
661 value_assign (ep,*Lcm(ex,ey));
662 p=(int)mpz_get_si(ep);
663 ne= (evalue *) malloc (sizeof(evalue));
664 value_init(ne->d);
665 value_set_si( ne->d,0);
666
667 ne->x.p=new_enode(res->x.p->type,p, res->x.p->pos);
668 for(i=0;i<p;i++) {
669 value_assign(ne->x.p->arr[i].d, res->x.p->arr[i%y].d);
670 if (value_notzero_p(ne->x.p->arr[i].d)) {
671 value_init(ne->x.p->arr[i].x.n);
672 value_assign(ne->x.p->arr[i].x.n, res->x.p->arr[i%y].x.n);
673 }
674 else {
675 ne->x.p->arr[i].x.p =ecopy(res->x.p->arr[i%y].x.p);
676 }
677 }
678 for(i=0;i<p;i++) {
679 eadd(&e1->x.p->arr[i%x], &ne->x.p->arr[i]);
680 }
681
682 value_assign(res->d, ne->d);
683 res->x.p=ne->x.p;
684
685 return ;
686 }
687 case evector:
688 fprintf(stderr, "eadd: ?cannot add vectors of different length\n");
689 return ;
690 }
691 }
692 return ;
693 }/* eadd */
694
695 static void emul_rev (evalue *e1, evalue *res)
696 {
697 evalue ev;
698 value_init(ev.d);
699 evalue_copy(&ev, e1);
700 emul(res, &ev);
701 free_evalue_refs(res);
702 *res = ev;
703 }
704
705 static void emul_poly (evalue *e1, evalue *res)
706 {
707 int i, j, o = res->x.p->type == modulo;
708 evalue tmp;
709 int size=(e1->x.p->size + res->x.p->size - o - 1);
710 value_init(tmp.d);
711 value_set_si(tmp.d,0);
712 tmp.x.p=new_enode(res->x.p->type, size, res->x.p->pos);
713 if (o)
714 evalue_copy(&tmp.x.p->arr[0], &e1->x.p->arr[0]);
715 for (i=o; i < e1->x.p->size; i++) {
716 evalue_copy(&tmp.x.p->arr[i], &e1->x.p->arr[i]);
717 emul(&res->x.p->arr[o], &tmp.x.p->arr[i]);
718 }
719 for (; i<size; i++)
720 evalue_set_si(&tmp.x.p->arr[i], 0, 1);
721 for (i=o+1; i<res->x.p->size; i++)
722 for (j=o; j<e1->x.p->size; j++) {
723 evalue ev;
724 value_init(ev.d);
725 evalue_copy(&ev, &e1->x.p->arr[j]);
726 emul(&res->x.p->arr[i], &ev);
727 eadd(&ev, &tmp.x.p->arr[i+j-o]);
728 free_evalue_refs(&ev);
729 }
730 free_evalue_refs(res);
731 *res = tmp;
732 }
733
734 void emul_partitions (evalue *e1,evalue *res)
735 {
736 int n, i, j, k;
737 Polyhedron *d;
738 struct section *s;
739 s = (struct section *)
740 malloc((e1->x.p->size/2) * (res->x.p->size/2) *
741 sizeof(struct section));
742 assert(s);
743
744 n = 0;
745 for (i = 0; i < res->x.p->size/2; ++i) {
746 for (j = 0; j < e1->x.p->size/2; ++j) {
747 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
748 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
749 if (emptyQ(d)) {
750 Domain_Free(d);
751 continue;
752 }
753
754 /* This code is only needed because the partitions
755 are not true partitions.
756 */
757 for (k = 0; k < n; ++k) {
758 if (DomainIncludes(s[k].D, d))
759 break;
760 if (DomainIncludes(d, s[k].D)) {
761 Domain_Free(s[k].D);
762 free_evalue_refs(&s[k].E);
763 if (n > k)
764 s[k] = s[--n];
765 --k;
766 }
767 }
768 if (k < n) {
769 Domain_Free(d);
770 continue;
771 }
772
773 value_init(s[n].E.d);
774 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
775 emul(&e1->x.p->arr[2*j+1], &s[n].E);
776 s[n].D = d;
777 ++n;
778 }
779 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
780 value_clear(res->x.p->arr[2*i].d);
781 free_evalue_refs(&res->x.p->arr[2*i+1]);
782 }
783
784 free(res->x.p);
785 res->x.p = new_enode(partition, 2*n, -1);
786 for (j = 0; j < n; ++j) {
787 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
788 value_clear(res->x.p->arr[2*j+1].d);
789 res->x.p->arr[2*j+1] = s[j].E;
790 }
791
792 free(s);
793 }
794
795 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
796 * do a copy of "res" befor calling this function if you nead it after. The vector type of
797 * evalues is not treated here */
798
799 void emul (evalue *e1, evalue *res ){
800 int i,j;
801
802 if((value_zero_p(e1->d)&&e1->x.p->type==evector)||(value_zero_p(res->d)&&(res->x.p->type==evector))) {
803 fprintf(stderr, "emul: do not proced on evector type !\n");
804 return;
805 }
806
807 if (value_zero_p(e1->d) && e1->x.p->type == partition) {
808 if (value_zero_p(res->d) && res->x.p->type == partition)
809 emul_partitions(e1, res);
810 else
811 emul_rev(e1, res);
812 } else if (value_zero_p(res->d) && res->x.p->type == partition) {
813 for (i = 0; i < res->x.p->size/2; ++i)
814 emul(e1, &res->x.p->arr[2*i+1]);
815 } else
816 if (value_zero_p(res->d) && res->x.p->type == relation) {
817 for (i = 1; i < res->x.p->size; ++i)
818 emul(e1, &res->x.p->arr[i]);
819 } else
820 if(value_zero_p(e1->d)&& value_zero_p(res->d)) {
821 switch(e1->x.p->type) {
822 case polynomial:
823 switch(res->x.p->type) {
824 case polynomial:
825 if(e1->x.p->pos == res->x.p->pos) {
826 /* Product of two polynomials of the same variable */
827 emul_poly(e1, res);
828 return;
829 }
830 else {
831 /* Product of two polynomials of different variables */
832
833 if(res->x.p->pos < e1->x.p->pos)
834 for( i=0; i<res->x.p->size ; i++)
835 emul(e1, &res->x.p->arr[i]);
836 else
837 emul_rev(e1, res);
838
839 return ;
840 }
841 case periodic:
842 case modulo:
843 /* Product of a polynomial and a periodic or modulo */
844 emul_rev(e1, res);
845 return;
846 }
847 case periodic:
848 switch(res->x.p->type) {
849 case periodic:
850 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size==res->x.p->size) {
851 /* Product of two periodics of the same parameter and period */
852
853 for(i=0; i<res->x.p->size;i++)
854 emul(&(e1->x.p->arr[i]), &(res->x.p->arr[i]));
855
856 return;
857 }
858 else{
859 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size!=res->x.p->size) {
860 /* Product of two periodics of the same parameter and different periods */
861 evalue *newp;
862 Value x,y,z;
863 int ix,iy,lcm;
864 value_init(x); value_init(y);value_init(z);
865 ix=e1->x.p->size;
866 iy=res->x.p->size;
867 value_set_si(x,e1->x.p->size);
868 value_set_si(y,res->x.p->size);
869 value_assign (z,*Lcm(x,y));
870 lcm=(int)mpz_get_si(z);
871 newp= (evalue *) malloc (sizeof(evalue));
872 value_init(newp->d);
873 value_set_si( newp->d,0);
874 newp->x.p=new_enode(periodic,lcm, e1->x.p->pos);
875 for(i=0;i<lcm;i++) {
876 value_assign(newp->x.p->arr[i].d, res->x.p->arr[i%iy].d);
877 if (value_notzero_p(newp->x.p->arr[i].d)) {
878 value_assign(newp->x.p->arr[i].x.n, res->x.p->arr[i%iy].x.n);
879 }
880 else {
881 newp->x.p->arr[i].x.p =ecopy(res->x.p->arr[i%iy].x.p);
882 }
883
884 }
885 for(i=0;i<lcm;i++)
886 emul(&e1->x.p->arr[i%ix], &newp->x.p->arr[i]);
887
888 value_assign(res->d,newp->d);
889 res->x.p=newp->x.p;
890
891 value_clear(x); value_clear(y);value_clear(z);
892 return ;
893 }
894 else {
895 /* Product of two periodics of different parameters */
896
897 for(i=0; i<res->x.p->size; i++)
898 emul(e1, &(res->x.p->arr[i]));
899
900 return;
901 }
902 }
903 case polynomial:
904 /* Product of a periodic and a polynomial */
905
906 for(i=0; i<res->x.p->size ; i++)
907 emul(e1, &(res->x.p->arr[i]));
908
909 return;
910
911 }
912 case modulo:
913 switch(res->x.p->type) {
914 case polynomial:
915 for(i=0; i<res->x.p->size ; i++)
916 emul(e1, &(res->x.p->arr[i]));
917 return;
918 case periodic:
919 assert(0);
920 case modulo:
921 if (e1->x.p->pos == res->x.p->pos &&
922 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
923 if (e1->x.p->pos != 2)
924 emul_poly(e1, res);
925 else {
926 evalue tmp;
927 /* x mod 2 == (x mod 2)^2 */
928 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1) (x mod 2) */
929 assert(e1->x.p->size == 3);
930 assert(res->x.p->size == 3);
931 value_init(tmp.d);
932 evalue_copy(&tmp, &res->x.p->arr[1]);
933 eadd(&res->x.p->arr[2], &tmp);
934 emul(&e1->x.p->arr[2], &tmp);
935 emul(&e1->x.p->arr[1], res);
936 eadd(&tmp, &res->x.p->arr[2]);
937 free_evalue_refs(&tmp);
938 }
939 } else {
940 if(mod_term_smaller(res, e1))
941 for(i=1; i<res->x.p->size ; i++)
942 emul(e1, &(res->x.p->arr[i]));
943 else
944 emul_rev(e1, res);
945 return;
946 }
947 }
948 break;
949 case relation:
950 emul_rev(e1, res);
951 break;
952 default:
953 assert(0);
954 }
955 }
956 else {
957 if (value_notzero_p(e1->d)&& value_notzero_p(res->d)) {
958 /* Product of two rational numbers */
959
960 Value g;
961 value_init(g);
962 value_multiply(res->d,e1->d,res->d);
963 value_multiply(res->x.n,e1->x.n,res->x.n );
964 Gcd(res->x.n, res->d,&g);
965 if (value_notone_p(g)) {
966 value_division(res->d,res->d,g);
967 value_division(res->x.n,res->x.n,g);
968 }
969 value_clear(g);
970 return ;
971 }
972 else {
973 if(value_zero_p(e1->d)&& value_notzero_p(res->d)) {
974 /* Product of an expression (polynomial or peririodic) and a rational number */
975
976 emul_rev(e1, res);
977 return ;
978 }
979 else {
980 /* Product of a rationel number and an expression (polynomial or peririodic) */
981
982 i = res->x.p->type == modulo ? 1 : 0;
983 for (; i<res->x.p->size; i++)
984 emul(e1, &res->x.p->arr[i]);
985
986 return ;
987 }
988 }
989 }
990
991 return ;
992 }
993
994 /* Frees mask ! */
995 void emask(evalue *mask, evalue *res) {
996 int n, i, j;
997 Polyhedron *d, *fd;
998 struct section *s;
999 evalue mone;
1000
1001 assert(mask->x.p->type == partition);
1002 assert(res->x.p->type == partition);
1003
1004 s = (struct section *)
1005 malloc((mask->x.p->size/2+1) * (res->x.p->size/2) *
1006 sizeof(struct section));
1007 assert(s);
1008
1009 value_init(mone.d);
1010 evalue_set_si(&mone, -1, 1);
1011
1012 n = 0;
1013 for (j = 0; j < res->x.p->size/2; ++j) {
1014 assert(mask->x.p->size >= 2);
1015 fd = DomainDifference(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1016 EVALUE_DOMAIN(mask->x.p->arr[0]), 0);
1017 if (!emptyQ(fd))
1018 for (i = 1; i < mask->x.p->size/2; ++i) {
1019 Polyhedron *t = fd;
1020 fd = DomainDifference(fd, EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1021 Domain_Free(t);
1022 if (emptyQ(fd))
1023 break;
1024 }
1025 if (emptyQ(fd)) {
1026 Domain_Free(fd);
1027 continue;
1028 }
1029 value_init(s[n].E.d);
1030 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1031 s[n].D = fd;
1032 ++n;
1033 }
1034 for (i = 0; i < mask->x.p->size/2; ++i) {
1035 fd = EVALUE_DOMAIN(mask->x.p->arr[2*i]);
1036 eadd(&mone, &mask->x.p->arr[2*i+1]);
1037 emul(&mone, &mask->x.p->arr[2*i+1]);
1038 for (j = 0; j < res->x.p->size/2; ++j) {
1039 Polyhedron *t;
1040 d = DomainIntersection(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1041 EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1042 if (emptyQ(d)) {
1043 Domain_Free(d);
1044 continue;
1045 }
1046 t = fd;
1047 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*j]), 0);
1048 if (t != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1049 Domain_Free(t);
1050 value_init(s[n].E.d);
1051 evalue_copy(&s[n].E, &mask->x.p->arr[2*i+1]);
1052 emul(&res->x.p->arr[2*j+1], &s[n].E);
1053 s[n].D = d;
1054 ++n;
1055 }
1056 if (!emptyQ(fd)) {
1057 assert(0); // We don't allow this.
1058 }
1059 }
1060
1061 free_evalue_refs(&mone);
1062 free_evalue_refs(mask);
1063 free_evalue_refs(res);
1064 value_init(res->d);
1065 res->x.p = new_enode(partition, 2*n, -1);
1066 for (j = 0; j < n; ++j) {
1067 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1068 value_clear(res->x.p->arr[2*j+1].d);
1069 res->x.p->arr[2*j+1] = s[j].E;
1070 }
1071
1072 free(s);
1073 }
1074
1075 void evalue_copy(evalue *dst, evalue *src)
1076 {
1077 value_assign(dst->d, src->d);
1078 if(value_notzero_p(src->d)) {
1079 value_init(dst->x.n);
1080 value_assign(dst->x.n, src->x.n);
1081 } else
1082 dst->x.p = ecopy(src->x.p);
1083 }
1084
1085 enode *new_enode(enode_type type,int size,int pos) {
1086
1087 enode *res;
1088 int i;
1089
1090 if(size == 0) {
1091 fprintf(stderr, "Allocating enode of size 0 !\n" );
1092 return NULL;
1093 }
1094 res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
1095 res->type = type;
1096 res->size = size;
1097 res->pos = pos;
1098 for(i=0; i<size; i++) {
1099 value_init(res->arr[i].d);
1100 value_set_si(res->arr[i].d,0);
1101 res->arr[i].x.p = 0;
1102 }
1103 return res;
1104 } /* new_enode */
1105
1106 enode *ecopy(enode *e) {
1107
1108 enode *res;
1109 int i;
1110
1111 res = new_enode(e->type,e->size,e->pos);
1112 for(i=0;i<e->size;++i) {
1113 value_assign(res->arr[i].d,e->arr[i].d);
1114 if(value_zero_p(res->arr[i].d))
1115 res->arr[i].x.p = ecopy(e->arr[i].x.p);
1116 else if (EVALUE_IS_DOMAIN(res->arr[i]))
1117 EVALUE_SET_DOMAIN(res->arr[i], Domain_Copy(EVALUE_DOMAIN(e->arr[i])));
1118 else {
1119 value_init(res->arr[i].x.n);
1120 value_assign(res->arr[i].x.n,e->arr[i].x.n);
1121 }
1122 }
1123 return(res);
1124 } /* ecopy */
1125
1126 int eequal(evalue *e1,evalue *e2) {
1127
1128 int i;
1129 enode *p1, *p2;
1130
1131 if (value_ne(e1->d,e2->d))
1132 return 0;
1133
1134 /* e1->d == e2->d */
1135 if (value_notzero_p(e1->d)) {
1136 if (value_ne(e1->x.n,e2->x.n))
1137 return 0;
1138
1139 /* e1->d == e2->d != 0 AND e1->n == e2->n */
1140 return 1;
1141 }
1142
1143 /* e1->d == e2->d == 0 */
1144 p1 = e1->x.p;
1145 p2 = e2->x.p;
1146 if (p1->type != p2->type) return 0;
1147 if (p1->size != p2->size) return 0;
1148 if (p1->pos != p2->pos) return 0;
1149 for (i=0; i<p1->size; i++)
1150 if (!eequal(&p1->arr[i], &p2->arr[i]) )
1151 return 0;
1152 return 1;
1153 } /* eequal */
1154
1155 void free_evalue_refs(evalue *e) {
1156
1157 enode *p;
1158 int i;
1159
1160 if (EVALUE_IS_DOMAIN(*e)) {
1161 Domain_Free(EVALUE_DOMAIN(*e));
1162 value_clear(e->d);
1163 return;
1164 } else if (value_pos_p(e->d)) {
1165
1166 /* 'e' stores a constant */
1167 value_clear(e->d);
1168 value_clear(e->x.n);
1169 return;
1170 }
1171 assert(value_zero_p(e->d));
1172 value_clear(e->d);
1173 p = e->x.p;
1174 if (!p) return; /* null pointer */
1175 for (i=0; i<p->size; i++) {
1176 free_evalue_refs(&(p->arr[i]));
1177 }
1178 free(p);
1179 return;
1180 } /* free_evalue_refs */
1181
1182 static void mod2table_r(evalue *e, Vector *periods, Value m, int p,
1183 Vector * val, evalue *res)
1184 {
1185 unsigned nparam = periods->Size;
1186
1187 if (p == nparam) {
1188 double d = compute_evalue(e, val->p);
1189 if (d > 0)
1190 d += .25;
1191 else
1192 d -= .25;
1193 value_set_si(res->d, 1);
1194 value_init(res->x.n);
1195 value_set_double(res->x.n, d);
1196 mpz_fdiv_r(res->x.n, res->x.n, m);
1197 return;
1198 }
1199 if (value_one_p(periods->p[p]))
1200 mod2table_r(e, periods, m, p+1, val, res);
1201 else {
1202 Value tmp;
1203 value_init(tmp);
1204
1205 value_assign(tmp, periods->p[p]);
1206 value_set_si(res->d, 0);
1207 res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
1208 do {
1209 value_decrement(tmp, tmp);
1210 value_assign(val->p[p], tmp);
1211 mod2table_r(e, periods, m, p+1, val,
1212 &res->x.p->arr[VALUE_TO_INT(tmp)]);
1213 } while (value_pos_p(tmp));
1214
1215 value_clear(tmp);
1216 }
1217 }
1218
1219 void evalue_mod2table(evalue *e, int nparam)
1220 {
1221 enode *p;
1222 int i;
1223
1224 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
1225 return;
1226 p = e->x.p;
1227 for (i=0; i<p->size; i++) {
1228 evalue_mod2table(&(p->arr[i]), nparam);
1229 }
1230 if (p->type == modulo) {
1231 Vector *periods = Vector_Alloc(nparam);
1232 Vector *val = Vector_Alloc(nparam);
1233 Value tmp;
1234 evalue *ev;
1235 evalue EP, res;
1236
1237 value_init(tmp);
1238 value_set_si(tmp, p->pos);
1239 Vector_Set(periods->p, 1, nparam);
1240 Vector_Set(val->p, 0, nparam);
1241 for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
1242 enode *p = ev->x.p;
1243
1244 assert(p->type == polynomial);
1245 assert(p->size == 2);
1246 assert(value_one_p(p->arr[1].d));
1247 Gcd(tmp, p->arr[1].x.n, &periods->p[p->pos-1]);
1248 value_division(periods->p[p->pos-1], tmp, periods->p[p->pos-1]);
1249 }
1250 value_set_si(tmp, p->pos);
1251 value_init(EP.d);
1252 mod2table_r(&p->arr[0], periods, tmp, 0, val, &EP);
1253
1254 value_init(res.d);
1255 evalue_set_si(&res, 0, 1);
1256 /* Compute the polynomial using Horner's rule */
1257 for (i=p->size-1;i>1;i--) {
1258 eadd(&p->arr[i], &res);
1259 emul(&EP, &res);
1260 }
1261 eadd(&p->arr[1], &res);
1262
1263 free_evalue_refs(e);
1264 free_evalue_refs(&EP);
1265 *e = res;
1266
1267 value_clear(tmp);
1268 Vector_Free(val);
1269 Vector_Free(periods);
1270 }
1271 } /* evalue_mod2table */
1272
1273 /********************************************************/
1274 /* function in domain */
1275 /* check if the parameters in list_args */
1276 /* verifies the constraints of Domain P */
1277 /********************************************************/
1278 int in_domain(Polyhedron *P, Value *list_args) {
1279
1280 int col,row;
1281 Value v; /* value of the constraint of a row when
1282 parameters are instanciated*/
1283 Value tmp;
1284
1285 value_init(v);
1286 value_init(tmp);
1287
1288 /*P->Constraint constraint matrice of polyhedron P*/
1289 for(row=0;row<P->NbConstraints;row++) {
1290 value_assign(v,P->Constraint[row][P->Dimension+1]); /*constant part*/
1291 for(col=1;col<P->Dimension+1;col++) {
1292 value_multiply(tmp,P->Constraint[row][col],list_args[col-1]);
1293 value_addto(v,v,tmp);
1294 }
1295 if (value_notzero_p(P->Constraint[row][0])) {
1296
1297 /*if v is not >=0 then this constraint is not respected */
1298 if (value_neg_p(v)) {
1299 next:
1300 value_clear(v);
1301 value_clear(tmp);
1302 return P->next ? in_domain(P->next, list_args) : 0;
1303 }
1304 }
1305 else {
1306
1307 /*if v is not = 0 then this constraint is not respected */
1308 if (value_notzero_p(v))
1309 goto next;
1310 }
1311 }
1312
1313 /*if not return before this point => all
1314 the constraints are respected */
1315 value_clear(v);
1316 value_clear(tmp);
1317 return 1;
1318 } /* in_domain */
1319
1320 /****************************************************/
1321 /* function compute enode */
1322 /* compute the value of enode p with parameters */
1323 /* list "list_args */
1324 /* compute the polynomial or the periodic */
1325 /****************************************************/
1326
1327 static double compute_enode(enode *p, Value *list_args) {
1328
1329 int i;
1330 Value m, param;
1331 double res=0.0;
1332
1333 if (!p)
1334 return(0.);
1335
1336 value_init(m);
1337 value_init(param);
1338
1339 if (p->type == polynomial) {
1340 if (p->size > 1)
1341 value_assign(param,list_args[p->pos-1]);
1342
1343 /* Compute the polynomial using Horner's rule */
1344 for (i=p->size-1;i>0;i--) {
1345 res +=compute_evalue(&p->arr[i],list_args);
1346 res *=VALUE_TO_DOUBLE(param);
1347 }
1348 res +=compute_evalue(&p->arr[0],list_args);
1349 }
1350 else if (p->type == modulo) {
1351 double d = compute_evalue(&p->arr[0], list_args);
1352 if (d > 0)
1353 d += .25;
1354 else
1355 d -= .25;
1356 value_set_double(param, d);
1357 value_set_si(m, p->pos);
1358 mpz_fdiv_r(param, param, m);
1359
1360 /* Compute the polynomial using Horner's rule */
1361 for (i=p->size-1;i>1;i--) {
1362 res +=compute_evalue(&p->arr[i],list_args);
1363 res *=VALUE_TO_DOUBLE(param);
1364 }
1365 res +=compute_evalue(&p->arr[1],list_args);
1366 }
1367 else if (p->type == periodic) {
1368 value_assign(param,list_args[p->pos-1]);
1369
1370 /* Choose the right element of the periodic */
1371 value_absolute(m,param);
1372 value_set_si(param,p->size);
1373 value_modulus(m,m,param);
1374 res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
1375 }
1376 else if (p->type == relation) {
1377 if (fabs(compute_evalue(&p->arr[0], list_args)) < 0.5)
1378 res = compute_evalue(&p->arr[1], list_args);
1379 else if (p->size > 2)
1380 res = compute_evalue(&p->arr[2], list_args);
1381 }
1382 else if (p->type == partition) {
1383 for (i = 0; i < p->size/2; ++i)
1384 if (in_domain(EVALUE_DOMAIN(p->arr[2*i]), list_args)) {
1385 res = compute_evalue(&p->arr[2*i+1], list_args);
1386 break;
1387 }
1388 }
1389 value_clear(m);
1390 value_clear(param);
1391 return res;
1392 } /* compute_enode */
1393
1394 /*************************************************/
1395 /* return the value of Ehrhart Polynomial */
1396 /* It returns a double, because since it is */
1397 /* a recursive function, some intermediate value */
1398 /* might not be integral */
1399 /*************************************************/
1400
1401 double compute_evalue(evalue *e,Value *list_args) {
1402
1403 double res;
1404
1405 if (value_notzero_p(e->d)) {
1406 if (value_notone_p(e->d))
1407 res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
1408 else
1409 res = VALUE_TO_DOUBLE(e->x.n);
1410 }
1411 else
1412 res = compute_enode(e->x.p,list_args);
1413 return res;
1414 } /* compute_evalue */
1415
1416
1417 /****************************************************/
1418 /* function compute_poly : */
1419 /* Check for the good validity domain */
1420 /* return the number of point in the Polyhedron */
1421 /* in allocated memory */
1422 /* Using the Ehrhart pseudo-polynomial */
1423 /****************************************************/
1424 Value *compute_poly(Enumeration *en,Value *list_args) {
1425
1426 Value *tmp;
1427 /* double d; int i; */
1428
1429 tmp = (Value *) malloc (sizeof(Value));
1430 assert(tmp != NULL);
1431 value_init(*tmp);
1432 value_set_si(*tmp,0);
1433
1434 if(!en)
1435 return(tmp); /* no ehrhart polynomial */
1436 if(en->ValidityDomain) {
1437 if(!en->ValidityDomain->Dimension) { /* no parameters */
1438 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
1439 return(tmp);
1440 }
1441 }
1442 else
1443 return(tmp); /* no Validity Domain */
1444 while(en) {
1445 if(in_domain(en->ValidityDomain,list_args)) {
1446
1447 #ifdef EVAL_EHRHART_DEBUG
1448 Print_Domain(stdout,en->ValidityDomain);
1449 print_evalue(stdout,&en->EP);
1450 #endif
1451
1452 /* d = compute_evalue(&en->EP,list_args);
1453 i = d;
1454 printf("(double)%lf = %d\n", d, i ); */
1455 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
1456 return(tmp);
1457 }
1458 else
1459 en=en->next;
1460 }
1461 value_set_si(*tmp,0);
1462 return(tmp); /* no compatible domain with the arguments */
1463 } /* compute_poly */
1464
lattice point enumeration library