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