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