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