?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
adapt to deal with existential variables
[barvinok.git] / ev_operations.c
blob5512dcfb311f00f120b6e61fa71efab5bbcca190
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 value_clear(s.fixed[j].m);
692 free_evalue_refs(&s.fixed[j].s);
693 }
694 }
695 }
696 if (e->x.p->size == 0) {
697 free(e->x.p);
698 evalue_set_si(e, 0, 1);
699 }
700 } else
701 _reduce_evalue(e, &s, 0);
702 if (s.max != 0)
703 free(s.fixed);
704 }
705
706 void print_evalue(FILE *DST,evalue *e,char **pname) {
707
708 if(value_notzero_p(e->d)) {
709 if(value_notone_p(e->d)) {
710 value_print(DST,VALUE_FMT,e->x.n);
711 fprintf(DST,"/");
712 value_print(DST,VALUE_FMT,e->d);
713 }
714 else {
715 value_print(DST,VALUE_FMT,e->x.n);
716 }
717 }
718 else
719 print_enode(DST,e->x.p,pname);
720 return;
721 } /* print_evalue */
722
723 void print_enode(FILE *DST,enode *p,char **pname) {
724
725 int i;
726
727 if (!p) {
728 fprintf(DST, "NULL");
729 return;
730 }
731 switch (p->type) {
732 case evector:
733 fprintf(DST, "{ ");
734 for (i=0; i<p->size; i++) {
735 print_evalue(DST, &p->arr[i], pname);
736 if (i!=(p->size-1))
737 fprintf(DST, ", ");
738 }
739 fprintf(DST, " }\n");
740 break;
741 case polynomial:
742 fprintf(DST, "( ");
743 for (i=p->size-1; i>=0; i--) {
744 print_evalue(DST, &p->arr[i], pname);
745 if (i==1) fprintf(DST, " * %s + ", pname[p->pos-1]);
746 else if (i>1)
747 fprintf(DST, " * %s^%d + ", pname[p->pos-1], i);
748 }
749 fprintf(DST, " )\n");
750 break;
751 case periodic:
752 fprintf(DST, "[ ");
753 for (i=0; i<p->size; i++) {
754 print_evalue(DST, &p->arr[i], pname);
755 if (i!=(p->size-1)) fprintf(DST, ", ");
756 }
757 fprintf(DST," ]_%s", pname[p->pos-1]);
758 break;
759 case fractional:
760 fprintf(DST, "( ");
761 for (i=p->size-1; i>=1; i--) {
762 print_evalue(DST, &p->arr[i], pname);
763 if (i >= 2) {
764 fprintf(DST, " * ");
765 fprintf(DST, "{");
766 print_evalue(DST, &p->arr[0], pname);
767 fprintf(DST, "}");
768 if (i>2)
769 fprintf(DST, "^%d + ", i-1);
770 else
771 fprintf(DST, " + ");
772 }
773 }
774 fprintf(DST, " )\n");
775 break;
776 case relation:
777 fprintf(DST, "[ ");
778 print_evalue(DST, &p->arr[0], pname);
779 fprintf(DST, "= 0 ] * \n");
780 print_evalue(DST, &p->arr[1], pname);
781 if (p->size > 2) {
782 fprintf(DST, " +\n [ ");
783 print_evalue(DST, &p->arr[0], pname);
784 fprintf(DST, "!= 0 ] * \n");
785 print_evalue(DST, &p->arr[2], pname);
786 }
787 break;
788 case partition:
789 for (i=0; i<p->size/2; i++) {
790 Print_Domain(DST, EVALUE_DOMAIN(p->arr[2*i]), pname);
791 print_evalue(DST, &p->arr[2*i+1], pname);
792 }
793 break;
794 default:
795 assert(0);
796 }
797 return;
798 } /* print_enode */
799
800 static void eadd_rev(evalue *e1, evalue *res)
801 {
802 evalue ev;
803 value_init(ev.d);
804 evalue_copy(&ev, e1);
805 eadd(res, &ev);
806 free_evalue_refs(res);
807 *res = ev;
808 }
809
810 static void eadd_rev_cst (evalue *e1, evalue *res)
811 {
812 evalue ev;
813 value_init(ev.d);
814 evalue_copy(&ev, e1);
815 eadd(res, &ev.x.p->arr[ev.x.p->type==fractional]);
816 free_evalue_refs(res);
817 *res = ev;
818 }
819
820 struct section { Polyhedron * D; evalue E; };
821
822 void eadd_partitions (evalue *e1,evalue *res)
823 {
824 int n, i, j;
825 Polyhedron *d, *fd;
826 struct section *s;
827 s = (struct section *)
828 malloc((e1->x.p->size/2+1) * (res->x.p->size/2+1) *
829 sizeof(struct section));
830 assert(s);
831
832 n = 0;
833 for (j = 0; j < e1->x.p->size/2; ++j) {
834 assert(res->x.p->size >= 2);
835 fd = DomainDifference(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
836 EVALUE_DOMAIN(res->x.p->arr[0]), 0);
837 if (!emptyQ(fd))
838 for (i = 1; i < res->x.p->size/2; ++i) {
839 Polyhedron *t = fd;
840 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
841 Domain_Free(t);
842 if (emptyQ(fd))
843 break;
844 }
845 if (emptyQ(fd)) {
846 Domain_Free(fd);
847 continue;
848 }
849 value_init(s[n].E.d);
850 evalue_copy(&s[n].E, &e1->x.p->arr[2*j+1]);
851 s[n].D = fd;
852 ++n;
853 }
854 for (i = 0; i < res->x.p->size/2; ++i) {
855 fd = EVALUE_DOMAIN(res->x.p->arr[2*i]);
856 for (j = 0; j < e1->x.p->size/2; ++j) {
857 Polyhedron *t;
858 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
859 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
860 if (emptyQ(d)) {
861 Domain_Free(d);
862 continue;
863 }
864 t = fd;
865 fd = DomainDifference(fd, EVALUE_DOMAIN(e1->x.p->arr[2*j]), 0);
866 if (t != EVALUE_DOMAIN(res->x.p->arr[2*i]))
867 Domain_Free(t);
868 value_init(s[n].E.d);
869 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
870 eadd(&e1->x.p->arr[2*j+1], &s[n].E);
871 s[n].D = d;
872 ++n;
873 }
874 if (!emptyQ(fd)) {
875 s[n].E = res->x.p->arr[2*i+1];
876 s[n].D = fd;
877 ++n;
878 } else {
879 free_evalue_refs(&res->x.p->arr[2*i+1]);
880 Domain_Free(fd);
881 }
882 if (fd != EVALUE_DOMAIN(res->x.p->arr[2*i]))
883 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
884 value_clear(res->x.p->arr[2*i].d);
885 }
886
887 free(res->x.p);
888 res->x.p = new_enode(partition, 2*n, -1);
889 for (j = 0; j < n; ++j) {
890 s[j].D = DomainConstraintSimplify(s[j].D, 0);
891 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
892 value_clear(res->x.p->arr[2*j+1].d);
893 res->x.p->arr[2*j+1] = s[j].E;
894 }
895
896 free(s);
897 }
898
899 static void explicit_complement(evalue *res)
900 {
901 enode *rel = new_enode(relation, 3, 0);
902 assert(rel);
903 value_clear(rel->arr[0].d);
904 rel->arr[0] = res->x.p->arr[0];
905 value_clear(rel->arr[1].d);
906 rel->arr[1] = res->x.p->arr[1];
907 value_set_si(rel->arr[2].d, 1);
908 value_init(rel->arr[2].x.n);
909 value_set_si(rel->arr[2].x.n, 0);
910 free(res->x.p);
911 res->x.p = rel;
912 }
913
914 void eadd(evalue *e1,evalue *res) {
915
916 int i;
917 if (value_notzero_p(e1->d) && value_notzero_p(res->d)) {
918 /* Add two rational numbers */
919 Value g,m1,m2;
920 value_init(g);
921 value_init(m1);
922 value_init(m2);
923
924 value_multiply(m1,e1->x.n,res->d);
925 value_multiply(m2,res->x.n,e1->d);
926 value_addto(res->x.n,m1,m2);
927 value_multiply(res->d,e1->d,res->d);
928 Gcd(res->x.n,res->d,&g);
929 if (value_notone_p(g)) {
930 value_division(res->d,res->d,g);
931 value_division(res->x.n,res->x.n,g);
932 }
933 value_clear(g); value_clear(m1); value_clear(m2);
934 return ;
935 }
936 else if (value_notzero_p(e1->d) && value_zero_p(res->d)) {
937 switch (res->x.p->type) {
938 case polynomial:
939 /* Add the constant to the constant term of a polynomial*/
940 eadd(e1, &res->x.p->arr[0]);
941 return ;
942 case periodic:
943 /* Add the constant to all elements of a periodic number */
944 for (i=0; i<res->x.p->size; i++) {
945 eadd(e1, &res->x.p->arr[i]);
946 }
947 return ;
948 case evector:
949 fprintf(stderr, "eadd: cannot add const with vector\n");
950 return;
951 case fractional:
952 eadd(e1, &res->x.p->arr[1]);
953 return ;
954 case partition:
955 assert(EVALUE_IS_ZERO(*e1));
956 break; /* Do nothing */
957 case relation:
958 /* Create (zero) complement if needed */
959 if (res->x.p->size < 3 && !EVALUE_IS_ZERO(*e1))
960 explicit_complement(res);
961 for (i = 1; i < res->x.p->size; ++i)
962 eadd(e1, &res->x.p->arr[i]);
963 break;
964 default:
965 assert(0);
966 }
967 }
968 /* add polynomial or periodic to constant
969 * you have to exchange e1 and res, before doing addition */
970
971 else if (value_zero_p(e1->d) && value_notzero_p(res->d)) {
972 eadd_rev(e1, res);
973 return;
974 }
975 else { // ((e1->d==0) && (res->d==0))
976 assert(!((e1->x.p->type == partition) ^
977 (res->x.p->type == partition)));
978 if (e1->x.p->type == partition) {
979 eadd_partitions(e1, res);
980 return;
981 }
982 if (e1->x.p->type == relation &&
983 (res->x.p->type != relation ||
984 mod_term_smaller(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
985 eadd_rev(e1, res);
986 return;
987 }
988 if (res->x.p->type == relation) {
989 if (e1->x.p->type == relation &&
990 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
991 if (res->x.p->size < 3 && e1->x.p->size == 3)
992 explicit_complement(res);
993 if (e1->x.p->size < 3 && res->x.p->size == 3)
994 explicit_complement(e1);
995 for (i = 1; i < res->x.p->size; ++i)
996 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
997 return;
998 }
999 if (res->x.p->size < 3)
1000 explicit_complement(res);
1001 for (i = 1; i < res->x.p->size; ++i)
1002 eadd(e1, &res->x.p->arr[i]);
1003 return;
1004 }
1005 if ((e1->x.p->type != res->x.p->type) ) {
1006 /* adding to evalues of different type. two cases are possible
1007 * res is periodic and e1 is polynomial, you have to exchange
1008 * e1 and res then to add e1 to the constant term of res */
1009 if (e1->x.p->type == polynomial) {
1010 eadd_rev_cst(e1, res);
1011 }
1012 else if (res->x.p->type == polynomial) {
1013 /* res is polynomial and e1 is periodic,
1014 add e1 to the constant term of res */
1015
1016 eadd(e1,&res->x.p->arr[0]);
1017 } else
1018 assert(0);
1019
1020 return;
1021 }
1022 else if (e1->x.p->pos != res->x.p->pos ||
1023 (res->x.p->type == fractional &&
1024 !eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1025 /* adding evalues of different position (i.e function of different unknowns
1026 * to case are possible */
1027
1028 switch (res->x.p->type) {
1029 case fractional:
1030 if(mod_term_smaller(res, e1))
1031 eadd(e1,&res->x.p->arr[1]);
1032 else
1033 eadd_rev_cst(e1, res);
1034 return;
1035 case polynomial: // res and e1 are polynomials
1036 // add e1 to the constant term of res
1037
1038 if(res->x.p->pos < e1->x.p->pos)
1039 eadd(e1,&res->x.p->arr[0]);
1040 else
1041 eadd_rev_cst(e1, res);
1042 // value_clear(g); value_clear(m1); value_clear(m2);
1043 return;
1044 case periodic: // res and e1 are pointers to periodic numbers
1045 //add e1 to all elements of res
1046
1047 if(res->x.p->pos < e1->x.p->pos)
1048 for (i=0;i<res->x.p->size;i++) {
1049 eadd(e1,&res->x.p->arr[i]);
1050 }
1051 else
1052 eadd_rev(e1, res);
1053 return;
1054 }
1055 }
1056
1057
1058 //same type , same pos and same size
1059 if (e1->x.p->size == res->x.p->size) {
1060 // add any element in e1 to the corresponding element in res
1061 if (res->x.p->type == fractional)
1062 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
1063 i = res->x.p->type == fractional ? 1 : 0;
1064 for (; i<res->x.p->size; i++) {
1065 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1066 }
1067 return ;
1068 }
1069
1070 /* Sizes are different */
1071 switch(res->x.p->type) {
1072 case polynomial:
1073 case fractional:
1074 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1075 /* new enode and add res to that new node. If you do not do */
1076 /* that, you lose the the upper weight part of e1 ! */
1077
1078 if(e1->x.p->size > res->x.p->size)
1079 eadd_rev(e1, res);
1080 else {
1081
1082 if (res->x.p->type == fractional)
1083 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
1084 i = res->x.p->type == fractional ? 1 : 0;
1085 for (; i<e1->x.p->size ; i++) {
1086 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1087 }
1088
1089 return ;
1090 }
1091 break;
1092
1093 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1094 case periodic:
1095 {
1096 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1097 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1098 to add periodicaly elements of e1 to elements of ne, and finaly to
1099 return ne. */
1100 int x,y,p;
1101 Value ex, ey ,ep;
1102 evalue *ne;
1103 value_init(ex); value_init(ey);value_init(ep);
1104 x=e1->x.p->size;
1105 y= res->x.p->size;
1106 value_set_si(ex,e1->x.p->size);
1107 value_set_si(ey,res->x.p->size);
1108 value_assign (ep,*Lcm(ex,ey));
1109 p=(int)mpz_get_si(ep);
1110 ne= (evalue *) malloc (sizeof(evalue));
1111 value_init(ne->d);
1112 value_set_si( ne->d,0);
1113
1114 ne->x.p=new_enode(res->x.p->type,p, res->x.p->pos);
1115 for(i=0;i<p;i++) {
1116 evalue_copy(&ne->x.p->arr[i], &res->x.p->arr[i%y]);
1117 }
1118 for(i=0;i<p;i++) {
1119 eadd(&e1->x.p->arr[i%x], &ne->x.p->arr[i]);
1120 }
1121
1122 value_assign(res->d, ne->d);
1123 res->x.p=ne->x.p;
1124
1125 return ;
1126 }
1127 case evector:
1128 fprintf(stderr, "eadd: ?cannot add vectors of different length\n");
1129 return ;
1130 }
1131 }
1132 return ;
1133 }/* eadd */
1134
1135 static void emul_rev (evalue *e1, evalue *res)
1136 {
1137 evalue ev;
1138 value_init(ev.d);
1139 evalue_copy(&ev, e1);
1140 emul(res, &ev);
1141 free_evalue_refs(res);
1142 *res = ev;
1143 }
1144
1145 static void emul_poly (evalue *e1, evalue *res)
1146 {
1147 int i, j, o = res->x.p->type == fractional;
1148 evalue tmp;
1149 int size=(e1->x.p->size + res->x.p->size - o - 1);
1150 value_init(tmp.d);
1151 value_set_si(tmp.d,0);
1152 tmp.x.p=new_enode(res->x.p->type, size, res->x.p->pos);
1153 if (o)
1154 evalue_copy(&tmp.x.p->arr[0], &e1->x.p->arr[0]);
1155 for (i=o; i < e1->x.p->size; i++) {
1156 evalue_copy(&tmp.x.p->arr[i], &e1->x.p->arr[i]);
1157 emul(&res->x.p->arr[o], &tmp.x.p->arr[i]);
1158 }
1159 for (; i<size; i++)
1160 evalue_set_si(&tmp.x.p->arr[i], 0, 1);
1161 for (i=o+1; i<res->x.p->size; i++)
1162 for (j=o; j<e1->x.p->size; j++) {
1163 evalue ev;
1164 value_init(ev.d);
1165 evalue_copy(&ev, &e1->x.p->arr[j]);
1166 emul(&res->x.p->arr[i], &ev);
1167 eadd(&ev, &tmp.x.p->arr[i+j-o]);
1168 free_evalue_refs(&ev);
1169 }
1170 free_evalue_refs(res);
1171 *res = tmp;
1172 }
1173
1174 void emul_partitions (evalue *e1,evalue *res)
1175 {
1176 int n, i, j, k;
1177 Polyhedron *d;
1178 struct section *s;
1179 s = (struct section *)
1180 malloc((e1->x.p->size/2) * (res->x.p->size/2) *
1181 sizeof(struct section));
1182 assert(s);
1183
1184 n = 0;
1185 for (i = 0; i < res->x.p->size/2; ++i) {
1186 for (j = 0; j < e1->x.p->size/2; ++j) {
1187 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1188 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1189 if (emptyQ(d)) {
1190 Domain_Free(d);
1191 continue;
1192 }
1193
1194 /* This code is only needed because the partitions
1195 are not true partitions.
1196 */
1197 for (k = 0; k < n; ++k) {
1198 if (DomainIncludes(s[k].D, d))
1199 break;
1200 if (DomainIncludes(d, s[k].D)) {
1201 Domain_Free(s[k].D);
1202 free_evalue_refs(&s[k].E);
1203 if (n > k)
1204 s[k] = s[--n];
1205 --k;
1206 }
1207 }
1208 if (k < n) {
1209 Domain_Free(d);
1210 continue;
1211 }
1212
1213 value_init(s[n].E.d);
1214 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
1215 emul(&e1->x.p->arr[2*j+1], &s[n].E);
1216 s[n].D = d;
1217 ++n;
1218 }
1219 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
1220 value_clear(res->x.p->arr[2*i].d);
1221 free_evalue_refs(&res->x.p->arr[2*i+1]);
1222 }
1223
1224 free(res->x.p);
1225 res->x.p = new_enode(partition, 2*n, -1);
1226 for (j = 0; j < n; ++j) {
1227 s[j].D = DomainConstraintSimplify(s[j].D, 0);
1228 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1229 value_clear(res->x.p->arr[2*j+1].d);
1230 res->x.p->arr[2*j+1] = s[j].E;
1231 }
1232
1233 free(s);
1234 }
1235
1236 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1237
1238 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1239 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1240 * evalues is not treated here */
1241
1242 void emul (evalue *e1, evalue *res ){
1243 int i,j;
1244
1245 if((value_zero_p(e1->d)&&e1->x.p->type==evector)||(value_zero_p(res->d)&&(res->x.p->type==evector))) {
1246 fprintf(stderr, "emul: do not proced on evector type !\n");
1247 return;
1248 }
1249
1250 if (value_zero_p(e1->d) && e1->x.p->type == partition) {
1251 if (value_zero_p(res->d) && res->x.p->type == partition)
1252 emul_partitions(e1, res);
1253 else
1254 emul_rev(e1, res);
1255 } else if (value_zero_p(res->d) && res->x.p->type == partition) {
1256 for (i = 0; i < res->x.p->size/2; ++i)
1257 emul(e1, &res->x.p->arr[2*i+1]);
1258 } else
1259 if (value_zero_p(res->d) && res->x.p->type == relation) {
1260 if (value_zero_p(e1->d) && e1->x.p->type == relation &&
1261 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1262 if (res->x.p->size < 3 && e1->x.p->size == 3)
1263 explicit_complement(res);
1264 if (e1->x.p->size < 3 && res->x.p->size == 3)
1265 explicit_complement(e1);
1266 for (i = 1; i < res->x.p->size; ++i)
1267 emul(&e1->x.p->arr[i], &res->x.p->arr[i]);
1268 return;
1269 }
1270 for (i = 1; i < res->x.p->size; ++i)
1271 emul(e1, &res->x.p->arr[i]);
1272 } else
1273 if(value_zero_p(e1->d)&& value_zero_p(res->d)) {
1274 switch(e1->x.p->type) {
1275 case polynomial:
1276 switch(res->x.p->type) {
1277 case polynomial:
1278 if(e1->x.p->pos == res->x.p->pos) {
1279 /* Product of two polynomials of the same variable */
1280 emul_poly(e1, res);
1281 return;
1282 }
1283 else {
1284 /* Product of two polynomials of different variables */
1285
1286 if(res->x.p->pos < e1->x.p->pos)
1287 for( i=0; i<res->x.p->size ; i++)
1288 emul(e1, &res->x.p->arr[i]);
1289 else
1290 emul_rev(e1, res);
1291
1292 return ;
1293 }
1294 case periodic:
1295 case fractional:
1296 /* Product of a polynomial and a periodic or fractional */
1297 emul_rev(e1, res);
1298 return;
1299 }
1300 case periodic:
1301 switch(res->x.p->type) {
1302 case periodic:
1303 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size==res->x.p->size) {
1304 /* Product of two periodics of the same parameter and period */
1305
1306 for(i=0; i<res->x.p->size;i++)
1307 emul(&(e1->x.p->arr[i]), &(res->x.p->arr[i]));
1308
1309 return;
1310 }
1311 else{
1312 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size!=res->x.p->size) {
1313 /* Product of two periodics of the same parameter and different periods */
1314 evalue *newp;
1315 Value x,y,z;
1316 int ix,iy,lcm;
1317 value_init(x); value_init(y);value_init(z);
1318 ix=e1->x.p->size;
1319 iy=res->x.p->size;
1320 value_set_si(x,e1->x.p->size);
1321 value_set_si(y,res->x.p->size);
1322 value_assign (z,*Lcm(x,y));
1323 lcm=(int)mpz_get_si(z);
1324 newp= (evalue *) malloc (sizeof(evalue));
1325 value_init(newp->d);
1326 value_set_si( newp->d,0);
1327 newp->x.p=new_enode(periodic,lcm, e1->x.p->pos);
1328 for(i=0;i<lcm;i++) {
1329 evalue_copy(&newp->x.p->arr[i],
1330 &res->x.p->arr[i%iy]);
1331 }
1332 for(i=0;i<lcm;i++)
1333 emul(&e1->x.p->arr[i%ix], &newp->x.p->arr[i]);
1334
1335 value_assign(res->d,newp->d);
1336 res->x.p=newp->x.p;
1337
1338 value_clear(x); value_clear(y);value_clear(z);
1339 return ;
1340 }
1341 else {
1342 /* Product of two periodics of different parameters */
1343
1344 if(res->x.p->pos < e1->x.p->pos)
1345 for(i=0; i<res->x.p->size; i++)
1346 emul(e1, &(res->x.p->arr[i]));
1347 else
1348 emul_rev(e1, res);
1349
1350 return;
1351 }
1352 }
1353 case polynomial:
1354 /* Product of a periodic and a polynomial */
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 fractional:
1363 switch(res->x.p->type) {
1364 case polynomial:
1365 for(i=0; i<res->x.p->size ; i++)
1366 emul(e1, &(res->x.p->arr[i]));
1367 return;
1368 case periodic:
1369 assert(0);
1370 case fractional:
1371 if (e1->x.p->pos == res->x.p->pos &&
1372 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1373 evalue d;
1374 value_init(d.d);
1375 poly_denom(&e1->x.p->arr[0], &d.d);
1376 if (!value_two_p(d.d))
1377 emul_poly(e1, res);
1378 else {
1379 value_init(d.x.n);
1380 value_set_si(d.x.n, 1);
1381 evalue tmp;
1382 /* { x }^2 == { x }/2 */
1383 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1384 assert(e1->x.p->size == 3);
1385 assert(res->x.p->size == 3);
1386 value_init(tmp.d);
1387 evalue_copy(&tmp, &res->x.p->arr[2]);
1388 emul(&d, &tmp);
1389 eadd(&res->x.p->arr[1], &tmp);
1390 emul(&e1->x.p->arr[2], &tmp);
1391 emul(&e1->x.p->arr[1], res);
1392 eadd(&tmp, &res->x.p->arr[2]);
1393 free_evalue_refs(&tmp);
1394 value_clear(d.x.n);
1395 }
1396 value_clear(d.d);
1397 } else {
1398 if(mod_term_smaller(res, e1))
1399 for(i=1; i<res->x.p->size ; i++)
1400 emul(e1, &(res->x.p->arr[i]));
1401 else
1402 emul_rev(e1, res);
1403 return;
1404 }
1405 }
1406 break;
1407 case relation:
1408 emul_rev(e1, res);
1409 break;
1410 default:
1411 assert(0);
1412 }
1413 }
1414 else {
1415 if (value_notzero_p(e1->d)&& value_notzero_p(res->d)) {
1416 /* Product of two rational numbers */
1417
1418 Value g;
1419 value_init(g);
1420 value_multiply(res->d,e1->d,res->d);
1421 value_multiply(res->x.n,e1->x.n,res->x.n );
1422 Gcd(res->x.n, res->d,&g);
1423 if (value_notone_p(g)) {
1424 value_division(res->d,res->d,g);
1425 value_division(res->x.n,res->x.n,g);
1426 }
1427 value_clear(g);
1428 return ;
1429 }
1430 else {
1431 if(value_zero_p(e1->d)&& value_notzero_p(res->d)) {
1432 /* Product of an expression (polynomial or peririodic) and a rational number */
1433
1434 emul_rev(e1, res);
1435 return ;
1436 }
1437 else {
1438 /* Product of a rationel number and an expression (polynomial or peririodic) */
1439
1440 i = res->x.p->type == fractional ? 1 : 0;
1441 for (; i<res->x.p->size; i++)
1442 emul(e1, &res->x.p->arr[i]);
1443
1444 return ;
1445 }
1446 }
1447 }
1448
1449 return ;
1450 }
1451
1452 /* Frees mask content ! */
1453 void emask(evalue *mask, evalue *res) {
1454 int n, i, j;
1455 Polyhedron *d, *fd;
1456 struct section *s;
1457 evalue mone;
1458
1459 if (EVALUE_IS_ZERO(*res)) {
1460 free_evalue_refs(mask);
1461 return;
1462 }
1463
1464 assert(value_zero_p(mask->d));
1465 assert(mask->x.p->type == partition);
1466 assert(value_zero_p(res->d));
1467 assert(res->x.p->type == partition);
1468
1469 s = (struct section *)
1470 malloc((mask->x.p->size/2+1) * (res->x.p->size/2) *
1471 sizeof(struct section));
1472 assert(s);
1473
1474 value_init(mone.d);
1475 evalue_set_si(&mone, -1, 1);
1476
1477 n = 0;
1478 for (j = 0; j < res->x.p->size/2; ++j) {
1479 assert(mask->x.p->size >= 2);
1480 fd = DomainDifference(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1481 EVALUE_DOMAIN(mask->x.p->arr[0]), 0);
1482 if (!emptyQ(fd))
1483 for (i = 1; i < mask->x.p->size/2; ++i) {
1484 Polyhedron *t = fd;
1485 fd = DomainDifference(fd, EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1486 Domain_Free(t);
1487 if (emptyQ(fd))
1488 break;
1489 }
1490 if (emptyQ(fd)) {
1491 Domain_Free(fd);
1492 continue;
1493 }
1494 value_init(s[n].E.d);
1495 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1496 s[n].D = fd;
1497 ++n;
1498 }
1499 for (i = 0; i < mask->x.p->size/2; ++i) {
1500 if (EVALUE_IS_ONE(mask->x.p->arr[2*i+1]))
1501 continue;
1502
1503 fd = EVALUE_DOMAIN(mask->x.p->arr[2*i]);
1504 eadd(&mone, &mask->x.p->arr[2*i+1]);
1505 emul(&mone, &mask->x.p->arr[2*i+1]);
1506 for (j = 0; j < res->x.p->size/2; ++j) {
1507 Polyhedron *t;
1508 d = DomainIntersection(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1509 EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1510 if (emptyQ(d)) {
1511 Domain_Free(d);
1512 continue;
1513 }
1514 t = fd;
1515 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*j]), 0);
1516 if (t != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1517 Domain_Free(t);
1518 value_init(s[n].E.d);
1519 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1520 emul(&mask->x.p->arr[2*i+1], &s[n].E);
1521 s[n].D = d;
1522 ++n;
1523 }
1524
1525 if (!emptyQ(fd)) {
1526 /* Just ignore; this may have been previously masked off */
1527 }
1528 if (fd != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1529 Domain_Free(fd);
1530 }
1531
1532 free_evalue_refs(&mone);
1533 free_evalue_refs(mask);
1534 free_evalue_refs(res);
1535 value_init(res->d);
1536 if (n == 0)
1537 evalue_set_si(res, 0, 1);
1538 else {
1539 res->x.p = new_enode(partition, 2*n, -1);
1540 for (j = 0; j < n; ++j) {
1541 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1542 value_clear(res->x.p->arr[2*j+1].d);
1543 res->x.p->arr[2*j+1] = s[j].E;
1544 }
1545 }
1546
1547 free(s);
1548 }
1549
1550 void evalue_copy(evalue *dst, evalue *src)
1551 {
1552 value_assign(dst->d, src->d);
1553 if(value_notzero_p(src->d)) {
1554 value_init(dst->x.n);
1555 value_assign(dst->x.n, src->x.n);
1556 } else
1557 dst->x.p = ecopy(src->x.p);
1558 }
1559
1560 enode *new_enode(enode_type type,int size,int pos) {
1561
1562 enode *res;
1563 int i;
1564
1565 if(size == 0) {
1566 fprintf(stderr, "Allocating enode of size 0 !\n" );
1567 return NULL;
1568 }
1569 res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
1570 res->type = type;
1571 res->size = size;
1572 res->pos = pos;
1573 for(i=0; i<size; i++) {
1574 value_init(res->arr[i].d);
1575 value_set_si(res->arr[i].d,0);
1576 res->arr[i].x.p = 0;
1577 }
1578 return res;
1579 } /* new_enode */
1580
1581 enode *ecopy(enode *e) {
1582
1583 enode *res;
1584 int i;
1585
1586 res = new_enode(e->type,e->size,e->pos);
1587 for(i=0;i<e->size;++i) {
1588 value_assign(res->arr[i].d,e->arr[i].d);
1589 if(value_zero_p(res->arr[i].d))
1590 res->arr[i].x.p = ecopy(e->arr[i].x.p);
1591 else if (EVALUE_IS_DOMAIN(res->arr[i]))
1592 EVALUE_SET_DOMAIN(res->arr[i], Domain_Copy(EVALUE_DOMAIN(e->arr[i])));
1593 else {
1594 value_init(res->arr[i].x.n);
1595 value_assign(res->arr[i].x.n,e->arr[i].x.n);
1596 }
1597 }
1598 return(res);
1599 } /* ecopy */
1600
1601 int ecmp(const evalue *e1, const evalue *e2)
1602 {
1603 enode *p1, *p2;
1604 int i;
1605 int r;
1606
1607 if (value_notzero_p(e1->d) && value_notzero_p(e2->d)) {
1608 Value m, m2;
1609 value_init(m);
1610 value_init(m2);
1611 value_multiply(m, e1->x.n, e2->d);
1612 value_multiply(m2, e2->x.n, e1->d);
1613
1614 if (value_lt(m, m2))
1615 r = -1;
1616 else if (value_gt(m, m2))
1617 r = 1;
1618 else
1619 r = 0;
1620
1621 value_clear(m);
1622 value_clear(m2);
1623
1624 return r;
1625 }
1626 if (value_notzero_p(e1->d))
1627 return -1;
1628 if (value_notzero_p(e2->d))
1629 return 1;
1630
1631 p1 = e1->x.p;
1632 p2 = e2->x.p;
1633
1634 if (p1->type != p2->type)
1635 return p1->type - p2->type;
1636 if (p1->pos != p2->pos)
1637 return p1->pos - p2->pos;
1638 if (p1->size != p2->size)
1639 return p1->size - p2->size;
1640
1641 for (i = p1->size-1; i >= 0; --i)
1642 if ((r = ecmp(&p1->arr[i], &p2->arr[i])) != 0)
1643 return r;
1644
1645 return 0;
1646 }
1647
1648 int eequal(evalue *e1,evalue *e2) {
1649
1650 int i;
1651 enode *p1, *p2;
1652
1653 if (value_ne(e1->d,e2->d))
1654 return 0;
1655
1656 /* e1->d == e2->d */
1657 if (value_notzero_p(e1->d)) {
1658 if (value_ne(e1->x.n,e2->x.n))
1659 return 0;
1660
1661 /* e1->d == e2->d != 0 AND e1->n == e2->n */
1662 return 1;
1663 }
1664
1665 /* e1->d == e2->d == 0 */
1666 p1 = e1->x.p;
1667 p2 = e2->x.p;
1668 if (p1->type != p2->type) return 0;
1669 if (p1->size != p2->size) return 0;
1670 if (p1->pos != p2->pos) return 0;
1671 for (i=0; i<p1->size; i++)
1672 if (!eequal(&p1->arr[i], &p2->arr[i]) )
1673 return 0;
1674 return 1;
1675 } /* eequal */
1676
1677 void free_evalue_refs(evalue *e) {
1678
1679 enode *p;
1680 int i;
1681
1682 if (EVALUE_IS_DOMAIN(*e)) {
1683 Domain_Free(EVALUE_DOMAIN(*e));
1684 value_clear(e->d);
1685 return;
1686 } else if (value_pos_p(e->d)) {
1687
1688 /* 'e' stores a constant */
1689 value_clear(e->d);
1690 value_clear(e->x.n);
1691 return;
1692 }
1693 assert(value_zero_p(e->d));
1694 value_clear(e->d);
1695 p = e->x.p;
1696 if (!p) return; /* null pointer */
1697 for (i=0; i<p->size; i++) {
1698 free_evalue_refs(&(p->arr[i]));
1699 }
1700 free(p);
1701 return;
1702 } /* free_evalue_refs */
1703
1704 static void mod2table_r(evalue *e, Vector *periods, Value m, int p,
1705 Vector * val, evalue *res)
1706 {
1707 unsigned nparam = periods->Size;
1708
1709 if (p == nparam) {
1710 double d = compute_evalue(e, val->p);
1711 d *= VALUE_TO_DOUBLE(m);
1712 if (d > 0)
1713 d += .25;
1714 else
1715 d -= .25;
1716 value_assign(res->d, m);
1717 value_init(res->x.n);
1718 value_set_double(res->x.n, d);
1719 mpz_fdiv_r(res->x.n, res->x.n, m);
1720 return;
1721 }
1722 if (value_one_p(periods->p[p]))
1723 mod2table_r(e, periods, m, p+1, val, res);
1724 else {
1725 Value tmp;
1726 value_init(tmp);
1727
1728 value_assign(tmp, periods->p[p]);
1729 value_set_si(res->d, 0);
1730 res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
1731 do {
1732 value_decrement(tmp, tmp);
1733 value_assign(val->p[p], tmp);
1734 mod2table_r(e, periods, m, p+1, val,
1735 &res->x.p->arr[VALUE_TO_INT(tmp)]);
1736 } while (value_pos_p(tmp));
1737
1738 value_clear(tmp);
1739 }
1740 }
1741
1742 static void rel2table(evalue *e, int zero)
1743 {
1744 if (value_pos_p(e->d)) {
1745 if (value_zero_p(e->x.n) == zero)
1746 value_set_si(e->x.n, 1);
1747 else
1748 value_set_si(e->x.n, 0);
1749 value_set_si(e->d, 1);
1750 } else {
1751 int i;
1752 for (i = 0; i < e->x.p->size; ++i)
1753 rel2table(&e->x.p->arr[i], zero);
1754 }
1755 }
1756
1757 void evalue_mod2table(evalue *e, int nparam)
1758 {
1759 enode *p;
1760 int i;
1761
1762 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
1763 return;
1764 p = e->x.p;
1765 for (i=0; i<p->size; i++) {
1766 evalue_mod2table(&(p->arr[i]), nparam);
1767 }
1768 if (p->type == relation) {
1769 evalue copy;
1770
1771 if (p->size > 2) {
1772 value_init(copy.d);
1773 evalue_copy(&copy, &p->arr[0]);
1774 }
1775 rel2table(&p->arr[0], 1);
1776 emul(&p->arr[0], &p->arr[1]);
1777 if (p->size > 2) {
1778 rel2table(&copy, 0);
1779 emul(&copy, &p->arr[2]);
1780 eadd(&p->arr[2], &p->arr[1]);
1781 free_evalue_refs(&p->arr[2]);
1782 free_evalue_refs(&copy);
1783 }
1784 free_evalue_refs(&p->arr[0]);
1785 value_clear(e->d);
1786 *e = p->arr[1];
1787 free(p);
1788 } else if (p->type == fractional) {
1789 Vector *periods = Vector_Alloc(nparam);
1790 Vector *val = Vector_Alloc(nparam);
1791 Value tmp;
1792 evalue *ev;
1793 evalue EP, res;
1794
1795 value_init(tmp);
1796 value_set_si(tmp, 1);
1797 Vector_Set(periods->p, 1, nparam);
1798 Vector_Set(val->p, 0, nparam);
1799 for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
1800 enode *p = ev->x.p;
1801
1802 assert(p->type == polynomial);
1803 assert(p->size == 2);
1804 value_assign(periods->p[p->pos-1], p->arr[1].d);
1805 Lcm3(tmp, p->arr[1].d, &tmp);
1806 }
1807 value_init(EP.d);
1808 mod2table_r(&p->arr[0], periods, tmp, 0, val, &EP);
1809
1810 value_init(res.d);
1811 evalue_set_si(&res, 0, 1);
1812 /* Compute the polynomial using Horner's rule */
1813 for (i=p->size-1;i>1;i--) {
1814 eadd(&p->arr[i], &res);
1815 emul(&EP, &res);
1816 }
1817 eadd(&p->arr[1], &res);
1818
1819 free_evalue_refs(e);
1820 free_evalue_refs(&EP);
1821 *e = res;
1822
1823 value_clear(tmp);
1824 Vector_Free(val);
1825 Vector_Free(periods);
1826 }
1827 } /* evalue_mod2table */
1828
1829 /********************************************************/
1830 /* function in domain */
1831 /* check if the parameters in list_args */
1832 /* verifies the constraints of Domain P */
1833 /********************************************************/
1834 int in_domain(Polyhedron *P, Value *list_args) {
1835
1836 int col,row;
1837 Value v; /* value of the constraint of a row when
1838 parameters are instanciated*/
1839 Value tmp;
1840
1841 value_init(v);
1842 value_init(tmp);
1843
1844 /*P->Constraint constraint matrice of polyhedron P*/
1845 for(row=0;row<P->NbConstraints;row++) {
1846 value_assign(v,P->Constraint[row][P->Dimension+1]); /*constant part*/
1847 for(col=1;col<P->Dimension+1;col++) {
1848 value_multiply(tmp,P->Constraint[row][col],list_args[col-1]);
1849 value_addto(v,v,tmp);
1850 }
1851 if (value_notzero_p(P->Constraint[row][0])) {
1852
1853 /*if v is not >=0 then this constraint is not respected */
1854 if (value_neg_p(v)) {
1855 next:
1856 value_clear(v);
1857 value_clear(tmp);
1858 return P->next ? in_domain(P->next, list_args) : 0;
1859 }
1860 }
1861 else {
1862
1863 /*if v is not = 0 then this constraint is not respected */
1864 if (value_notzero_p(v))
1865 goto next;
1866 }
1867 }
1868
1869 /*if not return before this point => all
1870 the constraints are respected */
1871 value_clear(v);
1872 value_clear(tmp);
1873 return 1;
1874 } /* in_domain */
1875
1876 /****************************************************/
1877 /* function compute enode */
1878 /* compute the value of enode p with parameters */
1879 /* list "list_args */
1880 /* compute the polynomial or the periodic */
1881 /****************************************************/
1882
1883 static double compute_enode(enode *p, Value *list_args) {
1884
1885 int i;
1886 Value m, param;
1887 double res=0.0;
1888
1889 if (!p)
1890 return(0.);
1891
1892 value_init(m);
1893 value_init(param);
1894
1895 if (p->type == polynomial) {
1896 if (p->size > 1)
1897 value_assign(param,list_args[p->pos-1]);
1898
1899 /* Compute the polynomial using Horner's rule */
1900 for (i=p->size-1;i>0;i--) {
1901 res +=compute_evalue(&p->arr[i],list_args);
1902 res *=VALUE_TO_DOUBLE(param);
1903 }
1904 res +=compute_evalue(&p->arr[0],list_args);
1905 }
1906 else if (p->type == fractional) {
1907 double d = compute_evalue(&p->arr[0], list_args);
1908 d -= floor(d+1e-10);
1909
1910 /* Compute the polynomial using Horner's rule */
1911 for (i=p->size-1;i>1;i--) {
1912 res +=compute_evalue(&p->arr[i],list_args);
1913 res *=d;
1914 }
1915 res +=compute_evalue(&p->arr[1],list_args);
1916 }
1917 else if (p->type == periodic) {
1918 value_assign(param,list_args[p->pos-1]);
1919
1920 /* Choose the right element of the periodic */
1921 value_absolute(m,param);
1922 value_set_si(param,p->size);
1923 value_modulus(m,m,param);
1924 res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
1925 }
1926 else if (p->type == relation) {
1927 if (fabs(compute_evalue(&p->arr[0], list_args)) < 1e-10)
1928 res = compute_evalue(&p->arr[1], list_args);
1929 else if (p->size > 2)
1930 res = compute_evalue(&p->arr[2], list_args);
1931 }
1932 else if (p->type == partition) {
1933 for (i = 0; i < p->size/2; ++i)
1934 if (in_domain(EVALUE_DOMAIN(p->arr[2*i]), list_args)) {
1935 res = compute_evalue(&p->arr[2*i+1], list_args);
1936 break;
1937 }
1938 }
1939 value_clear(m);
1940 value_clear(param);
1941 return res;
1942 } /* compute_enode */
1943
1944 /*************************************************/
1945 /* return the value of Ehrhart Polynomial */
1946 /* It returns a double, because since it is */
1947 /* a recursive function, some intermediate value */
1948 /* might not be integral */
1949 /*************************************************/
1950
1951 double compute_evalue(evalue *e,Value *list_args) {
1952
1953 double res;
1954
1955 if (value_notzero_p(e->d)) {
1956 if (value_notone_p(e->d))
1957 res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
1958 else
1959 res = VALUE_TO_DOUBLE(e->x.n);
1960 }
1961 else
1962 res = compute_enode(e->x.p,list_args);
1963 return res;
1964 } /* compute_evalue */
1965
1966
1967 /****************************************************/
1968 /* function compute_poly : */
1969 /* Check for the good validity domain */
1970 /* return the number of point in the Polyhedron */
1971 /* in allocated memory */
1972 /* Using the Ehrhart pseudo-polynomial */
1973 /****************************************************/
1974 Value *compute_poly(Enumeration *en,Value *list_args) {
1975
1976 Value *tmp;
1977 /* double d; int i; */
1978
1979 tmp = (Value *) malloc (sizeof(Value));
1980 assert(tmp != NULL);
1981 value_init(*tmp);
1982 value_set_si(*tmp,0);
1983
1984 if(!en)
1985 return(tmp); /* no ehrhart polynomial */
1986 if(en->ValidityDomain) {
1987 if(!en->ValidityDomain->Dimension) { /* no parameters */
1988 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
1989 return(tmp);
1990 }
1991 }
1992 else
1993 return(tmp); /* no Validity Domain */
1994 while(en) {
1995 if(in_domain(en->ValidityDomain,list_args)) {
1996
1997 #ifdef EVAL_EHRHART_DEBUG
1998 Print_Domain(stdout,en->ValidityDomain);
1999 print_evalue(stdout,&en->EP);
2000 #endif
2001
2002 /* d = compute_evalue(&en->EP,list_args);
2003 i = d;
2004 printf("(double)%lf = %d\n", d, i ); */
2005 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2006 return(tmp);
2007 }
2008 else
2009 en=en->next;
2010 }
2011 value_set_si(*tmp,0);
2012 return(tmp); /* no compatible domain with the arguments */
2013 } /* compute_poly */
2014
2015 size_t value_size(Value v) {
2016 return (v[0]._mp_size > 0 ? v[0]._mp_size : -v[0]._mp_size)
2017 * sizeof(v[0]._mp_d[0]);
2018 }
2019
2020 size_t domain_size(Polyhedron *D)
2021 {
2022 int i, j;
2023 size_t s = sizeof(*D);
2024
2025 for (i = 0; i < D->NbConstraints; ++i)
2026 for (j = 0; j < D->Dimension+2; ++j)
2027 s += value_size(D->Constraint[i][j]);
2028
2029 /*
2030 for (i = 0; i < D->NbRays; ++i)
2031 for (j = 0; j < D->Dimension+2; ++j)
2032 s += value_size(D->Ray[i][j]);
2033 */
2034
2035 return D->next ? s+domain_size(D->next) : s;
2036 }
2037
2038 size_t enode_size(enode *p) {
2039 size_t s = sizeof(*p) - sizeof(p->arr[0]);
2040 int i;
2041
2042 if (p->type == partition)
2043 for (i = 0; i < p->size/2; ++i) {
2044 s += domain_size(EVALUE_DOMAIN(p->arr[2*i]));
2045 s += evalue_size(&p->arr[2*i+1]);
2046 }
2047 else
2048 for (i = 0; i < p->size; ++i) {
2049 s += evalue_size(&p->arr[i]);
2050 }
2051 return s;
2052 }
2053
2054 size_t evalue_size(evalue *e)
2055 {
2056 size_t s = sizeof(*e);
2057 s += value_size(e->d);
2058 if (value_notzero_p(e->d))
2059 s += value_size(e->x.n);
2060 else
2061 s += enode_size(e->x.p);
2062 return s;
2063 }
2064
2065 static evalue *find_second(evalue *base, evalue *cst, evalue *e, Value m)
2066 {
2067 evalue *found = NULL;
2068 evalue offset;
2069 evalue copy;
2070 int i;
2071
2072 if (value_pos_p(e->d) || e->x.p->type != fractional)
2073 return NULL;
2074
2075 value_init(offset.d);
2076 value_init(offset.x.n);
2077 poly_denom(&e->x.p->arr[0], &offset.d);
2078 Lcm3(m, offset.d, &offset.d);
2079 value_set_si(offset.x.n, 1);
2080
2081 value_init(copy.d);
2082 evalue_copy(&copy, cst);
2083
2084 eadd(&offset, cst);
2085 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2086
2087 if (eequal(base, &e->x.p->arr[0]))
2088 found = &e->x.p->arr[0];
2089 else {
2090 value_set_si(offset.x.n, -2);
2091
2092 eadd(&offset, cst);
2093 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2094
2095 if (eequal(base, &e->x.p->arr[0]))
2096 found = base;
2097 }
2098 free_evalue_refs(cst);
2099 free_evalue_refs(&offset);
2100 *cst = copy;
2101
2102 for (i = 1; !found && i < e->x.p->size; ++i)
2103 found = find_second(base, cst, &e->x.p->arr[i], m);
2104
2105 return found;
2106 }
2107
2108 static evalue *find_relation_pair(evalue *e)
2109 {
2110 int i;
2111 evalue *found = NULL;
2112
2113 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2114 return NULL;
2115
2116 if (e->x.p->type == fractional) {
2117 Value m;
2118 evalue *cst;
2119
2120 value_init(m);
2121 poly_denom(&e->x.p->arr[0], &m);
2122
2123 for (cst = &e->x.p->arr[0]; value_zero_p(cst->d);
2124 cst = &cst->x.p->arr[0])
2125 ;
2126
2127 for (i = 1; !found && i < e->x.p->size; ++i)
2128 found = find_second(&e->x.p->arr[0], cst, &e->x.p->arr[i], m);
2129
2130 value_clear(m);
2131 }
2132
2133 i = e->x.p->type == relation;
2134 for (; !found && i < e->x.p->size; ++i)
2135 found = find_relation_pair(&e->x.p->arr[i]);
2136
2137 return found;
2138 }
2139
2140 void evalue_mod2relation(evalue *e) {
2141 evalue *d;
2142
2143 if (value_zero_p(e->d) && e->x.p->type == partition) {
2144 int i;
2145
2146 for (i = 0; i < e->x.p->size/2; ++i) {
2147 evalue_mod2relation(&e->x.p->arr[2*i+1]);
2148 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2149 value_clear(e->x.p->arr[2*i].d);
2150 free_evalue_refs(&e->x.p->arr[2*i+1]);
2151 e->x.p->size -= 2;
2152 if (2*i < e->x.p->size) {
2153 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2154 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2155 }
2156 --i;
2157 }
2158 }
2159 if (e->x.p->size == 0) {
2160 free(e->x.p);
2161 evalue_set_si(e, 0, 1);
2162 }
2163
2164 return;
2165 }
2166
2167 while ((d = find_relation_pair(e)) != NULL) {
2168 evalue split;
2169 evalue *ev;
2170
2171 value_init(split.d);
2172 value_set_si(split.d, 0);
2173 split.x.p = new_enode(relation, 3, 0);
2174 evalue_set_si(&split.x.p->arr[1], 1, 1);
2175 evalue_set_si(&split.x.p->arr[2], 1, 1);
2176
2177 ev = &split.x.p->arr[0];
2178 value_set_si(ev->d, 0);
2179 ev->x.p = new_enode(fractional, 3, -1);
2180 evalue_set_si(&ev->x.p->arr[1], 0, 1);
2181 evalue_set_si(&ev->x.p->arr[2], 1, 1);
2182 evalue_copy(&ev->x.p->arr[0], d);
2183
2184 emul(&split, e);
2185
2186 reduce_evalue(e);
2187
2188 free_evalue_refs(&split);
2189 }
2190 }
2191
2192 static int evalue_comp(const void * a, const void * b)
2193 {
2194 const evalue *e1 = *(const evalue **)a;
2195 const evalue *e2 = *(const evalue **)b;
2196 return ecmp(e1, e2);
2197 }
2198
2199 void evalue_combine(evalue *e)
2200 {
2201 evalue **evs;
2202 int i, k;
2203 enode *p;
2204 evalue tmp;
2205
2206 if (value_notzero_p(e->d) || e->x.p->type != partition)
2207 return;
2208
2209 NALLOC(evs, e->x.p->size/2);
2210 for (i = 0; i < e->x.p->size/2; ++i)
2211 evs[i] = &e->x.p->arr[2*i+1];
2212 qsort(evs, e->x.p->size/2, sizeof(evs[0]), evalue_comp);
2213 p = new_enode(partition, e->x.p->size, -1);
2214 for (i = 0, k = 0; i < e->x.p->size/2; ++i) {
2215 if (k == 0 || ecmp(&p->arr[2*k-1], evs[i]) != 0) {
2216 value_clear(p->arr[2*k].d);
2217 value_clear(p->arr[2*k+1].d);
2218 p->arr[2*k] = *(evs[i]-1);
2219 p->arr[2*k+1] = *(evs[i]);
2220 ++k;
2221 } else {
2222 Polyhedron *D = EVALUE_DOMAIN(*(evs[i]-1));
2223 Polyhedron *L = D;
2224
2225 value_clear((evs[i]-1)->d);
2226
2227 while (L->next)
2228 L = L->next;
2229 L->next = EVALUE_DOMAIN(p->arr[2*k-2]);
2230 EVALUE_SET_DOMAIN(p->arr[2*k-2], D);
2231 free_evalue_refs(evs[i]);
2232 }
2233 }
2234
2235 for (i = 2*k ; i < p->size; ++i)
2236 value_clear(p->arr[i].d);
2237
2238 free(evs);
2239 free(e->x.p);
2240 p->size = 2*k;
2241 e->x.p = p;
2242
2243 for (i = 0; i < e->x.p->size/2; ++i) {
2244 Polyhedron *H;
2245 if (value_notzero_p(e->x.p->arr[2*i+1].d))
2246 continue;
2247 H = DomainConvex(EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
2248 if (H == NULL)
2249 continue;
2250 for (k = 0; k < e->x.p->size/2; ++k) {
2251 Polyhedron *D, *N, **P;
2252 if (i == k)
2253 continue;
2254 P = &EVALUE_DOMAIN(e->x.p->arr[2*k]);
2255 D = *P;
2256 if (D == NULL)
2257 continue;
2258 for (; D; D = N) {
2259 *P = D;
2260 N = D->next;
2261 if (D->NbEq <= H->NbEq) {
2262 P = &D->next;
2263 continue;
2264 }
2265
2266 value_init(tmp.d);
2267 tmp.x.p = new_enode(partition, 2, -1);
2268 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Copy(D));
2269 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
2270 reduce_evalue(&tmp);
2271 if (value_notzero_p(tmp.d) ||
2272 ecmp(&tmp.x.p->arr[1], &e->x.p->arr[2*k+1]) != 0)
2273 P = &D->next;
2274 else {
2275 D->next = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2276 EVALUE_DOMAIN(e->x.p->arr[2*i]) = D;
2277 *P = NULL;
2278 }
2279 free_evalue_refs(&tmp);
2280 }
2281 }
2282 Polyhedron_Free(H);
2283 }
2284
2285 for (i = 0; i < e->x.p->size/2; ++i) {
2286 Polyhedron *H, *E;
2287 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2288 if (!D) {
2289 value_clear(e->x.p->arr[2*i].d);
2290 free_evalue_refs(&e->x.p->arr[2*i+1]);
2291 e->x.p->size -= 2;
2292 if (2*i < e->x.p->size) {
2293 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2294 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2295 }
2296 --i;
2297 continue;
2298 }
2299 if (!D->next)
2300 continue;
2301 H = DomainConvex(D, 0);
2302 E = DomainDifference(H, D, 0);
2303 Domain_Free(D);
2304 D = DomainDifference(H, E, 0);
2305 Domain_Free(H);
2306 Domain_Free(E);
2307 EVALUE_SET_DOMAIN(p->arr[2*i], D);
2308 }
2309 }
2310
2311 /* May change coefficients to become non-standard
2312 * => reduce p afterwards to correct
2313 */
2314 static Polyhedron *polynomial_projection(enode *p, Polyhedron *D, Value *d,
2315 Matrix **R)
2316 {
2317 Polyhedron *I, *H;
2318 evalue *pp;
2319 unsigned dim = D->Dimension;
2320 Matrix *T = Matrix_Alloc(2, dim+1);
2321 Value twice;
2322 value_init(twice);
2323 assert(T);
2324
2325 assert(p->type == fractional);
2326 pp = &p->arr[0];
2327 value_set_si(T->p[1][dim], 1);
2328 poly_denom(pp, d);
2329 while (value_zero_p(pp->d)) {
2330 assert(pp->x.p->type == polynomial);
2331 assert(pp->x.p->size == 2);
2332 assert(value_notzero_p(pp->x.p->arr[1].d));
2333 value_division(T->p[0][pp->x.p->pos-1], *d, pp->x.p->arr[1].d);
2334 mpz_mul_ui(twice, pp->x.p->arr[1].x.n, 2);
2335 if (value_gt(twice, pp->x.p->arr[1].d))
2336 value_substract(pp->x.p->arr[1].x.n,
2337 pp->x.p->arr[1].x.n, pp->x.p->arr[1].d);
2338 value_multiply(T->p[0][pp->x.p->pos-1],
2339 T->p[0][pp->x.p->pos-1], pp->x.p->arr[1].x.n);
2340 pp = &pp->x.p->arr[0];
2341 }
2342 value_division(T->p[0][dim], *d, pp->d);
2343 value_multiply(T->p[0][dim], T->p[0][dim], pp->x.n);
2344 I = DomainImage(D, T, 0);
2345 H = DomainConvex(I, 0);
2346 Domain_Free(I);
2347 *R = T;
2348
2349 value_clear(twice);
2350
2351 return H;
2352 }
2353
2354 static int reduce_in_domain(evalue *e, Polyhedron *D)
2355 {
2356 int i;
2357 enode *p;
2358 Value d, min, max;
2359 int r = 0;
2360 Polyhedron *I;
2361 Matrix *T;
2362
2363 if (value_notzero_p(e->d))
2364 return r;
2365
2366 p = e->x.p;
2367
2368 if (p->type == relation) {
2369 int equal;
2370 value_init(d);
2371 value_init(min);
2372 value_init(max);
2373
2374 I = polynomial_projection(p->arr[0].x.p, D, &d, &T);
2375 line_minmax(I, &min, &max); /* frees I */
2376 equal = value_eq(min, max);
2377 mpz_cdiv_q(min, min, d);
2378 mpz_fdiv_q(max, max, d);
2379
2380 if (value_gt(min, max)) {
2381 /* Never zero */
2382 if (p->size == 3) {
2383 value_clear(e->d);
2384 *e = p->arr[2];
2385 } else {
2386 evalue_set_si(e, 0, 1);
2387 r = 1;
2388 }
2389 free_evalue_refs(&(p->arr[1]));
2390 free_evalue_refs(&(p->arr[0]));
2391 free(p);
2392 value_clear(d);
2393 value_clear(min);
2394 value_clear(max);
2395 Matrix_Free(T);
2396 return r ? r : reduce_in_domain(e, D);
2397 } else if (equal) {
2398 /* Always zero */
2399 if (p->size == 3)
2400 free_evalue_refs(&(p->arr[2]));
2401 value_clear(e->d);
2402 *e = p->arr[1];
2403 free_evalue_refs(&(p->arr[0]));
2404 free(p);
2405 value_clear(d);
2406 value_clear(min);
2407 value_clear(max);
2408 Matrix_Free(T);
2409 return reduce_in_domain(e, D);
2410 } else if (value_eq(min, max)) {
2411 /* zero for a single value */
2412 Polyhedron *E;
2413 Matrix *M = Matrix_Alloc(1, D->Dimension+2);
2414 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
2415 value_multiply(min, min, d);
2416 value_substract(M->p[0][D->Dimension+1],
2417 M->p[0][D->Dimension+1], min);
2418 E = DomainAddConstraints(D, M, 0);
2419 value_clear(d);
2420 value_clear(min);
2421 value_clear(max);
2422 Matrix_Free(T);
2423 Matrix_Free(M);
2424 r = reduce_in_domain(&p->arr[1], E);
2425 if (p->size == 3)
2426 r |= reduce_in_domain(&p->arr[2], D);
2427 Domain_Free(E);
2428 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2429 return r;
2430 }
2431
2432 value_clear(d);
2433 value_clear(min);
2434 value_clear(max);
2435 Matrix_Free(T);
2436 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2437 }
2438
2439 i = p->type == relation ? 1 :
2440 p->type == fractional ? 1 : 0;
2441 for (; i<p->size; i++)
2442 r |= reduce_in_domain(&p->arr[i], D);
2443
2444 if (p->type != fractional) {
2445 if (r && p->type == polynomial) {
2446 evalue f;
2447 value_init(f.d);
2448 value_set_si(f.d, 0);
2449 f.x.p = new_enode(polynomial, 2, p->pos);
2450 evalue_set_si(&f.x.p->arr[0], 0, 1);
2451 evalue_set_si(&f.x.p->arr[1], 1, 1);
2452 reorder_terms(p, &f);
2453 value_clear(e->d);
2454 *e = p->arr[0];
2455 free(p);
2456 }
2457 return r;
2458 }
2459
2460 value_init(d);
2461 value_init(min);
2462 value_init(max);
2463 I = polynomial_projection(p, D, &d, &T);
2464 Matrix_Free(T);
2465 line_minmax(I, &min, &max); /* frees I */
2466 mpz_fdiv_q(min, min, d);
2467 mpz_fdiv_q(max, max, d);
2468
2469 if (value_eq(min, max)) {
2470 evalue inc;
2471 value_init(inc.d);
2472 value_init(inc.x.n);
2473 value_set_si(inc.d, 1);
2474 value_oppose(inc.x.n, min);
2475 eadd(&inc, &p->arr[0]);
2476 reorder_terms(p, &p->arr[0]); /* frees arr[0] */
2477 value_clear(e->d);
2478 *e = p->arr[1];
2479 free(p);
2480 free_evalue_refs(&inc);
2481 r = 1;
2482 } else {
2483 _reduce_evalue(&p->arr[0], 0, 1);
2484 if (r == 1) {
2485 evalue f;
2486 value_init(f.d);
2487 value_set_si(f.d, 0);
2488 f.x.p = new_enode(fractional, 3, -1);
2489 value_clear(f.x.p->arr[0].d);
2490 f.x.p->arr[0] = p->arr[0];
2491 evalue_set_si(&f.x.p->arr[1], 0, 1);
2492 evalue_set_si(&f.x.p->arr[2], 1, 1);
2493 reorder_terms(p, &f);
2494 value_clear(e->d);
2495 *e = p->arr[1];
2496 free(p);
2497 }
2498 }
2499
2500 value_clear(d);
2501 value_clear(min);
2502 value_clear(max);
2503
2504 return r;
2505 }
2506
2507 void evalue_range_reduction(evalue *e)
2508 {
2509 int i;
2510 if (value_notzero_p(e->d) || e->x.p->type != partition)
2511 return;
2512
2513 for (i = 0; i < e->x.p->size/2; ++i)
2514 if (reduce_in_domain(&e->x.p->arr[2*i+1],
2515 EVALUE_DOMAIN(e->x.p->arr[2*i]))) {
2516 reduce_evalue(&e->x.p->arr[2*i+1]);
2517
2518 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2519 free_evalue_refs(&e->x.p->arr[2*i+1]);
2520 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
2521 value_clear(e->x.p->arr[2*i].d);
2522 e->x.p->size -= 2;
2523 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2524 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2525 --i;
2526 }
2527 }
2528 }
lattice point enumeration library