?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
regression test
[barvinok.git] / ev_operations.c
blob0beab61c53d445ee085980c72ab4329024eeb4c7
1 #include <assert.h>
2 #include <math.h>
3 #include <stdlib.h>
4
5 #include "ev_operations.h"
6 #include "util.h"
7
8 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
9
10 void evalue_set_si(evalue *ev, int n, int d) {
11 value_set_si(ev->d, d);
12 value_init(ev->x.n);
13 value_set_si(ev->x.n, n);
14 }
15
16 void evalue_set(evalue *ev, Value n, Value d) {
17 value_assign(ev->d, d);
18 value_init(ev->x.n);
19 value_assign(ev->x.n, n);
20 }
21
22 void aep_evalue(evalue *e, int *ref) {
23
24 enode *p;
25 int i;
26
27 if (value_notzero_p(e->d))
28 return; /* a rational number, its already reduced */
29 if(!(p = e->x.p))
30 return; /* hum... an overflow probably occured */
31
32 /* First check the components of p */
33 for (i=0;i<p->size;i++)
34 aep_evalue(&p->arr[i],ref);
35
36 /* Then p itself */
37 switch (p->type) {
38 case polynomial:
39 case periodic:
40 case evector:
41 p->pos = ref[p->pos-1]+1;
42 }
43 return;
44 } /* aep_evalue */
45
46 /** Comments */
47 void addeliminatedparams_evalue(evalue *e,Matrix *CT) {
48
49 enode *p;
50 int i, j;
51 int *ref;
52
53 if (value_notzero_p(e->d))
54 return; /* a rational number, its already reduced */
55 if(!(p = e->x.p))
56 return; /* hum... an overflow probably occured */
57
58 /* Compute ref */
59 ref = (int *)malloc(sizeof(int)*(CT->NbRows-1));
60 for(i=0;i<CT->NbRows-1;i++)
61 for(j=0;j<CT->NbColumns;j++)
62 if(value_notzero_p(CT->p[i][j])) {
63 ref[i] = j;
64 break;
65 }
66
67 /* Transform the references in e, using ref */
68 aep_evalue(e,ref);
69 free( ref );
70 return;
71 } /* addeliminatedparams_evalue */
72
73 static int mod_rational_smaller(evalue *e1, evalue *e2)
74 {
75 int r;
76 Value m;
77 value_init(m);
78
79 assert(value_notzero_p(e1->d));
80 assert(value_notzero_p(e2->d));
81 value_multiply(m, e1->x.n, e2->d);
82 value_division(m, m, e1->d);
83 if (value_lt(m, e2->x.n))
84 r = 1;
85 else if (value_gt(m, e2->x.n))
86 r = 0;
87 else
88 r = -1;
89 value_clear(m);
90
91 return r;
92 }
93
94 static int mod_term_smaller_r(evalue *e1, evalue *e2)
95 {
96 if (value_notzero_p(e1->d)) {
97 if (value_zero_p(e2->d))
98 return 1;
99 int r = mod_rational_smaller(e1, e2);
100 return r == -1 ? 0 : r;
101 }
102 if (value_notzero_p(e2->d))
103 return 0;
104 if (e1->x.p->pos < e2->x.p->pos)
105 return 1;
106 else if (e1->x.p->pos > e2->x.p->pos)
107 return 0;
108 else {
109 int r = mod_rational_smaller(&e1->x.p->arr[1], &e2->x.p->arr[1]);
110 return r == -1
111 ? mod_term_smaller_r(&e1->x.p->arr[0], &e2->x.p->arr[0])
112 : r;
113 }
114 }
115
116 static int mod_term_smaller(evalue *e1, evalue *e2)
117 {
118 assert(value_zero_p(e1->d));
119 assert(value_zero_p(e2->d));
120 assert(e1->x.p->type == fractional);
121 assert(e2->x.p->type == fractional);
122 return mod_term_smaller_r(&e1->x.p->arr[0], &e2->x.p->arr[0]);
123 }
124
125 /* Negative pos means inequality */
126 /* s is negative of substitution if m is not zero */
127 struct fixed_param {
128 int pos;
129 evalue s;
130 Value d;
131 Value m;
132 };
133
134 struct subst {
135 struct fixed_param *fixed;
136 int n;
137 int max;
138 };
139
140 static int relations_depth(evalue *e)
141 {
142 int d;
143
144 for (d = 0;
145 value_zero_p(e->d) && e->x.p->type == relation;
146 e = &e->x.p->arr[1], ++d);
147 return d;
148 }
149
150 static void Lcm3(Value i, Value j, Value *res)
151 {
152 Value aux;
153
154 value_init(aux);
155 Gcd(i,j,&aux);
156 value_multiply(*res,i,j);
157 value_division(*res, *res, aux);
158 value_clear(aux);
159 }
160
161 static void poly_denom(evalue *p, Value *d)
162 {
163 value_set_si(*d, 1);
164
165 while (value_zero_p(p->d)) {
166 assert(p->x.p->type == polynomial);
167 assert(p->x.p->size == 2);
168 assert(value_notzero_p(p->x.p->arr[1].d));
169 Lcm3(*d, p->x.p->arr[1].d, d);
170 p = &p->x.p->arr[0];
171 }
172 Lcm3(*d, p->d, d);
173 }
174
175 #define EVALUE_IS_ONE(ev) (value_pos_p((ev).d) && value_one_p((ev).x.n))
176
177 static void realloc_substitution(struct subst *s, int d)
178 {
179 struct fixed_param *n;
180 int i;
181 NALLOC(n, s->max+d);
182 for (i = 0; i < s->n; ++i)
183 n[i] = s->fixed[i];
184 free(s->fixed);
185 s->fixed = n;
186 s->max += d;
187 }
188
189 static int add_modulo_substitution(struct subst *s, evalue *r)
190 {
191 evalue *p;
192 evalue *f;
193 evalue *m;
194
195 assert(value_zero_p(r->d) && r->x.p->type == relation);
196 m = &r->x.p->arr[0];
197
198 /* May have been reduced already */
199 if (value_notzero_p(m->d))
200 return 0;
201
202 assert(value_zero_p(m->d) && m->x.p->type == fractional);
203 assert(m->x.p->size == 3);
204
205 /* fractional was inverted during reduction
206 * invert it back and move constant in
207 */
208 if (!EVALUE_IS_ONE(m->x.p->arr[2])) {
209 assert(value_pos_p(m->x.p->arr[2].d));
210 assert(value_mone_p(m->x.p->arr[2].x.n));
211 value_set_si(m->x.p->arr[2].x.n, 1);
212 value_increment(m->x.p->arr[1].x.n, m->x.p->arr[1].x.n);
213 assert(value_eq(m->x.p->arr[1].x.n, m->x.p->arr[1].d));
214 value_set_si(m->x.p->arr[1].x.n, 1);
215 eadd(&m->x.p->arr[1], &m->x.p->arr[0]);
216 value_set_si(m->x.p->arr[1].x.n, 0);
217 value_set_si(m->x.p->arr[1].d, 1);
218 }
219
220 /* Oops. Nested identical relations. */
221 if (!EVALUE_IS_ZERO(m->x.p->arr[1]))
222 return 0;
223
224 if (s->n >= s->max) {
225 int d = relations_depth(r);
226 realloc_substitution(s, d);
227 }
228
229 p = &m->x.p->arr[0];
230 assert(value_zero_p(p->d) && p->x.p->type == polynomial);
231 assert(p->x.p->size == 2);
232 f = &p->x.p->arr[1];
233
234 assert(value_pos_p(f->x.n));
235
236 value_init(s->fixed[s->n].m);
237 value_assign(s->fixed[s->n].m, f->d);
238 s->fixed[s->n].pos = p->x.p->pos;
239 value_init(s->fixed[s->n].d);
240 value_assign(s->fixed[s->n].d, f->x.n);
241 value_init(s->fixed[s->n].s.d);
242 evalue_copy(&s->fixed[s->n].s, &p->x.p->arr[0]);
243 ++s->n;
244
245 return 1;
246 }
247
248 static void reorder_terms(enode *p, evalue *v)
249 {
250 int i;
251 int offset = p->type == fractional;
252
253 for (i = p->size-1; i >= offset+1; i--) {
254 emul(v, &p->arr[i]);
255 eadd(&p->arr[i], &p->arr[i-1]);
256 free_evalue_refs(&(p->arr[i]));
257 }
258 p->size = offset+1;
259 free_evalue_refs(v);
260 }
261
262 void _reduce_evalue (evalue *e, struct subst *s, int fract) {
263
264 enode *p;
265 int i, j, k;
266 int add;
267
268 if (value_notzero_p(e->d)) {
269 if (fract)
270 mpz_fdiv_r(e->x.n, e->x.n, e->d);
271 return; /* a rational number, its already reduced */
272 }
273
274 if(!(p = e->x.p))
275 return; /* hum... an overflow probably occured */
276
277 /* First reduce the components of p */
278 add = p->type == relation;
279 for (i=0; i<p->size; i++) {
280 if (add && i == 1)
281 add = add_modulo_substitution(s, e);
282
283 if (i == 0 && p->type==fractional)
284 _reduce_evalue(&p->arr[i], s, 1);
285 else
286 _reduce_evalue(&p->arr[i], s, fract);
287
288 if (add && i == p->size-1) {
289 --s->n;
290 value_clear(s->fixed[s->n].m);
291 value_clear(s->fixed[s->n].d);
292 free_evalue_refs(&s->fixed[s->n].s);
293 } else if (add && i == 1)
294 s->fixed[s->n-1].pos *= -1;
295 }
296
297 if (p->type==periodic) {
298
299 /* Try to reduce the period */
300 for (i=1; i<=(p->size)/2; i++) {
301 if ((p->size % i)==0) {
302
303 /* Can we reduce the size to i ? */
304 for (j=0; j<i; j++)
305 for (k=j+i; k<e->x.p->size; k+=i)
306 if (!eequal(&p->arr[j], &p->arr[k])) goto you_lose;
307
308 /* OK, lets do it */
309 for (j=i; j<p->size; j++) free_evalue_refs(&p->arr[j]);
310 p->size = i;
311 break;
312
313 you_lose: /* OK, lets not do it */
314 continue;
315 }
316 }
317
318 /* Try to reduce its strength */
319 if (p->size == 1) {
320 value_clear(e->d);
321 memcpy(e,&p->arr[0],sizeof(evalue));
322 free(p);
323 }
324 }
325 else if (p->type==polynomial) {
326 for (k = 0; s && k < s->n; ++k) {
327 if (s->fixed[k].pos == p->pos) {
328 int divide = value_notone_p(s->fixed[k].d);
329 evalue d;
330
331 if (value_notzero_p(s->fixed[k].m)) {
332 if (!fract)
333 continue;
334 assert(p->size == 2);
335 if (divide && value_ne(s->fixed[k].d, p->arr[1].x.n))
336 continue;
337 if (!mpz_divisible_p(s->fixed[k].m, p->arr[1].d))
338 continue;
339 divide = 1;
340 }
341
342 if (divide) {
343 value_init(d.d);
344 value_assign(d.d, s->fixed[k].d);
345 value_init(d.x.n);
346 if (value_notzero_p(s->fixed[k].m))
347 value_oppose(d.x.n, s->fixed[k].m);
348 else
349 value_set_si(d.x.n, 1);
350 }
351
352 for (i=p->size-1;i>=1;i--) {
353 emul(&s->fixed[k].s, &p->arr[i]);
354 if (divide)
355 emul(&d, &p->arr[i]);
356 eadd(&p->arr[i], &p->arr[i-1]);
357 free_evalue_refs(&(p->arr[i]));
358 }
359 p->size = 1;
360 _reduce_evalue(&p->arr[0], s, fract);
361
362 if (divide)
363 free_evalue_refs(&d);
364
365 break;
366 }
367 }
368
369 /* Try to reduce the degree */
370 for (i=p->size-1;i>=1;i--) {
371 if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
372 break;
373 /* Zero coefficient */
374 free_evalue_refs(&(p->arr[i]));
375 }
376 if (i+1<p->size)
377 p->size = i+1;
378
379 /* Try to reduce its strength */
380 if (p->size == 1) {
381 value_clear(e->d);
382 memcpy(e,&p->arr[0],sizeof(evalue));
383 free(p);
384 }
385 }
386 else if (p->type==fractional) {
387 int reorder = 0;
388 evalue v;
389
390 if (value_notzero_p(p->arr[0].d)) {
391 value_init(v.d);
392 value_assign(v.d, p->arr[0].d);
393 value_init(v.x.n);
394 mpz_fdiv_r(v.x.n, p->arr[0].x.n, p->arr[0].d);
395
396 reorder = 1;
397 } else {
398 evalue *f, *base;
399 evalue *pp = &p->arr[0];
400 assert(value_zero_p(pp->d) && pp->x.p->type == polynomial);
401 assert(pp->x.p->size == 2);
402
403 /* search for exact duplicate among the modulo inequalities */
404 do {
405 f = &pp->x.p->arr[1];
406 for (k = 0; s && k < s->n; ++k) {
407 if (-s->fixed[k].pos == pp->x.p->pos &&
408 value_eq(s->fixed[k].d, f->x.n) &&
409 value_eq(s->fixed[k].m, f->d) &&
410 eequal(&s->fixed[k].s, &pp->x.p->arr[0]))
411 break;
412 }
413 if (k < s->n) {
414 Value g;
415 value_init(g);
416
417 /* replace { E/m } by { (E-1)/m } + 1/m */
418 poly_denom(pp, &g);
419 if (reorder) {
420 evalue extra;
421 value_init(extra.d);
422 evalue_set_si(&extra, 1, 1);
423 value_assign(extra.d, g);
424 eadd(&extra, &v.x.p->arr[1]);
425 free_evalue_refs(&extra);
426
427 /* We've been going in circles; stop now */
428 if (value_ge(v.x.p->arr[1].x.n, v.x.p->arr[1].d)) {
429 free_evalue_refs(&v);
430 value_init(v.d);
431 evalue_set_si(&v, 0, 1);
432 break;
433 }
434 } else {
435 value_init(v.d);
436 value_set_si(v.d, 0);
437 v.x.p = new_enode(fractional, 3, -1);
438 evalue_set_si(&v.x.p->arr[1], 1, 1);
439 value_assign(v.x.p->arr[1].d, g);
440 evalue_set_si(&v.x.p->arr[2], 1, 1);
441 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
442 }
443
444 for (f = &v.x.p->arr[0]; value_zero_p(f->d);
445 f = &f->x.p->arr[0])
446 ;
447 value_division(f->d, g, f->d);
448 value_multiply(f->x.n, f->x.n, f->d);
449 value_assign(f->d, g);
450 value_decrement(f->x.n, f->x.n);
451 mpz_fdiv_r(f->x.n, f->x.n, f->d);
452
453 Gcd(f->d, f->x.n, &g);
454 value_division(f->d, f->d, g);
455 value_division(f->x.n, f->x.n, g);
456
457 value_clear(g);
458 pp = &v.x.p->arr[0];
459
460 reorder = 1;
461 }
462 } while (k < s->n);
463
464 /* reduction may have made this fractional arg smaller */
465 i = reorder ? p->size : 1;
466 for ( ; i < p->size; ++i)
467 if (value_zero_p(p->arr[i].d) &&
468 p->arr[i].x.p->type == fractional &&
469 !mod_term_smaller(e, &p->arr[i]))
470 break;
471 if (i < p->size) {
472 value_init(v.d);
473 value_set_si(v.d, 0);
474 v.x.p = new_enode(fractional, 3, -1);
475 evalue_set_si(&v.x.p->arr[1], 0, 1);
476 evalue_set_si(&v.x.p->arr[2], 1, 1);
477 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
478
479 reorder = 1;
480 }
481
482 if (!reorder) {
483 int invert = 0;
484 Value twice;
485 value_init(twice);
486
487 for (pp = &p->arr[0]; value_zero_p(pp->d);
488 pp = &pp->x.p->arr[0]) {
489 f = &pp->x.p->arr[1];
490 assert(value_pos_p(f->d));
491 mpz_mul_ui(twice, f->x.n, 2);
492 if (value_lt(twice, f->d))
493 break;
494 if (value_eq(twice, f->d))
495 continue;
496 invert = 1;
497 break;
498 }
499
500 if (invert) {
501 value_init(v.d);
502 value_set_si(v.d, 0);
503 v.x.p = new_enode(fractional, 3, -1);
504 evalue_set_si(&v.x.p->arr[1], 0, 1);
505 poly_denom(&p->arr[0], &twice);
506 value_assign(v.x.p->arr[1].d, twice);
507 value_decrement(v.x.p->arr[1].x.n, twice);
508 evalue_set_si(&v.x.p->arr[2], -1, 1);
509 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
510
511 for (pp = &v.x.p->arr[0]; value_zero_p(pp->d);
512 pp = &pp->x.p->arr[0]) {
513 f = &pp->x.p->arr[1];
514 value_oppose(f->x.n, f->x.n);
515 mpz_fdiv_r(f->x.n, f->x.n, f->d);
516 }
517 value_division(pp->d, twice, pp->d);
518 value_multiply(pp->x.n, pp->x.n, pp->d);
519 value_assign(pp->d, twice);
520 value_oppose(pp->x.n, pp->x.n);
521 value_decrement(pp->x.n, pp->x.n);
522 mpz_fdiv_r(pp->x.n, pp->x.n, pp->d);
523
524 reorder = 1;
525 }
526
527 value_clear(twice);
528 }
529 }
530
531 if (reorder) {
532 reorder_terms(p, &v);
533 _reduce_evalue(&p->arr[1], s, fract);
534 }
535
536 /* Try to reduce the degree */
537 for (i=p->size-1;i>=2;i--) {
538 if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
539 break;
540 /* Zero coefficient */
541 free_evalue_refs(&(p->arr[i]));
542 }
543 if (i+1<p->size)
544 p->size = i+1;
545
546 /* Try to reduce its strength */
547 if (p->size == 2) {
548 value_clear(e->d);
549 memcpy(e,&p->arr[1],sizeof(evalue));
550 free_evalue_refs(&(p->arr[0]));
551 free(p);
552 }
553 }
554 else if (p->type == relation) {
555 if (p->size == 3 && eequal(&p->arr[1], &p->arr[2])) {
556 free_evalue_refs(&(p->arr[2]));
557 free_evalue_refs(&(p->arr[0]));
558 p->size = 2;
559 value_clear(e->d);
560 *e = p->arr[1];
561 free(p);
562 return;
563 }
564 if (p->size == 3 && EVALUE_IS_ZERO(p->arr[2])) {
565 free_evalue_refs(&(p->arr[2]));
566 p->size = 2;
567 }
568 if (p->size == 2 && EVALUE_IS_ZERO(p->arr[1])) {
569 free_evalue_refs(&(p->arr[1]));
570 free_evalue_refs(&(p->arr[0]));
571 evalue_set_si(e, 0, 1);
572 free(p);
573 } else {
574 int reduced = 0;
575 evalue *m;
576 m = &p->arr[0];
577
578 /* Relation was reduced by means of an identical
579 * inequality => remove
580 */
581 if (value_zero_p(m->d) && !EVALUE_IS_ZERO(m->x.p->arr[1]))
582 reduced = 1;
583
584 if (reduced || value_notzero_p(p->arr[0].d)) {
585 if (!reduced && value_zero_p(p->arr[0].x.n)) {
586 value_clear(e->d);
587 memcpy(e,&p->arr[1],sizeof(evalue));
588 if (p->size == 3)
589 free_evalue_refs(&(p->arr[2]));
590 } else {
591 if (p->size == 3) {
592 value_clear(e->d);
593 memcpy(e,&p->arr[2],sizeof(evalue));
594 } else
595 evalue_set_si(e, 0, 1);
596 free_evalue_refs(&(p->arr[1]));
597 }
598 free_evalue_refs(&(p->arr[0]));
599 free(p);
600 }
601 }
602 }
603 return;
604 } /* reduce_evalue */
605
606 static void add_substitution(struct subst *s, Value *row, unsigned dim)
607 {
608 int k, l;
609 evalue *r;
610
611 for (k = 0; k < dim; ++k)
612 if (value_notzero_p(row[k+1]))
613 break;
614
615 Vector_Normalize_Positive(row+1, dim+1, k);
616 assert(s->n < s->max);
617 value_init(s->fixed[s->n].d);
618 value_init(s->fixed[s->n].m);
619 value_assign(s->fixed[s->n].d, row[k+1]);
620 s->fixed[s->n].pos = k+1;
621 value_set_si(s->fixed[s->n].m, 0);
622 r = &s->fixed[s->n].s;
623 value_init(r->d);
624 for (l = k+1; l < dim; ++l)
625 if (value_notzero_p(row[l+1])) {
626 value_set_si(r->d, 0);
627 r->x.p = new_enode(polynomial, 2, l + 1);
628 value_init(r->x.p->arr[1].x.n);
629 value_oppose(r->x.p->arr[1].x.n, row[l+1]);
630 value_set_si(r->x.p->arr[1].d, 1);
631 r = &r->x.p->arr[0];
632 }
633 value_init(r->x.n);
634 value_oppose(r->x.n, row[dim+1]);
635 value_set_si(r->d, 1);
636 ++s->n;
637 }
638
639 void reduce_evalue (evalue *e) {
640 struct subst s = { NULL, 0, 0 };
641
642 if (value_notzero_p(e->d))
643 return; /* a rational number, its already reduced */
644
645 if (e->x.p->type == partition) {
646 int i;
647 unsigned dim = -1;
648 for (i = 0; i < e->x.p->size/2; ++i) {
649 s.n = 0;
650 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
651 /* This shouldn't really happen;
652 * Empty domains should not be added.
653 */
654 if (emptyQ(D))
655 goto discard;
656
657 dim = D->Dimension;
658 if (D->next)
659 D = DomainConvex(D, 0);
660 if (!D->next && D->NbEq) {
661 int j, k;
662 if (s.max < dim) {
663 if (s.max != 0)
664 realloc_substitution(&s, dim);
665 else {
666 int d = relations_depth(&e->x.p->arr[2*i+1]);
667 s.max = dim+d;
668 NALLOC(s.fixed, s.max);
669 }
670 }
671 for (j = 0; j < D->NbEq; ++j)
672 add_substitution(&s, D->Constraint[j], dim);
673 }
674 if (D != EVALUE_DOMAIN(e->x.p->arr[2*i]))
675 Domain_Free(D);
676 _reduce_evalue(&e->x.p->arr[2*i+1], &s, 0);
677 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
678 discard:
679 free_evalue_refs(&e->x.p->arr[2*i+1]);
680 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
681 value_clear(e->x.p->arr[2*i].d);
682 e->x.p->size -= 2;
683 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
684 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
685 --i;
686 }
687 if (s.n != 0) {
688 int j;
689 for (j = 0; j < s.n; ++j) {
690 value_clear(s.fixed[j].d);
691 free_evalue_refs(&s.fixed[j].s);
692 }
693 }
694 }
695 if (e->x.p->size == 0) {
696 free(e->x.p);
697 evalue_set_si(e, 0, 1);
698 }
699 if (s.max != 0)
700 free(s.fixed);
701 } else
702 _reduce_evalue(e, &s, 0);
703 }
704
705 void print_evalue(FILE *DST,evalue *e,char **pname) {
706
707 if(value_notzero_p(e->d)) {
708 if(value_notone_p(e->d)) {
709 value_print(DST,VALUE_FMT,e->x.n);
710 fprintf(DST,"/");
711 value_print(DST,VALUE_FMT,e->d);
712 }
713 else {
714 value_print(DST,VALUE_FMT,e->x.n);
715 }
716 }
717 else
718 print_enode(DST,e->x.p,pname);
719 return;
720 } /* print_evalue */
721
722 void print_enode(FILE *DST,enode *p,char **pname) {
723
724 int i;
725
726 if (!p) {
727 fprintf(DST, "NULL");
728 return;
729 }
730 switch (p->type) {
731 case evector:
732 fprintf(DST, "{ ");
733 for (i=0; i<p->size; i++) {
734 print_evalue(DST, &p->arr[i], pname);
735 if (i!=(p->size-1))
736 fprintf(DST, ", ");
737 }
738 fprintf(DST, " }\n");
739 break;
740 case polynomial:
741 fprintf(DST, "( ");
742 for (i=p->size-1; i>=0; i--) {
743 print_evalue(DST, &p->arr[i], pname);
744 if (i==1) fprintf(DST, " * %s + ", pname[p->pos-1]);
745 else if (i>1)
746 fprintf(DST, " * %s^%d + ", pname[p->pos-1], i);
747 }
748 fprintf(DST, " )\n");
749 break;
750 case periodic:
751 fprintf(DST, "[ ");
752 for (i=0; i<p->size; i++) {
753 print_evalue(DST, &p->arr[i], pname);
754 if (i!=(p->size-1)) fprintf(DST, ", ");
755 }
756 fprintf(DST," ]_%s", pname[p->pos-1]);
757 break;
758 case fractional:
759 fprintf(DST, "( ");
760 for (i=p->size-1; i>=1; i--) {
761 print_evalue(DST, &p->arr[i], pname);
762 if (i >= 2) {
763 fprintf(DST, " * ");
764 fprintf(DST, "{");
765 print_evalue(DST, &p->arr[0], pname);
766 fprintf(DST, "}");
767 if (i>2)
768 fprintf(DST, "^%d + ", i-1);
769 else
770 fprintf(DST, " + ");
771 }
772 }
773 fprintf(DST, " )\n");
774 break;
775 case relation:
776 fprintf(DST, "[ ");
777 print_evalue(DST, &p->arr[0], pname);
778 fprintf(DST, "= 0 ] * \n");
779 print_evalue(DST, &p->arr[1], pname);
780 if (p->size > 2) {
781 fprintf(DST, " +\n [ ");
782 print_evalue(DST, &p->arr[0], pname);
783 fprintf(DST, "!= 0 ] * \n");
784 print_evalue(DST, &p->arr[2], pname);
785 }
786 break;
787 case partition:
788 for (i=0; i<p->size/2; i++) {
789 Print_Domain(DST, EVALUE_DOMAIN(p->arr[2*i]), pname);
790 print_evalue(DST, &p->arr[2*i+1], pname);
791 }
792 break;
793 default:
794 assert(0);
795 }
796 return;
797 } /* print_enode */
798
799 static void eadd_rev(evalue *e1, evalue *res)
800 {
801 evalue ev;
802 value_init(ev.d);
803 evalue_copy(&ev, e1);
804 eadd(res, &ev);
805 free_evalue_refs(res);
806 *res = ev;
807 }
808
809 static void eadd_rev_cst (evalue *e1, evalue *res)
810 {
811 evalue ev;
812 value_init(ev.d);
813 evalue_copy(&ev, e1);
814 eadd(res, &ev.x.p->arr[ev.x.p->type==fractional]);
815 free_evalue_refs(res);
816 *res = ev;
817 }
818
819 struct section { Polyhedron * D; evalue E; };
820
821 void eadd_partitions (evalue *e1,evalue *res)
822 {
823 int n, i, j;
824 Polyhedron *d, *fd;
825 struct section *s;
826 s = (struct section *)
827 malloc((e1->x.p->size/2+1) * (res->x.p->size/2+1) *
828 sizeof(struct section));
829 assert(s);
830
831 n = 0;
832 for (j = 0; j < e1->x.p->size/2; ++j) {
833 assert(res->x.p->size >= 2);
834 fd = DomainDifference(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
835 EVALUE_DOMAIN(res->x.p->arr[0]), 0);
836 if (!emptyQ(fd))
837 for (i = 1; i < res->x.p->size/2; ++i) {
838 Polyhedron *t = fd;
839 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
840 Domain_Free(t);
841 if (emptyQ(fd))
842 break;
843 }
844 if (emptyQ(fd)) {
845 Domain_Free(fd);
846 continue;
847 }
848 value_init(s[n].E.d);
849 evalue_copy(&s[n].E, &e1->x.p->arr[2*j+1]);
850 s[n].D = fd;
851 ++n;
852 }
853 for (i = 0; i < res->x.p->size/2; ++i) {
854 fd = EVALUE_DOMAIN(res->x.p->arr[2*i]);
855 for (j = 0; j < e1->x.p->size/2; ++j) {
856 Polyhedron *t;
857 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
858 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
859 if (emptyQ(d)) {
860 Domain_Free(d);
861 continue;
862 }
863 t = fd;
864 fd = DomainDifference(fd, EVALUE_DOMAIN(e1->x.p->arr[2*j]), 0);
865 if (t != EVALUE_DOMAIN(res->x.p->arr[2*i]))
866 Domain_Free(t);
867 value_init(s[n].E.d);
868 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
869 eadd(&e1->x.p->arr[2*j+1], &s[n].E);
870 s[n].D = d;
871 ++n;
872 }
873 if (!emptyQ(fd)) {
874 s[n].E = res->x.p->arr[2*i+1];
875 s[n].D = fd;
876 ++n;
877 } else {
878 free_evalue_refs(&res->x.p->arr[2*i+1]);
879 Domain_Free(fd);
880 }
881 if (fd != EVALUE_DOMAIN(res->x.p->arr[2*i]))
882 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
883 value_clear(res->x.p->arr[2*i].d);
884 }
885
886 free(res->x.p);
887 res->x.p = new_enode(partition, 2*n, -1);
888 for (j = 0; j < n; ++j) {
889 s[j].D = DomainConstraintSimplify(s[j].D, 0);
890 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
891 value_clear(res->x.p->arr[2*j+1].d);
892 res->x.p->arr[2*j+1] = s[j].E;
893 }
894
895 free(s);
896 }
897
898 static void explicit_complement(evalue *res)
899 {
900 enode *rel = new_enode(relation, 3, 0);
901 assert(rel);
902 value_clear(rel->arr[0].d);
903 rel->arr[0] = res->x.p->arr[0];
904 value_clear(rel->arr[1].d);
905 rel->arr[1] = res->x.p->arr[1];
906 value_set_si(rel->arr[2].d, 1);
907 value_init(rel->arr[2].x.n);
908 value_set_si(rel->arr[2].x.n, 0);
909 free(res->x.p);
910 res->x.p = rel;
911 }
912
913 void eadd(evalue *e1,evalue *res) {
914
915 int i;
916 if (value_notzero_p(e1->d) && value_notzero_p(res->d)) {
917 /* Add two rational numbers */
918 Value g,m1,m2;
919 value_init(g);
920 value_init(m1);
921 value_init(m2);
922
923 value_multiply(m1,e1->x.n,res->d);
924 value_multiply(m2,res->x.n,e1->d);
925 value_addto(res->x.n,m1,m2);
926 value_multiply(res->d,e1->d,res->d);
927 Gcd(res->x.n,res->d,&g);
928 if (value_notone_p(g)) {
929 value_division(res->d,res->d,g);
930 value_division(res->x.n,res->x.n,g);
931 }
932 value_clear(g); value_clear(m1); value_clear(m2);
933 return ;
934 }
935 else if (value_notzero_p(e1->d) && value_zero_p(res->d)) {
936 switch (res->x.p->type) {
937 case polynomial:
938 /* Add the constant to the constant term of a polynomial*/
939 eadd(e1, &res->x.p->arr[0]);
940 return ;
941 case periodic:
942 /* Add the constant to all elements of a periodic number */
943 for (i=0; i<res->x.p->size; i++) {
944 eadd(e1, &res->x.p->arr[i]);
945 }
946 return ;
947 case evector:
948 fprintf(stderr, "eadd: cannot add const with vector\n");
949 return;
950 case fractional:
951 eadd(e1, &res->x.p->arr[1]);
952 return ;
953 case partition:
954 assert(EVALUE_IS_ZERO(*e1));
955 break; /* Do nothing */
956 case relation:
957 /* Create (zero) complement if needed */
958 if (res->x.p->size < 3 && !EVALUE_IS_ZERO(*e1))
959 explicit_complement(res);
960 for (i = 1; i < res->x.p->size; ++i)
961 eadd(e1, &res->x.p->arr[i]);
962 break;
963 default:
964 assert(0);
965 }
966 }
967 /* add polynomial or periodic to constant
968 * you have to exchange e1 and res, before doing addition */
969
970 else if (value_zero_p(e1->d) && value_notzero_p(res->d)) {
971 eadd_rev(e1, res);
972 return;
973 }
974 else { // ((e1->d==0) && (res->d==0))
975 assert(!((e1->x.p->type == partition) ^
976 (res->x.p->type == partition)));
977 if (e1->x.p->type == partition) {
978 eadd_partitions(e1, res);
979 return;
980 }
981 if (e1->x.p->type == relation &&
982 (res->x.p->type != relation ||
983 mod_term_smaller(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
984 eadd_rev(e1, res);
985 return;
986 }
987 if (res->x.p->type == relation) {
988 if (e1->x.p->type == relation &&
989 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
990 if (res->x.p->size < 3 && e1->x.p->size == 3)
991 explicit_complement(res);
992 if (e1->x.p->size < 3 && res->x.p->size == 3)
993 explicit_complement(e1);
994 for (i = 1; i < res->x.p->size; ++i)
995 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
996 return;
997 }
998 if (res->x.p->size < 3)
999 explicit_complement(res);
1000 for (i = 1; i < res->x.p->size; ++i)
1001 eadd(e1, &res->x.p->arr[i]);
1002 return;
1003 }
1004 if ((e1->x.p->type != res->x.p->type) ) {
1005 /* adding to evalues of different type. two cases are possible
1006 * res is periodic and e1 is polynomial, you have to exchange
1007 * e1 and res then to add e1 to the constant term of res */
1008 if (e1->x.p->type == polynomial) {
1009 eadd_rev_cst(e1, res);
1010 }
1011 else if (res->x.p->type == polynomial) {
1012 /* res is polynomial and e1 is periodic,
1013 add e1 to the constant term of res */
1014
1015 eadd(e1,&res->x.p->arr[0]);
1016 } else
1017 assert(0);
1018
1019 return;
1020 }
1021 else if (e1->x.p->pos != res->x.p->pos ||
1022 (res->x.p->type == fractional &&
1023 !eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1024 /* adding evalues of different position (i.e function of different unknowns
1025 * to case are possible */
1026
1027 switch (res->x.p->type) {
1028 case fractional:
1029 if(mod_term_smaller(res, e1))
1030 eadd(e1,&res->x.p->arr[1]);
1031 else
1032 eadd_rev_cst(e1, res);
1033 return;
1034 case polynomial: // res and e1 are polynomials
1035 // add e1 to the constant term of res
1036
1037 if(res->x.p->pos < e1->x.p->pos)
1038 eadd(e1,&res->x.p->arr[0]);
1039 else
1040 eadd_rev_cst(e1, res);
1041 // value_clear(g); value_clear(m1); value_clear(m2);
1042 return;
1043 case periodic: // res and e1 are pointers to periodic numbers
1044 //add e1 to all elements of res
1045
1046 if(res->x.p->pos < e1->x.p->pos)
1047 for (i=0;i<res->x.p->size;i++) {
1048 eadd(e1,&res->x.p->arr[i]);
1049 }
1050 else
1051 eadd_rev(e1, res);
1052 return;
1053 }
1054 }
1055
1056
1057 //same type , same pos and same size
1058 if (e1->x.p->size == res->x.p->size) {
1059 // add any element in e1 to the corresponding element in res
1060 if (res->x.p->type == fractional)
1061 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
1062 i = res->x.p->type == fractional ? 1 : 0;
1063 for (; i<res->x.p->size; i++) {
1064 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1065 }
1066 return ;
1067 }
1068
1069 /* Sizes are different */
1070 switch(res->x.p->type) {
1071 case polynomial:
1072 case fractional:
1073 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1074 /* new enode and add res to that new node. If you do not do */
1075 /* that, you lose the the upper weight part of e1 ! */
1076
1077 if(e1->x.p->size > res->x.p->size)
1078 eadd_rev(e1, res);
1079 else {
1080
1081 if (res->x.p->type == fractional)
1082 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
1083 i = res->x.p->type == fractional ? 1 : 0;
1084 for (; i<e1->x.p->size ; i++) {
1085 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1086 }
1087
1088 return ;
1089 }
1090 break;
1091
1092 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1093 case periodic:
1094 {
1095 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1096 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1097 to add periodicaly elements of e1 to elements of ne, and finaly to
1098 return ne. */
1099 int x,y,p;
1100 Value ex, ey ,ep;
1101 evalue *ne;
1102 value_init(ex); value_init(ey);value_init(ep);
1103 x=e1->x.p->size;
1104 y= res->x.p->size;
1105 value_set_si(ex,e1->x.p->size);
1106 value_set_si(ey,res->x.p->size);
1107 value_assign (ep,*Lcm(ex,ey));
1108 p=(int)mpz_get_si(ep);
1109 ne= (evalue *) malloc (sizeof(evalue));
1110 value_init(ne->d);
1111 value_set_si( ne->d,0);
1112
1113 ne->x.p=new_enode(res->x.p->type,p, res->x.p->pos);
1114 for(i=0;i<p;i++) {
1115 value_assign(ne->x.p->arr[i].d, res->x.p->arr[i%y].d);
1116 if (value_notzero_p(ne->x.p->arr[i].d)) {
1117 value_init(ne->x.p->arr[i].x.n);
1118 value_assign(ne->x.p->arr[i].x.n, res->x.p->arr[i%y].x.n);
1119 }
1120 else {
1121 ne->x.p->arr[i].x.p =ecopy(res->x.p->arr[i%y].x.p);
1122 }
1123 }
1124 for(i=0;i<p;i++) {
1125 eadd(&e1->x.p->arr[i%x], &ne->x.p->arr[i]);
1126 }
1127
1128 value_assign(res->d, ne->d);
1129 res->x.p=ne->x.p;
1130
1131 return ;
1132 }
1133 case evector:
1134 fprintf(stderr, "eadd: ?cannot add vectors of different length\n");
1135 return ;
1136 }
1137 }
1138 return ;
1139 }/* eadd */
1140
1141 static void emul_rev (evalue *e1, evalue *res)
1142 {
1143 evalue ev;
1144 value_init(ev.d);
1145 evalue_copy(&ev, e1);
1146 emul(res, &ev);
1147 free_evalue_refs(res);
1148 *res = ev;
1149 }
1150
1151 static void emul_poly (evalue *e1, evalue *res)
1152 {
1153 int i, j, o = res->x.p->type == fractional;
1154 evalue tmp;
1155 int size=(e1->x.p->size + res->x.p->size - o - 1);
1156 value_init(tmp.d);
1157 value_set_si(tmp.d,0);
1158 tmp.x.p=new_enode(res->x.p->type, size, res->x.p->pos);
1159 if (o)
1160 evalue_copy(&tmp.x.p->arr[0], &e1->x.p->arr[0]);
1161 for (i=o; i < e1->x.p->size; i++) {
1162 evalue_copy(&tmp.x.p->arr[i], &e1->x.p->arr[i]);
1163 emul(&res->x.p->arr[o], &tmp.x.p->arr[i]);
1164 }
1165 for (; i<size; i++)
1166 evalue_set_si(&tmp.x.p->arr[i], 0, 1);
1167 for (i=o+1; i<res->x.p->size; i++)
1168 for (j=o; j<e1->x.p->size; j++) {
1169 evalue ev;
1170 value_init(ev.d);
1171 evalue_copy(&ev, &e1->x.p->arr[j]);
1172 emul(&res->x.p->arr[i], &ev);
1173 eadd(&ev, &tmp.x.p->arr[i+j-o]);
1174 free_evalue_refs(&ev);
1175 }
1176 free_evalue_refs(res);
1177 *res = tmp;
1178 }
1179
1180 void emul_partitions (evalue *e1,evalue *res)
1181 {
1182 int n, i, j, k;
1183 Polyhedron *d;
1184 struct section *s;
1185 s = (struct section *)
1186 malloc((e1->x.p->size/2) * (res->x.p->size/2) *
1187 sizeof(struct section));
1188 assert(s);
1189
1190 n = 0;
1191 for (i = 0; i < res->x.p->size/2; ++i) {
1192 for (j = 0; j < e1->x.p->size/2; ++j) {
1193 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1194 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1195 if (emptyQ(d)) {
1196 Domain_Free(d);
1197 continue;
1198 }
1199
1200 /* This code is only needed because the partitions
1201 are not true partitions.
1202 */
1203 for (k = 0; k < n; ++k) {
1204 if (DomainIncludes(s[k].D, d))
1205 break;
1206 if (DomainIncludes(d, s[k].D)) {
1207 Domain_Free(s[k].D);
1208 free_evalue_refs(&s[k].E);
1209 if (n > k)
1210 s[k] = s[--n];
1211 --k;
1212 }
1213 }
1214 if (k < n) {
1215 Domain_Free(d);
1216 continue;
1217 }
1218
1219 value_init(s[n].E.d);
1220 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
1221 emul(&e1->x.p->arr[2*j+1], &s[n].E);
1222 s[n].D = d;
1223 ++n;
1224 }
1225 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
1226 value_clear(res->x.p->arr[2*i].d);
1227 free_evalue_refs(&res->x.p->arr[2*i+1]);
1228 }
1229
1230 free(res->x.p);
1231 res->x.p = new_enode(partition, 2*n, -1);
1232 for (j = 0; j < n; ++j) {
1233 s[j].D = DomainConstraintSimplify(s[j].D, 0);
1234 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1235 value_clear(res->x.p->arr[2*j+1].d);
1236 res->x.p->arr[2*j+1] = s[j].E;
1237 }
1238
1239 free(s);
1240 }
1241
1242 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1243
1244 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1245 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1246 * evalues is not treated here */
1247
1248 void emul (evalue *e1, evalue *res ){
1249 int i,j;
1250
1251 if((value_zero_p(e1->d)&&e1->x.p->type==evector)||(value_zero_p(res->d)&&(res->x.p->type==evector))) {
1252 fprintf(stderr, "emul: do not proced on evector type !\n");
1253 return;
1254 }
1255
1256 if (value_zero_p(e1->d) && e1->x.p->type == partition) {
1257 if (value_zero_p(res->d) && res->x.p->type == partition)
1258 emul_partitions(e1, res);
1259 else
1260 emul_rev(e1, res);
1261 } else if (value_zero_p(res->d) && res->x.p->type == partition) {
1262 for (i = 0; i < res->x.p->size/2; ++i)
1263 emul(e1, &res->x.p->arr[2*i+1]);
1264 } else
1265 if (value_zero_p(res->d) && res->x.p->type == relation) {
1266 if (value_zero_p(e1->d) && e1->x.p->type == relation &&
1267 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1268 if (res->x.p->size < 3 && e1->x.p->size == 3)
1269 explicit_complement(res);
1270 if (e1->x.p->size < 3 && res->x.p->size == 3)
1271 explicit_complement(e1);
1272 for (i = 1; i < res->x.p->size; ++i)
1273 emul(&e1->x.p->arr[i], &res->x.p->arr[i]);
1274 return;
1275 }
1276 for (i = 1; i < res->x.p->size; ++i)
1277 emul(e1, &res->x.p->arr[i]);
1278 } else
1279 if(value_zero_p(e1->d)&& value_zero_p(res->d)) {
1280 switch(e1->x.p->type) {
1281 case polynomial:
1282 switch(res->x.p->type) {
1283 case polynomial:
1284 if(e1->x.p->pos == res->x.p->pos) {
1285 /* Product of two polynomials of the same variable */
1286 emul_poly(e1, res);
1287 return;
1288 }
1289 else {
1290 /* Product of two polynomials of different variables */
1291
1292 if(res->x.p->pos < e1->x.p->pos)
1293 for( i=0; i<res->x.p->size ; i++)
1294 emul(e1, &res->x.p->arr[i]);
1295 else
1296 emul_rev(e1, res);
1297
1298 return ;
1299 }
1300 case periodic:
1301 case fractional:
1302 /* Product of a polynomial and a periodic or fractional */
1303 emul_rev(e1, res);
1304 return;
1305 }
1306 case periodic:
1307 switch(res->x.p->type) {
1308 case periodic:
1309 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size==res->x.p->size) {
1310 /* Product of two periodics of the same parameter and period */
1311
1312 for(i=0; i<res->x.p->size;i++)
1313 emul(&(e1->x.p->arr[i]), &(res->x.p->arr[i]));
1314
1315 return;
1316 }
1317 else{
1318 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size!=res->x.p->size) {
1319 /* Product of two periodics of the same parameter and different periods */
1320 evalue *newp;
1321 Value x,y,z;
1322 int ix,iy,lcm;
1323 value_init(x); value_init(y);value_init(z);
1324 ix=e1->x.p->size;
1325 iy=res->x.p->size;
1326 value_set_si(x,e1->x.p->size);
1327 value_set_si(y,res->x.p->size);
1328 value_assign (z,*Lcm(x,y));
1329 lcm=(int)mpz_get_si(z);
1330 newp= (evalue *) malloc (sizeof(evalue));
1331 value_init(newp->d);
1332 value_set_si( newp->d,0);
1333 newp->x.p=new_enode(periodic,lcm, e1->x.p->pos);
1334 for(i=0;i<lcm;i++) {
1335 value_assign(newp->x.p->arr[i].d, res->x.p->arr[i%iy].d);
1336 if (value_notzero_p(newp->x.p->arr[i].d)) {
1337 value_assign(newp->x.p->arr[i].x.n, res->x.p->arr[i%iy].x.n);
1338 }
1339 else {
1340 newp->x.p->arr[i].x.p =ecopy(res->x.p->arr[i%iy].x.p);
1341 }
1342
1343 }
1344 for(i=0;i<lcm;i++)
1345 emul(&e1->x.p->arr[i%ix], &newp->x.p->arr[i]);
1346
1347 value_assign(res->d,newp->d);
1348 res->x.p=newp->x.p;
1349
1350 value_clear(x); value_clear(y);value_clear(z);
1351 return ;
1352 }
1353 else {
1354 /* Product of two periodics of different parameters */
1355
1356 for(i=0; i<res->x.p->size; i++)
1357 emul(e1, &(res->x.p->arr[i]));
1358
1359 return;
1360 }
1361 }
1362 case polynomial:
1363 /* Product of a periodic and a polynomial */
1364
1365 for(i=0; i<res->x.p->size ; i++)
1366 emul(e1, &(res->x.p->arr[i]));
1367
1368 return;
1369
1370 }
1371 case fractional:
1372 switch(res->x.p->type) {
1373 case polynomial:
1374 for(i=0; i<res->x.p->size ; i++)
1375 emul(e1, &(res->x.p->arr[i]));
1376 return;
1377 case periodic:
1378 assert(0);
1379 case fractional:
1380 if (e1->x.p->pos == res->x.p->pos &&
1381 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1382 evalue d;
1383 value_init(d.d);
1384 poly_denom(&e1->x.p->arr[0], &d.d);
1385 if (!value_two_p(d.d))
1386 emul_poly(e1, res);
1387 else {
1388 value_init(d.x.n);
1389 value_set_si(d.x.n, 1);
1390 evalue tmp;
1391 /* { x }^2 == { x }/2 */
1392 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1393 assert(e1->x.p->size == 3);
1394 assert(res->x.p->size == 3);
1395 value_init(tmp.d);
1396 evalue_copy(&tmp, &res->x.p->arr[2]);
1397 emul(&d, &tmp);
1398 eadd(&res->x.p->arr[1], &tmp);
1399 emul(&e1->x.p->arr[2], &tmp);
1400 emul(&e1->x.p->arr[1], res);
1401 eadd(&tmp, &res->x.p->arr[2]);
1402 free_evalue_refs(&tmp);
1403 value_clear(d.x.n);
1404 }
1405 value_clear(d.d);
1406 } else {
1407 if(mod_term_smaller(res, e1))
1408 for(i=1; i<res->x.p->size ; i++)
1409 emul(e1, &(res->x.p->arr[i]));
1410 else
1411 emul_rev(e1, res);
1412 return;
1413 }
1414 }
1415 break;
1416 case relation:
1417 emul_rev(e1, res);
1418 break;
1419 default:
1420 assert(0);
1421 }
1422 }
1423 else {
1424 if (value_notzero_p(e1->d)&& value_notzero_p(res->d)) {
1425 /* Product of two rational numbers */
1426
1427 Value g;
1428 value_init(g);
1429 value_multiply(res->d,e1->d,res->d);
1430 value_multiply(res->x.n,e1->x.n,res->x.n );
1431 Gcd(res->x.n, res->d,&g);
1432 if (value_notone_p(g)) {
1433 value_division(res->d,res->d,g);
1434 value_division(res->x.n,res->x.n,g);
1435 }
1436 value_clear(g);
1437 return ;
1438 }
1439 else {
1440 if(value_zero_p(e1->d)&& value_notzero_p(res->d)) {
1441 /* Product of an expression (polynomial or peririodic) and a rational number */
1442
1443 emul_rev(e1, res);
1444 return ;
1445 }
1446 else {
1447 /* Product of a rationel number and an expression (polynomial or peririodic) */
1448
1449 i = res->x.p->type == fractional ? 1 : 0;
1450 for (; i<res->x.p->size; i++)
1451 emul(e1, &res->x.p->arr[i]);
1452
1453 return ;
1454 }
1455 }
1456 }
1457
1458 return ;
1459 }
1460
1461 /* Frees mask ! */
1462 void emask(evalue *mask, evalue *res) {
1463 int n, i, j;
1464 Polyhedron *d, *fd;
1465 struct section *s;
1466 evalue mone;
1467
1468 if (EVALUE_IS_ZERO(*res))
1469 return;
1470
1471 assert(value_zero_p(mask->d));
1472 assert(mask->x.p->type == partition);
1473 assert(value_zero_p(res->d));
1474 assert(res->x.p->type == partition);
1475
1476 s = (struct section *)
1477 malloc((mask->x.p->size/2+1) * (res->x.p->size/2) *
1478 sizeof(struct section));
1479 assert(s);
1480
1481 value_init(mone.d);
1482 evalue_set_si(&mone, -1, 1);
1483
1484 n = 0;
1485 for (j = 0; j < res->x.p->size/2; ++j) {
1486 assert(mask->x.p->size >= 2);
1487 fd = DomainDifference(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1488 EVALUE_DOMAIN(mask->x.p->arr[0]), 0);
1489 if (!emptyQ(fd))
1490 for (i = 1; i < mask->x.p->size/2; ++i) {
1491 Polyhedron *t = fd;
1492 fd = DomainDifference(fd, EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1493 Domain_Free(t);
1494 if (emptyQ(fd))
1495 break;
1496 }
1497 if (emptyQ(fd)) {
1498 Domain_Free(fd);
1499 continue;
1500 }
1501 value_init(s[n].E.d);
1502 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1503 s[n].D = fd;
1504 ++n;
1505 }
1506 for (i = 0; i < mask->x.p->size/2; ++i) {
1507 if (EVALUE_IS_ONE(mask->x.p->arr[2*i+1]))
1508 continue;
1509
1510 fd = EVALUE_DOMAIN(mask->x.p->arr[2*i]);
1511 eadd(&mone, &mask->x.p->arr[2*i+1]);
1512 emul(&mone, &mask->x.p->arr[2*i+1]);
1513 for (j = 0; j < res->x.p->size/2; ++j) {
1514 Polyhedron *t;
1515 d = DomainIntersection(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1516 EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1517 if (emptyQ(d)) {
1518 Domain_Free(d);
1519 continue;
1520 }
1521 t = fd;
1522 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*j]), 0);
1523 if (t != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1524 Domain_Free(t);
1525 value_init(s[n].E.d);
1526 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1527 emul(&mask->x.p->arr[2*i+1], &s[n].E);
1528 s[n].D = d;
1529 ++n;
1530 }
1531
1532 if (!emptyQ(fd)) {
1533 /* Just ignore; this may have been previously masked off */
1534 }
1535 }
1536
1537 free_evalue_refs(&mone);
1538 free_evalue_refs(mask);
1539 free_evalue_refs(res);
1540 value_init(res->d);
1541 if (n == 0)
1542 evalue_set_si(res, 0, 1);
1543 else {
1544 res->x.p = new_enode(partition, 2*n, -1);
1545 for (j = 0; j < n; ++j) {
1546 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1547 value_clear(res->x.p->arr[2*j+1].d);
1548 res->x.p->arr[2*j+1] = s[j].E;
1549 }
1550 }
1551
1552 free(s);
1553 }
1554
1555 void evalue_copy(evalue *dst, evalue *src)
1556 {
1557 value_assign(dst->d, src->d);
1558 if(value_notzero_p(src->d)) {
1559 value_init(dst->x.n);
1560 value_assign(dst->x.n, src->x.n);
1561 } else
1562 dst->x.p = ecopy(src->x.p);
1563 }
1564
1565 enode *new_enode(enode_type type,int size,int pos) {
1566
1567 enode *res;
1568 int i;
1569
1570 if(size == 0) {
1571 fprintf(stderr, "Allocating enode of size 0 !\n" );
1572 return NULL;
1573 }
1574 res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
1575 res->type = type;
1576 res->size = size;
1577 res->pos = pos;
1578 for(i=0; i<size; i++) {
1579 value_init(res->arr[i].d);
1580 value_set_si(res->arr[i].d,0);
1581 res->arr[i].x.p = 0;
1582 }
1583 return res;
1584 } /* new_enode */
1585
1586 enode *ecopy(enode *e) {
1587
1588 enode *res;
1589 int i;
1590
1591 res = new_enode(e->type,e->size,e->pos);
1592 for(i=0;i<e->size;++i) {
1593 value_assign(res->arr[i].d,e->arr[i].d);
1594 if(value_zero_p(res->arr[i].d))
1595 res->arr[i].x.p = ecopy(e->arr[i].x.p);
1596 else if (EVALUE_IS_DOMAIN(res->arr[i]))
1597 EVALUE_SET_DOMAIN(res->arr[i], Domain_Copy(EVALUE_DOMAIN(e->arr[i])));
1598 else {
1599 value_init(res->arr[i].x.n);
1600 value_assign(res->arr[i].x.n,e->arr[i].x.n);
1601 }
1602 }
1603 return(res);
1604 } /* ecopy */
1605
1606 int ecmp(const evalue *e1, const evalue *e2)
1607 {
1608 enode *p1, *p2;
1609 int i;
1610 int r;
1611
1612 if (value_notzero_p(e1->d) && value_notzero_p(e2->d)) {
1613 Value m, m2;
1614 value_init(m);
1615 value_init(m2);
1616 value_multiply(m, e1->x.n, e2->d);
1617 value_multiply(m2, e2->x.n, e1->d);
1618
1619 if (value_lt(m, m2))
1620 r = -1;
1621 else if (value_gt(m, m2))
1622 r = 1;
1623 else
1624 r = 0;
1625
1626 value_clear(m);
1627 value_clear(m2);
1628
1629 return r;
1630 }
1631 if (value_notzero_p(e1->d))
1632 return -1;
1633 if (value_notzero_p(e2->d))
1634 return 1;
1635
1636 p1 = e1->x.p;
1637 p2 = e2->x.p;
1638
1639 if (p1->type != p2->type)
1640 return p1->type - p2->type;
1641 if (p1->pos != p2->pos)
1642 return p1->pos - p2->pos;
1643 if (p1->size != p2->size)
1644 return p1->size - p2->size;
1645
1646 for (i = p1->size-1; i >= 0; --i)
1647 if ((r = ecmp(&p1->arr[i], &p2->arr[i])) != 0)
1648 return r;
1649
1650 return 0;
1651 }
1652
1653 int eequal(evalue *e1,evalue *e2) {
1654
1655 int i;
1656 enode *p1, *p2;
1657
1658 if (value_ne(e1->d,e2->d))
1659 return 0;
1660
1661 /* e1->d == e2->d */
1662 if (value_notzero_p(e1->d)) {
1663 if (value_ne(e1->x.n,e2->x.n))
1664 return 0;
1665
1666 /* e1->d == e2->d != 0 AND e1->n == e2->n */
1667 return 1;
1668 }
1669
1670 /* e1->d == e2->d == 0 */
1671 p1 = e1->x.p;
1672 p2 = e2->x.p;
1673 if (p1->type != p2->type) return 0;
1674 if (p1->size != p2->size) return 0;
1675 if (p1->pos != p2->pos) return 0;
1676 for (i=0; i<p1->size; i++)
1677 if (!eequal(&p1->arr[i], &p2->arr[i]) )
1678 return 0;
1679 return 1;
1680 } /* eequal */
1681
1682 void free_evalue_refs(evalue *e) {
1683
1684 enode *p;
1685 int i;
1686
1687 if (EVALUE_IS_DOMAIN(*e)) {
1688 Domain_Free(EVALUE_DOMAIN(*e));
1689 value_clear(e->d);
1690 return;
1691 } else if (value_pos_p(e->d)) {
1692
1693 /* 'e' stores a constant */
1694 value_clear(e->d);
1695 value_clear(e->x.n);
1696 return;
1697 }
1698 assert(value_zero_p(e->d));
1699 value_clear(e->d);
1700 p = e->x.p;
1701 if (!p) return; /* null pointer */
1702 for (i=0; i<p->size; i++) {
1703 free_evalue_refs(&(p->arr[i]));
1704 }
1705 free(p);
1706 return;
1707 } /* free_evalue_refs */
1708
1709 static void mod2table_r(evalue *e, Vector *periods, Value m, int p,
1710 Vector * val, evalue *res)
1711 {
1712 unsigned nparam = periods->Size;
1713
1714 if (p == nparam) {
1715 double d = compute_evalue(e, val->p);
1716 d *= VALUE_TO_DOUBLE(m);
1717 if (d > 0)
1718 d += .25;
1719 else
1720 d -= .25;
1721 value_assign(res->d, m);
1722 value_init(res->x.n);
1723 value_set_double(res->x.n, d);
1724 mpz_fdiv_r(res->x.n, res->x.n, m);
1725 return;
1726 }
1727 if (value_one_p(periods->p[p]))
1728 mod2table_r(e, periods, m, p+1, val, res);
1729 else {
1730 Value tmp;
1731 value_init(tmp);
1732
1733 value_assign(tmp, periods->p[p]);
1734 value_set_si(res->d, 0);
1735 res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
1736 do {
1737 value_decrement(tmp, tmp);
1738 value_assign(val->p[p], tmp);
1739 mod2table_r(e, periods, m, p+1, val,
1740 &res->x.p->arr[VALUE_TO_INT(tmp)]);
1741 } while (value_pos_p(tmp));
1742
1743 value_clear(tmp);
1744 }
1745 }
1746
1747 static void rel2table(evalue *e, int zero)
1748 {
1749 if (value_pos_p(e->d)) {
1750 if (value_zero_p(e->x.n) == zero)
1751 value_set_si(e->x.n, 1);
1752 else
1753 value_set_si(e->x.n, 0);
1754 value_set_si(e->d, 1);
1755 } else {
1756 int i;
1757 for (i = 0; i < e->x.p->size; ++i)
1758 rel2table(&e->x.p->arr[i], zero);
1759 }
1760 }
1761
1762 void evalue_mod2table(evalue *e, int nparam)
1763 {
1764 enode *p;
1765 int i;
1766
1767 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
1768 return;
1769 p = e->x.p;
1770 for (i=0; i<p->size; i++) {
1771 evalue_mod2table(&(p->arr[i]), nparam);
1772 }
1773 if (p->type == relation) {
1774 evalue copy;
1775 evalue *ev;
1776
1777 if (p->size > 2) {
1778 value_init(copy.d);
1779 evalue_copy(&copy, &p->arr[0]);
1780 }
1781 rel2table(&p->arr[0], 1);
1782 emul(&p->arr[0], &p->arr[1]);
1783 if (p->size > 2) {
1784 rel2table(&copy, 0);
1785 emul(&copy, &p->arr[2]);
1786 eadd(&p->arr[2], &p->arr[1]);
1787 free_evalue_refs(&p->arr[2]);
1788 free_evalue_refs(&copy);
1789 }
1790 free_evalue_refs(&p->arr[0]);
1791 ev = &p->arr[1];
1792 free(p);
1793 value_clear(e->d);
1794 *e = *ev;
1795 } else if (p->type == fractional) {
1796 Vector *periods = Vector_Alloc(nparam);
1797 Vector *val = Vector_Alloc(nparam);
1798 Value tmp;
1799 evalue *ev;
1800 evalue EP, res;
1801
1802 value_init(tmp);
1803 value_set_si(tmp, 1);
1804 Vector_Set(periods->p, 1, nparam);
1805 Vector_Set(val->p, 0, nparam);
1806 for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
1807 enode *p = ev->x.p;
1808
1809 assert(p->type == polynomial);
1810 assert(p->size == 2);
1811 value_assign(periods->p[p->pos-1], p->arr[1].d);
1812 Lcm3(tmp, p->arr[1].d, &tmp);
1813 }
1814 value_init(EP.d);
1815 mod2table_r(&p->arr[0], periods, tmp, 0, val, &EP);
1816
1817 value_init(res.d);
1818 evalue_set_si(&res, 0, 1);
1819 /* Compute the polynomial using Horner's rule */
1820 for (i=p->size-1;i>1;i--) {
1821 eadd(&p->arr[i], &res);
1822 emul(&EP, &res);
1823 }
1824 eadd(&p->arr[1], &res);
1825
1826 free_evalue_refs(e);
1827 free_evalue_refs(&EP);
1828 *e = res;
1829
1830 value_clear(tmp);
1831 Vector_Free(val);
1832 Vector_Free(periods);
1833 }
1834 } /* evalue_mod2table */
1835
1836 /********************************************************/
1837 /* function in domain */
1838 /* check if the parameters in list_args */
1839 /* verifies the constraints of Domain P */
1840 /********************************************************/
1841 int in_domain(Polyhedron *P, Value *list_args) {
1842
1843 int col,row;
1844 Value v; /* value of the constraint of a row when
1845 parameters are instanciated*/
1846 Value tmp;
1847
1848 value_init(v);
1849 value_init(tmp);
1850
1851 /*P->Constraint constraint matrice of polyhedron P*/
1852 for(row=0;row<P->NbConstraints;row++) {
1853 value_assign(v,P->Constraint[row][P->Dimension+1]); /*constant part*/
1854 for(col=1;col<P->Dimension+1;col++) {
1855 value_multiply(tmp,P->Constraint[row][col],list_args[col-1]);
1856 value_addto(v,v,tmp);
1857 }
1858 if (value_notzero_p(P->Constraint[row][0])) {
1859
1860 /*if v is not >=0 then this constraint is not respected */
1861 if (value_neg_p(v)) {
1862 next:
1863 value_clear(v);
1864 value_clear(tmp);
1865 return P->next ? in_domain(P->next, list_args) : 0;
1866 }
1867 }
1868 else {
1869
1870 /*if v is not = 0 then this constraint is not respected */
1871 if (value_notzero_p(v))
1872 goto next;
1873 }
1874 }
1875
1876 /*if not return before this point => all
1877 the constraints are respected */
1878 value_clear(v);
1879 value_clear(tmp);
1880 return 1;
1881 } /* in_domain */
1882
1883 /****************************************************/
1884 /* function compute enode */
1885 /* compute the value of enode p with parameters */
1886 /* list "list_args */
1887 /* compute the polynomial or the periodic */
1888 /****************************************************/
1889
1890 static double compute_enode(enode *p, Value *list_args) {
1891
1892 int i;
1893 Value m, param;
1894 double res=0.0;
1895
1896 if (!p)
1897 return(0.);
1898
1899 value_init(m);
1900 value_init(param);
1901
1902 if (p->type == polynomial) {
1903 if (p->size > 1)
1904 value_assign(param,list_args[p->pos-1]);
1905
1906 /* Compute the polynomial using Horner's rule */
1907 for (i=p->size-1;i>0;i--) {
1908 res +=compute_evalue(&p->arr[i],list_args);
1909 res *=VALUE_TO_DOUBLE(param);
1910 }
1911 res +=compute_evalue(&p->arr[0],list_args);
1912 }
1913 else if (p->type == fractional) {
1914 double d = compute_evalue(&p->arr[0], list_args);
1915 d -= floor(d+1e-10);
1916
1917 /* Compute the polynomial using Horner's rule */
1918 for (i=p->size-1;i>1;i--) {
1919 res +=compute_evalue(&p->arr[i],list_args);
1920 res *=d;
1921 }
1922 res +=compute_evalue(&p->arr[1],list_args);
1923 }
1924 else if (p->type == periodic) {
1925 value_assign(param,list_args[p->pos-1]);
1926
1927 /* Choose the right element of the periodic */
1928 value_absolute(m,param);
1929 value_set_si(param,p->size);
1930 value_modulus(m,m,param);
1931 res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
1932 }
1933 else if (p->type == relation) {
1934 if (fabs(compute_evalue(&p->arr[0], list_args)) < 1e-10)
1935 res = compute_evalue(&p->arr[1], list_args);
1936 else if (p->size > 2)
1937 res = compute_evalue(&p->arr[2], list_args);
1938 }
1939 else if (p->type == partition) {
1940 for (i = 0; i < p->size/2; ++i)
1941 if (in_domain(EVALUE_DOMAIN(p->arr[2*i]), list_args)) {
1942 res = compute_evalue(&p->arr[2*i+1], list_args);
1943 break;
1944 }
1945 }
1946 value_clear(m);
1947 value_clear(param);
1948 return res;
1949 } /* compute_enode */
1950
1951 /*************************************************/
1952 /* return the value of Ehrhart Polynomial */
1953 /* It returns a double, because since it is */
1954 /* a recursive function, some intermediate value */
1955 /* might not be integral */
1956 /*************************************************/
1957
1958 double compute_evalue(evalue *e,Value *list_args) {
1959
1960 double res;
1961
1962 if (value_notzero_p(e->d)) {
1963 if (value_notone_p(e->d))
1964 res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
1965 else
1966 res = VALUE_TO_DOUBLE(e->x.n);
1967 }
1968 else
1969 res = compute_enode(e->x.p,list_args);
1970 return res;
1971 } /* compute_evalue */
1972
1973
1974 /****************************************************/
1975 /* function compute_poly : */
1976 /* Check for the good validity domain */
1977 /* return the number of point in the Polyhedron */
1978 /* in allocated memory */
1979 /* Using the Ehrhart pseudo-polynomial */
1980 /****************************************************/
1981 Value *compute_poly(Enumeration *en,Value *list_args) {
1982
1983 Value *tmp;
1984 /* double d; int i; */
1985
1986 tmp = (Value *) malloc (sizeof(Value));
1987 assert(tmp != NULL);
1988 value_init(*tmp);
1989 value_set_si(*tmp,0);
1990
1991 if(!en)
1992 return(tmp); /* no ehrhart polynomial */
1993 if(en->ValidityDomain) {
1994 if(!en->ValidityDomain->Dimension) { /* no parameters */
1995 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
1996 return(tmp);
1997 }
1998 }
1999 else
2000 return(tmp); /* no Validity Domain */
2001 while(en) {
2002 if(in_domain(en->ValidityDomain,list_args)) {
2003
2004 #ifdef EVAL_EHRHART_DEBUG
2005 Print_Domain(stdout,en->ValidityDomain);
2006 print_evalue(stdout,&en->EP);
2007 #endif
2008
2009 /* d = compute_evalue(&en->EP,list_args);
2010 i = d;
2011 printf("(double)%lf = %d\n", d, i ); */
2012 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2013 return(tmp);
2014 }
2015 else
2016 en=en->next;
2017 }
2018 value_set_si(*tmp,0);
2019 return(tmp); /* no compatible domain with the arguments */
2020 } /* compute_poly */
2021
2022 size_t value_size(Value v) {
2023 return (v[0]._mp_size > 0 ? v[0]._mp_size : -v[0]._mp_size)
2024 * sizeof(v[0]._mp_d[0]);
2025 }
2026
2027 size_t domain_size(Polyhedron *D)
2028 {
2029 int i, j;
2030 size_t s = sizeof(*D);
2031
2032 for (i = 0; i < D->NbConstraints; ++i)
2033 for (j = 0; j < D->Dimension+2; ++j)
2034 s += value_size(D->Constraint[i][j]);
2035
2036 for (i = 0; i < D->NbRays; ++i)
2037 for (j = 0; j < D->Dimension+2; ++j)
2038 s += value_size(D->Ray[i][j]);
2039
2040 return D->next ? s+domain_size(D->next) : s;
2041 }
2042
2043 size_t enode_size(enode *p) {
2044 size_t s = sizeof(*p) - sizeof(p->arr[0]);
2045 int i;
2046
2047 if (p->type == partition)
2048 for (i = 0; i < p->size/2; ++i) {
2049 s += domain_size(EVALUE_DOMAIN(p->arr[2*i]));
2050 s += evalue_size(&p->arr[2*i+1]);
2051 }
2052 else
2053 for (i = 0; i < p->size; ++i) {
2054 s += evalue_size(&p->arr[i]);
2055 }
2056 return s;
2057 }
2058
2059 size_t evalue_size(evalue *e)
2060 {
2061 size_t s = sizeof(*e);
2062 s += value_size(e->d);
2063 if (value_notzero_p(e->d))
2064 s += value_size(e->x.n);
2065 else
2066 s += enode_size(e->x.p);
2067 return s;
2068 }
2069
2070 static evalue *find_second(evalue *base, evalue *cst, evalue *e, Value m)
2071 {
2072 evalue *found = NULL;
2073 evalue offset;
2074 evalue copy;
2075 int i;
2076
2077 if (value_pos_p(e->d) || e->x.p->type != fractional)
2078 return NULL;
2079
2080 value_init(offset.d);
2081 value_init(offset.x.n);
2082 poly_denom(&e->x.p->arr[0], &offset.d);
2083 Lcm3(m, offset.d, &offset.d);
2084 value_set_si(offset.x.n, 1);
2085
2086 value_init(copy.d);
2087 evalue_copy(&copy, cst);
2088
2089 eadd(&offset, cst);
2090 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2091
2092 if (eequal(base, &e->x.p->arr[0]))
2093 found = &e->x.p->arr[0];
2094 else {
2095 value_set_si(offset.x.n, -2);
2096
2097 eadd(&offset, cst);
2098 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2099
2100 if (eequal(base, &e->x.p->arr[0]))
2101 found = base;
2102 }
2103 free_evalue_refs(cst);
2104 free_evalue_refs(&offset);
2105 *cst = copy;
2106
2107 for (i = 1; !found && i < e->x.p->size; ++i)
2108 found = find_second(base, cst, &e->x.p->arr[i], m);
2109
2110 return found;
2111 }
2112
2113 static evalue *find_relation_pair(evalue *e)
2114 {
2115 int i;
2116 evalue *found = NULL;
2117
2118 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2119 return NULL;
2120
2121 if (e->x.p->type == fractional) {
2122 Value m;
2123 evalue *cst;
2124
2125 value_init(m);
2126 poly_denom(&e->x.p->arr[0], &m);
2127
2128 for (cst = &e->x.p->arr[0]; value_zero_p(cst->d);
2129 cst = &cst->x.p->arr[0])
2130 ;
2131
2132 for (i = 1; !found && i < e->x.p->size; ++i)
2133 found = find_second(&e->x.p->arr[0], cst, &e->x.p->arr[i], m);
2134
2135 value_clear(m);
2136 }
2137
2138 i = e->x.p->type == relation;
2139 for (; !found && i < e->x.p->size; ++i)
2140 found = find_relation_pair(&e->x.p->arr[i]);
2141
2142 return found;
2143 }
2144
2145 void evalue_mod2relation(evalue *e) {
2146 evalue *d;
2147
2148 if (value_zero_p(e->d) && e->x.p->type == partition) {
2149 int i;
2150
2151 for (i = 0; i < e->x.p->size/2; ++i)
2152 evalue_mod2relation(&e->x.p->arr[2*i+1]);
2153
2154 return;
2155 }
2156
2157 while ((d = find_relation_pair(e)) != NULL) {
2158 evalue split;
2159 evalue *ev;
2160
2161 value_init(split.d);
2162 value_set_si(split.d, 0);
2163 split.x.p = new_enode(relation, 3, 0);
2164 evalue_set_si(&split.x.p->arr[1], 1, 1);
2165 evalue_set_si(&split.x.p->arr[2], 1, 1);
2166
2167 ev = &split.x.p->arr[0];
2168 value_set_si(ev->d, 0);
2169 ev->x.p = new_enode(fractional, 3, -1);
2170 evalue_set_si(&ev->x.p->arr[1], 0, 1);
2171 evalue_set_si(&ev->x.p->arr[2], 1, 1);
2172 evalue_copy(&ev->x.p->arr[0], d);
2173
2174 emul(&split, e);
2175
2176 reduce_evalue(e);
2177
2178 free_evalue_refs(&split);
2179 }
2180 }
2181
2182 static int evalue_comp(const void * a, const void * b)
2183 {
2184 const evalue *e1 = *(const evalue **)a;
2185 const evalue *e2 = *(const evalue **)b;
2186 return ecmp(e1, e2);
2187 }
2188
2189 void evalue_combine(evalue *e)
2190 {
2191 evalue **evs;
2192 int i, k;
2193 enode *p;
2194 evalue tmp;
2195
2196 if (value_notzero_p(e->d) || e->x.p->type != partition)
2197 return;
2198
2199 NALLOC(evs, e->x.p->size/2);
2200 for (i = 0; i < e->x.p->size/2; ++i)
2201 evs[i] = &e->x.p->arr[2*i+1];
2202 qsort(evs, e->x.p->size/2, sizeof(evs[0]), evalue_comp);
2203 p = new_enode(partition, e->x.p->size, -1);
2204 for (i = 0, k = 0; i < e->x.p->size/2; ++i) {
2205 if (k == 0 || ecmp(&p->arr[2*k-1], evs[i]) != 0) {
2206 p->arr[2*k] = *(evs[i]-1);
2207 p->arr[2*k+1] = *(evs[i]);
2208 ++k;
2209 } else {
2210 Polyhedron *D = EVALUE_DOMAIN(*(evs[i]-1));
2211 Polyhedron *L = D;
2212
2213 while (L->next)
2214 L = L->next;
2215 L->next = EVALUE_DOMAIN(p->arr[2*k-2]);
2216 EVALUE_SET_DOMAIN(p->arr[2*k-2], D);
2217 free_evalue_refs(evs[i]);
2218 }
2219 }
2220 free(e->x.p);
2221 p->size = 2*k;
2222 e->x.p = p;
2223
2224 for (i = 0; i < e->x.p->size/2; ++i) {
2225 Polyhedron *H;
2226 if (value_notzero_p(e->x.p->arr[2*i+1].d))
2227 continue;
2228 H = DomainConvex(EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
2229 if (H == NULL)
2230 continue;
2231 for (k = 0; k < e->x.p->size/2; ++k) {
2232 Polyhedron *D, *N, **P;
2233 if (i == k)
2234 continue;
2235 P = &EVALUE_DOMAIN(e->x.p->arr[2*k]);
2236 D = *P;
2237 if (D == NULL)
2238 continue;
2239 for (; D; D = N) {
2240 *P = D;
2241 N = D->next;
2242 if (D->NbEq <= H->NbEq) {
2243 P = &D->next;
2244 continue;
2245 }
2246
2247 value_init(tmp.d);
2248 tmp.x.p = new_enode(partition, 2, -1);
2249 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Copy(D));
2250 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
2251 reduce_evalue(&tmp);
2252 if (value_notzero_p(tmp.d) ||
2253 ecmp(&tmp.x.p->arr[1], &e->x.p->arr[2*k+1]) != 0)
2254 P = &D->next;
2255 else {
2256 D->next = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2257 EVALUE_DOMAIN(e->x.p->arr[2*i]) = D;
2258 *P = NULL;
2259 }
2260 free_evalue_refs(&tmp);
2261 }
2262 }
2263 }
2264
2265 for (i = 0; i < e->x.p->size/2; ++i) {
2266 Polyhedron *H, *E;
2267 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2268 if (!D) {
2269 free_evalue_refs(&e->x.p->arr[2*i+1]);
2270 e->x.p->size -= 2;
2271 if (2*i < e->x.p->size) {
2272 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2273 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2274 }
2275 --i;
2276 continue;
2277 }
2278 if (!D->next)
2279 continue;
2280 H = DomainConvex(D, 0);
2281 E = DomainDifference(H, D, 0);
2282 Domain_Free(D);
2283 D = DomainDifference(H, E, 0);
2284 Domain_Free(H);
2285 Domain_Free(E);
2286 EVALUE_SET_DOMAIN(p->arr[2*i], D);
2287 }
2288 }
2289
2290 /* May change coefficients to become non-standard
2291 * => reduce p afterwards to correct
2292 */
2293 static Polyhedron *polynomial_projection(enode *p, Polyhedron *D, Value *d,
2294 Matrix **R)
2295 {
2296 Polyhedron *I, *H;
2297 evalue *pp;
2298 unsigned dim = D->Dimension;
2299 Matrix *T = Matrix_Alloc(2, dim+1);
2300 Value twice;
2301 value_init(twice);
2302 assert(T);
2303
2304 assert(p->type == fractional);
2305 pp = &p->arr[0];
2306 value_set_si(T->p[1][dim], 1);
2307 poly_denom(pp, d);
2308 while (value_zero_p(pp->d)) {
2309 assert(pp->x.p->type == polynomial);
2310 assert(pp->x.p->size == 2);
2311 assert(value_notzero_p(pp->x.p->arr[1].d));
2312 value_division(T->p[0][pp->x.p->pos-1], *d, pp->x.p->arr[1].d);
2313 mpz_mul_ui(twice, pp->x.p->arr[1].x.n, 2);
2314 if (value_gt(twice, pp->x.p->arr[1].d))
2315 value_substract(pp->x.p->arr[1].x.n,
2316 pp->x.p->arr[1].x.n, pp->x.p->arr[1].d);
2317 value_multiply(T->p[0][pp->x.p->pos-1],
2318 T->p[0][pp->x.p->pos-1], pp->x.p->arr[1].x.n);
2319 pp = &pp->x.p->arr[0];
2320 }
2321 value_division(T->p[0][dim], *d, pp->d);
2322 value_multiply(T->p[0][dim], T->p[0][dim], pp->x.n);
2323 I = DomainImage(D, T, 0);
2324 H = DomainConvex(I, 0);
2325 Domain_Free(I);
2326 *R = T;
2327
2328 value_clear(twice);
2329
2330 return H;
2331 }
2332
2333 static int reduce_in_domain(evalue *e, Polyhedron *D)
2334 {
2335 int i;
2336 enode *p;
2337 Value d, min, max;
2338 int r = 0;
2339 Polyhedron *I;
2340 Matrix *T;
2341
2342 if (value_notzero_p(e->d))
2343 return r;
2344
2345 p = e->x.p;
2346
2347 if (p->type == relation) {
2348 int equal;
2349 value_init(d);
2350 value_init(min);
2351 value_init(max);
2352
2353 I = polynomial_projection(p->arr[0].x.p, D, &d, &T);
2354 line_minmax(I, &min, &max); /* frees I */
2355 equal = value_eq(min, max);
2356 mpz_cdiv_q(min, min, d);
2357 mpz_fdiv_q(max, max, d);
2358
2359 if (value_gt(min, max)) {
2360 /* Never zero */
2361 if (p->size == 3) {
2362 value_clear(e->d);
2363 *e = p->arr[2];
2364 } else {
2365 evalue_set_si(e, 0, 1);
2366 r = 1;
2367 }
2368 free_evalue_refs(&(p->arr[1]));
2369 free_evalue_refs(&(p->arr[0]));
2370 free(p);
2371 value_clear(d);
2372 value_clear(min);
2373 value_clear(max);
2374 Matrix_Free(T);
2375 return r ? r : reduce_in_domain(e, D);
2376 } else if (equal) {
2377 /* Always zero */
2378 if (p->size == 3)
2379 free_evalue_refs(&(p->arr[2]));
2380 *e = p->arr[1];
2381 free_evalue_refs(&(p->arr[0]));
2382 free(p);
2383 value_clear(d);
2384 value_clear(min);
2385 value_clear(max);
2386 Matrix_Free(T);
2387 return reduce_in_domain(e, D);
2388 } else if (value_eq(min, max)) {
2389 /* zero for a single value */
2390 Polyhedron *E;
2391 Matrix *M = Matrix_Alloc(1, D->Dimension+2);
2392 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
2393 value_multiply(min, min, d);
2394 value_substract(M->p[0][D->Dimension+1],
2395 M->p[0][D->Dimension+1], min);
2396 E = DomainAddConstraints(D, M, 0);
2397 value_clear(d);
2398 value_clear(min);
2399 value_clear(max);
2400 Matrix_Free(T);
2401 Matrix_Free(M);
2402 r = reduce_in_domain(&p->arr[1], E);
2403 if (p->size == 3)
2404 r |= reduce_in_domain(&p->arr[2], D);
2405 Domain_Free(E);
2406 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2407 return r;
2408 }
2409
2410 value_clear(d);
2411 value_clear(min);
2412 value_clear(max);
2413 Matrix_Free(T);
2414 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2415 }
2416
2417 i = p->type == relation ? 1 :
2418 p->type == fractional ? 1 : 0;
2419 for (; i<p->size; i++)
2420 r |= reduce_in_domain(&p->arr[i], D);
2421
2422 if (p->type != fractional) {
2423 if (r && p->type == polynomial) {
2424 evalue f;
2425 value_init(f.d);
2426 value_set_si(f.d, 0);
2427 f.x.p = new_enode(polynomial, 2, p->pos);
2428 evalue_set_si(&f.x.p->arr[0], 0, 1);
2429 evalue_set_si(&f.x.p->arr[1], 1, 1);
2430 reorder_terms(p, &f);
2431 *e = p->arr[0];
2432 free(p);
2433 }
2434 return r;
2435 }
2436
2437 value_init(d);
2438 value_init(min);
2439 value_init(max);
2440 I = polynomial_projection(p, D, &d, &T);
2441 Matrix_Free(T);
2442 line_minmax(I, &min, &max); /* frees I */
2443 mpz_fdiv_q(min, min, d);
2444 mpz_fdiv_q(max, max, d);
2445
2446 if (value_eq(min, max)) {
2447 evalue inc;
2448 value_init(inc.d);
2449 value_init(inc.x.n);
2450 value_set_si(inc.d, 1);
2451 value_oppose(inc.x.n, min);
2452 eadd(&inc, &p->arr[0]);
2453 reorder_terms(p, &p->arr[0]); /* frees arr[0] */
2454 value_clear(e->d);
2455 *e = p->arr[1];
2456 free(p);
2457 free_evalue_refs(&inc);
2458 r = 1;
2459 } else {
2460 _reduce_evalue(&p->arr[0], 0, 1);
2461 if (r == 1) {
2462 evalue f;
2463 value_init(f.d);
2464 value_set_si(f.d, 0);
2465 f.x.p = new_enode(fractional, 3, -1);
2466 value_clear(f.x.p->arr[0].d);
2467 f.x.p->arr[0] = p->arr[0];
2468 evalue_set_si(&f.x.p->arr[1], 0, 1);
2469 evalue_set_si(&f.x.p->arr[2], 1, 1);
2470 reorder_terms(p, &f);
2471 *e = p->arr[1];
2472 free(p);
2473 }
2474 }
2475
2476 value_clear(d);
2477 value_clear(min);
2478 value_clear(max);
2479
2480 return r;
2481 }
2482
2483 void evalue_range_reduction(evalue *e)
2484 {
2485 int i;
2486 if (value_notzero_p(e->d) || e->x.p->type != partition)
2487 return;
2488
2489 for (i = 0; i < e->x.p->size/2; ++i)
2490 if (reduce_in_domain(&e->x.p->arr[2*i+1],
2491 EVALUE_DOMAIN(e->x.p->arr[2*i]))) {
2492 reduce_evalue(&e->x.p->arr[2*i+1]);
2493
2494 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2495 free_evalue_refs(&e->x.p->arr[2*i+1]);
2496 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
2497 value_clear(e->x.p->arr[2*i].d);
2498 e->x.p->size -= 2;
2499 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2500 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2501 --i;
2502 }
2503 }
2504 }
lattice point enumeration library