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