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