?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
evalue.c: emul: handle another special case
[barvinok.git] / evalue.c
blob52adff07235fd653c915e6f606a5403de7a27bc2
1 /***********************************************************************/
2 /* copyright 1997, Doran Wilde */
3 /* copyright 1997-2000, Vincent Loechner */
4 /* copyright 2003-2006, Sven Verdoolaege */
5 /* Permission is granted to copy, use, and distribute */
6 /* for any commercial or noncommercial purpose under the terms */
7 /* of the GNU General Public license, version 2, June 1991 */
8 /* (see file : LICENSE). */
9 /***********************************************************************/
10
11 #include <assert.h>
12 #include <math.h>
13 #include <stdlib.h>
14 #include <string.h>
15 #include <barvinok/evalue.h>
16 #include <barvinok/barvinok.h>
17 #include <barvinok/util.h>
18
19 #ifndef value_pmodulus
20 #define value_pmodulus(ref,val1,val2) (mpz_fdiv_r((ref),(val1),(val2)))
21 #endif
22
23 #define ALLOC(type) (type*)malloc(sizeof(type))
24 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
25
26 #ifdef __GNUC__
27 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
28 #else
29 #define NALLOC(p,n) p = (void *)malloc((n) * sizeof(*p))
30 #endif
31
32 void evalue_set_si(evalue *ev, int n, int d) {
33 value_set_si(ev->d, d);
34 value_init(ev->x.n);
35 value_set_si(ev->x.n, n);
36 }
37
38 void evalue_set(evalue *ev, Value n, Value d) {
39 value_assign(ev->d, d);
40 value_init(ev->x.n);
41 value_assign(ev->x.n, n);
42 }
43
44 evalue* evalue_zero()
45 {
46 evalue *EP = ALLOC(evalue);
47 value_init(EP->d);
48 evalue_set_si(EP, 0, 1);
49 return EP;
50 }
51
52 void aep_evalue(evalue *e, int *ref) {
53
54 enode *p;
55 int i;
56
57 if (value_notzero_p(e->d))
58 return; /* a rational number, its already reduced */
59 if(!(p = e->x.p))
60 return; /* hum... an overflow probably occured */
61
62 /* First check the components of p */
63 for (i=0;i<p->size;i++)
64 aep_evalue(&p->arr[i],ref);
65
66 /* Then p itself */
67 switch (p->type) {
68 case polynomial:
69 case periodic:
70 case evector:
71 p->pos = ref[p->pos-1]+1;
72 }
73 return;
74 } /* aep_evalue */
75
76 /** Comments */
77 void addeliminatedparams_evalue(evalue *e,Matrix *CT) {
78
79 enode *p;
80 int i, j;
81 int *ref;
82
83 if (value_notzero_p(e->d))
84 return; /* a rational number, its already reduced */
85 if(!(p = e->x.p))
86 return; /* hum... an overflow probably occured */
87
88 /* Compute ref */
89 ref = (int *)malloc(sizeof(int)*(CT->NbRows-1));
90 for(i=0;i<CT->NbRows-1;i++)
91 for(j=0;j<CT->NbColumns;j++)
92 if(value_notzero_p(CT->p[i][j])) {
93 ref[i] = j;
94 break;
95 }
96
97 /* Transform the references in e, using ref */
98 aep_evalue(e,ref);
99 free( ref );
100 return;
101 } /* addeliminatedparams_evalue */
102
103 static void addeliminatedparams_partition(enode *p, Matrix *CT, Polyhedron *CEq,
104 unsigned nparam, unsigned MaxRays)
105 {
106 int i;
107 assert(p->type == partition);
108 p->pos = nparam;
109
110 for (i = 0; i < p->size/2; i++) {
111 Polyhedron *D = EVALUE_DOMAIN(p->arr[2*i]);
112 Polyhedron *T = DomainPreimage(D, CT, MaxRays);
113 Domain_Free(D);
114 if (CEq) {
115 D = T;
116 T = DomainIntersection(D, CEq, MaxRays);
117 Domain_Free(D);
118 }
119 EVALUE_SET_DOMAIN(p->arr[2*i], T);
120 }
121 }
122
123 void addeliminatedparams_enum(evalue *e, Matrix *CT, Polyhedron *CEq,
124 unsigned MaxRays, unsigned nparam)
125 {
126 enode *p;
127 int i;
128
129 if (CT->NbRows == CT->NbColumns)
130 return;
131
132 if (EVALUE_IS_ZERO(*e))
133 return;
134
135 if (value_notzero_p(e->d)) {
136 evalue res;
137 value_init(res.d);
138 value_set_si(res.d, 0);
139 res.x.p = new_enode(partition, 2, nparam);
140 EVALUE_SET_DOMAIN(res.x.p->arr[0],
141 DomainConstraintSimplify(Polyhedron_Copy(CEq), MaxRays));
142 value_clear(res.x.p->arr[1].d);
143 res.x.p->arr[1] = *e;
144 *e = res;
145 return;
146 }
147
148 p = e->x.p;
149 assert(p);
150
151 addeliminatedparams_partition(p, CT, CEq, nparam, MaxRays);
152 for (i = 0; i < p->size/2; i++)
153 addeliminatedparams_evalue(&p->arr[2*i+1], CT);
154 }
155
156 static int mod_rational_smaller(evalue *e1, evalue *e2)
157 {
158 int r;
159 Value m;
160 Value m2;
161 value_init(m);
162 value_init(m2);
163
164 assert(value_notzero_p(e1->d));
165 assert(value_notzero_p(e2->d));
166 value_multiply(m, e1->x.n, e2->d);
167 value_multiply(m2, e2->x.n, e1->d);
168 if (value_lt(m, m2))
169 r = 1;
170 else if (value_gt(m, m2))
171 r = 0;
172 else
173 r = -1;
174 value_clear(m);
175 value_clear(m2);
176
177 return r;
178 }
179
180 static int mod_term_smaller_r(evalue *e1, evalue *e2)
181 {
182 if (value_notzero_p(e1->d)) {
183 int r;
184 if (value_zero_p(e2->d))
185 return 1;
186 r = mod_rational_smaller(e1, e2);
187 return r == -1 ? 0 : r;
188 }
189 if (value_notzero_p(e2->d))
190 return 0;
191 if (e1->x.p->pos < e2->x.p->pos)
192 return 1;
193 else if (e1->x.p->pos > e2->x.p->pos)
194 return 0;
195 else {
196 int r = mod_rational_smaller(&e1->x.p->arr[1], &e2->x.p->arr[1]);
197 return r == -1
198 ? mod_term_smaller_r(&e1->x.p->arr[0], &e2->x.p->arr[0])
199 : r;
200 }
201 }
202
203 static int mod_term_smaller(const evalue *e1, const evalue *e2)
204 {
205 assert(value_zero_p(e1->d));
206 assert(value_zero_p(e2->d));
207 assert(e1->x.p->type == fractional || e1->x.p->type == flooring);
208 assert(e2->x.p->type == fractional || e2->x.p->type == flooring);
209 return mod_term_smaller_r(&e1->x.p->arr[0], &e2->x.p->arr[0]);
210 }
211
212 static void check_order(const evalue *e)
213 {
214 int i;
215 evalue *a;
216
217 if (value_notzero_p(e->d))
218 return;
219
220 switch (e->x.p->type) {
221 case partition:
222 for (i = 0; i < e->x.p->size/2; ++i)
223 check_order(&e->x.p->arr[2*i+1]);
224 break;
225 case relation:
226 for (i = 1; i < e->x.p->size; ++i) {
227 a = &e->x.p->arr[i];
228 if (value_notzero_p(a->d))
229 continue;
230 switch (a->x.p->type) {
231 case relation:
232 assert(mod_term_smaller(&e->x.p->arr[0], &a->x.p->arr[0]));
233 break;
234 case partition:
235 assert(0);
236 }
237 check_order(a);
238 }
239 break;
240 case polynomial:
241 for (i = 0; i < e->x.p->size; ++i) {
242 a = &e->x.p->arr[i];
243 if (value_notzero_p(a->d))
244 continue;
245 switch (a->x.p->type) {
246 case polynomial:
247 assert(e->x.p->pos < a->x.p->pos);
248 break;
249 case relation:
250 case partition:
251 assert(0);
252 }
253 check_order(a);
254 }
255 break;
256 case fractional:
257 case flooring:
258 for (i = 1; i < e->x.p->size; ++i) {
259 a = &e->x.p->arr[i];
260 if (value_notzero_p(a->d))
261 continue;
262 switch (a->x.p->type) {
263 case polynomial:
264 case relation:
265 case partition:
266 assert(0);
267 }
268 }
269 break;
270 }
271 }
272
273 /* Negative pos means inequality */
274 /* s is negative of substitution if m is not zero */
275 struct fixed_param {
276 int pos;
277 evalue s;
278 Value d;
279 Value m;
280 };
281
282 struct subst {
283 struct fixed_param *fixed;
284 int n;
285 int max;
286 };
287
288 static int relations_depth(evalue *e)
289 {
290 int d;
291
292 for (d = 0;
293 value_zero_p(e->d) && e->x.p->type == relation;
294 e = &e->x.p->arr[1], ++d);
295 return d;
296 }
297
298 static void poly_denom_not_constant(evalue **pp, Value *d)
299 {
300 evalue *p = *pp;
301 value_set_si(*d, 1);
302
303 while (value_zero_p(p->d)) {
304 assert(p->x.p->type == polynomial);
305 assert(p->x.p->size == 2);
306 assert(value_notzero_p(p->x.p->arr[1].d));
307 value_lcm(*d, *d, p->x.p->arr[1].d);
308 p = &p->x.p->arr[0];
309 }
310 *pp = p;
311 }
312
313 static void poly_denom(evalue *p, Value *d)
314 {
315 poly_denom_not_constant(&p, d);
316 value_lcm(*d, *d, p->d);
317 }
318
319 static void realloc_substitution(struct subst *s, int d)
320 {
321 struct fixed_param *n;
322 int i;
323 NALLOC(n, s->max+d);
324 for (i = 0; i < s->n; ++i)
325 n[i] = s->fixed[i];
326 free(s->fixed);
327 s->fixed = n;
328 s->max += d;
329 }
330
331 static int add_modulo_substitution(struct subst *s, evalue *r)
332 {
333 evalue *p;
334 evalue *f;
335 evalue *m;
336
337 assert(value_zero_p(r->d) && r->x.p->type == relation);
338 m = &r->x.p->arr[0];
339
340 /* May have been reduced already */
341 if (value_notzero_p(m->d))
342 return 0;
343
344 assert(value_zero_p(m->d) && m->x.p->type == fractional);
345 assert(m->x.p->size == 3);
346
347 /* fractional was inverted during reduction
348 * invert it back and move constant in
349 */
350 if (!EVALUE_IS_ONE(m->x.p->arr[2])) {
351 assert(value_pos_p(m->x.p->arr[2].d));
352 assert(value_mone_p(m->x.p->arr[2].x.n));
353 value_set_si(m->x.p->arr[2].x.n, 1);
354 value_increment(m->x.p->arr[1].x.n, m->x.p->arr[1].x.n);
355 assert(value_eq(m->x.p->arr[1].x.n, m->x.p->arr[1].d));
356 value_set_si(m->x.p->arr[1].x.n, 1);
357 eadd(&m->x.p->arr[1], &m->x.p->arr[0]);
358 value_set_si(m->x.p->arr[1].x.n, 0);
359 value_set_si(m->x.p->arr[1].d, 1);
360 }
361
362 /* Oops. Nested identical relations. */
363 if (!EVALUE_IS_ZERO(m->x.p->arr[1]))
364 return 0;
365
366 if (s->n >= s->max) {
367 int d = relations_depth(r);
368 realloc_substitution(s, d);
369 }
370
371 p = &m->x.p->arr[0];
372 assert(value_zero_p(p->d) && p->x.p->type == polynomial);
373 assert(p->x.p->size == 2);
374 f = &p->x.p->arr[1];
375
376 assert(value_pos_p(f->x.n));
377
378 value_init(s->fixed[s->n].m);
379 value_assign(s->fixed[s->n].m, f->d);
380 s->fixed[s->n].pos = p->x.p->pos;
381 value_init(s->fixed[s->n].d);
382 value_assign(s->fixed[s->n].d, f->x.n);
383 value_init(s->fixed[s->n].s.d);
384 evalue_copy(&s->fixed[s->n].s, &p->x.p->arr[0]);
385 ++s->n;
386
387 return 1;
388 }
389
390 static int type_offset(enode *p)
391 {
392 return p->type == fractional ? 1 :
393 p->type == flooring ? 1 : 0;
394 }
395
396 static void reorder_terms_about(enode *p, evalue *v)
397 {
398 int i;
399 int offset = type_offset(p);
400
401 for (i = p->size-1; i >= offset+1; i--) {
402 emul(v, &p->arr[i]);
403 eadd(&p->arr[i], &p->arr[i-1]);
404 free_evalue_refs(&(p->arr[i]));
405 }
406 p->size = offset+1;
407 free_evalue_refs(v);
408 }
409
410 static void reorder_terms(evalue *e)
411 {
412 enode *p;
413 evalue f;
414
415 assert(value_zero_p(e->d));
416 p = e->x.p;
417 assert(p->type == fractional); /* for now */
418
419 value_init(f.d);
420 value_set_si(f.d, 0);
421 f.x.p = new_enode(fractional, 3, -1);
422 value_clear(f.x.p->arr[0].d);
423 f.x.p->arr[0] = p->arr[0];
424 evalue_set_si(&f.x.p->arr[1], 0, 1);
425 evalue_set_si(&f.x.p->arr[2], 1, 1);
426 reorder_terms_about(p, &f);
427 value_clear(e->d);
428 *e = p->arr[1];
429 free(p);
430 }
431
432 void _reduce_evalue (evalue *e, struct subst *s, int fract) {
433
434 enode *p;
435 int i, j, k;
436 int add;
437
438 if (value_notzero_p(e->d)) {
439 if (fract)
440 mpz_fdiv_r(e->x.n, e->x.n, e->d);
441 return; /* a rational number, its already reduced */
442 }
443
444 if(!(p = e->x.p))
445 return; /* hum... an overflow probably occured */
446
447 /* First reduce the components of p */
448 add = p->type == relation;
449 for (i=0; i<p->size; i++) {
450 if (add && i == 1)
451 add = add_modulo_substitution(s, e);
452
453 if (i == 0 && p->type==fractional)
454 _reduce_evalue(&p->arr[i], s, 1);
455 else
456 _reduce_evalue(&p->arr[i], s, fract);
457
458 if (add && i == p->size-1) {
459 --s->n;
460 value_clear(s->fixed[s->n].m);
461 value_clear(s->fixed[s->n].d);
462 free_evalue_refs(&s->fixed[s->n].s);
463 } else if (add && i == 1)
464 s->fixed[s->n-1].pos *= -1;
465 }
466
467 if (p->type==periodic) {
468
469 /* Try to reduce the period */
470 for (i=1; i<=(p->size)/2; i++) {
471 if ((p->size % i)==0) {
472
473 /* Can we reduce the size to i ? */
474 for (j=0; j<i; j++)
475 for (k=j+i; k<e->x.p->size; k+=i)
476 if (!eequal(&p->arr[j], &p->arr[k])) goto you_lose;
477
478 /* OK, lets do it */
479 for (j=i; j<p->size; j++) free_evalue_refs(&p->arr[j]);
480 p->size = i;
481 break;
482
483 you_lose: /* OK, lets not do it */
484 continue;
485 }
486 }
487
488 /* Try to reduce its strength */
489 if (p->size == 1) {
490 value_clear(e->d);
491 memcpy(e,&p->arr[0],sizeof(evalue));
492 free(p);
493 }
494 }
495 else if (p->type==polynomial) {
496 for (k = 0; s && k < s->n; ++k) {
497 if (s->fixed[k].pos == p->pos) {
498 int divide = value_notone_p(s->fixed[k].d);
499 evalue d;
500
501 if (value_notzero_p(s->fixed[k].m)) {
502 if (!fract)
503 continue;
504 assert(p->size == 2);
505 if (divide && value_ne(s->fixed[k].d, p->arr[1].x.n))
506 continue;
507 if (!mpz_divisible_p(s->fixed[k].m, p->arr[1].d))
508 continue;
509 divide = 1;
510 }
511
512 if (divide) {
513 value_init(d.d);
514 value_assign(d.d, s->fixed[k].d);
515 value_init(d.x.n);
516 if (value_notzero_p(s->fixed[k].m))
517 value_oppose(d.x.n, s->fixed[k].m);
518 else
519 value_set_si(d.x.n, 1);
520 }
521
522 for (i=p->size-1;i>=1;i--) {
523 emul(&s->fixed[k].s, &p->arr[i]);
524 if (divide)
525 emul(&d, &p->arr[i]);
526 eadd(&p->arr[i], &p->arr[i-1]);
527 free_evalue_refs(&(p->arr[i]));
528 }
529 p->size = 1;
530 _reduce_evalue(&p->arr[0], s, fract);
531
532 if (divide)
533 free_evalue_refs(&d);
534
535 break;
536 }
537 }
538
539 /* Try to reduce the degree */
540 for (i=p->size-1;i>=1;i--) {
541 if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
542 break;
543 /* Zero coefficient */
544 free_evalue_refs(&(p->arr[i]));
545 }
546 if (i+1<p->size)
547 p->size = i+1;
548
549 /* Try to reduce its strength */
550 if (p->size == 1) {
551 value_clear(e->d);
552 memcpy(e,&p->arr[0],sizeof(evalue));
553 free(p);
554 }
555 }
556 else if (p->type==fractional) {
557 int reorder = 0;
558 evalue v;
559
560 if (value_notzero_p(p->arr[0].d)) {
561 value_init(v.d);
562 value_assign(v.d, p->arr[0].d);
563 value_init(v.x.n);
564 mpz_fdiv_r(v.x.n, p->arr[0].x.n, p->arr[0].d);
565
566 reorder = 1;
567 } else {
568 evalue *f, *base;
569 evalue *pp = &p->arr[0];
570 assert(value_zero_p(pp->d) && pp->x.p->type == polynomial);
571 assert(pp->x.p->size == 2);
572
573 /* search for exact duplicate among the modulo inequalities */
574 do {
575 f = &pp->x.p->arr[1];
576 for (k = 0; s && k < s->n; ++k) {
577 if (-s->fixed[k].pos == pp->x.p->pos &&
578 value_eq(s->fixed[k].d, f->x.n) &&
579 value_eq(s->fixed[k].m, f->d) &&
580 eequal(&s->fixed[k].s, &pp->x.p->arr[0]))
581 break;
582 }
583 if (k < s->n) {
584 Value g;
585 value_init(g);
586
587 /* replace { E/m } by { (E-1)/m } + 1/m */
588 poly_denom(pp, &g);
589 if (reorder) {
590 evalue extra;
591 value_init(extra.d);
592 evalue_set_si(&extra, 1, 1);
593 value_assign(extra.d, g);
594 eadd(&extra, &v.x.p->arr[1]);
595 free_evalue_refs(&extra);
596
597 /* We've been going in circles; stop now */
598 if (value_ge(v.x.p->arr[1].x.n, v.x.p->arr[1].d)) {
599 free_evalue_refs(&v);
600 value_init(v.d);
601 evalue_set_si(&v, 0, 1);
602 break;
603 }
604 } else {
605 value_init(v.d);
606 value_set_si(v.d, 0);
607 v.x.p = new_enode(fractional, 3, -1);
608 evalue_set_si(&v.x.p->arr[1], 1, 1);
609 value_assign(v.x.p->arr[1].d, g);
610 evalue_set_si(&v.x.p->arr[2], 1, 1);
611 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
612 }
613
614 for (f = &v.x.p->arr[0]; value_zero_p(f->d);
615 f = &f->x.p->arr[0])
616 ;
617 value_division(f->d, g, f->d);
618 value_multiply(f->x.n, f->x.n, f->d);
619 value_assign(f->d, g);
620 value_decrement(f->x.n, f->x.n);
621 mpz_fdiv_r(f->x.n, f->x.n, f->d);
622
623 value_gcd(g, f->d, f->x.n);
624 value_division(f->d, f->d, g);
625 value_division(f->x.n, f->x.n, g);
626
627 value_clear(g);
628 pp = &v.x.p->arr[0];
629
630 reorder = 1;
631 }
632 } while (k < s->n);
633
634 /* reduction may have made this fractional arg smaller */
635 i = reorder ? p->size : 1;
636 for ( ; i < p->size; ++i)
637 if (value_zero_p(p->arr[i].d) &&
638 p->arr[i].x.p->type == fractional &&
639 !mod_term_smaller(e, &p->arr[i]))
640 break;
641 if (i < p->size) {
642 value_init(v.d);
643 value_set_si(v.d, 0);
644 v.x.p = new_enode(fractional, 3, -1);
645 evalue_set_si(&v.x.p->arr[1], 0, 1);
646 evalue_set_si(&v.x.p->arr[2], 1, 1);
647 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
648
649 reorder = 1;
650 }
651
652 if (!reorder) {
653 Value m;
654 Value r;
655 evalue *pp = &p->arr[0];
656 value_init(m);
657 value_init(r);
658 poly_denom_not_constant(&pp, &m);
659 mpz_fdiv_r(r, m, pp->d);
660 if (value_notzero_p(r)) {
661 value_init(v.d);
662 value_set_si(v.d, 0);
663 v.x.p = new_enode(fractional, 3, -1);
664
665 value_multiply(r, m, pp->x.n);
666 value_multiply(v.x.p->arr[1].d, m, pp->d);
667 value_init(v.x.p->arr[1].x.n);
668 mpz_fdiv_r(v.x.p->arr[1].x.n, r, pp->d);
669 mpz_fdiv_q(r, r, pp->d);
670
671 evalue_set_si(&v.x.p->arr[2], 1, 1);
672 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
673 pp = &v.x.p->arr[0];
674 while (value_zero_p(pp->d))
675 pp = &pp->x.p->arr[0];
676
677 value_assign(pp->d, m);
678 value_assign(pp->x.n, r);
679
680 value_gcd(r, pp->d, pp->x.n);
681 value_division(pp->d, pp->d, r);
682 value_division(pp->x.n, pp->x.n, r);
683
684 reorder = 1;
685 }
686 value_clear(m);
687 value_clear(r);
688 }
689
690 if (!reorder) {
691 int invert = 0;
692 Value twice;
693 value_init(twice);
694
695 for (pp = &p->arr[0]; value_zero_p(pp->d);
696 pp = &pp->x.p->arr[0]) {
697 f = &pp->x.p->arr[1];
698 assert(value_pos_p(f->d));
699 mpz_mul_ui(twice, f->x.n, 2);
700 if (value_lt(twice, f->d))
701 break;
702 if (value_eq(twice, f->d))
703 continue;
704 invert = 1;
705 break;
706 }
707
708 if (invert) {
709 value_init(v.d);
710 value_set_si(v.d, 0);
711 v.x.p = new_enode(fractional, 3, -1);
712 evalue_set_si(&v.x.p->arr[1], 0, 1);
713 poly_denom(&p->arr[0], &twice);
714 value_assign(v.x.p->arr[1].d, twice);
715 value_decrement(v.x.p->arr[1].x.n, twice);
716 evalue_set_si(&v.x.p->arr[2], -1, 1);
717 evalue_copy(&v.x.p->arr[0], &p->arr[0]);
718
719 for (pp = &v.x.p->arr[0]; value_zero_p(pp->d);
720 pp = &pp->x.p->arr[0]) {
721 f = &pp->x.p->arr[1];
722 value_oppose(f->x.n, f->x.n);
723 mpz_fdiv_r(f->x.n, f->x.n, f->d);
724 }
725 value_division(pp->d, twice, pp->d);
726 value_multiply(pp->x.n, pp->x.n, pp->d);
727 value_assign(pp->d, twice);
728 value_oppose(pp->x.n, pp->x.n);
729 value_decrement(pp->x.n, pp->x.n);
730 mpz_fdiv_r(pp->x.n, pp->x.n, pp->d);
731
732 /* Maybe we should do this during reduction of
733 * the constant.
734 */
735 value_gcd(twice, pp->d, pp->x.n);
736 value_division(pp->d, pp->d, twice);
737 value_division(pp->x.n, pp->x.n, twice);
738
739 reorder = 1;
740 }
741
742 value_clear(twice);
743 }
744 }
745
746 if (reorder) {
747 reorder_terms_about(p, &v);
748 _reduce_evalue(&p->arr[1], s, fract);
749 }
750
751 /* Try to reduce the degree */
752 for (i=p->size-1;i>=2;i--) {
753 if (!(value_one_p(p->arr[i].d) && value_zero_p(p->arr[i].x.n)))
754 break;
755 /* Zero coefficient */
756 free_evalue_refs(&(p->arr[i]));
757 }
758 if (i+1<p->size)
759 p->size = i+1;
760
761 /* Try to reduce its strength */
762 if (p->size == 2) {
763 value_clear(e->d);
764 memcpy(e,&p->arr[1],sizeof(evalue));
765 free_evalue_refs(&(p->arr[0]));
766 free(p);
767 }
768 }
769 else if (p->type == flooring) {
770 /* Try to reduce the degree */
771 for (i=p->size-1;i>=2;i--) {
772 if (!EVALUE_IS_ZERO(p->arr[i]))
773 break;
774 /* Zero coefficient */
775 free_evalue_refs(&(p->arr[i]));
776 }
777 if (i+1<p->size)
778 p->size = i+1;
779
780 /* Try to reduce its strength */
781 if (p->size == 2) {
782 value_clear(e->d);
783 memcpy(e,&p->arr[1],sizeof(evalue));
784 free_evalue_refs(&(p->arr[0]));
785 free(p);
786 }
787 }
788 else if (p->type == relation) {
789 if (p->size == 3 && eequal(&p->arr[1], &p->arr[2])) {
790 free_evalue_refs(&(p->arr[2]));
791 free_evalue_refs(&(p->arr[0]));
792 p->size = 2;
793 value_clear(e->d);
794 *e = p->arr[1];
795 free(p);
796 return;
797 }
798 if (p->size == 3 && EVALUE_IS_ZERO(p->arr[2])) {
799 free_evalue_refs(&(p->arr[2]));
800 p->size = 2;
801 }
802 if (p->size == 2 && EVALUE_IS_ZERO(p->arr[1])) {
803 free_evalue_refs(&(p->arr[1]));
804 free_evalue_refs(&(p->arr[0]));
805 evalue_set_si(e, 0, 1);
806 free(p);
807 } else {
808 int reduced = 0;
809 evalue *m;
810 m = &p->arr[0];
811
812 /* Relation was reduced by means of an identical
813 * inequality => remove
814 */
815 if (value_zero_p(m->d) && !EVALUE_IS_ZERO(m->x.p->arr[1]))
816 reduced = 1;
817
818 if (reduced || value_notzero_p(p->arr[0].d)) {
819 if (!reduced && value_zero_p(p->arr[0].x.n)) {
820 value_clear(e->d);
821 memcpy(e,&p->arr[1],sizeof(evalue));
822 if (p->size == 3)
823 free_evalue_refs(&(p->arr[2]));
824 } else {
825 if (p->size == 3) {
826 value_clear(e->d);
827 memcpy(e,&p->arr[2],sizeof(evalue));
828 } else
829 evalue_set_si(e, 0, 1);
830 free_evalue_refs(&(p->arr[1]));
831 }
832 free_evalue_refs(&(p->arr[0]));
833 free(p);
834 }
835 }
836 }
837 return;
838 } /* reduce_evalue */
839
840 static void add_substitution(struct subst *s, Value *row, unsigned dim)
841 {
842 int k, l;
843 evalue *r;
844
845 for (k = 0; k < dim; ++k)
846 if (value_notzero_p(row[k+1]))
847 break;
848
849 Vector_Normalize_Positive(row+1, dim+1, k);
850 assert(s->n < s->max);
851 value_init(s->fixed[s->n].d);
852 value_init(s->fixed[s->n].m);
853 value_assign(s->fixed[s->n].d, row[k+1]);
854 s->fixed[s->n].pos = k+1;
855 value_set_si(s->fixed[s->n].m, 0);
856 r = &s->fixed[s->n].s;
857 value_init(r->d);
858 for (l = k+1; l < dim; ++l)
859 if (value_notzero_p(row[l+1])) {
860 value_set_si(r->d, 0);
861 r->x.p = new_enode(polynomial, 2, l + 1);
862 value_init(r->x.p->arr[1].x.n);
863 value_oppose(r->x.p->arr[1].x.n, row[l+1]);
864 value_set_si(r->x.p->arr[1].d, 1);
865 r = &r->x.p->arr[0];
866 }
867 value_init(r->x.n);
868 value_oppose(r->x.n, row[dim+1]);
869 value_set_si(r->d, 1);
870 ++s->n;
871 }
872
873 static void _reduce_evalue_in_domain(evalue *e, Polyhedron *D, struct subst *s)
874 {
875 unsigned dim;
876 Polyhedron *orig = D;
877
878 s->n = 0;
879 dim = D->Dimension;
880 if (D->next)
881 D = DomainConvex(D, 0);
882 if (!D->next && D->NbEq) {
883 int j, k;
884 if (s->max < dim) {
885 if (s->max != 0)
886 realloc_substitution(s, dim);
887 else {
888 int d = relations_depth(e);
889 s->max = dim+d;
890 NALLOC(s->fixed, s->max);
891 }
892 }
893 for (j = 0; j < D->NbEq; ++j)
894 add_substitution(s, D->Constraint[j], dim);
895 }
896 if (D != orig)
897 Domain_Free(D);
898 _reduce_evalue(e, s, 0);
899 if (s->n != 0) {
900 int j;
901 for (j = 0; j < s->n; ++j) {
902 value_clear(s->fixed[j].d);
903 value_clear(s->fixed[j].m);
904 free_evalue_refs(&s->fixed[j].s);
905 }
906 }
907 }
908
909 void reduce_evalue_in_domain(evalue *e, Polyhedron *D)
910 {
911 struct subst s = { NULL, 0, 0 };
912 if (emptyQ2(D)) {
913 if (EVALUE_IS_ZERO(*e))
914 return;
915 free_evalue_refs(e);
916 value_init(e->d);
917 evalue_set_si(e, 0, 1);
918 return;
919 }
920 _reduce_evalue_in_domain(e, D, &s);
921 if (s.max != 0)
922 free(s.fixed);
923 }
924
925 void reduce_evalue (evalue *e) {
926 struct subst s = { NULL, 0, 0 };
927
928 if (value_notzero_p(e->d))
929 return; /* a rational number, its already reduced */
930
931 if (e->x.p->type == partition) {
932 int i;
933 unsigned dim = -1;
934 for (i = 0; i < e->x.p->size/2; ++i) {
935 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
936
937 /* This shouldn't really happen;
938 * Empty domains should not be added.
939 */
940 POL_ENSURE_VERTICES(D);
941 if (!emptyQ(D))
942 _reduce_evalue_in_domain(&e->x.p->arr[2*i+1], D, &s);
943
944 if (emptyQ(D) || EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
945 free_evalue_refs(&e->x.p->arr[2*i+1]);
946 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
947 value_clear(e->x.p->arr[2*i].d);
948 e->x.p->size -= 2;
949 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
950 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
951 --i;
952 }
953 }
954 if (e->x.p->size == 0) {
955 free(e->x.p);
956 evalue_set_si(e, 0, 1);
957 }
958 } else
959 _reduce_evalue(e, &s, 0);
960 if (s.max != 0)
961 free(s.fixed);
962 }
963
964 static void print_evalue_r(FILE *DST, const evalue *e, const char *const *pname)
965 {
966 if(value_notzero_p(e->d)) {
967 if(value_notone_p(e->d)) {
968 value_print(DST,VALUE_FMT,e->x.n);
969 fprintf(DST,"/");
970 value_print(DST,VALUE_FMT,e->d);
971 }
972 else {
973 value_print(DST,VALUE_FMT,e->x.n);
974 }
975 }
976 else
977 print_enode(DST,e->x.p,pname);
978 return;
979 } /* print_evalue */
980
981 void print_evalue(FILE *DST, const evalue *e, const char * const *pname)
982 {
983 print_evalue_r(DST, e, pname);
984 if (value_notzero_p(e->d))
985 fprintf(DST, "\n");
986 }
987
988 void print_enode(FILE *DST, enode *p, const char *const *pname)
989 {
990 int i;
991
992 if (!p) {
993 fprintf(DST, "NULL");
994 return;
995 }
996 switch (p->type) {
997 case evector:
998 fprintf(DST, "{ ");
999 for (i=0; i<p->size; i++) {
1000 print_evalue_r(DST, &p->arr[i], pname);
1001 if (i!=(p->size-1))
1002 fprintf(DST, ", ");
1003 }
1004 fprintf(DST, " }\n");
1005 break;
1006 case polynomial:
1007 fprintf(DST, "( ");
1008 for (i=p->size-1; i>=0; i--) {
1009 print_evalue_r(DST, &p->arr[i], pname);
1010 if (i==1) fprintf(DST, " * %s + ", pname[p->pos-1]);
1011 else if (i>1)
1012 fprintf(DST, " * %s^%d + ", pname[p->pos-1], i);
1013 }
1014 fprintf(DST, " )\n");
1015 break;
1016 case periodic:
1017 fprintf(DST, "[ ");
1018 for (i=0; i<p->size; i++) {
1019 print_evalue_r(DST, &p->arr[i], pname);
1020 if (i!=(p->size-1)) fprintf(DST, ", ");
1021 }
1022 fprintf(DST," ]_%s", pname[p->pos-1]);
1023 break;
1024 case flooring:
1025 case fractional:
1026 fprintf(DST, "( ");
1027 for (i=p->size-1; i>=1; i--) {
1028 print_evalue_r(DST, &p->arr[i], pname);
1029 if (i >= 2) {
1030 fprintf(DST, " * ");
1031 fprintf(DST, p->type == flooring ? "[" : "{");
1032 print_evalue_r(DST, &p->arr[0], pname);
1033 fprintf(DST, p->type == flooring ? "]" : "}");
1034 if (i>2)
1035 fprintf(DST, "^%d + ", i-1);
1036 else
1037 fprintf(DST, " + ");
1038 }
1039 }
1040 fprintf(DST, " )\n");
1041 break;
1042 case relation:
1043 fprintf(DST, "[ ");
1044 print_evalue_r(DST, &p->arr[0], pname);
1045 fprintf(DST, "= 0 ] * \n");
1046 print_evalue_r(DST, &p->arr[1], pname);
1047 if (p->size > 2) {
1048 fprintf(DST, " +\n [ ");
1049 print_evalue_r(DST, &p->arr[0], pname);
1050 fprintf(DST, "!= 0 ] * \n");
1051 print_evalue_r(DST, &p->arr[2], pname);
1052 }
1053 break;
1054 case partition: {
1055 char **new_names = NULL;
1056 const char *const *names = pname;
1057 int maxdim = EVALUE_DOMAIN(p->arr[0])->Dimension;
1058 if (!pname || p->pos < maxdim) {
1059 new_names = ALLOCN(char *, maxdim);
1060 for (i = 0; i < p->pos; ++i) {
1061 if (pname)
1062 new_names[i] = (char *)pname[i];
1063 else {
1064 new_names[i] = ALLOCN(char, 10);
1065 snprintf(new_names[i], 10, "%c", 'P'+i);
1066 }
1067 }
1068 for ( ; i < maxdim; ++i) {
1069 new_names[i] = ALLOCN(char, 10);
1070 snprintf(new_names[i], 10, "_p%d", i);
1071 }
1072 names = (const char**)new_names;
1073 }
1074
1075 for (i=0; i<p->size/2; i++) {
1076 Print_Domain(DST, EVALUE_DOMAIN(p->arr[2*i]), names);
1077 print_evalue_r(DST, &p->arr[2*i+1], names);
1078 if (value_notzero_p(p->arr[2*i+1].d))
1079 fprintf(DST, "\n");
1080 }
1081
1082 if (!pname || p->pos < maxdim) {
1083 for (i = pname ? p->pos : 0; i < maxdim; ++i)
1084 free(new_names[i]);
1085 free(new_names);
1086 }
1087
1088 break;
1089 }
1090 default:
1091 assert(0);
1092 }
1093 return;
1094 } /* print_enode */
1095
1096 static void eadd_rev(const evalue *e1, evalue *res)
1097 {
1098 evalue ev;
1099 value_init(ev.d);
1100 evalue_copy(&ev, e1);
1101 eadd(res, &ev);
1102 free_evalue_refs(res);
1103 *res = ev;
1104 }
1105
1106 static void eadd_rev_cst(const evalue *e1, evalue *res)
1107 {
1108 evalue ev;
1109 value_init(ev.d);
1110 evalue_copy(&ev, e1);
1111 eadd(res, &ev.x.p->arr[type_offset(ev.x.p)]);
1112 free_evalue_refs(res);
1113 *res = ev;
1114 }
1115
1116 static int is_zero_on(evalue *e, Polyhedron *D)
1117 {
1118 int is_zero;
1119 evalue tmp;
1120 value_init(tmp.d);
1121 tmp.x.p = new_enode(partition, 2, D->Dimension);
1122 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Domain_Copy(D));
1123 evalue_copy(&tmp.x.p->arr[1], e);
1124 reduce_evalue(&tmp);
1125 is_zero = EVALUE_IS_ZERO(tmp);
1126 free_evalue_refs(&tmp);
1127 return is_zero;
1128 }
1129
1130 struct section { Polyhedron * D; evalue E; };
1131
1132 void eadd_partitions(const evalue *e1, evalue *res)
1133 {
1134 int n, i, j;
1135 Polyhedron *d, *fd;
1136 struct section *s;
1137 s = (struct section *)
1138 malloc((e1->x.p->size/2+1) * (res->x.p->size/2+1) *
1139 sizeof(struct section));
1140 assert(s);
1141 assert(e1->x.p->pos == res->x.p->pos);
1142 assert(e1->x.p->pos == EVALUE_DOMAIN(e1->x.p->arr[0])->Dimension);
1143 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1144
1145 n = 0;
1146 for (j = 0; j < e1->x.p->size/2; ++j) {
1147 assert(res->x.p->size >= 2);
1148 fd = DomainDifference(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1149 EVALUE_DOMAIN(res->x.p->arr[0]), 0);
1150 if (!emptyQ(fd))
1151 for (i = 1; i < res->x.p->size/2; ++i) {
1152 Polyhedron *t = fd;
1153 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1154 Domain_Free(t);
1155 if (emptyQ(fd))
1156 break;
1157 }
1158 if (emptyQ(fd)) {
1159 Domain_Free(fd);
1160 continue;
1161 }
1162 /* See if we can extend one of the domains in res to cover fd */
1163 for (i = 0; i < res->x.p->size/2; ++i)
1164 if (is_zero_on(&res->x.p->arr[2*i+1], fd))
1165 break;
1166 if (i < res->x.p->size/2) {
1167 EVALUE_SET_DOMAIN(res->x.p->arr[2*i],
1168 DomainConcat(fd, EVALUE_DOMAIN(res->x.p->arr[2*i])));
1169 continue;
1170 }
1171 value_init(s[n].E.d);
1172 evalue_copy(&s[n].E, &e1->x.p->arr[2*j+1]);
1173 s[n].D = fd;
1174 ++n;
1175 }
1176 for (i = 0; i < res->x.p->size/2; ++i) {
1177 fd = EVALUE_DOMAIN(res->x.p->arr[2*i]);
1178 for (j = 0; j < e1->x.p->size/2; ++j) {
1179 Polyhedron *t;
1180 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1181 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1182 if (emptyQ(d)) {
1183 Domain_Free(d);
1184 continue;
1185 }
1186 t = fd;
1187 fd = DomainDifference(fd, EVALUE_DOMAIN(e1->x.p->arr[2*j]), 0);
1188 if (t != EVALUE_DOMAIN(res->x.p->arr[2*i]))
1189 Domain_Free(t);
1190 value_init(s[n].E.d);
1191 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
1192 eadd(&e1->x.p->arr[2*j+1], &s[n].E);
1193 if (!emptyQ(fd) && is_zero_on(&e1->x.p->arr[2*j+1], fd)) {
1194 d = DomainConcat(fd, d);
1195 fd = Empty_Polyhedron(fd->Dimension);
1196 }
1197 s[n].D = d;
1198 ++n;
1199 }
1200 if (!emptyQ(fd)) {
1201 s[n].E = res->x.p->arr[2*i+1];
1202 s[n].D = fd;
1203 ++n;
1204 } else {
1205 free_evalue_refs(&res->x.p->arr[2*i+1]);
1206 Domain_Free(fd);
1207 }
1208 if (fd != EVALUE_DOMAIN(res->x.p->arr[2*i]))
1209 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
1210 value_clear(res->x.p->arr[2*i].d);
1211 }
1212
1213 free(res->x.p);
1214 assert(n > 0);
1215 res->x.p = new_enode(partition, 2*n, e1->x.p->pos);
1216 for (j = 0; j < n; ++j) {
1217 s[j].D = DomainConstraintSimplify(s[j].D, 0);
1218 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1219 value_clear(res->x.p->arr[2*j+1].d);
1220 res->x.p->arr[2*j+1] = s[j].E;
1221 }
1222
1223 free(s);
1224 }
1225
1226 static void explicit_complement(evalue *res)
1227 {
1228 enode *rel = new_enode(relation, 3, 0);
1229 assert(rel);
1230 value_clear(rel->arr[0].d);
1231 rel->arr[0] = res->x.p->arr[0];
1232 value_clear(rel->arr[1].d);
1233 rel->arr[1] = res->x.p->arr[1];
1234 value_set_si(rel->arr[2].d, 1);
1235 value_init(rel->arr[2].x.n);
1236 value_set_si(rel->arr[2].x.n, 0);
1237 free(res->x.p);
1238 res->x.p = rel;
1239 }
1240
1241 void eadd(const evalue *e1, evalue *res)
1242 {
1243 int i;
1244
1245 if (EVALUE_IS_ZERO(*e1))
1246 return;
1247
1248 if (EVALUE_IS_ZERO(*res)) {
1249 if (value_notzero_p(e1->d)) {
1250 value_assign(res->d, e1->d);
1251 value_assign(res->x.n, e1->x.n);
1252 } else {
1253 value_clear(res->x.n);
1254 value_set_si(res->d, 0);
1255 res->x.p = ecopy(e1->x.p);
1256 }
1257 return;
1258 }
1259
1260 if (value_notzero_p(e1->d) && value_notzero_p(res->d)) {
1261 /* Add two rational numbers */
1262 Value g;
1263 value_init(g);
1264
1265 if (value_eq(e1->d, res->d))
1266 value_addto(res->x.n, res->x.n, e1->x.n);
1267 else {
1268 value_multiply(res->x.n, res->x.n, e1->d);
1269 value_addmul(res->x.n, e1->x.n, res->d);
1270 value_multiply(res->d,e1->d,res->d);
1271 }
1272
1273 value_gcd(g, res->x.n, res->d);
1274 if (value_notone_p(g)) {
1275 value_division(res->d,res->d,g);
1276 value_division(res->x.n,res->x.n,g);
1277 }
1278 value_clear(g);
1279 return ;
1280 }
1281 else if (value_notzero_p(e1->d) && value_zero_p(res->d)) {
1282 switch (res->x.p->type) {
1283 case polynomial:
1284 /* Add the constant to the constant term of a polynomial*/
1285 eadd(e1, &res->x.p->arr[0]);
1286 return ;
1287 case periodic:
1288 /* Add the constant to all elements of a periodic number */
1289 for (i=0; i<res->x.p->size; i++) {
1290 eadd(e1, &res->x.p->arr[i]);
1291 }
1292 return ;
1293 case evector:
1294 fprintf(stderr, "eadd: cannot add const with vector\n");
1295 return;
1296 case flooring:
1297 case fractional:
1298 eadd(e1, &res->x.p->arr[1]);
1299 return ;
1300 case partition:
1301 assert(EVALUE_IS_ZERO(*e1));
1302 break; /* Do nothing */
1303 case relation:
1304 /* Create (zero) complement if needed */
1305 if (res->x.p->size < 3 && !EVALUE_IS_ZERO(*e1))
1306 explicit_complement(res);
1307 for (i = 1; i < res->x.p->size; ++i)
1308 eadd(e1, &res->x.p->arr[i]);
1309 break;
1310 default:
1311 assert(0);
1312 }
1313 }
1314 /* add polynomial or periodic to constant
1315 * you have to exchange e1 and res, before doing addition */
1316
1317 else if (value_zero_p(e1->d) && value_notzero_p(res->d)) {
1318 eadd_rev(e1, res);
1319 return;
1320 }
1321 else { // ((e1->d==0) && (res->d==0))
1322 assert(!((e1->x.p->type == partition) ^
1323 (res->x.p->type == partition)));
1324 if (e1->x.p->type == partition) {
1325 eadd_partitions(e1, res);
1326 return;
1327 }
1328 if (e1->x.p->type == relation &&
1329 (res->x.p->type != relation ||
1330 mod_term_smaller(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1331 eadd_rev(e1, res);
1332 return;
1333 }
1334 if (res->x.p->type == relation) {
1335 if (e1->x.p->type == relation &&
1336 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1337 if (res->x.p->size < 3 && e1->x.p->size == 3)
1338 explicit_complement(res);
1339 for (i = 1; i < e1->x.p->size; ++i)
1340 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1341 return;
1342 }
1343 if (res->x.p->size < 3)
1344 explicit_complement(res);
1345 for (i = 1; i < res->x.p->size; ++i)
1346 eadd(e1, &res->x.p->arr[i]);
1347 return;
1348 }
1349 if ((e1->x.p->type != res->x.p->type) ) {
1350 /* adding to evalues of different type. two cases are possible
1351 * res is periodic and e1 is polynomial, you have to exchange
1352 * e1 and res then to add e1 to the constant term of res */
1353 if (e1->x.p->type == polynomial) {
1354 eadd_rev_cst(e1, res);
1355 }
1356 else if (res->x.p->type == polynomial) {
1357 /* res is polynomial and e1 is periodic,
1358 add e1 to the constant term of res */
1359
1360 eadd(e1,&res->x.p->arr[0]);
1361 } else
1362 assert(0);
1363
1364 return;
1365 }
1366 else if (e1->x.p->pos != res->x.p->pos ||
1367 ((res->x.p->type == fractional ||
1368 res->x.p->type == flooring) &&
1369 !eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1370 /* adding evalues of different position (i.e function of different unknowns
1371 * to case are possible */
1372
1373 switch (res->x.p->type) {
1374 case flooring:
1375 case fractional:
1376 if (mod_term_smaller(res, e1))
1377 eadd(e1,&res->x.p->arr[1]);
1378 else
1379 eadd_rev_cst(e1, res);
1380 return;
1381 case polynomial: // res and e1 are polynomials
1382 // add e1 to the constant term of res
1383
1384 if(res->x.p->pos < e1->x.p->pos)
1385 eadd(e1,&res->x.p->arr[0]);
1386 else
1387 eadd_rev_cst(e1, res);
1388 // value_clear(g); value_clear(m1); value_clear(m2);
1389 return;
1390 case periodic: // res and e1 are pointers to periodic numbers
1391 //add e1 to all elements of res
1392
1393 if(res->x.p->pos < e1->x.p->pos)
1394 for (i=0;i<res->x.p->size;i++) {
1395 eadd(e1,&res->x.p->arr[i]);
1396 }
1397 else
1398 eadd_rev(e1, res);
1399 return;
1400 default:
1401 assert(0);
1402 }
1403 }
1404
1405
1406 //same type , same pos and same size
1407 if (e1->x.p->size == res->x.p->size) {
1408 // add any element in e1 to the corresponding element in res
1409 i = type_offset(res->x.p);
1410 if (i == 1)
1411 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
1412 for (; i<res->x.p->size; i++) {
1413 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1414 }
1415 return ;
1416 }
1417
1418 /* Sizes are different */
1419 switch(res->x.p->type) {
1420 case polynomial:
1421 case flooring:
1422 case fractional:
1423 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1424 /* new enode and add res to that new node. If you do not do */
1425 /* that, you lose the the upper weight part of e1 ! */
1426
1427 if(e1->x.p->size > res->x.p->size)
1428 eadd_rev(e1, res);
1429 else {
1430 i = type_offset(res->x.p);
1431 if (i == 1)
1432 assert(eequal(&e1->x.p->arr[0],
1433 &res->x.p->arr[0]));
1434 for (; i<e1->x.p->size ; i++) {
1435 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1436 }
1437
1438 return ;
1439 }
1440 break;
1441
1442 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1443 case periodic:
1444 {
1445 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1446 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1447 to add periodicaly elements of e1 to elements of ne, and finaly to
1448 return ne. */
1449 int x,y,p;
1450 Value ex, ey ,ep;
1451 evalue *ne;
1452 value_init(ex); value_init(ey);value_init(ep);
1453 x=e1->x.p->size;
1454 y= res->x.p->size;
1455 value_set_si(ex,e1->x.p->size);
1456 value_set_si(ey,res->x.p->size);
1457 value_assign (ep,*Lcm(ex,ey));
1458 p=(int)mpz_get_si(ep);
1459 ne= (evalue *) malloc (sizeof(evalue));
1460 value_init(ne->d);
1461 value_set_si( ne->d,0);
1462
1463 ne->x.p=new_enode(res->x.p->type,p, res->x.p->pos);
1464 for(i=0;i<p;i++) {
1465 evalue_copy(&ne->x.p->arr[i], &res->x.p->arr[i%y]);
1466 }
1467 for(i=0;i<p;i++) {
1468 eadd(&e1->x.p->arr[i%x], &ne->x.p->arr[i]);
1469 }
1470
1471 value_assign(res->d, ne->d);
1472 res->x.p=ne->x.p;
1473
1474 return ;
1475 }
1476 case evector:
1477 fprintf(stderr, "eadd: ?cannot add vectors of different length\n");
1478 return ;
1479 default:
1480 assert(0);
1481 }
1482 }
1483 return ;
1484 }/* eadd */
1485
1486 static void emul_rev(const evalue *e1, evalue *res)
1487 {
1488 evalue ev;
1489 value_init(ev.d);
1490 evalue_copy(&ev, e1);
1491 emul(res, &ev);
1492 free_evalue_refs(res);
1493 *res = ev;
1494 }
1495
1496 static void emul_poly(const evalue *e1, evalue *res)
1497 {
1498 int i, j, offset = type_offset(res->x.p);
1499 evalue tmp;
1500 enode *p;
1501 int size = (e1->x.p->size + res->x.p->size - offset - 1);
1502
1503 p = new_enode(res->x.p->type, size, res->x.p->pos);
1504
1505 for (i = offset; i < e1->x.p->size-1; ++i)
1506 if (!EVALUE_IS_ZERO(e1->x.p->arr[i]))
1507 break;
1508
1509 /* special case pure power */
1510 if (i == e1->x.p->size-1) {
1511 if (offset) {
1512 value_clear(p->arr[0].d);
1513 p->arr[0] = res->x.p->arr[0];
1514 }
1515 for (i = offset; i < e1->x.p->size-1; ++i)
1516 evalue_set_si(&p->arr[i], 0, 1);
1517 for (i = offset; i < res->x.p->size; ++i) {
1518 value_clear(p->arr[i+e1->x.p->size-offset-1].d);
1519 p->arr[i+e1->x.p->size-offset-1] = res->x.p->arr[i];
1520 emul(&e1->x.p->arr[e1->x.p->size-1],
1521 &p->arr[i+e1->x.p->size-offset-1]);
1522 }
1523 free(res->x.p);
1524 res->x.p = p;
1525 return;
1526 }
1527
1528 value_init(tmp.d);
1529 value_set_si(tmp.d,0);
1530 tmp.x.p = p;
1531 if (offset)
1532 evalue_copy(&p->arr[0], &e1->x.p->arr[0]);
1533 for (i = offset; i < e1->x.p->size; i++) {
1534 evalue_copy(&tmp.x.p->arr[i], &e1->x.p->arr[i]);
1535 emul(&res->x.p->arr[offset], &tmp.x.p->arr[i]);
1536 }
1537 for (; i<size; i++)
1538 evalue_set_si(&tmp.x.p->arr[i], 0, 1);
1539 for (i = offset+1; i<res->x.p->size; i++)
1540 for (j = offset; j<e1->x.p->size; j++) {
1541 evalue ev;
1542 value_init(ev.d);
1543 evalue_copy(&ev, &e1->x.p->arr[j]);
1544 emul(&res->x.p->arr[i], &ev);
1545 eadd(&ev, &tmp.x.p->arr[i+j-offset]);
1546 free_evalue_refs(&ev);
1547 }
1548 free_evalue_refs(res);
1549 *res = tmp;
1550 }
1551
1552 void emul_partitions(const evalue *e1, evalue *res)
1553 {
1554 int n, i, j, k;
1555 Polyhedron *d;
1556 struct section *s;
1557 s = (struct section *)
1558 malloc((e1->x.p->size/2) * (res->x.p->size/2) *
1559 sizeof(struct section));
1560 assert(s);
1561 assert(e1->x.p->pos == res->x.p->pos);
1562 assert(e1->x.p->pos == EVALUE_DOMAIN(e1->x.p->arr[0])->Dimension);
1563 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1564
1565 n = 0;
1566 for (i = 0; i < res->x.p->size/2; ++i) {
1567 for (j = 0; j < e1->x.p->size/2; ++j) {
1568 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1569 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1570 if (emptyQ(d)) {
1571 Domain_Free(d);
1572 continue;
1573 }
1574
1575 /* This code is only needed because the partitions
1576 are not true partitions.
1577 */
1578 for (k = 0; k < n; ++k) {
1579 if (DomainIncludes(s[k].D, d))
1580 break;
1581 if (DomainIncludes(d, s[k].D)) {
1582 Domain_Free(s[k].D);
1583 free_evalue_refs(&s[k].E);
1584 if (n > k)
1585 s[k] = s[--n];
1586 --k;
1587 }
1588 }
1589 if (k < n) {
1590 Domain_Free(d);
1591 continue;
1592 }
1593
1594 value_init(s[n].E.d);
1595 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
1596 emul(&e1->x.p->arr[2*j+1], &s[n].E);
1597 s[n].D = d;
1598 ++n;
1599 }
1600 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
1601 value_clear(res->x.p->arr[2*i].d);
1602 free_evalue_refs(&res->x.p->arr[2*i+1]);
1603 }
1604
1605 free(res->x.p);
1606 if (n == 0)
1607 evalue_set_si(res, 0, 1);
1608 else {
1609 res->x.p = new_enode(partition, 2*n, e1->x.p->pos);
1610 for (j = 0; j < n; ++j) {
1611 s[j].D = DomainConstraintSimplify(s[j].D, 0);
1612 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1613 value_clear(res->x.p->arr[2*j+1].d);
1614 res->x.p->arr[2*j+1] = s[j].E;
1615 }
1616 }
1617
1618 free(s);
1619 }
1620
1621 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1622
1623 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1624 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1625 * evalues is not treated here */
1626
1627 void emul(const evalue *e1, evalue *res)
1628 {
1629 int i,j;
1630
1631 if((value_zero_p(e1->d)&&e1->x.p->type==evector)||(value_zero_p(res->d)&&(res->x.p->type==evector))) {
1632 fprintf(stderr, "emul: do not proced on evector type !\n");
1633 return;
1634 }
1635
1636 if (EVALUE_IS_ZERO(*res))
1637 return;
1638
1639 if (EVALUE_IS_ONE(*e1))
1640 return;
1641
1642 if (EVALUE_IS_ZERO(*e1)) {
1643 if (value_notzero_p(res->d)) {
1644 value_assign(res->d, e1->d);
1645 value_assign(res->x.n, e1->x.n);
1646 } else {
1647 free_evalue_refs(res);
1648 value_init(res->d);
1649 evalue_set_si(res, 0, 1);
1650 }
1651 return;
1652 }
1653
1654 if (value_zero_p(e1->d) && e1->x.p->type == partition) {
1655 if (value_zero_p(res->d) && res->x.p->type == partition)
1656 emul_partitions(e1, res);
1657 else
1658 emul_rev(e1, res);
1659 } else if (value_zero_p(res->d) && res->x.p->type == partition) {
1660 for (i = 0; i < res->x.p->size/2; ++i)
1661 emul(e1, &res->x.p->arr[2*i+1]);
1662 } else
1663 if (value_zero_p(res->d) && res->x.p->type == relation) {
1664 if (value_zero_p(e1->d) && e1->x.p->type == relation &&
1665 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1666 if (e1->x.p->size < 3 && res->x.p->size == 3) {
1667 free_evalue_refs(&res->x.p->arr[2]);
1668 res->x.p->size = 2;
1669 }
1670 for (i = 1; i < res->x.p->size; ++i)
1671 emul(&e1->x.p->arr[i], &res->x.p->arr[i]);
1672 return;
1673 }
1674 for (i = 1; i < res->x.p->size; ++i)
1675 emul(e1, &res->x.p->arr[i]);
1676 } else
1677 if(value_zero_p(e1->d)&& value_zero_p(res->d)) {
1678 switch(e1->x.p->type) {
1679 case polynomial:
1680 switch(res->x.p->type) {
1681 case polynomial:
1682 if(e1->x.p->pos == res->x.p->pos) {
1683 /* Product of two polynomials of the same variable */
1684 emul_poly(e1, res);
1685 return;
1686 }
1687 else {
1688 /* Product of two polynomials of different variables */
1689
1690 if(res->x.p->pos < e1->x.p->pos)
1691 for( i=0; i<res->x.p->size ; i++)
1692 emul(e1, &res->x.p->arr[i]);
1693 else
1694 emul_rev(e1, res);
1695
1696 return ;
1697 }
1698 case periodic:
1699 case flooring:
1700 case fractional:
1701 /* Product of a polynomial and a periodic or fractional */
1702 emul_rev(e1, res);
1703 return;
1704 default:
1705 assert(0);
1706 }
1707 case periodic:
1708 switch(res->x.p->type) {
1709 case periodic:
1710 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size==res->x.p->size) {
1711 /* Product of two periodics of the same parameter and period */
1712
1713 for(i=0; i<res->x.p->size;i++)
1714 emul(&(e1->x.p->arr[i]), &(res->x.p->arr[i]));
1715
1716 return;
1717 }
1718 else{
1719 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size!=res->x.p->size) {
1720 /* Product of two periodics of the same parameter and different periods */
1721 evalue *newp;
1722 Value x,y,z;
1723 int ix,iy,lcm;
1724 value_init(x); value_init(y);value_init(z);
1725 ix=e1->x.p->size;
1726 iy=res->x.p->size;
1727 value_set_si(x,e1->x.p->size);
1728 value_set_si(y,res->x.p->size);
1729 value_assign (z,*Lcm(x,y));
1730 lcm=(int)mpz_get_si(z);
1731 newp= (evalue *) malloc (sizeof(evalue));
1732 value_init(newp->d);
1733 value_set_si( newp->d,0);
1734 newp->x.p=new_enode(periodic,lcm, e1->x.p->pos);
1735 for(i=0;i<lcm;i++) {
1736 evalue_copy(&newp->x.p->arr[i],
1737 &res->x.p->arr[i%iy]);
1738 }
1739 for(i=0;i<lcm;i++)
1740 emul(&e1->x.p->arr[i%ix], &newp->x.p->arr[i]);
1741
1742 value_assign(res->d,newp->d);
1743 res->x.p=newp->x.p;
1744
1745 value_clear(x); value_clear(y);value_clear(z);
1746 return ;
1747 }
1748 else {
1749 /* Product of two periodics of different parameters */
1750
1751 if(res->x.p->pos < e1->x.p->pos)
1752 for(i=0; i<res->x.p->size; i++)
1753 emul(e1, &(res->x.p->arr[i]));
1754 else
1755 emul_rev(e1, res);
1756
1757 return;
1758 }
1759 }
1760 case polynomial:
1761 /* Product of a periodic and a polynomial */
1762
1763 for(i=0; i<res->x.p->size ; i++)
1764 emul(e1, &(res->x.p->arr[i]));
1765
1766 return;
1767
1768 }
1769 case flooring:
1770 case fractional:
1771 switch(res->x.p->type) {
1772 case polynomial:
1773 for(i=0; i<res->x.p->size ; i++)
1774 emul(e1, &(res->x.p->arr[i]));
1775 return;
1776 default:
1777 case periodic:
1778 assert(0);
1779 case flooring:
1780 case fractional:
1781 assert(e1->x.p->type == res->x.p->type);
1782 if (e1->x.p->pos == res->x.p->pos &&
1783 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1784 evalue d;
1785 value_init(d.d);
1786 poly_denom(&e1->x.p->arr[0], &d.d);
1787 if (e1->x.p->type != fractional || !value_two_p(d.d))
1788 emul_poly(e1, res);
1789 else {
1790 evalue tmp;
1791 value_init(d.x.n);
1792 value_set_si(d.x.n, 1);
1793 /* { x }^2 == { x }/2 */
1794 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1795 assert(e1->x.p->size == 3);
1796 assert(res->x.p->size == 3);
1797 value_init(tmp.d);
1798 evalue_copy(&tmp, &res->x.p->arr[2]);
1799 emul(&d, &tmp);
1800 eadd(&res->x.p->arr[1], &tmp);
1801 emul(&e1->x.p->arr[2], &tmp);
1802 emul(&e1->x.p->arr[1], &res->x.p->arr[1]);
1803 emul(&e1->x.p->arr[1], &res->x.p->arr[2]);
1804 eadd(&tmp, &res->x.p->arr[2]);
1805 free_evalue_refs(&tmp);
1806 value_clear(d.x.n);
1807 }
1808 value_clear(d.d);
1809 } else {
1810 if(mod_term_smaller(res, e1))
1811 for(i=1; i<res->x.p->size ; i++)
1812 emul(e1, &(res->x.p->arr[i]));
1813 else
1814 emul_rev(e1, res);
1815 return;
1816 }
1817 }
1818 break;
1819 case relation:
1820 emul_rev(e1, res);
1821 break;
1822 default:
1823 assert(0);
1824 }
1825 }
1826 else {
1827 if (value_notzero_p(e1->d)&& value_notzero_p(res->d)) {
1828 /* Product of two rational numbers */
1829
1830 Value g;
1831 value_init(g);
1832 value_multiply(res->d,e1->d,res->d);
1833 value_multiply(res->x.n,e1->x.n,res->x.n );
1834 value_gcd(g, res->x.n, res->d);
1835 if (value_notone_p(g)) {
1836 value_division(res->d,res->d,g);
1837 value_division(res->x.n,res->x.n,g);
1838 }
1839 value_clear(g);
1840 return ;
1841 }
1842 else {
1843 if(value_zero_p(e1->d)&& value_notzero_p(res->d)) {
1844 /* Product of an expression (polynomial or peririodic) and a rational number */
1845
1846 emul_rev(e1, res);
1847 return ;
1848 }
1849 else {
1850 /* Product of a rationel number and an expression (polynomial or peririodic) */
1851
1852 i = type_offset(res->x.p);
1853 for (; i<res->x.p->size; i++)
1854 emul(e1, &res->x.p->arr[i]);
1855
1856 return ;
1857 }
1858 }
1859 }
1860
1861 return ;
1862 }
1863
1864 /* Frees mask content ! */
1865 void emask(evalue *mask, evalue *res) {
1866 int n, i, j;
1867 Polyhedron *d, *fd;
1868 struct section *s;
1869 evalue mone;
1870 int pos;
1871
1872 if (EVALUE_IS_ZERO(*res)) {
1873 free_evalue_refs(mask);
1874 return;
1875 }
1876
1877 assert(value_zero_p(mask->d));
1878 assert(mask->x.p->type == partition);
1879 assert(value_zero_p(res->d));
1880 assert(res->x.p->type == partition);
1881 assert(mask->x.p->pos == res->x.p->pos);
1882 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1883 assert(mask->x.p->pos == EVALUE_DOMAIN(mask->x.p->arr[0])->Dimension);
1884 pos = res->x.p->pos;
1885
1886 s = (struct section *)
1887 malloc((mask->x.p->size/2+1) * (res->x.p->size/2) *
1888 sizeof(struct section));
1889 assert(s);
1890
1891 value_init(mone.d);
1892 evalue_set_si(&mone, -1, 1);
1893
1894 n = 0;
1895 for (j = 0; j < res->x.p->size/2; ++j) {
1896 assert(mask->x.p->size >= 2);
1897 fd = DomainDifference(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1898 EVALUE_DOMAIN(mask->x.p->arr[0]), 0);
1899 if (!emptyQ(fd))
1900 for (i = 1; i < mask->x.p->size/2; ++i) {
1901 Polyhedron *t = fd;
1902 fd = DomainDifference(fd, EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1903 Domain_Free(t);
1904 if (emptyQ(fd))
1905 break;
1906 }
1907 if (emptyQ(fd)) {
1908 Domain_Free(fd);
1909 continue;
1910 }
1911 value_init(s[n].E.d);
1912 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1913 s[n].D = fd;
1914 ++n;
1915 }
1916 for (i = 0; i < mask->x.p->size/2; ++i) {
1917 if (EVALUE_IS_ONE(mask->x.p->arr[2*i+1]))
1918 continue;
1919
1920 fd = EVALUE_DOMAIN(mask->x.p->arr[2*i]);
1921 eadd(&mone, &mask->x.p->arr[2*i+1]);
1922 emul(&mone, &mask->x.p->arr[2*i+1]);
1923 for (j = 0; j < res->x.p->size/2; ++j) {
1924 Polyhedron *t;
1925 d = DomainIntersection(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1926 EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1927 if (emptyQ(d)) {
1928 Domain_Free(d);
1929 continue;
1930 }
1931 t = fd;
1932 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*j]), 0);
1933 if (t != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1934 Domain_Free(t);
1935 value_init(s[n].E.d);
1936 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1937 emul(&mask->x.p->arr[2*i+1], &s[n].E);
1938 s[n].D = d;
1939 ++n;
1940 }
1941
1942 if (!emptyQ(fd)) {
1943 /* Just ignore; this may have been previously masked off */
1944 }
1945 if (fd != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1946 Domain_Free(fd);
1947 }
1948
1949 free_evalue_refs(&mone);
1950 free_evalue_refs(mask);
1951 free_evalue_refs(res);
1952 value_init(res->d);
1953 if (n == 0)
1954 evalue_set_si(res, 0, 1);
1955 else {
1956 res->x.p = new_enode(partition, 2*n, pos);
1957 for (j = 0; j < n; ++j) {
1958 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1959 value_clear(res->x.p->arr[2*j+1].d);
1960 res->x.p->arr[2*j+1] = s[j].E;
1961 }
1962 }
1963
1964 free(s);
1965 }
1966
1967 void evalue_copy(evalue *dst, const evalue *src)
1968 {
1969 value_assign(dst->d, src->d);
1970 if(value_notzero_p(src->d)) {
1971 value_init(dst->x.n);
1972 value_assign(dst->x.n, src->x.n);
1973 } else
1974 dst->x.p = ecopy(src->x.p);
1975 }
1976
1977 evalue *evalue_dup(const evalue *e)
1978 {
1979 evalue *res = ALLOC(evalue);
1980 value_init(res->d);
1981 evalue_copy(res, e);
1982 return res;
1983 }
1984
1985 enode *new_enode(enode_type type,int size,int pos) {
1986
1987 enode *res;
1988 int i;
1989
1990 if(size == 0) {
1991 fprintf(stderr, "Allocating enode of size 0 !\n" );
1992 return NULL;
1993 }
1994 res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
1995 res->type = type;
1996 res->size = size;
1997 res->pos = pos;
1998 for(i=0; i<size; i++) {
1999 value_init(res->arr[i].d);
2000 value_set_si(res->arr[i].d,0);
2001 res->arr[i].x.p = 0;
2002 }
2003 return res;
2004 } /* new_enode */
2005
2006 enode *ecopy(enode *e) {
2007
2008 enode *res;
2009 int i;
2010
2011 res = new_enode(e->type,e->size,e->pos);
2012 for(i=0;i<e->size;++i) {
2013 value_assign(res->arr[i].d,e->arr[i].d);
2014 if(value_zero_p(res->arr[i].d))
2015 res->arr[i].x.p = ecopy(e->arr[i].x.p);
2016 else if (EVALUE_IS_DOMAIN(res->arr[i]))
2017 EVALUE_SET_DOMAIN(res->arr[i], Domain_Copy(EVALUE_DOMAIN(e->arr[i])));
2018 else {
2019 value_init(res->arr[i].x.n);
2020 value_assign(res->arr[i].x.n,e->arr[i].x.n);
2021 }
2022 }
2023 return(res);
2024 } /* ecopy */
2025
2026 int ecmp(const evalue *e1, const evalue *e2)
2027 {
2028 enode *p1, *p2;
2029 int i;
2030 int r;
2031
2032 if (value_notzero_p(e1->d) && value_notzero_p(e2->d)) {
2033 Value m, m2;
2034 value_init(m);
2035 value_init(m2);
2036 value_multiply(m, e1->x.n, e2->d);
2037 value_multiply(m2, e2->x.n, e1->d);
2038
2039 if (value_lt(m, m2))
2040 r = -1;
2041 else if (value_gt(m, m2))
2042 r = 1;
2043 else
2044 r = 0;
2045
2046 value_clear(m);
2047 value_clear(m2);
2048
2049 return r;
2050 }
2051 if (value_notzero_p(e1->d))
2052 return -1;
2053 if (value_notzero_p(e2->d))
2054 return 1;
2055
2056 p1 = e1->x.p;
2057 p2 = e2->x.p;
2058
2059 if (p1->type != p2->type)
2060 return p1->type - p2->type;
2061 if (p1->pos != p2->pos)
2062 return p1->pos - p2->pos;
2063 if (p1->size != p2->size)
2064 return p1->size - p2->size;
2065
2066 for (i = p1->size-1; i >= 0; --i)
2067 if ((r = ecmp(&p1->arr[i], &p2->arr[i])) != 0)
2068 return r;
2069
2070 return 0;
2071 }
2072
2073 int eequal(const evalue *e1, const evalue *e2)
2074 {
2075 int i;
2076 enode *p1, *p2;
2077
2078 if (value_ne(e1->d,e2->d))
2079 return 0;
2080
2081 /* e1->d == e2->d */
2082 if (value_notzero_p(e1->d)) {
2083 if (value_ne(e1->x.n,e2->x.n))
2084 return 0;
2085
2086 /* e1->d == e2->d != 0 AND e1->n == e2->n */
2087 return 1;
2088 }
2089
2090 /* e1->d == e2->d == 0 */
2091 p1 = e1->x.p;
2092 p2 = e2->x.p;
2093 if (p1->type != p2->type) return 0;
2094 if (p1->size != p2->size) return 0;
2095 if (p1->pos != p2->pos) return 0;
2096 for (i=0; i<p1->size; i++)
2097 if (!eequal(&p1->arr[i], &p2->arr[i]) )
2098 return 0;
2099 return 1;
2100 } /* eequal */
2101
2102 void free_evalue_refs(evalue *e) {
2103
2104 enode *p;
2105 int i;
2106
2107 if (EVALUE_IS_DOMAIN(*e)) {
2108 Domain_Free(EVALUE_DOMAIN(*e));
2109 value_clear(e->d);
2110 return;
2111 } else if (value_pos_p(e->d)) {
2112
2113 /* 'e' stores a constant */
2114 value_clear(e->d);
2115 value_clear(e->x.n);
2116 return;
2117 }
2118 assert(value_zero_p(e->d));
2119 value_clear(e->d);
2120 p = e->x.p;
2121 if (!p) return; /* null pointer */
2122 for (i=0; i<p->size; i++) {
2123 free_evalue_refs(&(p->arr[i]));
2124 }
2125 free(p);
2126 return;
2127 } /* free_evalue_refs */
2128
2129 void evalue_free(evalue *e)
2130 {
2131 free_evalue_refs(e);
2132 free(e);
2133 }
2134
2135 static void mod2table_r(evalue *e, Vector *periods, Value m, int p,
2136 Vector * val, evalue *res)
2137 {
2138 unsigned nparam = periods->Size;
2139
2140 if (p == nparam) {
2141 double d = compute_evalue(e, val->p);
2142 d *= VALUE_TO_DOUBLE(m);
2143 if (d > 0)
2144 d += .25;
2145 else
2146 d -= .25;
2147 value_assign(res->d, m);
2148 value_init(res->x.n);
2149 value_set_double(res->x.n, d);
2150 mpz_fdiv_r(res->x.n, res->x.n, m);
2151 return;
2152 }
2153 if (value_one_p(periods->p[p]))
2154 mod2table_r(e, periods, m, p+1, val, res);
2155 else {
2156 Value tmp;
2157 value_init(tmp);
2158
2159 value_assign(tmp, periods->p[p]);
2160 value_set_si(res->d, 0);
2161 res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
2162 do {
2163 value_decrement(tmp, tmp);
2164 value_assign(val->p[p], tmp);
2165 mod2table_r(e, periods, m, p+1, val,
2166 &res->x.p->arr[VALUE_TO_INT(tmp)]);
2167 } while (value_pos_p(tmp));
2168
2169 value_clear(tmp);
2170 }
2171 }
2172
2173 static void rel2table(evalue *e, int zero)
2174 {
2175 if (value_pos_p(e->d)) {
2176 if (value_zero_p(e->x.n) == zero)
2177 value_set_si(e->x.n, 1);
2178 else
2179 value_set_si(e->x.n, 0);
2180 value_set_si(e->d, 1);
2181 } else {
2182 int i;
2183 for (i = 0; i < e->x.p->size; ++i)
2184 rel2table(&e->x.p->arr[i], zero);
2185 }
2186 }
2187
2188 void evalue_mod2table(evalue *e, int nparam)
2189 {
2190 enode *p;
2191 int i;
2192
2193 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2194 return;
2195 p = e->x.p;
2196 for (i=0; i<p->size; i++) {
2197 evalue_mod2table(&(p->arr[i]), nparam);
2198 }
2199 if (p->type == relation) {
2200 evalue copy;
2201
2202 if (p->size > 2) {
2203 value_init(copy.d);
2204 evalue_copy(&copy, &p->arr[0]);
2205 }
2206 rel2table(&p->arr[0], 1);
2207 emul(&p->arr[0], &p->arr[1]);
2208 if (p->size > 2) {
2209 rel2table(&copy, 0);
2210 emul(&copy, &p->arr[2]);
2211 eadd(&p->arr[2], &p->arr[1]);
2212 free_evalue_refs(&p->arr[2]);
2213 free_evalue_refs(&copy);
2214 }
2215 free_evalue_refs(&p->arr[0]);
2216 value_clear(e->d);
2217 *e = p->arr[1];
2218 free(p);
2219 } else if (p->type == fractional) {
2220 Vector *periods = Vector_Alloc(nparam);
2221 Vector *val = Vector_Alloc(nparam);
2222 Value tmp;
2223 evalue *ev;
2224 evalue EP, res;
2225
2226 value_init(tmp);
2227 value_set_si(tmp, 1);
2228 Vector_Set(periods->p, 1, nparam);
2229 Vector_Set(val->p, 0, nparam);
2230 for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
2231 enode *p = ev->x.p;
2232
2233 assert(p->type == polynomial);
2234 assert(p->size == 2);
2235 value_assign(periods->p[p->pos-1], p->arr[1].d);
2236 value_lcm(tmp, tmp, p->arr[1].d);
2237 }
2238 value_lcm(tmp, tmp, ev->d);
2239 value_init(EP.d);
2240 mod2table_r(&p->arr[0], periods, tmp, 0, val, &EP);
2241
2242 value_init(res.d);
2243 evalue_set_si(&res, 0, 1);
2244 /* Compute the polynomial using Horner's rule */
2245 for (i=p->size-1;i>1;i--) {
2246 eadd(&p->arr[i], &res);
2247 emul(&EP, &res);
2248 }
2249 eadd(&p->arr[1], &res);
2250
2251 free_evalue_refs(e);
2252 free_evalue_refs(&EP);
2253 *e = res;
2254
2255 value_clear(tmp);
2256 Vector_Free(val);
2257 Vector_Free(periods);
2258 }
2259 } /* evalue_mod2table */
2260
2261 /********************************************************/
2262 /* function in domain */
2263 /* check if the parameters in list_args */
2264 /* verifies the constraints of Domain P */
2265 /********************************************************/
2266 int in_domain(Polyhedron *P, Value *list_args)
2267 {
2268 int row, in = 1;
2269 Value v; /* value of the constraint of a row when
2270 parameters are instantiated*/
2271
2272 value_init(v);
2273
2274 for (row = 0; row < P->NbConstraints; row++) {
2275 Inner_Product(P->Constraint[row]+1, list_args, P->Dimension, &v);
2276 value_addto(v, v, P->Constraint[row][P->Dimension+1]); /*constant part*/
2277 if (value_neg_p(v) ||
2278 value_zero_p(P->Constraint[row][0]) && value_notzero_p(v)) {
2279 in = 0;
2280 break;
2281 }
2282 }
2283
2284 value_clear(v);
2285 return in || (P->next && in_domain(P->next, list_args));
2286 } /* in_domain */
2287
2288 /****************************************************/
2289 /* function compute enode */
2290 /* compute the value of enode p with parameters */
2291 /* list "list_args */
2292 /* compute the polynomial or the periodic */
2293 /****************************************************/
2294
2295 static double compute_enode(enode *p, Value *list_args) {
2296
2297 int i;
2298 Value m, param;
2299 double res=0.0;
2300
2301 if (!p)
2302 return(0.);
2303
2304 value_init(m);
2305 value_init(param);
2306
2307 if (p->type == polynomial) {
2308 if (p->size > 1)
2309 value_assign(param,list_args[p->pos-1]);
2310
2311 /* Compute the polynomial using Horner's rule */
2312 for (i=p->size-1;i>0;i--) {
2313 res +=compute_evalue(&p->arr[i],list_args);
2314 res *=VALUE_TO_DOUBLE(param);
2315 }
2316 res +=compute_evalue(&p->arr[0],list_args);
2317 }
2318 else if (p->type == fractional) {
2319 double d = compute_evalue(&p->arr[0], list_args);
2320 d -= floor(d+1e-10);
2321
2322 /* Compute the polynomial using Horner's rule */
2323 for (i=p->size-1;i>1;i--) {
2324 res +=compute_evalue(&p->arr[i],list_args);
2325 res *=d;
2326 }
2327 res +=compute_evalue(&p->arr[1],list_args);
2328 }
2329 else if (p->type == flooring) {
2330 double d = compute_evalue(&p->arr[0], list_args);
2331 d = floor(d+1e-10);
2332
2333 /* Compute the polynomial using Horner's rule */
2334 for (i=p->size-1;i>1;i--) {
2335 res +=compute_evalue(&p->arr[i],list_args);
2336 res *=d;
2337 }
2338 res +=compute_evalue(&p->arr[1],list_args);
2339 }
2340 else if (p->type == periodic) {
2341 value_assign(m,list_args[p->pos-1]);
2342
2343 /* Choose the right element of the periodic */
2344 value_set_si(param,p->size);
2345 value_pmodulus(m,m,param);
2346 res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
2347 }
2348 else if (p->type == relation) {
2349 if (fabs(compute_evalue(&p->arr[0], list_args)) < 1e-10)
2350 res = compute_evalue(&p->arr[1], list_args);
2351 else if (p->size > 2)
2352 res = compute_evalue(&p->arr[2], list_args);
2353 }
2354 else if (p->type == partition) {
2355 int dim = EVALUE_DOMAIN(p->arr[0])->Dimension;
2356 Value *vals = list_args;
2357 if (p->pos < dim) {
2358 NALLOC(vals, dim);
2359 for (i = 0; i < dim; ++i) {
2360 value_init(vals[i]);
2361 if (i < p->pos)
2362 value_assign(vals[i], list_args[i]);
2363 }
2364 }
2365 for (i = 0; i < p->size/2; ++i)
2366 if (DomainContains(EVALUE_DOMAIN(p->arr[2*i]), vals, p->pos, 0, 1)) {
2367 res = compute_evalue(&p->arr[2*i+1], vals);
2368 break;
2369 }
2370 if (p->pos < dim) {
2371 for (i = 0; i < dim; ++i)
2372 value_clear(vals[i]);
2373 free(vals);
2374 }
2375 }
2376 else
2377 assert(0);
2378 value_clear(m);
2379 value_clear(param);
2380 return res;
2381 } /* compute_enode */
2382
2383 /*************************************************/
2384 /* return the value of Ehrhart Polynomial */
2385 /* It returns a double, because since it is */
2386 /* a recursive function, some intermediate value */
2387 /* might not be integral */
2388 /*************************************************/
2389
2390 double compute_evalue(const evalue *e, Value *list_args)
2391 {
2392 double res;
2393
2394 if (value_notzero_p(e->d)) {
2395 if (value_notone_p(e->d))
2396 res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
2397 else
2398 res = VALUE_TO_DOUBLE(e->x.n);
2399 }
2400 else
2401 res = compute_enode(e->x.p,list_args);
2402 return res;
2403 } /* compute_evalue */
2404
2405
2406 /****************************************************/
2407 /* function compute_poly : */
2408 /* Check for the good validity domain */
2409 /* return the number of point in the Polyhedron */
2410 /* in allocated memory */
2411 /* Using the Ehrhart pseudo-polynomial */
2412 /****************************************************/
2413 Value *compute_poly(Enumeration *en,Value *list_args) {
2414
2415 Value *tmp;
2416 /* double d; int i; */
2417
2418 tmp = (Value *) malloc (sizeof(Value));
2419 assert(tmp != NULL);
2420 value_init(*tmp);
2421 value_set_si(*tmp,0);
2422
2423 if(!en)
2424 return(tmp); /* no ehrhart polynomial */
2425 if(en->ValidityDomain) {
2426 if(!en->ValidityDomain->Dimension) { /* no parameters */
2427 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2428 return(tmp);
2429 }
2430 }
2431 else
2432 return(tmp); /* no Validity Domain */
2433 while(en) {
2434 if(in_domain(en->ValidityDomain,list_args)) {
2435
2436 #ifdef EVAL_EHRHART_DEBUG
2437 Print_Domain(stdout,en->ValidityDomain);
2438 print_evalue(stdout,&en->EP);
2439 #endif
2440
2441 /* d = compute_evalue(&en->EP,list_args);
2442 i = d;
2443 printf("(double)%lf = %d\n", d, i ); */
2444 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2445 return(tmp);
2446 }
2447 else
2448 en=en->next;
2449 }
2450 value_set_si(*tmp,0);
2451 return(tmp); /* no compatible domain with the arguments */
2452 } /* compute_poly */
2453
2454 static evalue *eval_polynomial(const enode *p, int offset,
2455 evalue *base, Value *values)
2456 {
2457 int i;
2458 evalue *res, *c;
2459
2460 res = evalue_zero();
2461 for (i = p->size-1; i > offset; --i) {
2462 c = evalue_eval(&p->arr[i], values);
2463 eadd(c, res);
2464 evalue_free(c);
2465 emul(base, res);
2466 }
2467 c = evalue_eval(&p->arr[offset], values);
2468 eadd(c, res);
2469 evalue_free(c);
2470
2471 return res;
2472 }
2473
2474 evalue *evalue_eval(const evalue *e, Value *values)
2475 {
2476 evalue *res = NULL;
2477 evalue param;
2478 evalue *param2;
2479 int i;
2480
2481 if (value_notzero_p(e->d)) {
2482 res = ALLOC(evalue);
2483 value_init(res->d);
2484 evalue_copy(res, e);
2485 return res;
2486 }
2487 switch (e->x.p->type) {
2488 case polynomial:
2489 value_init(param.x.n);
2490 value_assign(param.x.n, values[e->x.p->pos-1]);
2491 value_init(param.d);
2492 value_set_si(param.d, 1);
2493
2494 res = eval_polynomial(e->x.p, 0, &param, values);
2495 free_evalue_refs(&param);
2496 break;
2497 case fractional:
2498 param2 = evalue_eval(&e->x.p->arr[0], values);
2499 mpz_fdiv_r(param2->x.n, param2->x.n, param2->d);
2500
2501 res = eval_polynomial(e->x.p, 1, param2, values);
2502 evalue_free(param2);
2503 break;
2504 case flooring:
2505 param2 = evalue_eval(&e->x.p->arr[0], values);
2506 mpz_fdiv_q(param2->x.n, param2->x.n, param2->d);
2507 value_set_si(param2->d, 1);
2508
2509 res = eval_polynomial(e->x.p, 1, param2, values);
2510 evalue_free(param2);
2511 break;
2512 case relation:
2513 param2 = evalue_eval(&e->x.p->arr[0], values);
2514 if (value_zero_p(param2->x.n))
2515 res = evalue_eval(&e->x.p->arr[1], values);
2516 else if (e->x.p->size > 2)
2517 res = evalue_eval(&e->x.p->arr[2], values);
2518 else
2519 res = evalue_zero();
2520 evalue_free(param2);
2521 break;
2522 case partition:
2523 assert(e->x.p->pos == EVALUE_DOMAIN(e->x.p->arr[0])->Dimension);
2524 for (i = 0; i < e->x.p->size/2; ++i)
2525 if (in_domain(EVALUE_DOMAIN(e->x.p->arr[2*i]), values)) {
2526 res = evalue_eval(&e->x.p->arr[2*i+1], values);
2527 break;
2528 }
2529 if (!res)
2530 res = evalue_zero();
2531 break;
2532 default:
2533 assert(0);
2534 }
2535 return res;
2536 }
2537
2538 size_t value_size(Value v) {
2539 return (v[0]._mp_size > 0 ? v[0]._mp_size : -v[0]._mp_size)
2540 * sizeof(v[0]._mp_d[0]);
2541 }
2542
2543 size_t domain_size(Polyhedron *D)
2544 {
2545 int i, j;
2546 size_t s = sizeof(*D);
2547
2548 for (i = 0; i < D->NbConstraints; ++i)
2549 for (j = 0; j < D->Dimension+2; ++j)
2550 s += value_size(D->Constraint[i][j]);
2551
2552 /*
2553 for (i = 0; i < D->NbRays; ++i)
2554 for (j = 0; j < D->Dimension+2; ++j)
2555 s += value_size(D->Ray[i][j]);
2556 */
2557
2558 return D->next ? s+domain_size(D->next) : s;
2559 }
2560
2561 size_t enode_size(enode *p) {
2562 size_t s = sizeof(*p) - sizeof(p->arr[0]);
2563 int i;
2564
2565 if (p->type == partition)
2566 for (i = 0; i < p->size/2; ++i) {
2567 s += domain_size(EVALUE_DOMAIN(p->arr[2*i]));
2568 s += evalue_size(&p->arr[2*i+1]);
2569 }
2570 else
2571 for (i = 0; i < p->size; ++i) {
2572 s += evalue_size(&p->arr[i]);
2573 }
2574 return s;
2575 }
2576
2577 size_t evalue_size(evalue *e)
2578 {
2579 size_t s = sizeof(*e);
2580 s += value_size(e->d);
2581 if (value_notzero_p(e->d))
2582 s += value_size(e->x.n);
2583 else
2584 s += enode_size(e->x.p);
2585 return s;
2586 }
2587
2588 static evalue *find_second(evalue *base, evalue *cst, evalue *e, Value m)
2589 {
2590 evalue *found = NULL;
2591 evalue offset;
2592 evalue copy;
2593 int i;
2594
2595 if (value_pos_p(e->d) || e->x.p->type != fractional)
2596 return NULL;
2597
2598 value_init(offset.d);
2599 value_init(offset.x.n);
2600 poly_denom(&e->x.p->arr[0], &offset.d);
2601 value_lcm(offset.d, m, offset.d);
2602 value_set_si(offset.x.n, 1);
2603
2604 value_init(copy.d);
2605 evalue_copy(&copy, cst);
2606
2607 eadd(&offset, cst);
2608 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2609
2610 if (eequal(base, &e->x.p->arr[0]))
2611 found = &e->x.p->arr[0];
2612 else {
2613 value_set_si(offset.x.n, -2);
2614
2615 eadd(&offset, cst);
2616 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2617
2618 if (eequal(base, &e->x.p->arr[0]))
2619 found = base;
2620 }
2621 free_evalue_refs(cst);
2622 free_evalue_refs(&offset);
2623 *cst = copy;
2624
2625 for (i = 1; !found && i < e->x.p->size; ++i)
2626 found = find_second(base, cst, &e->x.p->arr[i], m);
2627
2628 return found;
2629 }
2630
2631 static evalue *find_relation_pair(evalue *e)
2632 {
2633 int i;
2634 evalue *found = NULL;
2635
2636 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2637 return NULL;
2638
2639 if (e->x.p->type == fractional) {
2640 Value m;
2641 evalue *cst;
2642
2643 value_init(m);
2644 poly_denom(&e->x.p->arr[0], &m);
2645
2646 for (cst = &e->x.p->arr[0]; value_zero_p(cst->d);
2647 cst = &cst->x.p->arr[0])
2648 ;
2649
2650 for (i = 1; !found && i < e->x.p->size; ++i)
2651 found = find_second(&e->x.p->arr[0], cst, &e->x.p->arr[i], m);
2652
2653 value_clear(m);
2654 }
2655
2656 i = e->x.p->type == relation;
2657 for (; !found && i < e->x.p->size; ++i)
2658 found = find_relation_pair(&e->x.p->arr[i]);
2659
2660 return found;
2661 }
2662
2663 void evalue_mod2relation(evalue *e) {
2664 evalue *d;
2665
2666 if (value_zero_p(e->d) && e->x.p->type == partition) {
2667 int i;
2668
2669 for (i = 0; i < e->x.p->size/2; ++i) {
2670 evalue_mod2relation(&e->x.p->arr[2*i+1]);
2671 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2672 value_clear(e->x.p->arr[2*i].d);
2673 free_evalue_refs(&e->x.p->arr[2*i+1]);
2674 e->x.p->size -= 2;
2675 if (2*i < e->x.p->size) {
2676 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2677 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2678 }
2679 --i;
2680 }
2681 }
2682 if (e->x.p->size == 0) {
2683 free(e->x.p);
2684 evalue_set_si(e, 0, 1);
2685 }
2686
2687 return;
2688 }
2689
2690 while ((d = find_relation_pair(e)) != NULL) {
2691 evalue split;
2692 evalue *ev;
2693
2694 value_init(split.d);
2695 value_set_si(split.d, 0);
2696 split.x.p = new_enode(relation, 3, 0);
2697 evalue_set_si(&split.x.p->arr[1], 1, 1);
2698 evalue_set_si(&split.x.p->arr[2], 1, 1);
2699
2700 ev = &split.x.p->arr[0];
2701 value_set_si(ev->d, 0);
2702 ev->x.p = new_enode(fractional, 3, -1);
2703 evalue_set_si(&ev->x.p->arr[1], 0, 1);
2704 evalue_set_si(&ev->x.p->arr[2], 1, 1);
2705 evalue_copy(&ev->x.p->arr[0], d);
2706
2707 emul(&split, e);
2708
2709 reduce_evalue(e);
2710
2711 free_evalue_refs(&split);
2712 }
2713 }
2714
2715 static int evalue_comp(const void * a, const void * b)
2716 {
2717 const evalue *e1 = *(const evalue **)a;
2718 const evalue *e2 = *(const evalue **)b;
2719 return ecmp(e1, e2);
2720 }
2721
2722 void evalue_combine(evalue *e)
2723 {
2724 evalue **evs;
2725 int i, k;
2726 enode *p;
2727 evalue tmp;
2728
2729 if (value_notzero_p(e->d) || e->x.p->type != partition)
2730 return;
2731
2732 NALLOC(evs, e->x.p->size/2);
2733 for (i = 0; i < e->x.p->size/2; ++i)
2734 evs[i] = &e->x.p->arr[2*i+1];
2735 qsort(evs, e->x.p->size/2, sizeof(evs[0]), evalue_comp);
2736 p = new_enode(partition, e->x.p->size, e->x.p->pos);
2737 for (i = 0, k = 0; i < e->x.p->size/2; ++i) {
2738 if (k == 0 || ecmp(&p->arr[2*k-1], evs[i]) != 0) {
2739 value_clear(p->arr[2*k].d);
2740 value_clear(p->arr[2*k+1].d);
2741 p->arr[2*k] = *(evs[i]-1);
2742 p->arr[2*k+1] = *(evs[i]);
2743 ++k;
2744 } else {
2745 Polyhedron *D = EVALUE_DOMAIN(*(evs[i]-1));
2746 Polyhedron *L = D;
2747
2748 value_clear((evs[i]-1)->d);
2749
2750 while (L->next)
2751 L = L->next;
2752 L->next = EVALUE_DOMAIN(p->arr[2*k-2]);
2753 EVALUE_SET_DOMAIN(p->arr[2*k-2], D);
2754 free_evalue_refs(evs[i]);
2755 }
2756 }
2757
2758 for (i = 2*k ; i < p->size; ++i)
2759 value_clear(p->arr[i].d);
2760
2761 free(evs);
2762 free(e->x.p);
2763 p->size = 2*k;
2764 e->x.p = p;
2765
2766 for (i = 0; i < e->x.p->size/2; ++i) {
2767 Polyhedron *H;
2768 if (value_notzero_p(e->x.p->arr[2*i+1].d))
2769 continue;
2770 H = DomainConvex(EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
2771 if (H == NULL)
2772 continue;
2773 for (k = 0; k < e->x.p->size/2; ++k) {
2774 Polyhedron *D, *N, **P;
2775 if (i == k)
2776 continue;
2777 P = &EVALUE_DOMAIN(e->x.p->arr[2*k]);
2778 D = *P;
2779 if (D == NULL)
2780 continue;
2781 for (; D; D = N) {
2782 *P = D;
2783 N = D->next;
2784 if (D->NbEq <= H->NbEq) {
2785 P = &D->next;
2786 continue;
2787 }
2788
2789 value_init(tmp.d);
2790 tmp.x.p = new_enode(partition, 2, e->x.p->pos);
2791 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Copy(D));
2792 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
2793 reduce_evalue(&tmp);
2794 if (value_notzero_p(tmp.d) ||
2795 ecmp(&tmp.x.p->arr[1], &e->x.p->arr[2*k+1]) != 0)
2796 P = &D->next;
2797 else {
2798 D->next = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2799 EVALUE_DOMAIN(e->x.p->arr[2*i]) = D;
2800 *P = NULL;
2801 }
2802 free_evalue_refs(&tmp);
2803 }
2804 }
2805 Polyhedron_Free(H);
2806 }
2807
2808 for (i = 0; i < e->x.p->size/2; ++i) {
2809 Polyhedron *H, *E;
2810 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2811 if (!D) {
2812 value_clear(e->x.p->arr[2*i].d);
2813 free_evalue_refs(&e->x.p->arr[2*i+1]);
2814 e->x.p->size -= 2;
2815 if (2*i < e->x.p->size) {
2816 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2817 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2818 }
2819 --i;
2820 continue;
2821 }
2822 if (!D->next)
2823 continue;
2824 H = DomainConvex(D, 0);
2825 E = DomainDifference(H, D, 0);
2826 Domain_Free(D);
2827 D = DomainDifference(H, E, 0);
2828 Domain_Free(H);
2829 Domain_Free(E);
2830 EVALUE_SET_DOMAIN(p->arr[2*i], D);
2831 }
2832 }
2833
2834 /* Use smallest representative for coefficients in affine form in
2835 * argument of fractional.
2836 * Since any change will make the argument non-standard,
2837 * the containing evalue will have to be reduced again afterward.
2838 */
2839 static void fractional_minimal_coefficients(enode *p)
2840 {
2841 evalue *pp;
2842 Value twice;
2843 value_init(twice);
2844
2845 assert(p->type == fractional);
2846 pp = &p->arr[0];
2847 while (value_zero_p(pp->d)) {
2848 assert(pp->x.p->type == polynomial);
2849 assert(pp->x.p->size == 2);
2850 assert(value_notzero_p(pp->x.p->arr[1].d));
2851 mpz_mul_ui(twice, pp->x.p->arr[1].x.n, 2);
2852 if (value_gt(twice, pp->x.p->arr[1].d))
2853 value_subtract(pp->x.p->arr[1].x.n,
2854 pp->x.p->arr[1].x.n, pp->x.p->arr[1].d);
2855 pp = &pp->x.p->arr[0];
2856 }
2857
2858 value_clear(twice);
2859 }
2860
2861 static Polyhedron *polynomial_projection(enode *p, Polyhedron *D, Value *d,
2862 Matrix **R)
2863 {
2864 Polyhedron *I, *H;
2865 evalue *pp;
2866 unsigned dim = D->Dimension;
2867 Matrix *T = Matrix_Alloc(2, dim+1);
2868 assert(T);
2869
2870 assert(p->type == fractional || p->type == flooring);
2871 value_set_si(T->p[1][dim], 1);
2872 evalue_extract_affine(&p->arr[0], T->p[0], &T->p[0][dim], d);
2873 I = DomainImage(D, T, 0);
2874 H = DomainConvex(I, 0);
2875 Domain_Free(I);
2876 if (R)
2877 *R = T;
2878 else
2879 Matrix_Free(T);
2880
2881 return H;
2882 }
2883
2884 static void replace_by_affine(evalue *e, Value offset)
2885 {
2886 enode *p;
2887 evalue inc;
2888
2889 p = e->x.p;
2890 value_init(inc.d);
2891 value_init(inc.x.n);
2892 value_set_si(inc.d, 1);
2893 value_oppose(inc.x.n, offset);
2894 eadd(&inc, &p->arr[0]);
2895 reorder_terms_about(p, &p->arr[0]); /* frees arr[0] */
2896 value_clear(e->d);
2897 *e = p->arr[1];
2898 free(p);
2899 free_evalue_refs(&inc);
2900 }
2901
2902 int evalue_range_reduction_in_domain(evalue *e, Polyhedron *D)
2903 {
2904 int i;
2905 enode *p;
2906 Value d, min, max;
2907 int r = 0;
2908 Polyhedron *I;
2909 int bounded;
2910
2911 if (value_notzero_p(e->d))
2912 return r;
2913
2914 p = e->x.p;
2915
2916 if (p->type == relation) {
2917 Matrix *T;
2918 int equal;
2919 value_init(d);
2920 value_init(min);
2921 value_init(max);
2922
2923 fractional_minimal_coefficients(p->arr[0].x.p);
2924 I = polynomial_projection(p->arr[0].x.p, D, &d, &T);
2925 bounded = line_minmax(I, &min, &max); /* frees I */
2926 equal = value_eq(min, max);
2927 mpz_cdiv_q(min, min, d);
2928 mpz_fdiv_q(max, max, d);
2929
2930 if (bounded && value_gt(min, max)) {
2931 /* Never zero */
2932 if (p->size == 3) {
2933 value_clear(e->d);
2934 *e = p->arr[2];
2935 } else {
2936 evalue_set_si(e, 0, 1);
2937 r = 1;
2938 }
2939 free_evalue_refs(&(p->arr[1]));
2940 free_evalue_refs(&(p->arr[0]));
2941 free(p);
2942 value_clear(d);
2943 value_clear(min);
2944 value_clear(max);
2945 Matrix_Free(T);
2946 return r ? r : evalue_range_reduction_in_domain(e, D);
2947 } else if (bounded && equal) {
2948 /* Always zero */
2949 if (p->size == 3)
2950 free_evalue_refs(&(p->arr[2]));
2951 value_clear(e->d);
2952 *e = p->arr[1];
2953 free_evalue_refs(&(p->arr[0]));
2954 free(p);
2955 value_clear(d);
2956 value_clear(min);
2957 value_clear(max);
2958 Matrix_Free(T);
2959 return evalue_range_reduction_in_domain(e, D);
2960 } else if (bounded && value_eq(min, max)) {
2961 /* zero for a single value */
2962 Polyhedron *E;
2963 Matrix *M = Matrix_Alloc(1, D->Dimension+2);
2964 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
2965 value_multiply(min, min, d);
2966 value_subtract(M->p[0][D->Dimension+1],
2967 M->p[0][D->Dimension+1], min);
2968 E = DomainAddConstraints(D, M, 0);
2969 value_clear(d);
2970 value_clear(min);
2971 value_clear(max);
2972 Matrix_Free(T);
2973 Matrix_Free(M);
2974 r = evalue_range_reduction_in_domain(&p->arr[1], E);
2975 if (p->size == 3)
2976 r |= evalue_range_reduction_in_domain(&p->arr[2], D);
2977 Domain_Free(E);
2978 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2979 return r;
2980 }
2981
2982 value_clear(d);
2983 value_clear(min);
2984 value_clear(max);
2985 Matrix_Free(T);
2986 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2987 }
2988
2989 i = p->type == relation ? 1 :
2990 p->type == fractional ? 1 : 0;
2991 for (; i<p->size; i++)
2992 r |= evalue_range_reduction_in_domain(&p->arr[i], D);
2993
2994 if (p->type != fractional) {
2995 if (r && p->type == polynomial) {
2996 evalue f;
2997 value_init(f.d);
2998 value_set_si(f.d, 0);
2999 f.x.p = new_enode(polynomial, 2, p->pos);
3000 evalue_set_si(&f.x.p->arr[0], 0, 1);
3001 evalue_set_si(&f.x.p->arr[1], 1, 1);
3002 reorder_terms_about(p, &f);
3003 value_clear(e->d);
3004 *e = p->arr[0];
3005 free(p);
3006 }
3007 return r;
3008 }
3009
3010 value_init(d);
3011 value_init(min);
3012 value_init(max);
3013 fractional_minimal_coefficients(p);
3014 I = polynomial_projection(p, D, &d, NULL);
3015 bounded = line_minmax(I, &min, &max); /* frees I */
3016 mpz_fdiv_q(min, min, d);
3017 mpz_fdiv_q(max, max, d);
3018 value_subtract(d, max, min);
3019
3020 if (bounded && value_eq(min, max)) {
3021 replace_by_affine(e, min);
3022 r = 1;
3023 } else if (bounded && value_one_p(d) && p->size > 3) {
3024 /* replace {g}^2 by -(g-min)^2 + (2{g}+1)*(g-min) - {g}
3025 * See pages 199-200 of PhD thesis.
3026 */
3027 evalue rem;
3028 evalue inc;
3029 evalue t;
3030 evalue f;
3031 evalue factor;
3032 value_init(rem.d);
3033 value_set_si(rem.d, 0);
3034 rem.x.p = new_enode(fractional, 3, -1);
3035 evalue_copy(&rem.x.p->arr[0], &p->arr[0]);
3036 value_clear(rem.x.p->arr[1].d);
3037 value_clear(rem.x.p->arr[2].d);
3038 rem.x.p->arr[1] = p->arr[1];
3039 rem.x.p->arr[2] = p->arr[2];
3040 for (i = 3; i < p->size; ++i)
3041 p->arr[i-2] = p->arr[i];
3042 p->size -= 2;
3043
3044 value_init(inc.d);
3045 value_init(inc.x.n);
3046 value_set_si(inc.d, 1);
3047 value_oppose(inc.x.n, min);
3048
3049 value_init(t.d);
3050 evalue_copy(&t, &p->arr[0]);
3051 eadd(&inc, &t);
3052
3053 value_init(f.d);
3054 value_set_si(f.d, 0);
3055 f.x.p = new_enode(fractional, 3, -1);
3056 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
3057 evalue_set_si(&f.x.p->arr[1], 1, 1);
3058 evalue_set_si(&f.x.p->arr[2], 2, 1);
3059
3060 value_init(factor.d);
3061 evalue_set_si(&factor, -1, 1);
3062 emul(&t, &factor);
3063
3064 eadd(&f, &factor);
3065 emul(&t, &factor);
3066
3067 value_clear(f.x.p->arr[1].x.n);
3068 value_clear(f.x.p->arr[2].x.n);
3069 evalue_set_si(&f.x.p->arr[1], 0, 1);
3070 evalue_set_si(&f.x.p->arr[2], -1, 1);
3071 eadd(&f, &factor);
3072
3073 if (r) {
3074 reorder_terms(&rem);
3075 reorder_terms(e);
3076 }
3077
3078 emul(&factor, e);
3079 eadd(&rem, e);
3080
3081 free_evalue_refs(&inc);
3082 free_evalue_refs(&t);
3083 free_evalue_refs(&f);
3084 free_evalue_refs(&factor);
3085 free_evalue_refs(&rem);
3086
3087 evalue_range_reduction_in_domain(e, D);
3088
3089 r = 1;
3090 } else {
3091 _reduce_evalue(&p->arr[0], 0, 1);
3092 if (r)
3093 reorder_terms(e);
3094 }
3095
3096 value_clear(d);
3097 value_clear(min);
3098 value_clear(max);
3099
3100 return r;
3101 }
3102
3103 void evalue_range_reduction(evalue *e)
3104 {
3105 int i;
3106 if (value_notzero_p(e->d) || e->x.p->type != partition)
3107 return;
3108
3109 for (i = 0; i < e->x.p->size/2; ++i)
3110 if (evalue_range_reduction_in_domain(&e->x.p->arr[2*i+1],
3111 EVALUE_DOMAIN(e->x.p->arr[2*i]))) {
3112 reduce_evalue(&e->x.p->arr[2*i+1]);
3113
3114 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
3115 free_evalue_refs(&e->x.p->arr[2*i+1]);
3116 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
3117 value_clear(e->x.p->arr[2*i].d);
3118 e->x.p->size -= 2;
3119 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
3120 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
3121 --i;
3122 }
3123 }
3124 }
3125
3126 /* Frees EP
3127 */
3128 Enumeration* partition2enumeration(evalue *EP)
3129 {
3130 int i;
3131 Enumeration *en, *res = NULL;
3132
3133 if (EVALUE_IS_ZERO(*EP)) {
3134 free(EP);
3135 return res;
3136 }
3137
3138 for (i = 0; i < EP->x.p->size/2; ++i) {
3139 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[2*i])->Dimension);
3140 en = (Enumeration *)malloc(sizeof(Enumeration));
3141 en->next = res;
3142 res = en;
3143 res->ValidityDomain = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
3144 value_clear(EP->x.p->arr[2*i].d);
3145 res->EP = EP->x.p->arr[2*i+1];
3146 }
3147 free(EP->x.p);
3148 value_clear(EP->d);
3149 free(EP);
3150 return res;
3151 }
3152
3153 int evalue_frac2floor_in_domain3(evalue *e, Polyhedron *D, int shift)
3154 {
3155 enode *p;
3156 int r = 0;
3157 int i;
3158 Polyhedron *I;
3159 Value d, min;
3160 evalue fl;
3161
3162 if (value_notzero_p(e->d))
3163 return r;
3164
3165 p = e->x.p;
3166
3167 i = p->type == relation ? 1 :
3168 p->type == fractional ? 1 : 0;
3169 for (; i<p->size; i++)
3170 r |= evalue_frac2floor_in_domain3(&p->arr[i], D, shift);
3171
3172 if (p->type != fractional) {
3173 if (r && p->type == polynomial) {
3174 evalue f;
3175 value_init(f.d);
3176 value_set_si(f.d, 0);
3177 f.x.p = new_enode(polynomial, 2, p->pos);
3178 evalue_set_si(&f.x.p->arr[0], 0, 1);
3179 evalue_set_si(&f.x.p->arr[1], 1, 1);
3180 reorder_terms_about(p, &f);
3181 value_clear(e->d);
3182 *e = p->arr[0];
3183 free(p);
3184 }
3185 return r;
3186 }
3187
3188 if (shift) {
3189 value_init(d);
3190 I = polynomial_projection(p, D, &d, NULL);
3191
3192 /*
3193 Polyhedron_Print(stderr, P_VALUE_FMT, I);
3194 */
3195
3196 assert(I->NbEq == 0); /* Should have been reduced */
3197
3198 /* Find minimum */
3199 for (i = 0; i < I->NbConstraints; ++i)
3200 if (value_pos_p(I->Constraint[i][1]))
3201 break;
3202
3203 if (i < I->NbConstraints) {
3204 value_init(min);
3205 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
3206 mpz_cdiv_q(min, I->Constraint[i][2], I->Constraint[i][1]);
3207 if (value_neg_p(min)) {
3208 evalue offset;
3209 mpz_fdiv_q(min, min, d);
3210 value_init(offset.d);
3211 value_set_si(offset.d, 1);
3212 value_init(offset.x.n);
3213 value_oppose(offset.x.n, min);
3214 eadd(&offset, &p->arr[0]);
3215 free_evalue_refs(&offset);
3216 }
3217 value_clear(min);
3218 }
3219
3220 Polyhedron_Free(I);
3221 value_clear(d);
3222 }
3223
3224 value_init(fl.d);
3225 value_set_si(fl.d, 0);
3226 fl.x.p = new_enode(flooring, 3, -1);
3227 evalue_set_si(&fl.x.p->arr[1], 0, 1);
3228 evalue_set_si(&fl.x.p->arr[2], -1, 1);
3229 evalue_copy(&fl.x.p->arr[0], &p->arr[0]);
3230
3231 eadd(&fl, &p->arr[0]);
3232 reorder_terms_about(p, &p->arr[0]);
3233 value_clear(e->d);
3234 *e = p->arr[1];
3235 free(p);
3236 free_evalue_refs(&fl);
3237
3238 return 1;
3239 }
3240
3241 int evalue_frac2floor_in_domain(evalue *e, Polyhedron *D)
3242 {
3243 return evalue_frac2floor_in_domain3(e, D, 1);
3244 }
3245
3246 void evalue_frac2floor2(evalue *e, int shift)
3247 {
3248 int i;
3249 if (value_notzero_p(e->d) || e->x.p->type != partition) {
3250 if (!shift) {
3251 if (evalue_frac2floor_in_domain3(e, NULL, 0))
3252 reduce_evalue(e);
3253 }
3254 return;
3255 }
3256
3257 for (i = 0; i < e->x.p->size/2; ++i)
3258 if (evalue_frac2floor_in_domain3(&e->x.p->arr[2*i+1],
3259 EVALUE_DOMAIN(e->x.p->arr[2*i]), shift))
3260 reduce_evalue(&e->x.p->arr[2*i+1]);
3261 }
3262
3263 void evalue_frac2floor(evalue *e)
3264 {
3265 evalue_frac2floor2(e, 1);
3266 }
3267
3268 /* Add a new variable with lower bound 1 and upper bound specified
3269 * by row. If negative is true, then the new variable has upper
3270 * bound -1 and lower bound specified by row.
3271 */
3272 static Matrix *esum_add_constraint(int nvar, Polyhedron *D, Matrix *C,
3273 Vector *row, int negative)
3274 {
3275 int nr, nc;
3276 int i;
3277 int nparam = D->Dimension - nvar;
3278
3279 if (C == 0) {
3280 nr = D->NbConstraints + 2;
3281 nc = D->Dimension + 2 + 1;
3282 C = Matrix_Alloc(nr, nc);
3283 for (i = 0; i < D->NbConstraints; ++i) {
3284 Vector_Copy(D->Constraint[i], C->p[i], 1 + nvar);
3285 Vector_Copy(D->Constraint[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
3286 D->Dimension + 1 - nvar);
3287 }
3288 } else {
3289 Matrix *oldC = C;
3290 nr = C->NbRows + 2;
3291 nc = C->NbColumns + 1;
3292 C = Matrix_Alloc(nr, nc);
3293 for (i = 0; i < oldC->NbRows; ++i) {
3294 Vector_Copy(oldC->p[i], C->p[i], 1 + nvar);
3295 Vector_Copy(oldC->p[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
3296 oldC->NbColumns - 1 - nvar);
3297 }
3298 }
3299 value_set_si(C->p[nr-2][0], 1);
3300 if (negative)
3301 value_set_si(C->p[nr-2][1 + nvar], -1);
3302 else
3303 value_set_si(C->p[nr-2][1 + nvar], 1);
3304 value_set_si(C->p[nr-2][nc - 1], -1);
3305
3306 Vector_Copy(row->p, C->p[nr-1], 1 + nvar + 1);
3307 Vector_Copy(row->p + 1 + nvar + 1, C->p[nr-1] + C->NbColumns - 1 - nparam,
3308 1 + nparam);
3309
3310 return C;
3311 }
3312
3313 static void floor2frac_r(evalue *e, int nvar)
3314 {
3315 enode *p;
3316 int i;
3317 evalue f;
3318 evalue *pp;
3319
3320 if (value_notzero_p(e->d))
3321 return;
3322
3323 p = e->x.p;
3324
3325 assert(p->type == flooring);
3326 for (i = 1; i < p->size; i++)
3327 floor2frac_r(&p->arr[i], nvar);
3328
3329 for (pp = &p->arr[0]; value_zero_p(pp->d); pp = &pp->x.p->arr[0]) {
3330 assert(pp->x.p->type == polynomial);
3331 pp->x.p->pos -= nvar;
3332 }
3333
3334 value_init(f.d);
3335 value_set_si(f.d, 0);
3336 f.x.p = new_enode(fractional, 3, -1);
3337 evalue_set_si(&f.x.p->arr[1], 0, 1);
3338 evalue_set_si(&f.x.p->arr[2], -1, 1);
3339 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
3340
3341 eadd(&f, &p->arr[0]);
3342 reorder_terms_about(p, &p->arr[0]);
3343 value_clear(e->d);
3344 *e = p->arr[1];
3345 free(p);
3346 free_evalue_refs(&f);
3347 }
3348
3349 /* Convert flooring back to fractional and shift position
3350 * of the parameters by nvar
3351 */
3352 static void floor2frac(evalue *e, int nvar)
3353 {
3354 floor2frac_r(e, nvar);
3355 reduce_evalue(e);
3356 }
3357
3358 evalue *esum_over_domain_cst(int nvar, Polyhedron *D, Matrix *C)
3359 {
3360 evalue *t;
3361 int nparam = D->Dimension - nvar;
3362
3363 if (C != 0) {
3364 C = Matrix_Copy(C);
3365 D = Constraints2Polyhedron(C, 0);
3366 Matrix_Free(C);
3367 }
3368
3369 t = barvinok_enumerate_e(D, 0, nparam, 0);
3370
3371 /* Double check that D was not unbounded. */
3372 assert(!(value_pos_p(t->d) && value_neg_p(t->x.n)));
3373
3374 if (C != 0)
3375 Polyhedron_Free(D);
3376
3377 return t;
3378 }
3379
3380 static evalue *esum_over_domain(evalue *e, int nvar, Polyhedron *D,
3381 int *signs, Matrix *C, unsigned MaxRays)
3382 {
3383 Vector *row = NULL;
3384 int i, offset;
3385 evalue *res;
3386 Matrix *origC;
3387 evalue *factor = NULL;
3388 evalue cum;
3389 int negative = 0;
3390
3391 if (EVALUE_IS_ZERO(*e))
3392 return 0;
3393
3394 if (D->next) {
3395 Polyhedron *DD = Disjoint_Domain(D, 0, MaxRays);
3396 Polyhedron *Q;
3397
3398 Q = DD;
3399 DD = Q->next;
3400 Q->next = 0;
3401
3402 res = esum_over_domain(e, nvar, Q, signs, C, MaxRays);
3403 Polyhedron_Free(Q);
3404
3405 for (Q = DD; Q; Q = DD) {
3406 evalue *t;
3407
3408 DD = Q->next;
3409 Q->next = 0;
3410
3411 t = esum_over_domain(e, nvar, Q, signs, C, MaxRays);
3412 Polyhedron_Free(Q);
3413
3414 if (!res)
3415 res = t;
3416 else if (t) {
3417 eadd(t, res);
3418 evalue_free(t);
3419 }
3420 }
3421 return res;
3422 }
3423
3424 if (value_notzero_p(e->d)) {
3425 evalue *t;
3426
3427 t = esum_over_domain_cst(nvar, D, C);
3428
3429 if (!EVALUE_IS_ONE(*e))
3430 emul(e, t);
3431
3432 return t;
3433 }
3434
3435 switch (e->x.p->type) {
3436 case flooring: {
3437 evalue *pp = &e->x.p->arr[0];
3438
3439 if (pp->x.p->pos > nvar) {
3440 /* remainder is independent of the summated vars */
3441 evalue f;
3442 evalue *t;
3443
3444 value_init(f.d);
3445 evalue_copy(&f, e);
3446 floor2frac(&f, nvar);
3447
3448 t = esum_over_domain_cst(nvar, D, C);
3449
3450 emul(&f, t);
3451
3452 free_evalue_refs(&f);
3453
3454 return t;
3455 }
3456
3457 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3458 poly_denom(pp, &row->p[1 + nvar]);
3459 value_set_si(row->p[0], 1);
3460 for (pp = &e->x.p->arr[0]; value_zero_p(pp->d);
3461 pp = &pp->x.p->arr[0]) {
3462 int pos;
3463 assert(pp->x.p->type == polynomial);
3464 pos = pp->x.p->pos;
3465 if (pos >= 1 + nvar)
3466 ++pos;
3467 value_assign(row->p[pos], row->p[1+nvar]);
3468 value_division(row->p[pos], row->p[pos], pp->x.p->arr[1].d);
3469 value_multiply(row->p[pos], row->p[pos], pp->x.p->arr[1].x.n);
3470 }
3471 value_assign(row->p[1 + D->Dimension + 1], row->p[1+nvar]);
3472 value_division(row->p[1 + D->Dimension + 1],
3473 row->p[1 + D->Dimension + 1],
3474 pp->d);
3475 value_multiply(row->p[1 + D->Dimension + 1],
3476 row->p[1 + D->Dimension + 1],
3477 pp->x.n);
3478 value_oppose(row->p[1 + nvar], row->p[1 + nvar]);
3479 break;
3480 }
3481 case polynomial: {
3482 int pos = e->x.p->pos;
3483
3484 if (pos > nvar) {
3485 factor = ALLOC(evalue);
3486 value_init(factor->d);
3487 value_set_si(factor->d, 0);
3488 factor->x.p = new_enode(polynomial, 2, pos - nvar);
3489 evalue_set_si(&factor->x.p->arr[0], 0, 1);
3490 evalue_set_si(&factor->x.p->arr[1], 1, 1);
3491 break;
3492 }
3493
3494 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3495 negative = signs[pos-1] < 0;
3496 value_set_si(row->p[0], 1);
3497 if (negative) {
3498 value_set_si(row->p[pos], -1);
3499 value_set_si(row->p[1 + nvar], 1);
3500 } else {
3501 value_set_si(row->p[pos], 1);
3502 value_set_si(row->p[1 + nvar], -1);
3503 }
3504 break;
3505 }
3506 default:
3507 assert(0);
3508 }
3509
3510 offset = type_offset(e->x.p);
3511
3512 res = esum_over_domain(&e->x.p->arr[offset], nvar, D, signs, C, MaxRays);
3513
3514 if (factor) {
3515 value_init(cum.d);
3516 evalue_copy(&cum, factor);
3517 }
3518
3519 origC = C;
3520 for (i = 1; offset+i < e->x.p->size; ++i) {
3521 evalue *t;
3522 if (row) {
3523 Matrix *prevC = C;
3524 C = esum_add_constraint(nvar, D, C, row, negative);
3525 if (prevC != origC)
3526 Matrix_Free(prevC);
3527 }
3528 /*
3529 if (row)
3530 Vector_Print(stderr, P_VALUE_FMT, row);
3531 if (C)
3532 Matrix_Print(stderr, P_VALUE_FMT, C);
3533 */
3534 t = esum_over_domain(&e->x.p->arr[offset+i], nvar, D, signs, C, MaxRays);
3535
3536 if (t) {
3537 if (factor)
3538 emul(&cum, t);
3539 if (negative && (i % 2))
3540 evalue_negate(t);
3541 }
3542
3543 if (!res)
3544 res = t;
3545 else if (t) {
3546 eadd(t, res);
3547 evalue_free(t);
3548 }
3549 if (factor && offset+i+1 < e->x.p->size)
3550 emul(factor, &cum);
3551 }
3552 if (C != origC)
3553 Matrix_Free(C);
3554
3555 if (factor) {
3556 free_evalue_refs(&cum);
3557 evalue_free(factor);
3558 }
3559
3560 if (row)
3561 Vector_Free(row);
3562
3563 reduce_evalue(res);
3564
3565 return res;
3566 }
3567
3568 static void domain_signs(Polyhedron *D, int *signs)
3569 {
3570 int j, k;
3571
3572 POL_ENSURE_VERTICES(D);
3573 for (j = 0; j < D->Dimension; ++j) {
3574 signs[j] = 0;
3575 for (k = 0; k < D->NbRays; ++k) {
3576 signs[j] = value_sign(D->Ray[k][1+j]);
3577 if (signs[j])
3578 break;
3579 }
3580 }
3581 }
3582
3583 static void shift_floor_in_domain(evalue *e, Polyhedron *D)
3584 {
3585 Value d, m;
3586 Polyhedron *I;
3587 int i;
3588 enode *p;
3589
3590 if (value_notzero_p(e->d))
3591 return;
3592
3593 p = e->x.p;
3594
3595 for (i = type_offset(p); i < p->size; ++i)
3596 shift_floor_in_domain(&p->arr[i], D);
3597
3598 if (p->type != flooring)
3599 return;
3600
3601 value_init(d);
3602 value_init(m);
3603
3604 I = polynomial_projection(p, D, &d, NULL);
3605 assert(I->NbEq == 0); /* Should have been reduced */
3606
3607 for (i = 0; i < I->NbConstraints; ++i)
3608 if (value_pos_p(I->Constraint[i][1]))
3609 break;
3610 assert(i < I->NbConstraints);
3611 if (i < I->NbConstraints) {
3612 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
3613 mpz_fdiv_q(m, I->Constraint[i][2], I->Constraint[i][1]);
3614 if (value_neg_p(m)) {
3615 /* replace [e] by [e-m]+m such that e-m >= 0 */
3616 evalue f;
3617
3618 value_init(f.d);
3619 value_init(f.x.n);
3620 value_set_si(f.d, 1);
3621 value_oppose(f.x.n, m);
3622 eadd(&f, &p->arr[0]);
3623 value_clear(f.x.n);
3624
3625 value_set_si(f.d, 0);
3626 f.x.p = new_enode(flooring, 3, -1);
3627 value_clear(f.x.p->arr[0].d);
3628 f.x.p->arr[0] = p->arr[0];
3629 evalue_set_si(&f.x.p->arr[2], 1, 1);
3630 value_set_si(f.x.p->arr[1].d, 1);
3631 value_init(f.x.p->arr[1].x.n);
3632 value_assign(f.x.p->arr[1].x.n, m);
3633 reorder_terms_about(p, &f);
3634 value_clear(e->d);
3635 *e = p->arr[1];
3636 free(p);
3637 }
3638 }
3639 Polyhedron_Free(I);
3640 value_clear(d);
3641 value_clear(m);
3642 }
3643
3644 /* Make arguments of all floors non-negative */
3645 static void shift_floor_arguments(evalue *e)
3646 {
3647 int i;
3648
3649 if (value_notzero_p(e->d) || e->x.p->type != partition)
3650 return;
3651
3652 for (i = 0; i < e->x.p->size/2; ++i)
3653 shift_floor_in_domain(&e->x.p->arr[2*i+1],
3654 EVALUE_DOMAIN(e->x.p->arr[2*i]));
3655 }
3656
3657 evalue *evalue_sum(evalue *e, int nvar, unsigned MaxRays)
3658 {
3659 int i, dim;
3660 int *signs;
3661 evalue *res = ALLOC(evalue);
3662 value_init(res->d);
3663
3664 assert(nvar >= 0);
3665 if (nvar == 0 || EVALUE_IS_ZERO(*e)) {
3666 evalue_copy(res, e);
3667 return res;
3668 }
3669
3670 evalue_split_domains_into_orthants(e, MaxRays);
3671 evalue_frac2floor2(e, 0);
3672 evalue_set_si(res, 0, 1);
3673
3674 assert(value_zero_p(e->d));
3675 assert(e->x.p->type == partition);
3676 shift_floor_arguments(e);
3677
3678 assert(e->x.p->size >= 2);
3679 dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
3680
3681 signs = alloca(sizeof(int) * dim);
3682
3683 for (i = 0; i < e->x.p->size/2; ++i) {
3684 evalue *t;
3685 domain_signs(EVALUE_DOMAIN(e->x.p->arr[2*i]), signs);
3686 t = esum_over_domain(&e->x.p->arr[2*i+1], nvar,
3687 EVALUE_DOMAIN(e->x.p->arr[2*i]), signs, 0,
3688 MaxRays);
3689 eadd(t, res);
3690 evalue_free(t);
3691 }
3692
3693 reduce_evalue(res);
3694
3695 return res;
3696 }
3697
3698 evalue *esum(evalue *e, int nvar)
3699 {
3700 return evalue_sum(e, nvar, 0);
3701 }
3702
3703 /* Initial silly implementation */
3704 void eor(evalue *e1, evalue *res)
3705 {
3706 evalue E;
3707 evalue mone;
3708 value_init(E.d);
3709 value_init(mone.d);
3710 evalue_set_si(&mone, -1, 1);
3711
3712 evalue_copy(&E, res);
3713 eadd(e1, &E);
3714 emul(e1, res);
3715 emul(&mone, res);
3716 eadd(&E, res);
3717
3718 free_evalue_refs(&E);
3719 free_evalue_refs(&mone);
3720 }
3721
3722 /* computes denominator of polynomial evalue
3723 * d should point to a value initialized to 1
3724 */
3725 void evalue_denom(const evalue *e, Value *d)
3726 {
3727 int i, offset;
3728
3729 if (value_notzero_p(e->d)) {
3730 value_lcm(*d, *d, e->d);
3731 return;
3732 }
3733 assert(e->x.p->type == polynomial ||
3734 e->x.p->type == fractional ||
3735 e->x.p->type == flooring);
3736 offset = type_offset(e->x.p);
3737 for (i = e->x.p->size-1; i >= offset; --i)
3738 evalue_denom(&e->x.p->arr[i], d);
3739 }
3740
3741 /* Divides the evalue e by the integer n */
3742 void evalue_div(evalue *e, Value n)
3743 {
3744 int i, offset;
3745
3746 if (value_one_p(n) || EVALUE_IS_ZERO(*e))
3747 return;
3748
3749 if (value_notzero_p(e->d)) {
3750 Value gc;
3751 value_init(gc);
3752 value_multiply(e->d, e->d, n);
3753 value_gcd(gc, e->x.n, e->d);
3754 if (value_notone_p(gc)) {
3755 value_division(e->d, e->d, gc);
3756 value_division(e->x.n, e->x.n, gc);
3757 }
3758 value_clear(gc);
3759 return;
3760 }
3761 if (e->x.p->type == partition) {
3762 for (i = 0; i < e->x.p->size/2; ++i)
3763 evalue_div(&e->x.p->arr[2*i+1], n);
3764 return;
3765 }
3766 offset = type_offset(e->x.p);
3767 for (i = e->x.p->size-1; i >= offset; --i)
3768 evalue_div(&e->x.p->arr[i], n);
3769 }
3770
3771 /* Multiplies the evalue e by the integer n */
3772 void evalue_mul(evalue *e, Value n)
3773 {
3774 int i, offset;
3775
3776 if (value_one_p(n) || EVALUE_IS_ZERO(*e))
3777 return;
3778
3779 if (value_notzero_p(e->d)) {
3780 Value gc;
3781 value_init(gc);
3782 value_multiply(e->x.n, e->x.n, n);
3783 value_gcd(gc, e->x.n, e->d);
3784 if (value_notone_p(gc)) {
3785 value_division(e->d, e->d, gc);
3786 value_division(e->x.n, e->x.n, gc);
3787 }
3788 value_clear(gc);
3789 return;
3790 }
3791 if (e->x.p->type == partition) {
3792 for (i = 0; i < e->x.p->size/2; ++i)
3793 evalue_mul(&e->x.p->arr[2*i+1], n);
3794 return;
3795 }
3796 offset = type_offset(e->x.p);
3797 for (i = e->x.p->size-1; i >= offset; --i)
3798 evalue_mul(&e->x.p->arr[i], n);
3799 }
3800
3801 /* Multiplies the evalue e by the n/d */
3802 void evalue_mul_div(evalue *e, Value n, Value d)
3803 {
3804 int i, offset;
3805
3806 if ((value_one_p(n) && value_one_p(d)) || EVALUE_IS_ZERO(*e))
3807 return;
3808
3809 if (value_notzero_p(e->d)) {
3810 Value gc;
3811 value_init(gc);
3812 value_multiply(e->x.n, e->x.n, n);
3813 value_multiply(e->d, e->d, d);
3814 value_gcd(gc, e->x.n, e->d);
3815 if (value_notone_p(gc)) {
3816 value_division(e->d, e->d, gc);
3817 value_division(e->x.n, e->x.n, gc);
3818 }
3819 value_clear(gc);
3820 return;
3821 }
3822 if (e->x.p->type == partition) {
3823 for (i = 0; i < e->x.p->size/2; ++i)
3824 evalue_mul_div(&e->x.p->arr[2*i+1], n, d);
3825 return;
3826 }
3827 offset = type_offset(e->x.p);
3828 for (i = e->x.p->size-1; i >= offset; --i)
3829 evalue_mul_div(&e->x.p->arr[i], n, d);
3830 }
3831
3832 void evalue_negate(evalue *e)
3833 {
3834 int i, offset;
3835
3836 if (value_notzero_p(e->d)) {
3837 value_oppose(e->x.n, e->x.n);
3838 return;
3839 }
3840 if (e->x.p->type == partition) {
3841 for (i = 0; i < e->x.p->size/2; ++i)
3842 evalue_negate(&e->x.p->arr[2*i+1]);
3843 return;
3844 }
3845 offset = type_offset(e->x.p);
3846 for (i = e->x.p->size-1; i >= offset; --i)
3847 evalue_negate(&e->x.p->arr[i]);
3848 }
3849
3850 void evalue_add_constant(evalue *e, const Value cst)
3851 {
3852 int i;
3853
3854 if (value_zero_p(e->d)) {
3855 if (e->x.p->type == partition) {
3856 for (i = 0; i < e->x.p->size/2; ++i)
3857 evalue_add_constant(&e->x.p->arr[2*i+1], cst);
3858 return;
3859 }
3860 if (e->x.p->type == relation) {
3861 for (i = 1; i < e->x.p->size; ++i)
3862 evalue_add_constant(&e->x.p->arr[i], cst);
3863 return;
3864 }
3865 do {
3866 e = &e->x.p->arr[type_offset(e->x.p)];
3867 } while (value_zero_p(e->d));
3868 }
3869 value_addmul(e->x.n, cst, e->d);
3870 }
3871
3872 static void evalue_frac2polynomial_r(evalue *e, int *signs, int sign, int in_frac)
3873 {
3874 int i, offset;
3875 Value d;
3876 enode *p;
3877 int sign_odd = sign;
3878
3879 if (value_notzero_p(e->d)) {
3880 if (in_frac && sign * value_sign(e->x.n) < 0) {
3881 value_set_si(e->x.n, 0);
3882 value_set_si(e->d, 1);
3883 }
3884 return;
3885 }
3886
3887 if (e->x.p->type == relation) {
3888 for (i = e->x.p->size-1; i >= 1; --i)
3889 evalue_frac2polynomial_r(&e->x.p->arr[i], signs, sign, in_frac);
3890 return;
3891 }
3892
3893 if (e->x.p->type == polynomial)
3894 sign_odd *= signs[e->x.p->pos-1];
3895 offset = type_offset(e->x.p);
3896 evalue_frac2polynomial_r(&e->x.p->arr[offset], signs, sign, in_frac);
3897 in_frac |= e->x.p->type == fractional;
3898 for (i = e->x.p->size-1; i > offset; --i)
3899 evalue_frac2polynomial_r(&e->x.p->arr[i], signs,
3900 (i - offset) % 2 ? sign_odd : sign, in_frac);
3901
3902 if (e->x.p->type != fractional)
3903 return;
3904
3905 /* replace { a/m } by (m-1)/m if sign != 0
3906 * and by (m-1)/(2m) if sign == 0
3907 */
3908 value_init(d);
3909 value_set_si(d, 1);
3910 evalue_denom(&e->x.p->arr[0], &d);
3911 free_evalue_refs(&e->x.p->arr[0]);
3912 value_init(e->x.p->arr[0].d);
3913 value_init(e->x.p->arr[0].x.n);
3914 if (sign == 0)
3915 value_addto(e->x.p->arr[0].d, d, d);
3916 else
3917 value_assign(e->x.p->arr[0].d, d);
3918 value_decrement(e->x.p->arr[0].x.n, d);
3919 value_clear(d);
3920
3921 p = e->x.p;
3922 reorder_terms_about(p, &p->arr[0]);
3923 value_clear(e->d);
3924 *e = p->arr[1];
3925 free(p);
3926 }
3927
3928 /* Approximate the evalue in fractional representation by a polynomial.
3929 * If sign > 0, the result is an upper bound;
3930 * if sign < 0, the result is a lower bound;
3931 * if sign = 0, the result is an intermediate approximation.
3932 */
3933 void evalue_frac2polynomial(evalue *e, int sign, unsigned MaxRays)
3934 {
3935 int i, dim;
3936 int *signs;
3937
3938 if (value_notzero_p(e->d))
3939 return;
3940 assert(e->x.p->type == partition);
3941 /* make sure all variables in the domains have a fixed sign */
3942 if (sign) {
3943 evalue_split_domains_into_orthants(e, MaxRays);
3944 if (EVALUE_IS_ZERO(*e))
3945 return;
3946 }
3947
3948 assert(e->x.p->size >= 2);
3949 dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
3950
3951 signs = alloca(sizeof(int) * dim);
3952
3953 if (!sign)
3954 for (i = 0; i < dim; ++i)
3955 signs[i] = 0;
3956 for (i = 0; i < e->x.p->size/2; ++i) {
3957 if (sign)
3958 domain_signs(EVALUE_DOMAIN(e->x.p->arr[2*i]), signs);
3959 evalue_frac2polynomial_r(&e->x.p->arr[2*i+1], signs, sign, 0);
3960 }
3961 }
3962
3963 /* Split the domains of e (which is assumed to be a partition)
3964 * such that each resulting domain lies entirely in one orthant.
3965 */
3966 void evalue_split_domains_into_orthants(evalue *e, unsigned MaxRays)
3967 {
3968 int i, dim;
3969 assert(value_zero_p(e->d));
3970 assert(e->x.p->type == partition);
3971 assert(e->x.p->size >= 2);
3972 dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
3973
3974 for (i = 0; i < dim; ++i) {
3975 evalue split;
3976 Matrix *C, *C2;
3977 C = Matrix_Alloc(1, 1 + dim + 1);
3978 value_set_si(C->p[0][0], 1);
3979 value_init(split.d);
3980 value_set_si(split.d, 0);
3981 split.x.p = new_enode(partition, 4, dim);
3982 value_set_si(C->p[0][1+i], 1);
3983 C2 = Matrix_Copy(C);
3984 EVALUE_SET_DOMAIN(split.x.p->arr[0], Constraints2Polyhedron(C2, MaxRays));
3985 Matrix_Free(C2);
3986 evalue_set_si(&split.x.p->arr[1], 1, 1);
3987 value_set_si(C->p[0][1+i], -1);
3988 value_set_si(C->p[0][1+dim], -1);
3989 EVALUE_SET_DOMAIN(split.x.p->arr[2], Constraints2Polyhedron(C, MaxRays));
3990 evalue_set_si(&split.x.p->arr[3], 1, 1);
3991 emul(&split, e);
3992 free_evalue_refs(&split);
3993 Matrix_Free(C);
3994 }
3995 reduce_evalue(e);
3996 }
3997
3998 static evalue *find_fractional_with_max_periods(evalue *e, Polyhedron *D,
3999 int max_periods,
4000 Matrix **TT,
4001 Value *min, Value *max)
4002 {
4003 Matrix *T;
4004 evalue *f = NULL;
4005 Value d;
4006 int i;
4007
4008 if (value_notzero_p(e->d))
4009 return NULL;
4010
4011 if (e->x.p->type == fractional) {
4012 Polyhedron *I;
4013 int bounded;
4014
4015 value_init(d);
4016 I = polynomial_projection(e->x.p, D, &d, &T);
4017 bounded = line_minmax(I, min, max); /* frees I */
4018 if (bounded) {
4019 Value mp;
4020 value_init(mp);
4021 value_set_si(mp, max_periods);
4022 mpz_fdiv_q(*min, *min, d);
4023 mpz_fdiv_q(*max, *max, d);
4024 value_assign(T->p[1][D->Dimension], d);
4025 value_subtract(d, *max, *min);
4026 if (value_ge(d, mp))
4027 Matrix_Free(T);
4028 else
4029 f = evalue_dup(&e->x.p->arr[0]);
4030 value_clear(mp);
4031 } else
4032 Matrix_Free(T);
4033 value_clear(d);
4034 if (f) {
4035 *TT = T;
4036 return f;
4037 }
4038 }
4039
4040 for (i = type_offset(e->x.p); i < e->x.p->size; ++i)
4041 if ((f = find_fractional_with_max_periods(&e->x.p->arr[i], D, max_periods,
4042 TT, min, max)))
4043 return f;
4044
4045 return NULL;
4046 }
4047
4048 static void replace_fract_by_affine(evalue *e, evalue *f, Value val)
4049 {
4050 int i, offset;
4051
4052 if (value_notzero_p(e->d))
4053 return;
4054
4055 offset = type_offset(e->x.p);
4056 for (i = e->x.p->size-1; i >= offset; --i)
4057 replace_fract_by_affine(&e->x.p->arr[i], f, val);
4058
4059 if (e->x.p->type != fractional)
4060 return;
4061
4062 if (!eequal(&e->x.p->arr[0], f))
4063 return;
4064
4065 replace_by_affine(e, val);
4066 }
4067
4068 /* Look for fractional parts that can be removed by splitting the corresponding
4069 * domain into at most max_periods parts.
4070 * We use a very simply strategy that looks for the first fractional part
4071 * that satisfies the condition, performs the split and then continues
4072 * looking for other fractional parts in the split domains until no
4073 * such fractional part can be found anymore.
4074 */
4075 void evalue_split_periods(evalue *e, int max_periods, unsigned int MaxRays)
4076 {
4077 int i, j, n;
4078 Value min;
4079 Value max;
4080 Value d;
4081
4082 if (EVALUE_IS_ZERO(*e))
4083 return;
4084 if (value_notzero_p(e->d) || e->x.p->type != partition) {
4085 fprintf(stderr,
4086 "WARNING: evalue_split_periods called on incorrect evalue type\n");
4087 return;
4088 }
4089
4090 value_init(min);
4091 value_init(max);
4092 value_init(d);
4093
4094 for (i = 0; i < e->x.p->size/2; ++i) {
4095 enode *p;
4096 evalue *f;
4097 Matrix *T;
4098 Matrix *M;
4099 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
4100 Polyhedron *E;
4101 f = find_fractional_with_max_periods(&e->x.p->arr[2*i+1], D, max_periods,
4102 &T, &min, &max);
4103 if (!f)
4104 continue;
4105
4106 M = Matrix_Alloc(2, 2+D->Dimension);
4107
4108 value_subtract(d, max, min);
4109 n = VALUE_TO_INT(d)+1;
4110
4111 value_set_si(M->p[0][0], 1);
4112 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
4113 value_multiply(d, max, T->p[1][D->Dimension]);
4114 value_subtract(M->p[0][1+D->Dimension], M->p[0][1+D->Dimension], d);
4115 value_set_si(d, -1);
4116 value_set_si(M->p[1][0], 1);
4117 Vector_Scale(T->p[0], M->p[1]+1, d, D->Dimension+1);
4118 value_addmul(M->p[1][1+D->Dimension], max, T->p[1][D->Dimension]);
4119 value_addto(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension],
4120 T->p[1][D->Dimension]);
4121 value_decrement(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension]);
4122
4123 p = new_enode(partition, e->x.p->size + (n-1)*2, e->x.p->pos);
4124 for (j = 0; j < 2*i; ++j) {
4125 value_clear(p->arr[j].d);
4126 p->arr[j] = e->x.p->arr[j];
4127 }
4128 for (j = 2*i+2; j < e->x.p->size; ++j) {
4129 value_clear(p->arr[j+2*(n-1)].d);
4130 p->arr[j+2*(n-1)] = e->x.p->arr[j];
4131 }
4132 for (j = n-1; j >= 0; --j) {
4133 if (j == 0) {
4134 value_clear(p->arr[2*i+1].d);
4135 p->arr[2*i+1] = e->x.p->arr[2*i+1];
4136 } else
4137 evalue_copy(&p->arr[2*(i+j)+1], &e->x.p->arr[2*i+1]);
4138 if (j != n-1) {
4139 value_subtract(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension],
4140 T->p[1][D->Dimension]);
4141 value_addto(M->p[0][1+D->Dimension], M->p[0][1+D->Dimension],
4142 T->p[1][D->Dimension]);
4143 }
4144 replace_fract_by_affine(&p->arr[2*(i+j)+1], f, max);
4145 E = DomainAddConstraints(D, M, MaxRays);
4146 EVALUE_SET_DOMAIN(p->arr[2*(i+j)], E);
4147 if (evalue_range_reduction_in_domain(&p->arr[2*(i+j)+1], E))
4148 reduce_evalue(&p->arr[2*(i+j)+1]);
4149 value_decrement(max, max);
4150 }
4151 value_clear(e->x.p->arr[2*i].d);
4152 Domain_Free(D);
4153 Matrix_Free(M);
4154 Matrix_Free(T);
4155 evalue_free(f);
4156 free(e->x.p);
4157 e->x.p = p;
4158 --i;
4159 }
4160
4161 value_clear(d);
4162 value_clear(min);
4163 value_clear(max);
4164 }
4165
4166 void evalue_extract_affine(const evalue *e, Value *coeff, Value *cst, Value *d)
4167 {
4168 value_set_si(*d, 1);
4169 evalue_denom(e, d);
4170 for ( ; value_zero_p(e->d); e = &e->x.p->arr[0]) {
4171 evalue *c;
4172 assert(e->x.p->type == polynomial);
4173 assert(e->x.p->size == 2);
4174 c = &e->x.p->arr[1];
4175 value_multiply(coeff[e->x.p->pos-1], *d, c->x.n);
4176 value_division(coeff[e->x.p->pos-1], coeff[e->x.p->pos-1], c->d);
4177 }
4178 value_multiply(*cst, *d, e->x.n);
4179 value_division(*cst, *cst, e->d);
4180 }
4181
4182 /* returns an evalue that corresponds to
4183 *
4184 * c/den x_param
4185 */
4186 static evalue *term(int param, Value c, Value den)
4187 {
4188 evalue *EP = ALLOC(evalue);
4189 value_init(EP->d);
4190 value_set_si(EP->d,0);
4191 EP->x.p = new_enode(polynomial, 2, param + 1);
4192 evalue_set_si(&EP->x.p->arr[0], 0, 1);
4193 value_init(EP->x.p->arr[1].x.n);
4194 value_assign(EP->x.p->arr[1].d, den);
4195 value_assign(EP->x.p->arr[1].x.n, c);
4196 return EP;
4197 }
4198
4199 evalue *affine2evalue(Value *coeff, Value denom, int nvar)
4200 {
4201 int i;
4202 evalue *E = ALLOC(evalue);
4203 value_init(E->d);
4204 evalue_set(E, coeff[nvar], denom);
4205 for (i = 0; i < nvar; ++i) {
4206 evalue *t;
4207 if (value_zero_p(coeff[i]))
4208 continue;
4209 t = term(i, coeff[i], denom);
4210 eadd(t, E);
4211 evalue_free(t);
4212 }
4213 return E;
4214 }
4215
4216 void evalue_substitute(evalue *e, evalue **subs)
4217 {
4218 int i, offset;
4219 evalue *v;
4220 enode *p;
4221
4222 if (value_notzero_p(e->d))
4223 return;
4224
4225 p = e->x.p;
4226 assert(p->type != partition);
4227
4228 for (i = 0; i < p->size; ++i)
4229 evalue_substitute(&p->arr[i], subs);
4230
4231 if (p->type == polynomial)
4232 v = subs[p->pos-1];
4233 else {
4234 v = ALLOC(evalue);
4235 value_init(v->d);
4236 value_set_si(v->d, 0);
4237 v->x.p = new_enode(p->type, 3, -1);
4238 value_clear(v->x.p->arr[0].d);
4239 v->x.p->arr[0] = p->arr[0];
4240 evalue_set_si(&v->x.p->arr[1], 0, 1);
4241 evalue_set_si(&v->x.p->arr[2], 1, 1);
4242 }
4243
4244 offset = type_offset(p);
4245
4246 for (i = p->size-1; i >= offset+1; i--) {
4247 emul(v, &p->arr[i]);
4248 eadd(&p->arr[i], &p->arr[i-1]);
4249 free_evalue_refs(&(p->arr[i]));
4250 }
4251
4252 if (p->type != polynomial)
4253 evalue_free(v);
4254
4255 value_clear(e->d);
4256 *e = p->arr[offset];
4257 free(p);
4258 }
4259
4260 /* evalue e is given in terms of "new" parameter; CP maps the new
4261 * parameters back to the old parameters.
4262 * Transforms e such that it refers back to the old parameters.
4263 */
4264 void evalue_backsubstitute(evalue *e, Matrix *CP, unsigned MaxRays)
4265 {
4266 Matrix *eq;
4267 Matrix *inv;
4268 evalue **subs;
4269 enode *p;
4270 int i;
4271 unsigned nparam = CP->NbColumns-1;
4272 Polyhedron *CEq;
4273
4274 if (EVALUE_IS_ZERO(*e))
4275 return;
4276
4277 assert(value_zero_p(e->d));
4278 p = e->x.p;
4279 assert(p->type == partition);
4280
4281 inv = left_inverse(CP, &eq);
4282 subs = ALLOCN(evalue *, nparam);
4283 for (i = 0; i < nparam; ++i)
4284 subs[i] = affine2evalue(inv->p[i], inv->p[nparam][inv->NbColumns-1],
4285 inv->NbColumns-1);
4286
4287 CEq = Constraints2Polyhedron(eq, MaxRays);
4288 addeliminatedparams_partition(p, inv, CEq, inv->NbColumns-1, MaxRays);
4289 Polyhedron_Free(CEq);
4290
4291 for (i = 0; i < p->size/2; ++i)
4292 evalue_substitute(&p->arr[2*i+1], subs);
4293
4294 for (i = 0; i < nparam; ++i)
4295 evalue_free(subs[i]);
4296 free(subs);
4297 Matrix_Free(eq);
4298 Matrix_Free(inv);
4299 }
4300
4301 /* Computes
4302 *
4303 * \sum_{i=0}^n c_i/d X^i
4304 *
4305 * where d is the last element in the vector c.
4306 */
4307 evalue *evalue_polynomial(Vector *c, const evalue* X)
4308 {
4309 unsigned dim = c->Size-2;
4310 evalue EC;
4311 evalue *EP = ALLOC(evalue);
4312 int i;
4313
4314 value_init(EC.d);
4315 evalue_set(&EC, c->p[dim], c->p[dim+1]);
4316
4317 value_init(EP->d);
4318 evalue_set(EP, c->p[dim], c->p[dim+1]);
4319
4320 for (i = dim-1; i >= 0; --i) {
4321 emul(X, EP);
4322 value_assign(EC.x.n, c->p[i]);
4323 eadd(&EC, EP);
4324 }
4325 free_evalue_refs(&EC);
4326 return EP;
4327 }
4328
4329 /* Create an evalue from an array of pairs of domains and evalues. */
4330 evalue *evalue_from_section_array(struct evalue_section *s, int n)
4331 {
4332 int i;
4333 evalue *res;
4334
4335 res = ALLOC(evalue);
4336 value_init(res->d);
4337
4338 if (n == 0)
4339 evalue_set_si(res, 0, 1);
4340 else {
4341 value_set_si(res->d, 0);
4342 res->x.p = new_enode(partition, 2*n, s[0].D->Dimension);
4343 for (i = 0; i < n; ++i) {
4344 EVALUE_SET_DOMAIN(res->x.p->arr[2*i], s[i].D);
4345 value_clear(res->x.p->arr[2*i+1].d);
4346 res->x.p->arr[2*i+1] = *s[i].E;
4347 free(s[i].E);
4348 }
4349 }
4350 return res;
4351 }
lattice point enumeration library