?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
barvinok_maximize: barf on unexpected first token
[barvinok.git] / evalue.c
blob3ad59c4a111c22aaf83c61f4842ea1d79266e21b
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 static void print_evalue_r(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_evalue(FILE *DST, const evalue *e, char **pname)
982 {
983 print_evalue_r(DST, e, pname);
984 if (value_notzero_p(e->d))
985 fprintf(DST, "\n");
986 }
987
988 void print_enode(FILE *DST,enode *p,char **pname) {
989
990 int i;
991
992 if (!p) {
993 fprintf(DST, "NULL");
994 return;
995 }
996 switch (p->type) {
997 case evector:
998 fprintf(DST, "{ ");
999 for (i=0; i<p->size; i++) {
1000 print_evalue_r(DST, &p->arr[i], pname);
1001 if (i!=(p->size-1))
1002 fprintf(DST, ", ");
1003 }
1004 fprintf(DST, " }\n");
1005 break;
1006 case polynomial:
1007 fprintf(DST, "( ");
1008 for (i=p->size-1; i>=0; i--) {
1009 print_evalue_r(DST, &p->arr[i], pname);
1010 if (i==1) fprintf(DST, " * %s + ", pname[p->pos-1]);
1011 else if (i>1)
1012 fprintf(DST, " * %s^%d + ", pname[p->pos-1], i);
1013 }
1014 fprintf(DST, " )\n");
1015 break;
1016 case periodic:
1017 fprintf(DST, "[ ");
1018 for (i=0; i<p->size; i++) {
1019 print_evalue_r(DST, &p->arr[i], pname);
1020 if (i!=(p->size-1)) fprintf(DST, ", ");
1021 }
1022 fprintf(DST," ]_%s", pname[p->pos-1]);
1023 break;
1024 case flooring:
1025 case fractional:
1026 fprintf(DST, "( ");
1027 for (i=p->size-1; i>=1; i--) {
1028 print_evalue_r(DST, &p->arr[i], pname);
1029 if (i >= 2) {
1030 fprintf(DST, " * ");
1031 fprintf(DST, p->type == flooring ? "[" : "{");
1032 print_evalue_r(DST, &p->arr[0], pname);
1033 fprintf(DST, p->type == flooring ? "]" : "}");
1034 if (i>2)
1035 fprintf(DST, "^%d + ", i-1);
1036 else
1037 fprintf(DST, " + ");
1038 }
1039 }
1040 fprintf(DST, " )\n");
1041 break;
1042 case relation:
1043 fprintf(DST, "[ ");
1044 print_evalue_r(DST, &p->arr[0], pname);
1045 fprintf(DST, "= 0 ] * \n");
1046 print_evalue_r(DST, &p->arr[1], pname);
1047 if (p->size > 2) {
1048 fprintf(DST, " +\n [ ");
1049 print_evalue_r(DST, &p->arr[0], pname);
1050 fprintf(DST, "!= 0 ] * \n");
1051 print_evalue_r(DST, &p->arr[2], pname);
1052 }
1053 break;
1054 case partition: {
1055 char **names = pname;
1056 int maxdim = EVALUE_DOMAIN(p->arr[0])->Dimension;
1057 if (!pname || p->pos < maxdim) {
1058 NALLOC(names, maxdim);
1059 for (i = 0; i < p->pos; ++i) {
1060 if (pname)
1061 names[i] = pname[i];
1062 else {
1063 NALLOC(names[i], 10);
1064 snprintf(names[i], 10, "%c", 'P'+i);
1065 }
1066 }
1067 for ( ; i < maxdim; ++i) {
1068 NALLOC(names[i], 10);
1069 snprintf(names[i], 10, "_p%d", i);
1070 }
1071 }
1072
1073 for (i=0; i<p->size/2; i++) {
1074 Print_Domain(DST, EVALUE_DOMAIN(p->arr[2*i]), names);
1075 print_evalue_r(DST, &p->arr[2*i+1], names);
1076 }
1077
1078 if (!pname || p->pos < maxdim) {
1079 for (i = pname ? p->pos : 0; i < maxdim; ++i)
1080 free(names[i]);
1081 free(names);
1082 }
1083
1084 break;
1085 }
1086 default:
1087 assert(0);
1088 }
1089 return;
1090 } /* print_enode */
1091
1092 static void eadd_rev(const evalue *e1, evalue *res)
1093 {
1094 evalue ev;
1095 value_init(ev.d);
1096 evalue_copy(&ev, e1);
1097 eadd(res, &ev);
1098 free_evalue_refs(res);
1099 *res = ev;
1100 }
1101
1102 static void eadd_rev_cst(const evalue *e1, evalue *res)
1103 {
1104 evalue ev;
1105 value_init(ev.d);
1106 evalue_copy(&ev, e1);
1107 eadd(res, &ev.x.p->arr[type_offset(ev.x.p)]);
1108 free_evalue_refs(res);
1109 *res = ev;
1110 }
1111
1112 static int is_zero_on(evalue *e, Polyhedron *D)
1113 {
1114 int is_zero;
1115 evalue tmp;
1116 value_init(tmp.d);
1117 tmp.x.p = new_enode(partition, 2, D->Dimension);
1118 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Domain_Copy(D));
1119 evalue_copy(&tmp.x.p->arr[1], e);
1120 reduce_evalue(&tmp);
1121 is_zero = EVALUE_IS_ZERO(tmp);
1122 free_evalue_refs(&tmp);
1123 return is_zero;
1124 }
1125
1126 struct section { Polyhedron * D; evalue E; };
1127
1128 void eadd_partitions(const evalue *e1, evalue *res)
1129 {
1130 int n, i, j;
1131 Polyhedron *d, *fd;
1132 struct section *s;
1133 s = (struct section *)
1134 malloc((e1->x.p->size/2+1) * (res->x.p->size/2+1) *
1135 sizeof(struct section));
1136 assert(s);
1137 assert(e1->x.p->pos == res->x.p->pos);
1138 assert(e1->x.p->pos == EVALUE_DOMAIN(e1->x.p->arr[0])->Dimension);
1139 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1140
1141 n = 0;
1142 for (j = 0; j < e1->x.p->size/2; ++j) {
1143 assert(res->x.p->size >= 2);
1144 fd = DomainDifference(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1145 EVALUE_DOMAIN(res->x.p->arr[0]), 0);
1146 if (!emptyQ(fd))
1147 for (i = 1; i < res->x.p->size/2; ++i) {
1148 Polyhedron *t = fd;
1149 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1150 Domain_Free(t);
1151 if (emptyQ(fd))
1152 break;
1153 }
1154 if (emptyQ(fd)) {
1155 Domain_Free(fd);
1156 continue;
1157 }
1158 /* See if we can extend one of the domains in res to cover fd */
1159 for (i = 0; i < res->x.p->size/2; ++i)
1160 if (is_zero_on(&res->x.p->arr[2*i+1], fd))
1161 break;
1162 if (i < res->x.p->size/2) {
1163 EVALUE_SET_DOMAIN(res->x.p->arr[2*i],
1164 DomainConcat(fd, EVALUE_DOMAIN(res->x.p->arr[2*i])));
1165 continue;
1166 }
1167 value_init(s[n].E.d);
1168 evalue_copy(&s[n].E, &e1->x.p->arr[2*j+1]);
1169 s[n].D = fd;
1170 ++n;
1171 }
1172 for (i = 0; i < res->x.p->size/2; ++i) {
1173 fd = EVALUE_DOMAIN(res->x.p->arr[2*i]);
1174 for (j = 0; j < e1->x.p->size/2; ++j) {
1175 Polyhedron *t;
1176 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1177 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1178 if (emptyQ(d)) {
1179 Domain_Free(d);
1180 continue;
1181 }
1182 t = fd;
1183 fd = DomainDifference(fd, EVALUE_DOMAIN(e1->x.p->arr[2*j]), 0);
1184 if (t != EVALUE_DOMAIN(res->x.p->arr[2*i]))
1185 Domain_Free(t);
1186 value_init(s[n].E.d);
1187 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
1188 eadd(&e1->x.p->arr[2*j+1], &s[n].E);
1189 if (!emptyQ(fd) && is_zero_on(&e1->x.p->arr[2*j+1], fd)) {
1190 d = DomainConcat(fd, d);
1191 fd = Empty_Polyhedron(fd->Dimension);
1192 }
1193 s[n].D = d;
1194 ++n;
1195 }
1196 if (!emptyQ(fd)) {
1197 s[n].E = res->x.p->arr[2*i+1];
1198 s[n].D = fd;
1199 ++n;
1200 } else {
1201 free_evalue_refs(&res->x.p->arr[2*i+1]);
1202 Domain_Free(fd);
1203 }
1204 if (fd != EVALUE_DOMAIN(res->x.p->arr[2*i]))
1205 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
1206 value_clear(res->x.p->arr[2*i].d);
1207 }
1208
1209 free(res->x.p);
1210 assert(n > 0);
1211 res->x.p = new_enode(partition, 2*n, e1->x.p->pos);
1212 for (j = 0; j < n; ++j) {
1213 s[j].D = DomainConstraintSimplify(s[j].D, 0);
1214 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1215 value_clear(res->x.p->arr[2*j+1].d);
1216 res->x.p->arr[2*j+1] = s[j].E;
1217 }
1218
1219 free(s);
1220 }
1221
1222 static void explicit_complement(evalue *res)
1223 {
1224 enode *rel = new_enode(relation, 3, 0);
1225 assert(rel);
1226 value_clear(rel->arr[0].d);
1227 rel->arr[0] = res->x.p->arr[0];
1228 value_clear(rel->arr[1].d);
1229 rel->arr[1] = res->x.p->arr[1];
1230 value_set_si(rel->arr[2].d, 1);
1231 value_init(rel->arr[2].x.n);
1232 value_set_si(rel->arr[2].x.n, 0);
1233 free(res->x.p);
1234 res->x.p = rel;
1235 }
1236
1237 void eadd(const evalue *e1, evalue *res)
1238 {
1239 int i;
1240 if (value_notzero_p(e1->d) && value_notzero_p(res->d)) {
1241 /* Add two rational numbers */
1242 Value g,m1,m2;
1243 value_init(g);
1244 value_init(m1);
1245 value_init(m2);
1246
1247 value_multiply(m1,e1->x.n,res->d);
1248 value_multiply(m2,res->x.n,e1->d);
1249 value_addto(res->x.n,m1,m2);
1250 value_multiply(res->d,e1->d,res->d);
1251 Gcd(res->x.n,res->d,&g);
1252 if (value_notone_p(g)) {
1253 value_division(res->d,res->d,g);
1254 value_division(res->x.n,res->x.n,g);
1255 }
1256 value_clear(g); value_clear(m1); value_clear(m2);
1257 return ;
1258 }
1259 else if (value_notzero_p(e1->d) && value_zero_p(res->d)) {
1260 switch (res->x.p->type) {
1261 case polynomial:
1262 /* Add the constant to the constant term of a polynomial*/
1263 eadd(e1, &res->x.p->arr[0]);
1264 return ;
1265 case periodic:
1266 /* Add the constant to all elements of a periodic number */
1267 for (i=0; i<res->x.p->size; i++) {
1268 eadd(e1, &res->x.p->arr[i]);
1269 }
1270 return ;
1271 case evector:
1272 fprintf(stderr, "eadd: cannot add const with vector\n");
1273 return;
1274 case flooring:
1275 case fractional:
1276 eadd(e1, &res->x.p->arr[1]);
1277 return ;
1278 case partition:
1279 assert(EVALUE_IS_ZERO(*e1));
1280 break; /* Do nothing */
1281 case relation:
1282 /* Create (zero) complement if needed */
1283 if (res->x.p->size < 3 && !EVALUE_IS_ZERO(*e1))
1284 explicit_complement(res);
1285 for (i = 1; i < res->x.p->size; ++i)
1286 eadd(e1, &res->x.p->arr[i]);
1287 break;
1288 default:
1289 assert(0);
1290 }
1291 }
1292 /* add polynomial or periodic to constant
1293 * you have to exchange e1 and res, before doing addition */
1294
1295 else if (value_zero_p(e1->d) && value_notzero_p(res->d)) {
1296 eadd_rev(e1, res);
1297 return;
1298 }
1299 else { // ((e1->d==0) && (res->d==0))
1300 assert(!((e1->x.p->type == partition) ^
1301 (res->x.p->type == partition)));
1302 if (e1->x.p->type == partition) {
1303 eadd_partitions(e1, res);
1304 return;
1305 }
1306 if (e1->x.p->type == relation &&
1307 (res->x.p->type != relation ||
1308 mod_term_smaller(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1309 eadd_rev(e1, res);
1310 return;
1311 }
1312 if (res->x.p->type == relation) {
1313 if (e1->x.p->type == relation &&
1314 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1315 if (res->x.p->size < 3 && e1->x.p->size == 3)
1316 explicit_complement(res);
1317 for (i = 1; i < e1->x.p->size; ++i)
1318 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1319 return;
1320 }
1321 if (res->x.p->size < 3)
1322 explicit_complement(res);
1323 for (i = 1; i < res->x.p->size; ++i)
1324 eadd(e1, &res->x.p->arr[i]);
1325 return;
1326 }
1327 if ((e1->x.p->type != res->x.p->type) ) {
1328 /* adding to evalues of different type. two cases are possible
1329 * res is periodic and e1 is polynomial, you have to exchange
1330 * e1 and res then to add e1 to the constant term of res */
1331 if (e1->x.p->type == polynomial) {
1332 eadd_rev_cst(e1, res);
1333 }
1334 else if (res->x.p->type == polynomial) {
1335 /* res is polynomial and e1 is periodic,
1336 add e1 to the constant term of res */
1337
1338 eadd(e1,&res->x.p->arr[0]);
1339 } else
1340 assert(0);
1341
1342 return;
1343 }
1344 else if (e1->x.p->pos != res->x.p->pos ||
1345 ((res->x.p->type == fractional ||
1346 res->x.p->type == flooring) &&
1347 !eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1348 /* adding evalues of different position (i.e function of different unknowns
1349 * to case are possible */
1350
1351 switch (res->x.p->type) {
1352 case flooring:
1353 case fractional:
1354 if (mod_term_smaller(res, e1))
1355 eadd(e1,&res->x.p->arr[1]);
1356 else
1357 eadd_rev_cst(e1, res);
1358 return;
1359 case polynomial: // res and e1 are polynomials
1360 // add e1 to the constant term of res
1361
1362 if(res->x.p->pos < e1->x.p->pos)
1363 eadd(e1,&res->x.p->arr[0]);
1364 else
1365 eadd_rev_cst(e1, res);
1366 // value_clear(g); value_clear(m1); value_clear(m2);
1367 return;
1368 case periodic: // res and e1 are pointers to periodic numbers
1369 //add e1 to all elements of res
1370
1371 if(res->x.p->pos < e1->x.p->pos)
1372 for (i=0;i<res->x.p->size;i++) {
1373 eadd(e1,&res->x.p->arr[i]);
1374 }
1375 else
1376 eadd_rev(e1, res);
1377 return;
1378 default:
1379 assert(0);
1380 }
1381 }
1382
1383
1384 //same type , same pos and same size
1385 if (e1->x.p->size == res->x.p->size) {
1386 // add any element in e1 to the corresponding element in res
1387 i = type_offset(res->x.p);
1388 if (i == 1)
1389 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
1390 for (; i<res->x.p->size; i++) {
1391 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1392 }
1393 return ;
1394 }
1395
1396 /* Sizes are different */
1397 switch(res->x.p->type) {
1398 case polynomial:
1399 case flooring:
1400 case fractional:
1401 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1402 /* new enode and add res to that new node. If you do not do */
1403 /* that, you lose the the upper weight part of e1 ! */
1404
1405 if(e1->x.p->size > res->x.p->size)
1406 eadd_rev(e1, res);
1407 else {
1408 i = type_offset(res->x.p);
1409 if (i == 1)
1410 assert(eequal(&e1->x.p->arr[0],
1411 &res->x.p->arr[0]));
1412 for (; i<e1->x.p->size ; i++) {
1413 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1414 }
1415
1416 return ;
1417 }
1418 break;
1419
1420 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1421 case periodic:
1422 {
1423 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1424 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1425 to add periodicaly elements of e1 to elements of ne, and finaly to
1426 return ne. */
1427 int x,y,p;
1428 Value ex, ey ,ep;
1429 evalue *ne;
1430 value_init(ex); value_init(ey);value_init(ep);
1431 x=e1->x.p->size;
1432 y= res->x.p->size;
1433 value_set_si(ex,e1->x.p->size);
1434 value_set_si(ey,res->x.p->size);
1435 value_assign (ep,*Lcm(ex,ey));
1436 p=(int)mpz_get_si(ep);
1437 ne= (evalue *) malloc (sizeof(evalue));
1438 value_init(ne->d);
1439 value_set_si( ne->d,0);
1440
1441 ne->x.p=new_enode(res->x.p->type,p, res->x.p->pos);
1442 for(i=0;i<p;i++) {
1443 evalue_copy(&ne->x.p->arr[i], &res->x.p->arr[i%y]);
1444 }
1445 for(i=0;i<p;i++) {
1446 eadd(&e1->x.p->arr[i%x], &ne->x.p->arr[i]);
1447 }
1448
1449 value_assign(res->d, ne->d);
1450 res->x.p=ne->x.p;
1451
1452 return ;
1453 }
1454 case evector:
1455 fprintf(stderr, "eadd: ?cannot add vectors of different length\n");
1456 return ;
1457 default:
1458 assert(0);
1459 }
1460 }
1461 return ;
1462 }/* eadd */
1463
1464 static void emul_rev (evalue *e1, evalue *res)
1465 {
1466 evalue ev;
1467 value_init(ev.d);
1468 evalue_copy(&ev, e1);
1469 emul(res, &ev);
1470 free_evalue_refs(res);
1471 *res = ev;
1472 }
1473
1474 static void emul_poly (evalue *e1, evalue *res)
1475 {
1476 int i, j, o = type_offset(res->x.p);
1477 evalue tmp;
1478 int size=(e1->x.p->size + res->x.p->size - o - 1);
1479 value_init(tmp.d);
1480 value_set_si(tmp.d,0);
1481 tmp.x.p=new_enode(res->x.p->type, size, res->x.p->pos);
1482 if (o)
1483 evalue_copy(&tmp.x.p->arr[0], &e1->x.p->arr[0]);
1484 for (i=o; i < e1->x.p->size; i++) {
1485 evalue_copy(&tmp.x.p->arr[i], &e1->x.p->arr[i]);
1486 emul(&res->x.p->arr[o], &tmp.x.p->arr[i]);
1487 }
1488 for (; i<size; i++)
1489 evalue_set_si(&tmp.x.p->arr[i], 0, 1);
1490 for (i=o+1; i<res->x.p->size; i++)
1491 for (j=o; j<e1->x.p->size; j++) {
1492 evalue ev;
1493 value_init(ev.d);
1494 evalue_copy(&ev, &e1->x.p->arr[j]);
1495 emul(&res->x.p->arr[i], &ev);
1496 eadd(&ev, &tmp.x.p->arr[i+j-o]);
1497 free_evalue_refs(&ev);
1498 }
1499 free_evalue_refs(res);
1500 *res = tmp;
1501 }
1502
1503 void emul_partitions (evalue *e1,evalue *res)
1504 {
1505 int n, i, j, k;
1506 Polyhedron *d;
1507 struct section *s;
1508 s = (struct section *)
1509 malloc((e1->x.p->size/2) * (res->x.p->size/2) *
1510 sizeof(struct section));
1511 assert(s);
1512 assert(e1->x.p->pos == res->x.p->pos);
1513 assert(e1->x.p->pos == EVALUE_DOMAIN(e1->x.p->arr[0])->Dimension);
1514 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1515
1516 n = 0;
1517 for (i = 0; i < res->x.p->size/2; ++i) {
1518 for (j = 0; j < e1->x.p->size/2; ++j) {
1519 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1520 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1521 if (emptyQ(d)) {
1522 Domain_Free(d);
1523 continue;
1524 }
1525
1526 /* This code is only needed because the partitions
1527 are not true partitions.
1528 */
1529 for (k = 0; k < n; ++k) {
1530 if (DomainIncludes(s[k].D, d))
1531 break;
1532 if (DomainIncludes(d, s[k].D)) {
1533 Domain_Free(s[k].D);
1534 free_evalue_refs(&s[k].E);
1535 if (n > k)
1536 s[k] = s[--n];
1537 --k;
1538 }
1539 }
1540 if (k < n) {
1541 Domain_Free(d);
1542 continue;
1543 }
1544
1545 value_init(s[n].E.d);
1546 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
1547 emul(&e1->x.p->arr[2*j+1], &s[n].E);
1548 s[n].D = d;
1549 ++n;
1550 }
1551 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
1552 value_clear(res->x.p->arr[2*i].d);
1553 free_evalue_refs(&res->x.p->arr[2*i+1]);
1554 }
1555
1556 free(res->x.p);
1557 if (n == 0)
1558 evalue_set_si(res, 0, 1);
1559 else {
1560 res->x.p = new_enode(partition, 2*n, e1->x.p->pos);
1561 for (j = 0; j < n; ++j) {
1562 s[j].D = DomainConstraintSimplify(s[j].D, 0);
1563 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1564 value_clear(res->x.p->arr[2*j+1].d);
1565 res->x.p->arr[2*j+1] = s[j].E;
1566 }
1567 }
1568
1569 free(s);
1570 }
1571
1572 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1573
1574 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1575 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1576 * evalues is not treated here */
1577
1578 void emul (evalue *e1, evalue *res ){
1579 int i,j;
1580
1581 if((value_zero_p(e1->d)&&e1->x.p->type==evector)||(value_zero_p(res->d)&&(res->x.p->type==evector))) {
1582 fprintf(stderr, "emul: do not proced on evector type !\n");
1583 return;
1584 }
1585
1586 if (EVALUE_IS_ZERO(*res))
1587 return;
1588
1589 if (value_zero_p(e1->d) && e1->x.p->type == partition) {
1590 if (value_zero_p(res->d) && res->x.p->type == partition)
1591 emul_partitions(e1, res);
1592 else
1593 emul_rev(e1, res);
1594 } else if (value_zero_p(res->d) && res->x.p->type == partition) {
1595 for (i = 0; i < res->x.p->size/2; ++i)
1596 emul(e1, &res->x.p->arr[2*i+1]);
1597 } else
1598 if (value_zero_p(res->d) && res->x.p->type == relation) {
1599 if (value_zero_p(e1->d) && e1->x.p->type == relation &&
1600 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1601 if (res->x.p->size < 3 && e1->x.p->size == 3)
1602 explicit_complement(res);
1603 if (e1->x.p->size < 3 && res->x.p->size == 3)
1604 explicit_complement(e1);
1605 for (i = 1; i < res->x.p->size; ++i)
1606 emul(&e1->x.p->arr[i], &res->x.p->arr[i]);
1607 return;
1608 }
1609 for (i = 1; i < res->x.p->size; ++i)
1610 emul(e1, &res->x.p->arr[i]);
1611 } else
1612 if(value_zero_p(e1->d)&& value_zero_p(res->d)) {
1613 switch(e1->x.p->type) {
1614 case polynomial:
1615 switch(res->x.p->type) {
1616 case polynomial:
1617 if(e1->x.p->pos == res->x.p->pos) {
1618 /* Product of two polynomials of the same variable */
1619 emul_poly(e1, res);
1620 return;
1621 }
1622 else {
1623 /* Product of two polynomials of different variables */
1624
1625 if(res->x.p->pos < e1->x.p->pos)
1626 for( i=0; i<res->x.p->size ; i++)
1627 emul(e1, &res->x.p->arr[i]);
1628 else
1629 emul_rev(e1, res);
1630
1631 return ;
1632 }
1633 case periodic:
1634 case flooring:
1635 case fractional:
1636 /* Product of a polynomial and a periodic or fractional */
1637 emul_rev(e1, res);
1638 return;
1639 default:
1640 assert(0);
1641 }
1642 case periodic:
1643 switch(res->x.p->type) {
1644 case periodic:
1645 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size==res->x.p->size) {
1646 /* Product of two periodics of the same parameter and period */
1647
1648 for(i=0; i<res->x.p->size;i++)
1649 emul(&(e1->x.p->arr[i]), &(res->x.p->arr[i]));
1650
1651 return;
1652 }
1653 else{
1654 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size!=res->x.p->size) {
1655 /* Product of two periodics of the same parameter and different periods */
1656 evalue *newp;
1657 Value x,y,z;
1658 int ix,iy,lcm;
1659 value_init(x); value_init(y);value_init(z);
1660 ix=e1->x.p->size;
1661 iy=res->x.p->size;
1662 value_set_si(x,e1->x.p->size);
1663 value_set_si(y,res->x.p->size);
1664 value_assign (z,*Lcm(x,y));
1665 lcm=(int)mpz_get_si(z);
1666 newp= (evalue *) malloc (sizeof(evalue));
1667 value_init(newp->d);
1668 value_set_si( newp->d,0);
1669 newp->x.p=new_enode(periodic,lcm, e1->x.p->pos);
1670 for(i=0;i<lcm;i++) {
1671 evalue_copy(&newp->x.p->arr[i],
1672 &res->x.p->arr[i%iy]);
1673 }
1674 for(i=0;i<lcm;i++)
1675 emul(&e1->x.p->arr[i%ix], &newp->x.p->arr[i]);
1676
1677 value_assign(res->d,newp->d);
1678 res->x.p=newp->x.p;
1679
1680 value_clear(x); value_clear(y);value_clear(z);
1681 return ;
1682 }
1683 else {
1684 /* Product of two periodics of different parameters */
1685
1686 if(res->x.p->pos < e1->x.p->pos)
1687 for(i=0; i<res->x.p->size; i++)
1688 emul(e1, &(res->x.p->arr[i]));
1689 else
1690 emul_rev(e1, res);
1691
1692 return;
1693 }
1694 }
1695 case polynomial:
1696 /* Product of a periodic and a polynomial */
1697
1698 for(i=0; i<res->x.p->size ; i++)
1699 emul(e1, &(res->x.p->arr[i]));
1700
1701 return;
1702
1703 }
1704 case flooring:
1705 case fractional:
1706 switch(res->x.p->type) {
1707 case polynomial:
1708 for(i=0; i<res->x.p->size ; i++)
1709 emul(e1, &(res->x.p->arr[i]));
1710 return;
1711 default:
1712 case periodic:
1713 assert(0);
1714 case flooring:
1715 case fractional:
1716 assert(e1->x.p->type == res->x.p->type);
1717 if (e1->x.p->pos == res->x.p->pos &&
1718 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1719 evalue d;
1720 value_init(d.d);
1721 poly_denom(&e1->x.p->arr[0], &d.d);
1722 if (e1->x.p->type != fractional || !value_two_p(d.d))
1723 emul_poly(e1, res);
1724 else {
1725 evalue tmp;
1726 value_init(d.x.n);
1727 value_set_si(d.x.n, 1);
1728 /* { x }^2 == { x }/2 */
1729 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1730 assert(e1->x.p->size == 3);
1731 assert(res->x.p->size == 3);
1732 value_init(tmp.d);
1733 evalue_copy(&tmp, &res->x.p->arr[2]);
1734 emul(&d, &tmp);
1735 eadd(&res->x.p->arr[1], &tmp);
1736 emul(&e1->x.p->arr[2], &tmp);
1737 emul(&e1->x.p->arr[1], res);
1738 eadd(&tmp, &res->x.p->arr[2]);
1739 free_evalue_refs(&tmp);
1740 value_clear(d.x.n);
1741 }
1742 value_clear(d.d);
1743 } else {
1744 if(mod_term_smaller(res, e1))
1745 for(i=1; i<res->x.p->size ; i++)
1746 emul(e1, &(res->x.p->arr[i]));
1747 else
1748 emul_rev(e1, res);
1749 return;
1750 }
1751 }
1752 break;
1753 case relation:
1754 emul_rev(e1, res);
1755 break;
1756 default:
1757 assert(0);
1758 }
1759 }
1760 else {
1761 if (value_notzero_p(e1->d)&& value_notzero_p(res->d)) {
1762 /* Product of two rational numbers */
1763
1764 Value g;
1765 value_init(g);
1766 value_multiply(res->d,e1->d,res->d);
1767 value_multiply(res->x.n,e1->x.n,res->x.n );
1768 Gcd(res->x.n, res->d,&g);
1769 if (value_notone_p(g)) {
1770 value_division(res->d,res->d,g);
1771 value_division(res->x.n,res->x.n,g);
1772 }
1773 value_clear(g);
1774 return ;
1775 }
1776 else {
1777 if(value_zero_p(e1->d)&& value_notzero_p(res->d)) {
1778 /* Product of an expression (polynomial or peririodic) and a rational number */
1779
1780 emul_rev(e1, res);
1781 return ;
1782 }
1783 else {
1784 /* Product of a rationel number and an expression (polynomial or peririodic) */
1785
1786 i = type_offset(res->x.p);
1787 for (; i<res->x.p->size; i++)
1788 emul(e1, &res->x.p->arr[i]);
1789
1790 return ;
1791 }
1792 }
1793 }
1794
1795 return ;
1796 }
1797
1798 /* Frees mask content ! */
1799 void emask(evalue *mask, evalue *res) {
1800 int n, i, j;
1801 Polyhedron *d, *fd;
1802 struct section *s;
1803 evalue mone;
1804 int pos;
1805
1806 if (EVALUE_IS_ZERO(*res)) {
1807 free_evalue_refs(mask);
1808 return;
1809 }
1810
1811 assert(value_zero_p(mask->d));
1812 assert(mask->x.p->type == partition);
1813 assert(value_zero_p(res->d));
1814 assert(res->x.p->type == partition);
1815 assert(mask->x.p->pos == res->x.p->pos);
1816 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1817 assert(mask->x.p->pos == EVALUE_DOMAIN(mask->x.p->arr[0])->Dimension);
1818 pos = res->x.p->pos;
1819
1820 s = (struct section *)
1821 malloc((mask->x.p->size/2+1) * (res->x.p->size/2) *
1822 sizeof(struct section));
1823 assert(s);
1824
1825 value_init(mone.d);
1826 evalue_set_si(&mone, -1, 1);
1827
1828 n = 0;
1829 for (j = 0; j < res->x.p->size/2; ++j) {
1830 assert(mask->x.p->size >= 2);
1831 fd = DomainDifference(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1832 EVALUE_DOMAIN(mask->x.p->arr[0]), 0);
1833 if (!emptyQ(fd))
1834 for (i = 1; i < mask->x.p->size/2; ++i) {
1835 Polyhedron *t = fd;
1836 fd = DomainDifference(fd, EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1837 Domain_Free(t);
1838 if (emptyQ(fd))
1839 break;
1840 }
1841 if (emptyQ(fd)) {
1842 Domain_Free(fd);
1843 continue;
1844 }
1845 value_init(s[n].E.d);
1846 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1847 s[n].D = fd;
1848 ++n;
1849 }
1850 for (i = 0; i < mask->x.p->size/2; ++i) {
1851 if (EVALUE_IS_ONE(mask->x.p->arr[2*i+1]))
1852 continue;
1853
1854 fd = EVALUE_DOMAIN(mask->x.p->arr[2*i]);
1855 eadd(&mone, &mask->x.p->arr[2*i+1]);
1856 emul(&mone, &mask->x.p->arr[2*i+1]);
1857 for (j = 0; j < res->x.p->size/2; ++j) {
1858 Polyhedron *t;
1859 d = DomainIntersection(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1860 EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1861 if (emptyQ(d)) {
1862 Domain_Free(d);
1863 continue;
1864 }
1865 t = fd;
1866 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*j]), 0);
1867 if (t != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1868 Domain_Free(t);
1869 value_init(s[n].E.d);
1870 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1871 emul(&mask->x.p->arr[2*i+1], &s[n].E);
1872 s[n].D = d;
1873 ++n;
1874 }
1875
1876 if (!emptyQ(fd)) {
1877 /* Just ignore; this may have been previously masked off */
1878 }
1879 if (fd != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1880 Domain_Free(fd);
1881 }
1882
1883 free_evalue_refs(&mone);
1884 free_evalue_refs(mask);
1885 free_evalue_refs(res);
1886 value_init(res->d);
1887 if (n == 0)
1888 evalue_set_si(res, 0, 1);
1889 else {
1890 res->x.p = new_enode(partition, 2*n, pos);
1891 for (j = 0; j < n; ++j) {
1892 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1893 value_clear(res->x.p->arr[2*j+1].d);
1894 res->x.p->arr[2*j+1] = s[j].E;
1895 }
1896 }
1897
1898 free(s);
1899 }
1900
1901 void evalue_copy(evalue *dst, const evalue *src)
1902 {
1903 value_assign(dst->d, src->d);
1904 if(value_notzero_p(src->d)) {
1905 value_init(dst->x.n);
1906 value_assign(dst->x.n, src->x.n);
1907 } else
1908 dst->x.p = ecopy(src->x.p);
1909 }
1910
1911 enode *new_enode(enode_type type,int size,int pos) {
1912
1913 enode *res;
1914 int i;
1915
1916 if(size == 0) {
1917 fprintf(stderr, "Allocating enode of size 0 !\n" );
1918 return NULL;
1919 }
1920 res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
1921 res->type = type;
1922 res->size = size;
1923 res->pos = pos;
1924 for(i=0; i<size; i++) {
1925 value_init(res->arr[i].d);
1926 value_set_si(res->arr[i].d,0);
1927 res->arr[i].x.p = 0;
1928 }
1929 return res;
1930 } /* new_enode */
1931
1932 enode *ecopy(enode *e) {
1933
1934 enode *res;
1935 int i;
1936
1937 res = new_enode(e->type,e->size,e->pos);
1938 for(i=0;i<e->size;++i) {
1939 value_assign(res->arr[i].d,e->arr[i].d);
1940 if(value_zero_p(res->arr[i].d))
1941 res->arr[i].x.p = ecopy(e->arr[i].x.p);
1942 else if (EVALUE_IS_DOMAIN(res->arr[i]))
1943 EVALUE_SET_DOMAIN(res->arr[i], Domain_Copy(EVALUE_DOMAIN(e->arr[i])));
1944 else {
1945 value_init(res->arr[i].x.n);
1946 value_assign(res->arr[i].x.n,e->arr[i].x.n);
1947 }
1948 }
1949 return(res);
1950 } /* ecopy */
1951
1952 int ecmp(const evalue *e1, const evalue *e2)
1953 {
1954 enode *p1, *p2;
1955 int i;
1956 int r;
1957
1958 if (value_notzero_p(e1->d) && value_notzero_p(e2->d)) {
1959 Value m, m2;
1960 value_init(m);
1961 value_init(m2);
1962 value_multiply(m, e1->x.n, e2->d);
1963 value_multiply(m2, e2->x.n, e1->d);
1964
1965 if (value_lt(m, m2))
1966 r = -1;
1967 else if (value_gt(m, m2))
1968 r = 1;
1969 else
1970 r = 0;
1971
1972 value_clear(m);
1973 value_clear(m2);
1974
1975 return r;
1976 }
1977 if (value_notzero_p(e1->d))
1978 return -1;
1979 if (value_notzero_p(e2->d))
1980 return 1;
1981
1982 p1 = e1->x.p;
1983 p2 = e2->x.p;
1984
1985 if (p1->type != p2->type)
1986 return p1->type - p2->type;
1987 if (p1->pos != p2->pos)
1988 return p1->pos - p2->pos;
1989 if (p1->size != p2->size)
1990 return p1->size - p2->size;
1991
1992 for (i = p1->size-1; i >= 0; --i)
1993 if ((r = ecmp(&p1->arr[i], &p2->arr[i])) != 0)
1994 return r;
1995
1996 return 0;
1997 }
1998
1999 int eequal(const evalue *e1, const evalue *e2)
2000 {
2001 int i;
2002 enode *p1, *p2;
2003
2004 if (value_ne(e1->d,e2->d))
2005 return 0;
2006
2007 /* e1->d == e2->d */
2008 if (value_notzero_p(e1->d)) {
2009 if (value_ne(e1->x.n,e2->x.n))
2010 return 0;
2011
2012 /* e1->d == e2->d != 0 AND e1->n == e2->n */
2013 return 1;
2014 }
2015
2016 /* e1->d == e2->d == 0 */
2017 p1 = e1->x.p;
2018 p2 = e2->x.p;
2019 if (p1->type != p2->type) return 0;
2020 if (p1->size != p2->size) return 0;
2021 if (p1->pos != p2->pos) return 0;
2022 for (i=0; i<p1->size; i++)
2023 if (!eequal(&p1->arr[i], &p2->arr[i]) )
2024 return 0;
2025 return 1;
2026 } /* eequal */
2027
2028 void free_evalue_refs(evalue *e) {
2029
2030 enode *p;
2031 int i;
2032
2033 if (EVALUE_IS_DOMAIN(*e)) {
2034 Domain_Free(EVALUE_DOMAIN(*e));
2035 value_clear(e->d);
2036 return;
2037 } else if (value_pos_p(e->d)) {
2038
2039 /* 'e' stores a constant */
2040 value_clear(e->d);
2041 value_clear(e->x.n);
2042 return;
2043 }
2044 assert(value_zero_p(e->d));
2045 value_clear(e->d);
2046 p = e->x.p;
2047 if (!p) return; /* null pointer */
2048 for (i=0; i<p->size; i++) {
2049 free_evalue_refs(&(p->arr[i]));
2050 }
2051 free(p);
2052 return;
2053 } /* free_evalue_refs */
2054
2055 static void mod2table_r(evalue *e, Vector *periods, Value m, int p,
2056 Vector * val, evalue *res)
2057 {
2058 unsigned nparam = periods->Size;
2059
2060 if (p == nparam) {
2061 double d = compute_evalue(e, val->p);
2062 d *= VALUE_TO_DOUBLE(m);
2063 if (d > 0)
2064 d += .25;
2065 else
2066 d -= .25;
2067 value_assign(res->d, m);
2068 value_init(res->x.n);
2069 value_set_double(res->x.n, d);
2070 mpz_fdiv_r(res->x.n, res->x.n, m);
2071 return;
2072 }
2073 if (value_one_p(periods->p[p]))
2074 mod2table_r(e, periods, m, p+1, val, res);
2075 else {
2076 Value tmp;
2077 value_init(tmp);
2078
2079 value_assign(tmp, periods->p[p]);
2080 value_set_si(res->d, 0);
2081 res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
2082 do {
2083 value_decrement(tmp, tmp);
2084 value_assign(val->p[p], tmp);
2085 mod2table_r(e, periods, m, p+1, val,
2086 &res->x.p->arr[VALUE_TO_INT(tmp)]);
2087 } while (value_pos_p(tmp));
2088
2089 value_clear(tmp);
2090 }
2091 }
2092
2093 static void rel2table(evalue *e, int zero)
2094 {
2095 if (value_pos_p(e->d)) {
2096 if (value_zero_p(e->x.n) == zero)
2097 value_set_si(e->x.n, 1);
2098 else
2099 value_set_si(e->x.n, 0);
2100 value_set_si(e->d, 1);
2101 } else {
2102 int i;
2103 for (i = 0; i < e->x.p->size; ++i)
2104 rel2table(&e->x.p->arr[i], zero);
2105 }
2106 }
2107
2108 void evalue_mod2table(evalue *e, int nparam)
2109 {
2110 enode *p;
2111 int i;
2112
2113 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2114 return;
2115 p = e->x.p;
2116 for (i=0; i<p->size; i++) {
2117 evalue_mod2table(&(p->arr[i]), nparam);
2118 }
2119 if (p->type == relation) {
2120 evalue copy;
2121
2122 if (p->size > 2) {
2123 value_init(copy.d);
2124 evalue_copy(&copy, &p->arr[0]);
2125 }
2126 rel2table(&p->arr[0], 1);
2127 emul(&p->arr[0], &p->arr[1]);
2128 if (p->size > 2) {
2129 rel2table(&copy, 0);
2130 emul(&copy, &p->arr[2]);
2131 eadd(&p->arr[2], &p->arr[1]);
2132 free_evalue_refs(&p->arr[2]);
2133 free_evalue_refs(&copy);
2134 }
2135 free_evalue_refs(&p->arr[0]);
2136 value_clear(e->d);
2137 *e = p->arr[1];
2138 free(p);
2139 } else if (p->type == fractional) {
2140 Vector *periods = Vector_Alloc(nparam);
2141 Vector *val = Vector_Alloc(nparam);
2142 Value tmp;
2143 evalue *ev;
2144 evalue EP, res;
2145
2146 value_init(tmp);
2147 value_set_si(tmp, 1);
2148 Vector_Set(periods->p, 1, nparam);
2149 Vector_Set(val->p, 0, nparam);
2150 for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
2151 enode *p = ev->x.p;
2152
2153 assert(p->type == polynomial);
2154 assert(p->size == 2);
2155 value_assign(periods->p[p->pos-1], p->arr[1].d);
2156 value_lcm(tmp, p->arr[1].d, &tmp);
2157 }
2158 value_lcm(tmp, ev->d, &tmp);
2159 value_init(EP.d);
2160 mod2table_r(&p->arr[0], periods, tmp, 0, val, &EP);
2161
2162 value_init(res.d);
2163 evalue_set_si(&res, 0, 1);
2164 /* Compute the polynomial using Horner's rule */
2165 for (i=p->size-1;i>1;i--) {
2166 eadd(&p->arr[i], &res);
2167 emul(&EP, &res);
2168 }
2169 eadd(&p->arr[1], &res);
2170
2171 free_evalue_refs(e);
2172 free_evalue_refs(&EP);
2173 *e = res;
2174
2175 value_clear(tmp);
2176 Vector_Free(val);
2177 Vector_Free(periods);
2178 }
2179 } /* evalue_mod2table */
2180
2181 /********************************************************/
2182 /* function in domain */
2183 /* check if the parameters in list_args */
2184 /* verifies the constraints of Domain P */
2185 /********************************************************/
2186 int in_domain(Polyhedron *P, Value *list_args)
2187 {
2188 int row, in = 1;
2189 Value v; /* value of the constraint of a row when
2190 parameters are instantiated*/
2191
2192 value_init(v);
2193
2194 for (row = 0; row < P->NbConstraints; row++) {
2195 Inner_Product(P->Constraint[row]+1, list_args, P->Dimension, &v);
2196 value_addto(v, v, P->Constraint[row][P->Dimension+1]); /*constant part*/
2197 if (value_neg_p(v) ||
2198 value_zero_p(P->Constraint[row][0]) && value_notzero_p(v)) {
2199 in = 0;
2200 break;
2201 }
2202 }
2203
2204 value_clear(v);
2205 return in || (P->next && in_domain(P->next, list_args));
2206 } /* in_domain */
2207
2208 /****************************************************/
2209 /* function compute enode */
2210 /* compute the value of enode p with parameters */
2211 /* list "list_args */
2212 /* compute the polynomial or the periodic */
2213 /****************************************************/
2214
2215 static double compute_enode(enode *p, Value *list_args) {
2216
2217 int i;
2218 Value m, param;
2219 double res=0.0;
2220
2221 if (!p)
2222 return(0.);
2223
2224 value_init(m);
2225 value_init(param);
2226
2227 if (p->type == polynomial) {
2228 if (p->size > 1)
2229 value_assign(param,list_args[p->pos-1]);
2230
2231 /* Compute the polynomial using Horner's rule */
2232 for (i=p->size-1;i>0;i--) {
2233 res +=compute_evalue(&p->arr[i],list_args);
2234 res *=VALUE_TO_DOUBLE(param);
2235 }
2236 res +=compute_evalue(&p->arr[0],list_args);
2237 }
2238 else if (p->type == fractional) {
2239 double d = compute_evalue(&p->arr[0], list_args);
2240 d -= floor(d+1e-10);
2241
2242 /* Compute the polynomial using Horner's rule */
2243 for (i=p->size-1;i>1;i--) {
2244 res +=compute_evalue(&p->arr[i],list_args);
2245 res *=d;
2246 }
2247 res +=compute_evalue(&p->arr[1],list_args);
2248 }
2249 else if (p->type == flooring) {
2250 double d = compute_evalue(&p->arr[0], list_args);
2251 d = floor(d+1e-10);
2252
2253 /* Compute the polynomial using Horner's rule */
2254 for (i=p->size-1;i>1;i--) {
2255 res +=compute_evalue(&p->arr[i],list_args);
2256 res *=d;
2257 }
2258 res +=compute_evalue(&p->arr[1],list_args);
2259 }
2260 else if (p->type == periodic) {
2261 value_assign(m,list_args[p->pos-1]);
2262
2263 /* Choose the right element of the periodic */
2264 value_set_si(param,p->size);
2265 value_pmodulus(m,m,param);
2266 res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
2267 }
2268 else if (p->type == relation) {
2269 if (fabs(compute_evalue(&p->arr[0], list_args)) < 1e-10)
2270 res = compute_evalue(&p->arr[1], list_args);
2271 else if (p->size > 2)
2272 res = compute_evalue(&p->arr[2], list_args);
2273 }
2274 else if (p->type == partition) {
2275 int dim = EVALUE_DOMAIN(p->arr[0])->Dimension;
2276 Value *vals = list_args;
2277 if (p->pos < dim) {
2278 NALLOC(vals, dim);
2279 for (i = 0; i < dim; ++i) {
2280 value_init(vals[i]);
2281 if (i < p->pos)
2282 value_assign(vals[i], list_args[i]);
2283 }
2284 }
2285 for (i = 0; i < p->size/2; ++i)
2286 if (DomainContains(EVALUE_DOMAIN(p->arr[2*i]), vals, p->pos, 0, 1)) {
2287 res = compute_evalue(&p->arr[2*i+1], vals);
2288 break;
2289 }
2290 if (p->pos < dim) {
2291 for (i = 0; i < dim; ++i)
2292 value_clear(vals[i]);
2293 free(vals);
2294 }
2295 }
2296 else
2297 assert(0);
2298 value_clear(m);
2299 value_clear(param);
2300 return res;
2301 } /* compute_enode */
2302
2303 /*************************************************/
2304 /* return the value of Ehrhart Polynomial */
2305 /* It returns a double, because since it is */
2306 /* a recursive function, some intermediate value */
2307 /* might not be integral */
2308 /*************************************************/
2309
2310 double compute_evalue(const evalue *e, Value *list_args)
2311 {
2312 double res;
2313
2314 if (value_notzero_p(e->d)) {
2315 if (value_notone_p(e->d))
2316 res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
2317 else
2318 res = VALUE_TO_DOUBLE(e->x.n);
2319 }
2320 else
2321 res = compute_enode(e->x.p,list_args);
2322 return res;
2323 } /* compute_evalue */
2324
2325
2326 /****************************************************/
2327 /* function compute_poly : */
2328 /* Check for the good validity domain */
2329 /* return the number of point in the Polyhedron */
2330 /* in allocated memory */
2331 /* Using the Ehrhart pseudo-polynomial */
2332 /****************************************************/
2333 Value *compute_poly(Enumeration *en,Value *list_args) {
2334
2335 Value *tmp;
2336 /* double d; int i; */
2337
2338 tmp = (Value *) malloc (sizeof(Value));
2339 assert(tmp != NULL);
2340 value_init(*tmp);
2341 value_set_si(*tmp,0);
2342
2343 if(!en)
2344 return(tmp); /* no ehrhart polynomial */
2345 if(en->ValidityDomain) {
2346 if(!en->ValidityDomain->Dimension) { /* no parameters */
2347 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2348 return(tmp);
2349 }
2350 }
2351 else
2352 return(tmp); /* no Validity Domain */
2353 while(en) {
2354 if(in_domain(en->ValidityDomain,list_args)) {
2355
2356 #ifdef EVAL_EHRHART_DEBUG
2357 Print_Domain(stdout,en->ValidityDomain);
2358 print_evalue(stdout,&en->EP);
2359 #endif
2360
2361 /* d = compute_evalue(&en->EP,list_args);
2362 i = d;
2363 printf("(double)%lf = %d\n", d, i ); */
2364 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2365 return(tmp);
2366 }
2367 else
2368 en=en->next;
2369 }
2370 value_set_si(*tmp,0);
2371 return(tmp); /* no compatible domain with the arguments */
2372 } /* compute_poly */
2373
2374 static evalue *eval_polynomial(const enode *p, int offset,
2375 evalue *base, Value *values)
2376 {
2377 int i;
2378 evalue *res, *c;
2379
2380 res = evalue_zero();
2381 for (i = p->size-1; i > offset; --i) {
2382 c = evalue_eval(&p->arr[i], values);
2383 eadd(c, res);
2384 free_evalue_refs(c);
2385 free(c);
2386 emul(base, res);
2387 }
2388 c = evalue_eval(&p->arr[offset], values);
2389 eadd(c, res);
2390 free_evalue_refs(c);
2391 free(c);
2392
2393 return res;
2394 }
2395
2396 evalue *evalue_eval(const evalue *e, Value *values)
2397 {
2398 evalue *res = NULL;
2399 evalue param;
2400 evalue *param2;
2401 int i;
2402
2403 if (value_notzero_p(e->d)) {
2404 res = ALLOC(evalue);
2405 value_init(res->d);
2406 evalue_copy(res, e);
2407 return res;
2408 }
2409 switch (e->x.p->type) {
2410 case polynomial:
2411 value_init(param.x.n);
2412 value_assign(param.x.n, values[e->x.p->pos-1]);
2413 value_init(param.d);
2414 value_set_si(param.d, 1);
2415
2416 res = eval_polynomial(e->x.p, 0, &param, values);
2417 free_evalue_refs(&param);
2418 break;
2419 case fractional:
2420 param2 = evalue_eval(&e->x.p->arr[0], values);
2421 mpz_fdiv_r(param2->x.n, param2->x.n, param2->d);
2422
2423 res = eval_polynomial(e->x.p, 1, param2, values);
2424 free_evalue_refs(param2);
2425 free(param2);
2426 break;
2427 case flooring:
2428 param2 = evalue_eval(&e->x.p->arr[0], values);
2429 mpz_fdiv_q(param2->x.n, param2->x.n, param2->d);
2430 value_set_si(param2->d, 1);
2431
2432 res = eval_polynomial(e->x.p, 1, param2, values);
2433 free_evalue_refs(param2);
2434 free(param2);
2435 break;
2436 case relation:
2437 param2 = evalue_eval(&e->x.p->arr[0], values);
2438 if (value_zero_p(param2->x.n))
2439 res = evalue_eval(&e->x.p->arr[1], values);
2440 else if (e->x.p->size > 2)
2441 res = evalue_eval(&e->x.p->arr[2], values);
2442 else
2443 res = evalue_zero();
2444 free_evalue_refs(param2);
2445 free(param2);
2446 break;
2447 case partition:
2448 assert(e->x.p->pos == EVALUE_DOMAIN(e->x.p->arr[0])->Dimension);
2449 for (i = 0; i < e->x.p->size/2; ++i)
2450 if (in_domain(EVALUE_DOMAIN(e->x.p->arr[2*i]), values)) {
2451 res = evalue_eval(&e->x.p->arr[2*i+1], values);
2452 break;
2453 }
2454 if (!res)
2455 res = evalue_zero();
2456 break;
2457 default:
2458 assert(0);
2459 }
2460 return res;
2461 }
2462
2463 size_t value_size(Value v) {
2464 return (v[0]._mp_size > 0 ? v[0]._mp_size : -v[0]._mp_size)
2465 * sizeof(v[0]._mp_d[0]);
2466 }
2467
2468 size_t domain_size(Polyhedron *D)
2469 {
2470 int i, j;
2471 size_t s = sizeof(*D);
2472
2473 for (i = 0; i < D->NbConstraints; ++i)
2474 for (j = 0; j < D->Dimension+2; ++j)
2475 s += value_size(D->Constraint[i][j]);
2476
2477 /*
2478 for (i = 0; i < D->NbRays; ++i)
2479 for (j = 0; j < D->Dimension+2; ++j)
2480 s += value_size(D->Ray[i][j]);
2481 */
2482
2483 return D->next ? s+domain_size(D->next) : s;
2484 }
2485
2486 size_t enode_size(enode *p) {
2487 size_t s = sizeof(*p) - sizeof(p->arr[0]);
2488 int i;
2489
2490 if (p->type == partition)
2491 for (i = 0; i < p->size/2; ++i) {
2492 s += domain_size(EVALUE_DOMAIN(p->arr[2*i]));
2493 s += evalue_size(&p->arr[2*i+1]);
2494 }
2495 else
2496 for (i = 0; i < p->size; ++i) {
2497 s += evalue_size(&p->arr[i]);
2498 }
2499 return s;
2500 }
2501
2502 size_t evalue_size(evalue *e)
2503 {
2504 size_t s = sizeof(*e);
2505 s += value_size(e->d);
2506 if (value_notzero_p(e->d))
2507 s += value_size(e->x.n);
2508 else
2509 s += enode_size(e->x.p);
2510 return s;
2511 }
2512
2513 static evalue *find_second(evalue *base, evalue *cst, evalue *e, Value m)
2514 {
2515 evalue *found = NULL;
2516 evalue offset;
2517 evalue copy;
2518 int i;
2519
2520 if (value_pos_p(e->d) || e->x.p->type != fractional)
2521 return NULL;
2522
2523 value_init(offset.d);
2524 value_init(offset.x.n);
2525 poly_denom(&e->x.p->arr[0], &offset.d);
2526 value_lcm(m, offset.d, &offset.d);
2527 value_set_si(offset.x.n, 1);
2528
2529 value_init(copy.d);
2530 evalue_copy(&copy, cst);
2531
2532 eadd(&offset, cst);
2533 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2534
2535 if (eequal(base, &e->x.p->arr[0]))
2536 found = &e->x.p->arr[0];
2537 else {
2538 value_set_si(offset.x.n, -2);
2539
2540 eadd(&offset, cst);
2541 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2542
2543 if (eequal(base, &e->x.p->arr[0]))
2544 found = base;
2545 }
2546 free_evalue_refs(cst);
2547 free_evalue_refs(&offset);
2548 *cst = copy;
2549
2550 for (i = 1; !found && i < e->x.p->size; ++i)
2551 found = find_second(base, cst, &e->x.p->arr[i], m);
2552
2553 return found;
2554 }
2555
2556 static evalue *find_relation_pair(evalue *e)
2557 {
2558 int i;
2559 evalue *found = NULL;
2560
2561 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2562 return NULL;
2563
2564 if (e->x.p->type == fractional) {
2565 Value m;
2566 evalue *cst;
2567
2568 value_init(m);
2569 poly_denom(&e->x.p->arr[0], &m);
2570
2571 for (cst = &e->x.p->arr[0]; value_zero_p(cst->d);
2572 cst = &cst->x.p->arr[0])
2573 ;
2574
2575 for (i = 1; !found && i < e->x.p->size; ++i)
2576 found = find_second(&e->x.p->arr[0], cst, &e->x.p->arr[i], m);
2577
2578 value_clear(m);
2579 }
2580
2581 i = e->x.p->type == relation;
2582 for (; !found && i < e->x.p->size; ++i)
2583 found = find_relation_pair(&e->x.p->arr[i]);
2584
2585 return found;
2586 }
2587
2588 void evalue_mod2relation(evalue *e) {
2589 evalue *d;
2590
2591 if (value_zero_p(e->d) && e->x.p->type == partition) {
2592 int i;
2593
2594 for (i = 0; i < e->x.p->size/2; ++i) {
2595 evalue_mod2relation(&e->x.p->arr[2*i+1]);
2596 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2597 value_clear(e->x.p->arr[2*i].d);
2598 free_evalue_refs(&e->x.p->arr[2*i+1]);
2599 e->x.p->size -= 2;
2600 if (2*i < e->x.p->size) {
2601 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2602 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2603 }
2604 --i;
2605 }
2606 }
2607 if (e->x.p->size == 0) {
2608 free(e->x.p);
2609 evalue_set_si(e, 0, 1);
2610 }
2611
2612 return;
2613 }
2614
2615 while ((d = find_relation_pair(e)) != NULL) {
2616 evalue split;
2617 evalue *ev;
2618
2619 value_init(split.d);
2620 value_set_si(split.d, 0);
2621 split.x.p = new_enode(relation, 3, 0);
2622 evalue_set_si(&split.x.p->arr[1], 1, 1);
2623 evalue_set_si(&split.x.p->arr[2], 1, 1);
2624
2625 ev = &split.x.p->arr[0];
2626 value_set_si(ev->d, 0);
2627 ev->x.p = new_enode(fractional, 3, -1);
2628 evalue_set_si(&ev->x.p->arr[1], 0, 1);
2629 evalue_set_si(&ev->x.p->arr[2], 1, 1);
2630 evalue_copy(&ev->x.p->arr[0], d);
2631
2632 emul(&split, e);
2633
2634 reduce_evalue(e);
2635
2636 free_evalue_refs(&split);
2637 }
2638 }
2639
2640 static int evalue_comp(const void * a, const void * b)
2641 {
2642 const evalue *e1 = *(const evalue **)a;
2643 const evalue *e2 = *(const evalue **)b;
2644 return ecmp(e1, e2);
2645 }
2646
2647 void evalue_combine(evalue *e)
2648 {
2649 evalue **evs;
2650 int i, k;
2651 enode *p;
2652 evalue tmp;
2653
2654 if (value_notzero_p(e->d) || e->x.p->type != partition)
2655 return;
2656
2657 NALLOC(evs, e->x.p->size/2);
2658 for (i = 0; i < e->x.p->size/2; ++i)
2659 evs[i] = &e->x.p->arr[2*i+1];
2660 qsort(evs, e->x.p->size/2, sizeof(evs[0]), evalue_comp);
2661 p = new_enode(partition, e->x.p->size, e->x.p->pos);
2662 for (i = 0, k = 0; i < e->x.p->size/2; ++i) {
2663 if (k == 0 || ecmp(&p->arr[2*k-1], evs[i]) != 0) {
2664 value_clear(p->arr[2*k].d);
2665 value_clear(p->arr[2*k+1].d);
2666 p->arr[2*k] = *(evs[i]-1);
2667 p->arr[2*k+1] = *(evs[i]);
2668 ++k;
2669 } else {
2670 Polyhedron *D = EVALUE_DOMAIN(*(evs[i]-1));
2671 Polyhedron *L = D;
2672
2673 value_clear((evs[i]-1)->d);
2674
2675 while (L->next)
2676 L = L->next;
2677 L->next = EVALUE_DOMAIN(p->arr[2*k-2]);
2678 EVALUE_SET_DOMAIN(p->arr[2*k-2], D);
2679 free_evalue_refs(evs[i]);
2680 }
2681 }
2682
2683 for (i = 2*k ; i < p->size; ++i)
2684 value_clear(p->arr[i].d);
2685
2686 free(evs);
2687 free(e->x.p);
2688 p->size = 2*k;
2689 e->x.p = p;
2690
2691 for (i = 0; i < e->x.p->size/2; ++i) {
2692 Polyhedron *H;
2693 if (value_notzero_p(e->x.p->arr[2*i+1].d))
2694 continue;
2695 H = DomainConvex(EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
2696 if (H == NULL)
2697 continue;
2698 for (k = 0; k < e->x.p->size/2; ++k) {
2699 Polyhedron *D, *N, **P;
2700 if (i == k)
2701 continue;
2702 P = &EVALUE_DOMAIN(e->x.p->arr[2*k]);
2703 D = *P;
2704 if (D == NULL)
2705 continue;
2706 for (; D; D = N) {
2707 *P = D;
2708 N = D->next;
2709 if (D->NbEq <= H->NbEq) {
2710 P = &D->next;
2711 continue;
2712 }
2713
2714 value_init(tmp.d);
2715 tmp.x.p = new_enode(partition, 2, e->x.p->pos);
2716 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Copy(D));
2717 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
2718 reduce_evalue(&tmp);
2719 if (value_notzero_p(tmp.d) ||
2720 ecmp(&tmp.x.p->arr[1], &e->x.p->arr[2*k+1]) != 0)
2721 P = &D->next;
2722 else {
2723 D->next = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2724 EVALUE_DOMAIN(e->x.p->arr[2*i]) = D;
2725 *P = NULL;
2726 }
2727 free_evalue_refs(&tmp);
2728 }
2729 }
2730 Polyhedron_Free(H);
2731 }
2732
2733 for (i = 0; i < e->x.p->size/2; ++i) {
2734 Polyhedron *H, *E;
2735 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2736 if (!D) {
2737 value_clear(e->x.p->arr[2*i].d);
2738 free_evalue_refs(&e->x.p->arr[2*i+1]);
2739 e->x.p->size -= 2;
2740 if (2*i < e->x.p->size) {
2741 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2742 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2743 }
2744 --i;
2745 continue;
2746 }
2747 if (!D->next)
2748 continue;
2749 H = DomainConvex(D, 0);
2750 E = DomainDifference(H, D, 0);
2751 Domain_Free(D);
2752 D = DomainDifference(H, E, 0);
2753 Domain_Free(H);
2754 Domain_Free(E);
2755 EVALUE_SET_DOMAIN(p->arr[2*i], D);
2756 }
2757 }
2758
2759 /* Use smallest representative for coefficients in affine form in
2760 * argument of fractional.
2761 * Since any change will make the argument non-standard,
2762 * the containing evalue will have to be reduced again afterward.
2763 */
2764 static void fractional_minimal_coefficients(enode *p)
2765 {
2766 evalue *pp;
2767 Value twice;
2768 value_init(twice);
2769
2770 assert(p->type == fractional);
2771 pp = &p->arr[0];
2772 while (value_zero_p(pp->d)) {
2773 assert(pp->x.p->type == polynomial);
2774 assert(pp->x.p->size == 2);
2775 assert(value_notzero_p(pp->x.p->arr[1].d));
2776 mpz_mul_ui(twice, pp->x.p->arr[1].x.n, 2);
2777 if (value_gt(twice, pp->x.p->arr[1].d))
2778 value_subtract(pp->x.p->arr[1].x.n,
2779 pp->x.p->arr[1].x.n, pp->x.p->arr[1].d);
2780 pp = &pp->x.p->arr[0];
2781 }
2782
2783 value_clear(twice);
2784 }
2785
2786 static Polyhedron *polynomial_projection(enode *p, Polyhedron *D, Value *d,
2787 Matrix **R)
2788 {
2789 Polyhedron *I, *H;
2790 evalue *pp;
2791 unsigned dim = D->Dimension;
2792 Matrix *T = Matrix_Alloc(2, dim+1);
2793 assert(T);
2794
2795 assert(p->type == fractional);
2796 pp = &p->arr[0];
2797 value_set_si(T->p[1][dim], 1);
2798 poly_denom(pp, d);
2799 while (value_zero_p(pp->d)) {
2800 assert(pp->x.p->type == polynomial);
2801 assert(pp->x.p->size == 2);
2802 assert(value_notzero_p(pp->x.p->arr[1].d));
2803 value_division(T->p[0][pp->x.p->pos-1], *d, pp->x.p->arr[1].d);
2804 value_multiply(T->p[0][pp->x.p->pos-1],
2805 T->p[0][pp->x.p->pos-1], pp->x.p->arr[1].x.n);
2806 pp = &pp->x.p->arr[0];
2807 }
2808 value_division(T->p[0][dim], *d, pp->d);
2809 value_multiply(T->p[0][dim], T->p[0][dim], pp->x.n);
2810 I = DomainImage(D, T, 0);
2811 H = DomainConvex(I, 0);
2812 Domain_Free(I);
2813 if (R)
2814 *R = T;
2815 else
2816 Matrix_Free(T);
2817
2818 return H;
2819 }
2820
2821 int evalue_range_reduction_in_domain(evalue *e, Polyhedron *D)
2822 {
2823 int i;
2824 enode *p;
2825 Value d, min, max;
2826 int r = 0;
2827 Polyhedron *I;
2828 int bounded;
2829
2830 if (value_notzero_p(e->d))
2831 return r;
2832
2833 p = e->x.p;
2834
2835 if (p->type == relation) {
2836 Matrix *T;
2837 int equal;
2838 value_init(d);
2839 value_init(min);
2840 value_init(max);
2841
2842 fractional_minimal_coefficients(p->arr[0].x.p);
2843 I = polynomial_projection(p->arr[0].x.p, D, &d, &T);
2844 bounded = line_minmax(I, &min, &max); /* frees I */
2845 equal = value_eq(min, max);
2846 mpz_cdiv_q(min, min, d);
2847 mpz_fdiv_q(max, max, d);
2848
2849 if (bounded && value_gt(min, max)) {
2850 /* Never zero */
2851 if (p->size == 3) {
2852 value_clear(e->d);
2853 *e = p->arr[2];
2854 } else {
2855 evalue_set_si(e, 0, 1);
2856 r = 1;
2857 }
2858 free_evalue_refs(&(p->arr[1]));
2859 free_evalue_refs(&(p->arr[0]));
2860 free(p);
2861 value_clear(d);
2862 value_clear(min);
2863 value_clear(max);
2864 Matrix_Free(T);
2865 return r ? r : evalue_range_reduction_in_domain(e, D);
2866 } else if (bounded && equal) {
2867 /* Always zero */
2868 if (p->size == 3)
2869 free_evalue_refs(&(p->arr[2]));
2870 value_clear(e->d);
2871 *e = p->arr[1];
2872 free_evalue_refs(&(p->arr[0]));
2873 free(p);
2874 value_clear(d);
2875 value_clear(min);
2876 value_clear(max);
2877 Matrix_Free(T);
2878 return evalue_range_reduction_in_domain(e, D);
2879 } else if (bounded && value_eq(min, max)) {
2880 /* zero for a single value */
2881 Polyhedron *E;
2882 Matrix *M = Matrix_Alloc(1, D->Dimension+2);
2883 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
2884 value_multiply(min, min, d);
2885 value_subtract(M->p[0][D->Dimension+1],
2886 M->p[0][D->Dimension+1], min);
2887 E = DomainAddConstraints(D, M, 0);
2888 value_clear(d);
2889 value_clear(min);
2890 value_clear(max);
2891 Matrix_Free(T);
2892 Matrix_Free(M);
2893 r = evalue_range_reduction_in_domain(&p->arr[1], E);
2894 if (p->size == 3)
2895 r |= evalue_range_reduction_in_domain(&p->arr[2], D);
2896 Domain_Free(E);
2897 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2898 return r;
2899 }
2900
2901 value_clear(d);
2902 value_clear(min);
2903 value_clear(max);
2904 Matrix_Free(T);
2905 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2906 }
2907
2908 i = p->type == relation ? 1 :
2909 p->type == fractional ? 1 : 0;
2910 for (; i<p->size; i++)
2911 r |= evalue_range_reduction_in_domain(&p->arr[i], D);
2912
2913 if (p->type != fractional) {
2914 if (r && p->type == polynomial) {
2915 evalue f;
2916 value_init(f.d);
2917 value_set_si(f.d, 0);
2918 f.x.p = new_enode(polynomial, 2, p->pos);
2919 evalue_set_si(&f.x.p->arr[0], 0, 1);
2920 evalue_set_si(&f.x.p->arr[1], 1, 1);
2921 reorder_terms_about(p, &f);
2922 value_clear(e->d);
2923 *e = p->arr[0];
2924 free(p);
2925 }
2926 return r;
2927 }
2928
2929 value_init(d);
2930 value_init(min);
2931 value_init(max);
2932 fractional_minimal_coefficients(p);
2933 I = polynomial_projection(p, D, &d, NULL);
2934 bounded = line_minmax(I, &min, &max); /* frees I */
2935 mpz_fdiv_q(min, min, d);
2936 mpz_fdiv_q(max, max, d);
2937 value_subtract(d, max, min);
2938
2939 if (bounded && value_eq(min, max)) {
2940 evalue inc;
2941 value_init(inc.d);
2942 value_init(inc.x.n);
2943 value_set_si(inc.d, 1);
2944 value_oppose(inc.x.n, min);
2945 eadd(&inc, &p->arr[0]);
2946 reorder_terms_about(p, &p->arr[0]); /* frees arr[0] */
2947 value_clear(e->d);
2948 *e = p->arr[1];
2949 free(p);
2950 free_evalue_refs(&inc);
2951 r = 1;
2952 } else if (bounded && value_one_p(d) && p->size > 3) {
2953 /* replace {g}^2 by -(g-min)^2 + (2{g}+1)*(g-min) - {g}
2954 * See pages 199-200 of PhD thesis.
2955 */
2956 evalue rem;
2957 evalue inc;
2958 evalue t;
2959 evalue f;
2960 evalue factor;
2961 value_init(rem.d);
2962 value_set_si(rem.d, 0);
2963 rem.x.p = new_enode(fractional, 3, -1);
2964 evalue_copy(&rem.x.p->arr[0], &p->arr[0]);
2965 value_clear(rem.x.p->arr[1].d);
2966 value_clear(rem.x.p->arr[2].d);
2967 rem.x.p->arr[1] = p->arr[1];
2968 rem.x.p->arr[2] = p->arr[2];
2969 for (i = 3; i < p->size; ++i)
2970 p->arr[i-2] = p->arr[i];
2971 p->size -= 2;
2972
2973 value_init(inc.d);
2974 value_init(inc.x.n);
2975 value_set_si(inc.d, 1);
2976 value_oppose(inc.x.n, min);
2977
2978 value_init(t.d);
2979 evalue_copy(&t, &p->arr[0]);
2980 eadd(&inc, &t);
2981
2982 value_init(f.d);
2983 value_set_si(f.d, 0);
2984 f.x.p = new_enode(fractional, 3, -1);
2985 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
2986 evalue_set_si(&f.x.p->arr[1], 1, 1);
2987 evalue_set_si(&f.x.p->arr[2], 2, 1);
2988
2989 value_init(factor.d);
2990 evalue_set_si(&factor, -1, 1);
2991 emul(&t, &factor);
2992
2993 eadd(&f, &factor);
2994 emul(&t, &factor);
2995
2996 value_clear(f.x.p->arr[1].x.n);
2997 value_clear(f.x.p->arr[2].x.n);
2998 evalue_set_si(&f.x.p->arr[1], 0, 1);
2999 evalue_set_si(&f.x.p->arr[2], -1, 1);
3000 eadd(&f, &factor);
3001
3002 if (r) {
3003 reorder_terms(&rem);
3004 reorder_terms(e);
3005 }
3006
3007 emul(&factor, e);
3008 eadd(&rem, e);
3009
3010 free_evalue_refs(&inc);
3011 free_evalue_refs(&t);
3012 free_evalue_refs(&f);
3013 free_evalue_refs(&factor);
3014 free_evalue_refs(&rem);
3015
3016 evalue_range_reduction_in_domain(e, D);
3017
3018 r = 1;
3019 } else {
3020 _reduce_evalue(&p->arr[0], 0, 1);
3021 if (r)
3022 reorder_terms(e);
3023 }
3024
3025 value_clear(d);
3026 value_clear(min);
3027 value_clear(max);
3028
3029 return r;
3030 }
3031
3032 void evalue_range_reduction(evalue *e)
3033 {
3034 int i;
3035 if (value_notzero_p(e->d) || e->x.p->type != partition)
3036 return;
3037
3038 for (i = 0; i < e->x.p->size/2; ++i)
3039 if (evalue_range_reduction_in_domain(&e->x.p->arr[2*i+1],
3040 EVALUE_DOMAIN(e->x.p->arr[2*i]))) {
3041 reduce_evalue(&e->x.p->arr[2*i+1]);
3042
3043 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
3044 free_evalue_refs(&e->x.p->arr[2*i+1]);
3045 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
3046 value_clear(e->x.p->arr[2*i].d);
3047 e->x.p->size -= 2;
3048 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
3049 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
3050 --i;
3051 }
3052 }
3053 }
3054
3055 /* Frees EP
3056 */
3057 Enumeration* partition2enumeration(evalue *EP)
3058 {
3059 int i;
3060 Enumeration *en, *res = NULL;
3061
3062 if (EVALUE_IS_ZERO(*EP)) {
3063 free(EP);
3064 return res;
3065 }
3066
3067 for (i = 0; i < EP->x.p->size/2; ++i) {
3068 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[2*i])->Dimension);
3069 en = (Enumeration *)malloc(sizeof(Enumeration));
3070 en->next = res;
3071 res = en;
3072 res->ValidityDomain = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
3073 value_clear(EP->x.p->arr[2*i].d);
3074 res->EP = EP->x.p->arr[2*i+1];
3075 }
3076 free(EP->x.p);
3077 value_clear(EP->d);
3078 free(EP);
3079 return res;
3080 }
3081
3082 int evalue_frac2floor_in_domain3(evalue *e, Polyhedron *D, int shift)
3083 {
3084 enode *p;
3085 int r = 0;
3086 int i;
3087 Polyhedron *I;
3088 Value d, min;
3089 evalue fl;
3090
3091 if (value_notzero_p(e->d))
3092 return r;
3093
3094 p = e->x.p;
3095
3096 i = p->type == relation ? 1 :
3097 p->type == fractional ? 1 : 0;
3098 for (; i<p->size; i++)
3099 r |= evalue_frac2floor_in_domain3(&p->arr[i], D, shift);
3100
3101 if (p->type != fractional) {
3102 if (r && p->type == polynomial) {
3103 evalue f;
3104 value_init(f.d);
3105 value_set_si(f.d, 0);
3106 f.x.p = new_enode(polynomial, 2, p->pos);
3107 evalue_set_si(&f.x.p->arr[0], 0, 1);
3108 evalue_set_si(&f.x.p->arr[1], 1, 1);
3109 reorder_terms_about(p, &f);
3110 value_clear(e->d);
3111 *e = p->arr[0];
3112 free(p);
3113 }
3114 return r;
3115 }
3116
3117 if (shift) {
3118 value_init(d);
3119 I = polynomial_projection(p, D, &d, NULL);
3120
3121 /*
3122 Polyhedron_Print(stderr, P_VALUE_FMT, I);
3123 */
3124
3125 assert(I->NbEq == 0); /* Should have been reduced */
3126
3127 /* Find minimum */
3128 for (i = 0; i < I->NbConstraints; ++i)
3129 if (value_pos_p(I->Constraint[i][1]))
3130 break;
3131
3132 if (i < I->NbConstraints) {
3133 value_init(min);
3134 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
3135 mpz_cdiv_q(min, I->Constraint[i][2], I->Constraint[i][1]);
3136 if (value_neg_p(min)) {
3137 evalue offset;
3138 mpz_fdiv_q(min, min, d);
3139 value_init(offset.d);
3140 value_set_si(offset.d, 1);
3141 value_init(offset.x.n);
3142 value_oppose(offset.x.n, min);
3143 eadd(&offset, &p->arr[0]);
3144 free_evalue_refs(&offset);
3145 }
3146 value_clear(min);
3147 }
3148
3149 Polyhedron_Free(I);
3150 value_clear(d);
3151 }
3152
3153 value_init(fl.d);
3154 value_set_si(fl.d, 0);
3155 fl.x.p = new_enode(flooring, 3, -1);
3156 evalue_set_si(&fl.x.p->arr[1], 0, 1);
3157 evalue_set_si(&fl.x.p->arr[2], -1, 1);
3158 evalue_copy(&fl.x.p->arr[0], &p->arr[0]);
3159
3160 eadd(&fl, &p->arr[0]);
3161 reorder_terms_about(p, &p->arr[0]);
3162 value_clear(e->d);
3163 *e = p->arr[1];
3164 free(p);
3165 free_evalue_refs(&fl);
3166
3167 return 1;
3168 }
3169
3170 int evalue_frac2floor_in_domain(evalue *e, Polyhedron *D)
3171 {
3172 return evalue_frac2floor_in_domain3(e, D, 1);
3173 }
3174
3175 void evalue_frac2floor2(evalue *e, int shift)
3176 {
3177 int i;
3178 if (value_notzero_p(e->d) || e->x.p->type != partition) {
3179 if (!shift) {
3180 if (evalue_frac2floor_in_domain3(e, NULL, 0))
3181 reduce_evalue(e);
3182 }
3183 return;
3184 }
3185
3186 for (i = 0; i < e->x.p->size/2; ++i)
3187 if (evalue_frac2floor_in_domain3(&e->x.p->arr[2*i+1],
3188 EVALUE_DOMAIN(e->x.p->arr[2*i]), shift))
3189 reduce_evalue(&e->x.p->arr[2*i+1]);
3190 }
3191
3192 void evalue_frac2floor(evalue *e)
3193 {
3194 evalue_frac2floor2(e, 1);
3195 }
3196
3197 static Matrix *esum_add_constraint(int nvar, Polyhedron *D, Matrix *C,
3198 Vector *row)
3199 {
3200 int nr, nc;
3201 int i;
3202 int nparam = D->Dimension - nvar;
3203
3204 if (C == 0) {
3205 nr = D->NbConstraints + 2;
3206 nc = D->Dimension + 2 + 1;
3207 C = Matrix_Alloc(nr, nc);
3208 for (i = 0; i < D->NbConstraints; ++i) {
3209 Vector_Copy(D->Constraint[i], C->p[i], 1 + nvar);
3210 Vector_Copy(D->Constraint[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
3211 D->Dimension + 1 - nvar);
3212 }
3213 } else {
3214 Matrix *oldC = C;
3215 nr = C->NbRows + 2;
3216 nc = C->NbColumns + 1;
3217 C = Matrix_Alloc(nr, nc);
3218 for (i = 0; i < oldC->NbRows; ++i) {
3219 Vector_Copy(oldC->p[i], C->p[i], 1 + nvar);
3220 Vector_Copy(oldC->p[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
3221 oldC->NbColumns - 1 - nvar);
3222 }
3223 }
3224 value_set_si(C->p[nr-2][0], 1);
3225 value_set_si(C->p[nr-2][1 + nvar], 1);
3226 value_set_si(C->p[nr-2][nc - 1], -1);
3227
3228 Vector_Copy(row->p, C->p[nr-1], 1 + nvar + 1);
3229 Vector_Copy(row->p + 1 + nvar + 1, C->p[nr-1] + C->NbColumns - 1 - nparam,
3230 1 + nparam);
3231
3232 return C;
3233 }
3234
3235 static void floor2frac_r(evalue *e, int nvar)
3236 {
3237 enode *p;
3238 int i;
3239 evalue f;
3240 evalue *pp;
3241
3242 if (value_notzero_p(e->d))
3243 return;
3244
3245 p = e->x.p;
3246
3247 assert(p->type == flooring);
3248 for (i = 1; i < p->size; i++)
3249 floor2frac_r(&p->arr[i], nvar);
3250
3251 for (pp = &p->arr[0]; value_zero_p(pp->d); pp = &pp->x.p->arr[0]) {
3252 assert(pp->x.p->type == polynomial);
3253 pp->x.p->pos -= nvar;
3254 }
3255
3256 value_init(f.d);
3257 value_set_si(f.d, 0);
3258 f.x.p = new_enode(fractional, 3, -1);
3259 evalue_set_si(&f.x.p->arr[1], 0, 1);
3260 evalue_set_si(&f.x.p->arr[2], -1, 1);
3261 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
3262
3263 eadd(&f, &p->arr[0]);
3264 reorder_terms_about(p, &p->arr[0]);
3265 value_clear(e->d);
3266 *e = p->arr[1];
3267 free(p);
3268 free_evalue_refs(&f);
3269 }
3270
3271 /* Convert flooring back to fractional and shift position
3272 * of the parameters by nvar
3273 */
3274 static void floor2frac(evalue *e, int nvar)
3275 {
3276 floor2frac_r(e, nvar);
3277 reduce_evalue(e);
3278 }
3279
3280 evalue *esum_over_domain_cst(int nvar, Polyhedron *D, Matrix *C)
3281 {
3282 evalue *t;
3283 int nparam = D->Dimension - nvar;
3284
3285 if (C != 0) {
3286 C = Matrix_Copy(C);
3287 D = Constraints2Polyhedron(C, 0);
3288 Matrix_Free(C);
3289 }
3290
3291 t = barvinok_enumerate_e(D, 0, nparam, 0);
3292
3293 /* Double check that D was not unbounded. */
3294 assert(!(value_pos_p(t->d) && value_neg_p(t->x.n)));
3295
3296 if (C != 0)
3297 Polyhedron_Free(D);
3298
3299 return t;
3300 }
3301
3302 evalue *esum_over_domain(evalue *e, int nvar, Polyhedron *D,
3303 Matrix *C)
3304 {
3305 Vector *row = NULL;
3306 int i;
3307 evalue *res;
3308 Matrix *origC;
3309 evalue *factor = NULL;
3310 evalue cum;
3311
3312 if (EVALUE_IS_ZERO(*e))
3313 return 0;
3314
3315 if (D->next) {
3316 Polyhedron *DD = Disjoint_Domain(D, 0, 0);
3317 Polyhedron *Q;
3318
3319 Q = DD;
3320 DD = Q->next;
3321 Q->next = 0;
3322
3323 res = esum_over_domain(e, nvar, Q, C);
3324 Polyhedron_Free(Q);
3325
3326 for (Q = DD; Q; Q = DD) {
3327 evalue *t;
3328
3329 DD = Q->next;
3330 Q->next = 0;
3331
3332 t = esum_over_domain(e, nvar, Q, C);
3333 Polyhedron_Free(Q);
3334
3335 if (!res)
3336 res = t;
3337 else if (t) {
3338 eadd(t, res);
3339 free_evalue_refs(t);
3340 free(t);
3341 }
3342 }
3343 return res;
3344 }
3345
3346 if (value_notzero_p(e->d)) {
3347 evalue *t;
3348
3349 t = esum_over_domain_cst(nvar, D, C);
3350
3351 if (!EVALUE_IS_ONE(*e))
3352 emul(e, t);
3353
3354 return t;
3355 }
3356
3357 switch (e->x.p->type) {
3358 case flooring: {
3359 evalue *pp = &e->x.p->arr[0];
3360
3361 if (pp->x.p->pos > nvar) {
3362 /* remainder is independent of the summated vars */
3363 evalue f;
3364 evalue *t;
3365
3366 value_init(f.d);
3367 evalue_copy(&f, e);
3368 floor2frac(&f, nvar);
3369
3370 t = esum_over_domain_cst(nvar, D, C);
3371
3372 emul(&f, t);
3373
3374 free_evalue_refs(&f);
3375
3376 return t;
3377 }
3378
3379 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3380 poly_denom(pp, &row->p[1 + nvar]);
3381 value_set_si(row->p[0], 1);
3382 for (pp = &e->x.p->arr[0]; value_zero_p(pp->d);
3383 pp = &pp->x.p->arr[0]) {
3384 int pos;
3385 assert(pp->x.p->type == polynomial);
3386 pos = pp->x.p->pos;
3387 if (pos >= 1 + nvar)
3388 ++pos;
3389 value_assign(row->p[pos], row->p[1+nvar]);
3390 value_division(row->p[pos], row->p[pos], pp->x.p->arr[1].d);
3391 value_multiply(row->p[pos], row->p[pos], pp->x.p->arr[1].x.n);
3392 }
3393 value_assign(row->p[1 + D->Dimension + 1], row->p[1+nvar]);
3394 value_division(row->p[1 + D->Dimension + 1],
3395 row->p[1 + D->Dimension + 1],
3396 pp->d);
3397 value_multiply(row->p[1 + D->Dimension + 1],
3398 row->p[1 + D->Dimension + 1],
3399 pp->x.n);
3400 value_oppose(row->p[1 + nvar], row->p[1 + nvar]);
3401 break;
3402 }
3403 case polynomial: {
3404 int pos = e->x.p->pos;
3405
3406 if (pos > nvar) {
3407 factor = ALLOC(evalue);
3408 value_init(factor->d);
3409 value_set_si(factor->d, 0);
3410 factor->x.p = new_enode(polynomial, 2, pos - nvar);
3411 evalue_set_si(&factor->x.p->arr[0], 0, 1);
3412 evalue_set_si(&factor->x.p->arr[1], 1, 1);
3413 break;
3414 }
3415
3416 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3417 for (i = 0; i < D->NbRays; ++i)
3418 if (value_notzero_p(D->Ray[i][pos]))
3419 break;
3420 assert(i < D->NbRays);
3421 if (value_neg_p(D->Ray[i][pos])) {
3422 factor = ALLOC(evalue);
3423 value_init(factor->d);
3424 evalue_set_si(factor, -1, 1);
3425 }
3426 value_set_si(row->p[0], 1);
3427 value_set_si(row->p[pos], 1);
3428 value_set_si(row->p[1 + nvar], -1);
3429 break;
3430 }
3431 default:
3432 assert(0);
3433 }
3434
3435 i = type_offset(e->x.p);
3436
3437 res = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3438 ++i;
3439
3440 if (factor) {
3441 value_init(cum.d);
3442 evalue_copy(&cum, factor);
3443 }
3444
3445 origC = C;
3446 for (; i < e->x.p->size; ++i) {
3447 evalue *t;
3448 if (row) {
3449 Matrix *prevC = C;
3450 C = esum_add_constraint(nvar, D, C, row);
3451 if (prevC != origC)
3452 Matrix_Free(prevC);
3453 }
3454 /*
3455 if (row)
3456 Vector_Print(stderr, P_VALUE_FMT, row);
3457 if (C)
3458 Matrix_Print(stderr, P_VALUE_FMT, C);
3459 */
3460 t = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3461
3462 if (t && factor)
3463 emul(&cum, t);
3464
3465 if (!res)
3466 res = t;
3467 else if (t) {
3468 eadd(t, res);
3469 free_evalue_refs(t);
3470 free(t);
3471 }
3472 if (factor && i+1 < e->x.p->size)
3473 emul(factor, &cum);
3474 }
3475 if (C != origC)
3476 Matrix_Free(C);
3477
3478 if (factor) {
3479 free_evalue_refs(factor);
3480 free_evalue_refs(&cum);
3481 free(factor);
3482 }
3483
3484 if (row)
3485 Vector_Free(row);
3486
3487 reduce_evalue(res);
3488
3489 return res;
3490 }
3491
3492 evalue *esum(evalue *e, int nvar)
3493 {
3494 int i;
3495 evalue *res = ALLOC(evalue);
3496 value_init(res->d);
3497
3498 assert(nvar >= 0);
3499 if (nvar == 0 || EVALUE_IS_ZERO(*e)) {
3500 evalue_copy(res, e);
3501 return res;
3502 }
3503
3504 evalue_set_si(res, 0, 1);
3505
3506 assert(value_zero_p(e->d));
3507 assert(e->x.p->type == partition);
3508
3509 for (i = 0; i < e->x.p->size/2; ++i) {
3510 evalue *t;
3511 t = esum_over_domain(&e->x.p->arr[2*i+1], nvar,
3512 EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
3513 eadd(t, res);
3514 free_evalue_refs(t);
3515 free(t);
3516 }
3517
3518 reduce_evalue(res);
3519
3520 return res;
3521 }
3522
3523 /* Initial silly implementation */
3524 void eor(evalue *e1, evalue *res)
3525 {
3526 evalue E;
3527 evalue mone;
3528 value_init(E.d);
3529 value_init(mone.d);
3530 evalue_set_si(&mone, -1, 1);
3531
3532 evalue_copy(&E, res);
3533 eadd(e1, &E);
3534 emul(e1, res);
3535 emul(&mone, res);
3536 eadd(&E, res);
3537
3538 free_evalue_refs(&E);
3539 free_evalue_refs(&mone);
3540 }
3541
3542 /* computes denominator of polynomial evalue
3543 * d should point to a value initialized to 1
3544 */
3545 void evalue_denom(const evalue *e, Value *d)
3546 {
3547 int i, offset;
3548
3549 if (value_notzero_p(e->d)) {
3550 value_lcm(*d, e->d, d);
3551 return;
3552 }
3553 assert(e->x.p->type == polynomial ||
3554 e->x.p->type == fractional ||
3555 e->x.p->type == flooring);
3556 offset = type_offset(e->x.p);
3557 for (i = e->x.p->size-1; i >= offset; --i)
3558 evalue_denom(&e->x.p->arr[i], d);
3559 }
3560
3561 /* Divides the evalue e by the integer n */
3562 void evalue_div(evalue * e, Value n)
3563 {
3564 int i, offset;
3565
3566 if (value_notzero_p(e->d)) {
3567 Value gc;
3568 value_init(gc);
3569 value_multiply(e->d, e->d, n);
3570 Gcd(e->x.n, e->d, &gc);
3571 if (value_notone_p(gc)) {
3572 value_division(e->d, e->d, gc);
3573 value_division(e->x.n, e->x.n, gc);
3574 }
3575 value_clear(gc);
3576 return;
3577 }
3578 if (e->x.p->type == partition) {
3579 for (i = 0; i < e->x.p->size/2; ++i)
3580 evalue_div(&e->x.p->arr[2*i+1], n);
3581 return;
3582 }
3583 offset = type_offset(e->x.p);
3584 for (i = e->x.p->size-1; i >= offset; --i)
3585 evalue_div(&e->x.p->arr[i], n);
3586 }
3587
3588 static void evalue_frac2polynomial_r(evalue *e, int *signs, int sign, int in_frac)
3589 {
3590 int i, offset;
3591 Value d;
3592 enode *p;
3593 int sign_odd = sign;
3594
3595 if (value_notzero_p(e->d)) {
3596 if (in_frac && sign * value_sign(e->x.n) < 0) {
3597 value_set_si(e->x.n, 0);
3598 value_set_si(e->d, 1);
3599 }
3600 return;
3601 }
3602
3603 if (e->x.p->type == relation) {
3604 for (i = e->x.p->size-1; i >= 1; --i)
3605 evalue_frac2polynomial_r(&e->x.p->arr[i], signs, sign, in_frac);
3606 return;
3607 }
3608
3609 if (e->x.p->type == polynomial)
3610 sign_odd *= signs[e->x.p->pos-1];
3611 offset = type_offset(e->x.p);
3612 evalue_frac2polynomial_r(&e->x.p->arr[offset], signs, sign, in_frac);
3613 in_frac |= e->x.p->type == fractional;
3614 for (i = e->x.p->size-1; i > offset; --i)
3615 evalue_frac2polynomial_r(&e->x.p->arr[i], signs,
3616 (i - offset) % 2 ? sign_odd : sign, in_frac);
3617
3618 if (e->x.p->type != fractional)
3619 return;
3620
3621 /* replace { a/m } by (m-1)/m if sign != 0
3622 * and by (m-1)/(2m) if sign == 0
3623 */
3624 value_init(d);
3625 value_set_si(d, 1);
3626 evalue_denom(&e->x.p->arr[0], &d);
3627 free_evalue_refs(&e->x.p->arr[0]);
3628 value_init(e->x.p->arr[0].d);
3629 value_init(e->x.p->arr[0].x.n);
3630 if (sign == 0)
3631 value_addto(e->x.p->arr[0].d, d, d);
3632 else
3633 value_assign(e->x.p->arr[0].d, d);
3634 value_decrement(e->x.p->arr[0].x.n, d);
3635 value_clear(d);
3636
3637 p = e->x.p;
3638 reorder_terms_about(p, &p->arr[0]);
3639 value_clear(e->d);
3640 *e = p->arr[1];
3641 free(p);
3642 }
3643
3644 /* Approximate the evalue in fractional representation by a polynomial.
3645 * If sign > 0, the result is an upper bound;
3646 * if sign < 0, the result is a lower bound;
3647 * if sign = 0, the result is an intermediate approximation.
3648 */
3649 void evalue_frac2polynomial(evalue *e, int sign, unsigned MaxRays)
3650 {
3651 int i, j, k, dim;
3652 int *signs;
3653
3654 if (value_notzero_p(e->d))
3655 return;
3656 assert(e->x.p->type == partition);
3657 /* make sure all variables in the domains have a fixed sign */
3658 if (sign)
3659 evalue_split_domains_into_orthants(e, MaxRays);
3660
3661 assert(e->x.p->size >= 2);
3662 dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
3663
3664 signs = alloca(sizeof(int) * dim);
3665
3666 for (i = 0; i < e->x.p->size/2; ++i) {
3667 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
3668 POL_ENSURE_VERTICES(D);
3669 for (j = 0; j < dim; ++j) {
3670 signs[j] = 0;
3671 if (!sign)
3672 continue;
3673 for (k = 0; k < D->NbRays; ++k) {
3674 signs[j] = value_sign(D->Ray[k][1+j]);
3675 if (signs[j])
3676 break;
3677 }
3678 }
3679 evalue_frac2polynomial_r(&e->x.p->arr[2*i+1], signs, sign, 0);
3680 }
3681 }
3682
3683 /* Split the domains of e (which is assumed to be a partition)
3684 * such that each resulting domain lies entirely in one orthant.
3685 */
3686 void evalue_split_domains_into_orthants(evalue *e, unsigned MaxRays)
3687 {
3688 int i, dim;
3689 assert(value_zero_p(e->d));
3690 assert(e->x.p->type == partition);
3691 assert(e->x.p->size >= 2);
3692 dim = EVALUE_DOMAIN(e->x.p->arr[0])->Dimension;
3693
3694 for (i = 0; i < dim; ++i) {
3695 evalue split;
3696 Matrix *C, *C2;
3697 C = Matrix_Alloc(1, 1 + dim + 1);
3698 value_set_si(C->p[0][0], 1);
3699 value_init(split.d);
3700 value_set_si(split.d, 0);
3701 split.x.p = new_enode(partition, 4, dim);
3702 value_set_si(C->p[0][1+i], 1);
3703 C2 = Matrix_Copy(C);
3704 EVALUE_SET_DOMAIN(split.x.p->arr[0], Constraints2Polyhedron(C2, MaxRays));
3705 Matrix_Free(C2);
3706 evalue_set_si(&split.x.p->arr[1], 1, 1);
3707 value_set_si(C->p[0][1+i], -1);
3708 value_set_si(C->p[0][1+dim], -1);
3709 EVALUE_SET_DOMAIN(split.x.p->arr[2], Constraints2Polyhedron(C, MaxRays));
3710 evalue_set_si(&split.x.p->arr[3], 1, 1);
3711 emul(&split, e);
3712 free_evalue_refs(&split);
3713 Matrix_Free(C);
3714 }
3715 }
3716
3717 static Matrix *find_fractional_with_max_periods(evalue *e, Polyhedron *D,
3718 int max_periods,
3719 Value *min, Value *max)
3720 {
3721 Matrix *T;
3722 Value d;
3723 int i;
3724
3725 if (value_notzero_p(e->d))
3726 return NULL;
3727
3728 if (e->x.p->type == fractional) {
3729 Polyhedron *I;
3730 int bounded;
3731
3732 value_init(d);
3733 I = polynomial_projection(e->x.p, D, &d, &T);
3734 bounded = line_minmax(I, min, max); /* frees I */
3735 if (bounded) {
3736 Value mp;
3737 value_init(mp);
3738 value_set_si(mp, max_periods);
3739 mpz_fdiv_q(*min, *min, d);
3740 mpz_fdiv_q(*max, *max, d);
3741 value_assign(T->p[1][D->Dimension], d);
3742 value_subtract(d, *max, *min);
3743 if (value_ge(d, mp)) {
3744 Matrix_Free(T);
3745 T = NULL;
3746 }
3747 value_clear(mp);
3748 } else {
3749 Matrix_Free(T);
3750 T = NULL;
3751 }
3752 value_clear(d);
3753 if (T)
3754 return T;
3755 }
3756
3757 for (i = type_offset(e->x.p); i < e->x.p->size; ++i)
3758 if ((T = find_fractional_with_max_periods(&e->x.p->arr[i], D, max_periods,
3759 min, max)))
3760 return T;
3761
3762 return NULL;
3763 }
3764
3765 /* Look for fractional parts that can be removed by splitting the corresponding
3766 * domain into at most max_periods parts.
3767 * We use a very simply strategy that looks for the first fractional part
3768 * that satisfies the condition, performs the split and then continues
3769 * looking for other fractional parts in the split domains until no
3770 * such fractional part can be found anymore.
3771 */
3772 void evalue_split_periods(evalue *e, int max_periods, unsigned int MaxRays)
3773 {
3774 int i, j, n;
3775 Value min;
3776 Value max;
3777 Value d;
3778
3779 if (EVALUE_IS_ZERO(*e))
3780 return;
3781 if (value_notzero_p(e->d) || e->x.p->type != partition) {
3782 fprintf(stderr,
3783 "WARNING: evalue_split_periods called on incorrect evalue type\n");
3784 return;
3785 }
3786
3787 value_init(min);
3788 value_init(max);
3789 value_init(d);
3790
3791 for (i = 0; i < e->x.p->size/2; ++i) {
3792 enode *p;
3793 Matrix *T = NULL;
3794 Matrix *M;
3795 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
3796 Polyhedron *E;
3797 T = find_fractional_with_max_periods(&e->x.p->arr[2*i+1], D, max_periods,
3798 &min, &max);
3799 if (!T)
3800 continue;
3801
3802 M = Matrix_Alloc(2, 2+D->Dimension);
3803
3804 value_subtract(d, max, min);
3805 n = VALUE_TO_INT(d)+1;
3806
3807 value_set_si(M->p[0][0], 1);
3808 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
3809 value_multiply(d, max, T->p[1][D->Dimension]);
3810 value_subtract(M->p[0][1+D->Dimension], M->p[0][1+D->Dimension], d);
3811 value_set_si(d, -1);
3812 value_set_si(M->p[1][0], 1);
3813 Vector_Scale(T->p[0], M->p[1]+1, d, D->Dimension+1);
3814 value_addmul(M->p[1][1+D->Dimension], max, T->p[1][D->Dimension]);
3815 value_addto(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension],
3816 T->p[1][D->Dimension]);
3817 value_decrement(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension]);
3818
3819 p = new_enode(partition, e->x.p->size + (n-1)*2, e->x.p->pos);
3820 for (j = 0; j < 2*i; ++j) {
3821 value_clear(p->arr[j].d);
3822 p->arr[j] = e->x.p->arr[j];
3823 }
3824 for (j = 2*i+2; j < e->x.p->size; ++j) {
3825 value_clear(p->arr[j+2*(n-1)].d);
3826 p->arr[j+2*(n-1)] = e->x.p->arr[j];
3827 }
3828 for (j = n-1; j >= 0; --j) {
3829 if (j == 0) {
3830 value_clear(p->arr[2*i+1].d);
3831 p->arr[2*i+1] = e->x.p->arr[2*i+1];
3832 } else
3833 evalue_copy(&p->arr[2*(i+j)+1], &e->x.p->arr[2*i+1]);
3834 if (j != n-1) {
3835 value_subtract(M->p[1][1+D->Dimension], M->p[1][1+D->Dimension],
3836 T->p[1][D->Dimension]);
3837 value_addto(M->p[0][1+D->Dimension], M->p[0][1+D->Dimension],
3838 T->p[1][D->Dimension]);
3839 }
3840 E = DomainAddConstraints(D, M, MaxRays);
3841 EVALUE_SET_DOMAIN(p->arr[2*(i+j)], E);
3842 if (evalue_range_reduction_in_domain(&p->arr[2*(i+j)+1], E))
3843 reduce_evalue(&p->arr[2*(i+j)+1]);
3844 }
3845 value_clear(e->x.p->arr[2*i].d);
3846 Domain_Free(D);
3847 Matrix_Free(M);
3848 Matrix_Free(T);
3849 free(e->x.p);
3850 e->x.p = p;
3851 --i;
3852 }
3853
3854 value_clear(d);
3855 value_clear(min);
3856 value_clear(max);
3857 }
3858
3859 void evalue_extract_affine(const evalue *e, Value *coeff, Value *cst, Value *d)
3860 {
3861 value_set_si(*d, 1);
3862 evalue_denom(e, d);
3863 for ( ; value_zero_p(e->d); e = &e->x.p->arr[0]) {
3864 evalue *c;
3865 assert(e->x.p->type == polynomial);
3866 assert(e->x.p->size == 2);
3867 c = &e->x.p->arr[1];
3868 value_multiply(coeff[e->x.p->pos-1], *d, c->x.n);
3869 value_division(coeff[e->x.p->pos-1], coeff[e->x.p->pos-1], c->d);
3870 }
3871 value_multiply(*cst, *d, e->x.n);
3872 value_division(*cst, *cst, e->d);
3873 }
3874
3875 /* returns an evalue that corresponds to
3876 *
3877 * c/den x_param
3878 */
3879 static evalue *term(int param, Value c, Value den)
3880 {
3881 evalue *EP = ALLOC(evalue);
3882 value_init(EP->d);
3883 value_set_si(EP->d,0);
3884 EP->x.p = new_enode(polynomial, 2, param + 1);
3885 evalue_set_si(&EP->x.p->arr[0], 0, 1);
3886 value_init(EP->x.p->arr[1].x.n);
3887 value_assign(EP->x.p->arr[1].d, den);
3888 value_assign(EP->x.p->arr[1].x.n, c);
3889 return EP;
3890 }
3891
3892 evalue *affine2evalue(Value *coeff, Value denom, int nvar)
3893 {
3894 int i;
3895 evalue *E = ALLOC(evalue);
3896 value_init(E->d);
3897 evalue_set(E, coeff[nvar], denom);
3898 for (i = 0; i < nvar; ++i) {
3899 evalue *t = term(i, coeff[i], denom);
3900 eadd(t, E);
3901 free_evalue_refs(t);
3902 free(t);
3903 }
3904 return E;
3905 }
3906
3907 void evalue_substitute(evalue *e, evalue **subs)
3908 {
3909 int i, offset;
3910 evalue *v;
3911 enode *p;
3912
3913 if (value_notzero_p(e->d))
3914 return;
3915
3916 p = e->x.p;
3917 assert(p->type != partition);
3918
3919 for (i = 0; i < p->size; ++i)
3920 evalue_substitute(&p->arr[i], subs);
3921
3922 if (p->type == polynomial)
3923 v = subs[p->pos-1];
3924 else {
3925 v = ALLOC(evalue);
3926 value_init(v->d);
3927 value_set_si(v->d, 0);
3928 v->x.p = new_enode(p->type, 3, -1);
3929 value_clear(v->x.p->arr[0].d);
3930 v->x.p->arr[0] = p->arr[0];
3931 evalue_set_si(&v->x.p->arr[1], 0, 1);
3932 evalue_set_si(&v->x.p->arr[2], 1, 1);
3933 }
3934
3935 offset = type_offset(p);
3936
3937 for (i = p->size-1; i >= offset+1; i--) {
3938 emul(v, &p->arr[i]);
3939 eadd(&p->arr[i], &p->arr[i-1]);
3940 free_evalue_refs(&(p->arr[i]));
3941 }
3942
3943 if (p->type != polynomial) {
3944 free_evalue_refs(v);
3945 free(v);
3946 }
3947
3948 value_clear(e->d);
3949 *e = p->arr[offset];
3950 free(p);
3951 }
3952
3953 /* evalue e is given in terms of "new" parameter; CP maps the new
3954 * parameters back to the old parameters.
3955 * Transforms e such that it refers back to the old parameters.
3956 */
3957 void evalue_backsubstitute(evalue *e, Matrix *CP, unsigned MaxRays)
3958 {
3959 Matrix *eq;
3960 Matrix *inv;
3961 evalue **subs;
3962 enode *p;
3963 int i;
3964 unsigned nparam = CP->NbColumns-1;
3965 Polyhedron *CEq;
3966
3967 if (EVALUE_IS_ZERO(*e))
3968 return;
3969
3970 assert(value_zero_p(e->d));
3971 p = e->x.p;
3972 assert(p->type == partition);
3973
3974 inv = left_inverse(CP, &eq);
3975 subs = ALLOCN(evalue *, nparam);
3976 for (i = 0; i < nparam; ++i)
3977 subs[i] = affine2evalue(inv->p[i], inv->p[nparam][inv->NbColumns-1],
3978 inv->NbColumns-1);
3979
3980 CEq = Constraints2Polyhedron(eq, MaxRays);
3981 addeliminatedparams_partition(p, inv, CEq, inv->NbColumns-1, MaxRays);
3982 Polyhedron_Free(CEq);
3983
3984 for (i = 0; i < p->size/2; ++i)
3985 evalue_substitute(&p->arr[2*i+1], subs);
3986
3987 for (i = 0; i < nparam; ++i) {
3988 free_evalue_refs(subs[i]);
3989 free(subs[i]);
3990 }
3991 free(subs);
3992 Matrix_Free(eq);
3993 Matrix_Free(inv);
3994 }
lattice point enumeration library