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