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