?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
evalue.c: clean up in_domain
[barvinok.git] / evalue.c
blob4600cbca03aedc1c5a6ad326834762a5ef433291
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 row, in = 1;
2078 Value v; /* value of the constraint of a row when
2079 parameters are instantiated*/
2080
2081 value_init(v);
2082
2083 for (row = 0; row < P->NbConstraints; row++) {
2084 Inner_Product(P->Constraint[row]+1, list_args, P->Dimension, &v);
2085 value_addto(v, v, P->Constraint[row][P->Dimension+1]); /*constant part*/
2086 if (value_neg_p(v) ||
2087 value_zero_p(P->Constraint[row][0]) && value_notzero_p(v)) {
2088 in = 0;
2089 break;
2090 }
2091 }
2092
2093 value_clear(v);
2094 return in || (P->next && in_domain(P->next, list_args));
2095 } /* in_domain */
2096
2097 /****************************************************/
2098 /* function compute enode */
2099 /* compute the value of enode p with parameters */
2100 /* list "list_args */
2101 /* compute the polynomial or the periodic */
2102 /****************************************************/
2103
2104 static double compute_enode(enode *p, Value *list_args) {
2105
2106 int i;
2107 Value m, param;
2108 double res=0.0;
2109
2110 if (!p)
2111 return(0.);
2112
2113 value_init(m);
2114 value_init(param);
2115
2116 if (p->type == polynomial) {
2117 if (p->size > 1)
2118 value_assign(param,list_args[p->pos-1]);
2119
2120 /* Compute the polynomial using Horner's rule */
2121 for (i=p->size-1;i>0;i--) {
2122 res +=compute_evalue(&p->arr[i],list_args);
2123 res *=VALUE_TO_DOUBLE(param);
2124 }
2125 res +=compute_evalue(&p->arr[0],list_args);
2126 }
2127 else if (p->type == fractional) {
2128 double d = compute_evalue(&p->arr[0], list_args);
2129 d -= floor(d+1e-10);
2130
2131 /* Compute the polynomial using Horner's rule */
2132 for (i=p->size-1;i>1;i--) {
2133 res +=compute_evalue(&p->arr[i],list_args);
2134 res *=d;
2135 }
2136 res +=compute_evalue(&p->arr[1],list_args);
2137 }
2138 else if (p->type == flooring) {
2139 double d = compute_evalue(&p->arr[0], list_args);
2140 d = floor(d+1e-10);
2141
2142 /* Compute the polynomial using Horner's rule */
2143 for (i=p->size-1;i>1;i--) {
2144 res +=compute_evalue(&p->arr[i],list_args);
2145 res *=d;
2146 }
2147 res +=compute_evalue(&p->arr[1],list_args);
2148 }
2149 else if (p->type == periodic) {
2150 value_assign(m,list_args[p->pos-1]);
2151
2152 /* Choose the right element of the periodic */
2153 value_set_si(param,p->size);
2154 value_pmodulus(m,m,param);
2155 res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
2156 }
2157 else if (p->type == relation) {
2158 if (fabs(compute_evalue(&p->arr[0], list_args)) < 1e-10)
2159 res = compute_evalue(&p->arr[1], list_args);
2160 else if (p->size > 2)
2161 res = compute_evalue(&p->arr[2], list_args);
2162 }
2163 else if (p->type == partition) {
2164 int dim = EVALUE_DOMAIN(p->arr[0])->Dimension;
2165 Value *vals = list_args;
2166 if (p->pos < dim) {
2167 NALLOC(vals, dim);
2168 for (i = 0; i < dim; ++i) {
2169 value_init(vals[i]);
2170 if (i < p->pos)
2171 value_assign(vals[i], list_args[i]);
2172 }
2173 }
2174 for (i = 0; i < p->size/2; ++i)
2175 if (DomainContains(EVALUE_DOMAIN(p->arr[2*i]), vals, p->pos, 0, 1)) {
2176 res = compute_evalue(&p->arr[2*i+1], vals);
2177 break;
2178 }
2179 if (p->pos < dim) {
2180 for (i = 0; i < dim; ++i)
2181 value_clear(vals[i]);
2182 free(vals);
2183 }
2184 }
2185 else
2186 assert(0);
2187 value_clear(m);
2188 value_clear(param);
2189 return res;
2190 } /* compute_enode */
2191
2192 /*************************************************/
2193 /* return the value of Ehrhart Polynomial */
2194 /* It returns a double, because since it is */
2195 /* a recursive function, some intermediate value */
2196 /* might not be integral */
2197 /*************************************************/
2198
2199 double compute_evalue(evalue *e,Value *list_args) {
2200
2201 double res;
2202
2203 if (value_notzero_p(e->d)) {
2204 if (value_notone_p(e->d))
2205 res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
2206 else
2207 res = VALUE_TO_DOUBLE(e->x.n);
2208 }
2209 else
2210 res = compute_enode(e->x.p,list_args);
2211 return res;
2212 } /* compute_evalue */
2213
2214
2215 /****************************************************/
2216 /* function compute_poly : */
2217 /* Check for the good validity domain */
2218 /* return the number of point in the Polyhedron */
2219 /* in allocated memory */
2220 /* Using the Ehrhart pseudo-polynomial */
2221 /****************************************************/
2222 Value *compute_poly(Enumeration *en,Value *list_args) {
2223
2224 Value *tmp;
2225 /* double d; int i; */
2226
2227 tmp = (Value *) malloc (sizeof(Value));
2228 assert(tmp != NULL);
2229 value_init(*tmp);
2230 value_set_si(*tmp,0);
2231
2232 if(!en)
2233 return(tmp); /* no ehrhart polynomial */
2234 if(en->ValidityDomain) {
2235 if(!en->ValidityDomain->Dimension) { /* no parameters */
2236 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2237 return(tmp);
2238 }
2239 }
2240 else
2241 return(tmp); /* no Validity Domain */
2242 while(en) {
2243 if(in_domain(en->ValidityDomain,list_args)) {
2244
2245 #ifdef EVAL_EHRHART_DEBUG
2246 Print_Domain(stdout,en->ValidityDomain);
2247 print_evalue(stdout,&en->EP);
2248 #endif
2249
2250 /* d = compute_evalue(&en->EP,list_args);
2251 i = d;
2252 printf("(double)%lf = %d\n", d, i ); */
2253 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2254 return(tmp);
2255 }
2256 else
2257 en=en->next;
2258 }
2259 value_set_si(*tmp,0);
2260 return(tmp); /* no compatible domain with the arguments */
2261 } /* compute_poly */
2262
2263 size_t value_size(Value v) {
2264 return (v[0]._mp_size > 0 ? v[0]._mp_size : -v[0]._mp_size)
2265 * sizeof(v[0]._mp_d[0]);
2266 }
2267
2268 size_t domain_size(Polyhedron *D)
2269 {
2270 int i, j;
2271 size_t s = sizeof(*D);
2272
2273 for (i = 0; i < D->NbConstraints; ++i)
2274 for (j = 0; j < D->Dimension+2; ++j)
2275 s += value_size(D->Constraint[i][j]);
2276
2277 /*
2278 for (i = 0; i < D->NbRays; ++i)
2279 for (j = 0; j < D->Dimension+2; ++j)
2280 s += value_size(D->Ray[i][j]);
2281 */
2282
2283 return D->next ? s+domain_size(D->next) : s;
2284 }
2285
2286 size_t enode_size(enode *p) {
2287 size_t s = sizeof(*p) - sizeof(p->arr[0]);
2288 int i;
2289
2290 if (p->type == partition)
2291 for (i = 0; i < p->size/2; ++i) {
2292 s += domain_size(EVALUE_DOMAIN(p->arr[2*i]));
2293 s += evalue_size(&p->arr[2*i+1]);
2294 }
2295 else
2296 for (i = 0; i < p->size; ++i) {
2297 s += evalue_size(&p->arr[i]);
2298 }
2299 return s;
2300 }
2301
2302 size_t evalue_size(evalue *e)
2303 {
2304 size_t s = sizeof(*e);
2305 s += value_size(e->d);
2306 if (value_notzero_p(e->d))
2307 s += value_size(e->x.n);
2308 else
2309 s += enode_size(e->x.p);
2310 return s;
2311 }
2312
2313 static evalue *find_second(evalue *base, evalue *cst, evalue *e, Value m)
2314 {
2315 evalue *found = NULL;
2316 evalue offset;
2317 evalue copy;
2318 int i;
2319
2320 if (value_pos_p(e->d) || e->x.p->type != fractional)
2321 return NULL;
2322
2323 value_init(offset.d);
2324 value_init(offset.x.n);
2325 poly_denom(&e->x.p->arr[0], &offset.d);
2326 value_lcm(m, offset.d, &offset.d);
2327 value_set_si(offset.x.n, 1);
2328
2329 value_init(copy.d);
2330 evalue_copy(&copy, cst);
2331
2332 eadd(&offset, cst);
2333 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2334
2335 if (eequal(base, &e->x.p->arr[0]))
2336 found = &e->x.p->arr[0];
2337 else {
2338 value_set_si(offset.x.n, -2);
2339
2340 eadd(&offset, cst);
2341 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2342
2343 if (eequal(base, &e->x.p->arr[0]))
2344 found = base;
2345 }
2346 free_evalue_refs(cst);
2347 free_evalue_refs(&offset);
2348 *cst = copy;
2349
2350 for (i = 1; !found && i < e->x.p->size; ++i)
2351 found = find_second(base, cst, &e->x.p->arr[i], m);
2352
2353 return found;
2354 }
2355
2356 static evalue *find_relation_pair(evalue *e)
2357 {
2358 int i;
2359 evalue *found = NULL;
2360
2361 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2362 return NULL;
2363
2364 if (e->x.p->type == fractional) {
2365 Value m;
2366 evalue *cst;
2367
2368 value_init(m);
2369 poly_denom(&e->x.p->arr[0], &m);
2370
2371 for (cst = &e->x.p->arr[0]; value_zero_p(cst->d);
2372 cst = &cst->x.p->arr[0])
2373 ;
2374
2375 for (i = 1; !found && i < e->x.p->size; ++i)
2376 found = find_second(&e->x.p->arr[0], cst, &e->x.p->arr[i], m);
2377
2378 value_clear(m);
2379 }
2380
2381 i = e->x.p->type == relation;
2382 for (; !found && i < e->x.p->size; ++i)
2383 found = find_relation_pair(&e->x.p->arr[i]);
2384
2385 return found;
2386 }
2387
2388 void evalue_mod2relation(evalue *e) {
2389 evalue *d;
2390
2391 if (value_zero_p(e->d) && e->x.p->type == partition) {
2392 int i;
2393
2394 for (i = 0; i < e->x.p->size/2; ++i) {
2395 evalue_mod2relation(&e->x.p->arr[2*i+1]);
2396 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2397 value_clear(e->x.p->arr[2*i].d);
2398 free_evalue_refs(&e->x.p->arr[2*i+1]);
2399 e->x.p->size -= 2;
2400 if (2*i < e->x.p->size) {
2401 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2402 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2403 }
2404 --i;
2405 }
2406 }
2407 if (e->x.p->size == 0) {
2408 free(e->x.p);
2409 evalue_set_si(e, 0, 1);
2410 }
2411
2412 return;
2413 }
2414
2415 while ((d = find_relation_pair(e)) != NULL) {
2416 evalue split;
2417 evalue *ev;
2418
2419 value_init(split.d);
2420 value_set_si(split.d, 0);
2421 split.x.p = new_enode(relation, 3, 0);
2422 evalue_set_si(&split.x.p->arr[1], 1, 1);
2423 evalue_set_si(&split.x.p->arr[2], 1, 1);
2424
2425 ev = &split.x.p->arr[0];
2426 value_set_si(ev->d, 0);
2427 ev->x.p = new_enode(fractional, 3, -1);
2428 evalue_set_si(&ev->x.p->arr[1], 0, 1);
2429 evalue_set_si(&ev->x.p->arr[2], 1, 1);
2430 evalue_copy(&ev->x.p->arr[0], d);
2431
2432 emul(&split, e);
2433
2434 reduce_evalue(e);
2435
2436 free_evalue_refs(&split);
2437 }
2438 }
2439
2440 static int evalue_comp(const void * a, const void * b)
2441 {
2442 const evalue *e1 = *(const evalue **)a;
2443 const evalue *e2 = *(const evalue **)b;
2444 return ecmp(e1, e2);
2445 }
2446
2447 void evalue_combine(evalue *e)
2448 {
2449 evalue **evs;
2450 int i, k;
2451 enode *p;
2452 evalue tmp;
2453
2454 if (value_notzero_p(e->d) || e->x.p->type != partition)
2455 return;
2456
2457 NALLOC(evs, e->x.p->size/2);
2458 for (i = 0; i < e->x.p->size/2; ++i)
2459 evs[i] = &e->x.p->arr[2*i+1];
2460 qsort(evs, e->x.p->size/2, sizeof(evs[0]), evalue_comp);
2461 p = new_enode(partition, e->x.p->size, e->x.p->pos);
2462 for (i = 0, k = 0; i < e->x.p->size/2; ++i) {
2463 if (k == 0 || ecmp(&p->arr[2*k-1], evs[i]) != 0) {
2464 value_clear(p->arr[2*k].d);
2465 value_clear(p->arr[2*k+1].d);
2466 p->arr[2*k] = *(evs[i]-1);
2467 p->arr[2*k+1] = *(evs[i]);
2468 ++k;
2469 } else {
2470 Polyhedron *D = EVALUE_DOMAIN(*(evs[i]-1));
2471 Polyhedron *L = D;
2472
2473 value_clear((evs[i]-1)->d);
2474
2475 while (L->next)
2476 L = L->next;
2477 L->next = EVALUE_DOMAIN(p->arr[2*k-2]);
2478 EVALUE_SET_DOMAIN(p->arr[2*k-2], D);
2479 free_evalue_refs(evs[i]);
2480 }
2481 }
2482
2483 for (i = 2*k ; i < p->size; ++i)
2484 value_clear(p->arr[i].d);
2485
2486 free(evs);
2487 free(e->x.p);
2488 p->size = 2*k;
2489 e->x.p = p;
2490
2491 for (i = 0; i < e->x.p->size/2; ++i) {
2492 Polyhedron *H;
2493 if (value_notzero_p(e->x.p->arr[2*i+1].d))
2494 continue;
2495 H = DomainConvex(EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
2496 if (H == NULL)
2497 continue;
2498 for (k = 0; k < e->x.p->size/2; ++k) {
2499 Polyhedron *D, *N, **P;
2500 if (i == k)
2501 continue;
2502 P = &EVALUE_DOMAIN(e->x.p->arr[2*k]);
2503 D = *P;
2504 if (D == NULL)
2505 continue;
2506 for (; D; D = N) {
2507 *P = D;
2508 N = D->next;
2509 if (D->NbEq <= H->NbEq) {
2510 P = &D->next;
2511 continue;
2512 }
2513
2514 value_init(tmp.d);
2515 tmp.x.p = new_enode(partition, 2, e->x.p->pos);
2516 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Copy(D));
2517 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
2518 reduce_evalue(&tmp);
2519 if (value_notzero_p(tmp.d) ||
2520 ecmp(&tmp.x.p->arr[1], &e->x.p->arr[2*k+1]) != 0)
2521 P = &D->next;
2522 else {
2523 D->next = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2524 EVALUE_DOMAIN(e->x.p->arr[2*i]) = D;
2525 *P = NULL;
2526 }
2527 free_evalue_refs(&tmp);
2528 }
2529 }
2530 Polyhedron_Free(H);
2531 }
2532
2533 for (i = 0; i < e->x.p->size/2; ++i) {
2534 Polyhedron *H, *E;
2535 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2536 if (!D) {
2537 value_clear(e->x.p->arr[2*i].d);
2538 free_evalue_refs(&e->x.p->arr[2*i+1]);
2539 e->x.p->size -= 2;
2540 if (2*i < e->x.p->size) {
2541 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2542 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2543 }
2544 --i;
2545 continue;
2546 }
2547 if (!D->next)
2548 continue;
2549 H = DomainConvex(D, 0);
2550 E = DomainDifference(H, D, 0);
2551 Domain_Free(D);
2552 D = DomainDifference(H, E, 0);
2553 Domain_Free(H);
2554 Domain_Free(E);
2555 EVALUE_SET_DOMAIN(p->arr[2*i], D);
2556 }
2557 }
2558
2559 /* May change coefficients to become non-standard if fiddle is set
2560 * => reduce p afterwards to correct
2561 */
2562 static Polyhedron *polynomial_projection(enode *p, Polyhedron *D, Value *d,
2563 Matrix **R, int fiddle)
2564 {
2565 Polyhedron *I, *H;
2566 evalue *pp;
2567 unsigned dim = D->Dimension;
2568 Matrix *T = Matrix_Alloc(2, dim+1);
2569 Value twice;
2570 value_init(twice);
2571 assert(T);
2572
2573 assert(p->type == fractional);
2574 pp = &p->arr[0];
2575 value_set_si(T->p[1][dim], 1);
2576 poly_denom(pp, d);
2577 while (value_zero_p(pp->d)) {
2578 assert(pp->x.p->type == polynomial);
2579 assert(pp->x.p->size == 2);
2580 assert(value_notzero_p(pp->x.p->arr[1].d));
2581 value_division(T->p[0][pp->x.p->pos-1], *d, pp->x.p->arr[1].d);
2582 mpz_mul_ui(twice, pp->x.p->arr[1].x.n, 2);
2583 if (fiddle && value_gt(twice, pp->x.p->arr[1].d))
2584 value_subtract(pp->x.p->arr[1].x.n,
2585 pp->x.p->arr[1].x.n, pp->x.p->arr[1].d);
2586 value_multiply(T->p[0][pp->x.p->pos-1],
2587 T->p[0][pp->x.p->pos-1], pp->x.p->arr[1].x.n);
2588 pp = &pp->x.p->arr[0];
2589 }
2590 value_division(T->p[0][dim], *d, pp->d);
2591 value_multiply(T->p[0][dim], T->p[0][dim], pp->x.n);
2592 I = DomainImage(D, T, 0);
2593 H = DomainConvex(I, 0);
2594 Domain_Free(I);
2595 *R = T;
2596
2597 value_clear(twice);
2598
2599 return H;
2600 }
2601
2602 int evalue_range_reduction_in_domain(evalue *e, Polyhedron *D)
2603 {
2604 int i;
2605 enode *p;
2606 Value d, min, max;
2607 int r = 0;
2608 Polyhedron *I;
2609 Matrix *T;
2610 int bounded;
2611
2612 if (value_notzero_p(e->d))
2613 return r;
2614
2615 p = e->x.p;
2616
2617 if (p->type == relation) {
2618 int equal;
2619 value_init(d);
2620 value_init(min);
2621 value_init(max);
2622
2623 I = polynomial_projection(p->arr[0].x.p, D, &d, &T, 1);
2624 bounded = line_minmax(I, &min, &max); /* frees I */
2625 equal = value_eq(min, max);
2626 mpz_cdiv_q(min, min, d);
2627 mpz_fdiv_q(max, max, d);
2628
2629 if (bounded && value_gt(min, max)) {
2630 /* Never zero */
2631 if (p->size == 3) {
2632 value_clear(e->d);
2633 *e = p->arr[2];
2634 } else {
2635 evalue_set_si(e, 0, 1);
2636 r = 1;
2637 }
2638 free_evalue_refs(&(p->arr[1]));
2639 free_evalue_refs(&(p->arr[0]));
2640 free(p);
2641 value_clear(d);
2642 value_clear(min);
2643 value_clear(max);
2644 Matrix_Free(T);
2645 return r ? r : evalue_range_reduction_in_domain(e, D);
2646 } else if (bounded && equal) {
2647 /* Always zero */
2648 if (p->size == 3)
2649 free_evalue_refs(&(p->arr[2]));
2650 value_clear(e->d);
2651 *e = p->arr[1];
2652 free_evalue_refs(&(p->arr[0]));
2653 free(p);
2654 value_clear(d);
2655 value_clear(min);
2656 value_clear(max);
2657 Matrix_Free(T);
2658 return evalue_range_reduction_in_domain(e, D);
2659 } else if (bounded && value_eq(min, max)) {
2660 /* zero for a single value */
2661 Polyhedron *E;
2662 Matrix *M = Matrix_Alloc(1, D->Dimension+2);
2663 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
2664 value_multiply(min, min, d);
2665 value_subtract(M->p[0][D->Dimension+1],
2666 M->p[0][D->Dimension+1], min);
2667 E = DomainAddConstraints(D, M, 0);
2668 value_clear(d);
2669 value_clear(min);
2670 value_clear(max);
2671 Matrix_Free(T);
2672 Matrix_Free(M);
2673 r = evalue_range_reduction_in_domain(&p->arr[1], E);
2674 if (p->size == 3)
2675 r |= evalue_range_reduction_in_domain(&p->arr[2], D);
2676 Domain_Free(E);
2677 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2678 return r;
2679 }
2680
2681 value_clear(d);
2682 value_clear(min);
2683 value_clear(max);
2684 Matrix_Free(T);
2685 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2686 }
2687
2688 i = p->type == relation ? 1 :
2689 p->type == fractional ? 1 : 0;
2690 for (; i<p->size; i++)
2691 r |= evalue_range_reduction_in_domain(&p->arr[i], D);
2692
2693 if (p->type != fractional) {
2694 if (r && p->type == polynomial) {
2695 evalue f;
2696 value_init(f.d);
2697 value_set_si(f.d, 0);
2698 f.x.p = new_enode(polynomial, 2, p->pos);
2699 evalue_set_si(&f.x.p->arr[0], 0, 1);
2700 evalue_set_si(&f.x.p->arr[1], 1, 1);
2701 reorder_terms(p, &f);
2702 value_clear(e->d);
2703 *e = p->arr[0];
2704 free(p);
2705 }
2706 return r;
2707 }
2708
2709 value_init(d);
2710 value_init(min);
2711 value_init(max);
2712 I = polynomial_projection(p, D, &d, &T, 1);
2713 Matrix_Free(T);
2714 bounded = line_minmax(I, &min, &max); /* frees I */
2715 mpz_fdiv_q(min, min, d);
2716 mpz_fdiv_q(max, max, d);
2717 value_subtract(d, max, min);
2718
2719 if (bounded && value_eq(min, max)) {
2720 evalue inc;
2721 value_init(inc.d);
2722 value_init(inc.x.n);
2723 value_set_si(inc.d, 1);
2724 value_oppose(inc.x.n, min);
2725 eadd(&inc, &p->arr[0]);
2726 reorder_terms(p, &p->arr[0]); /* frees arr[0] */
2727 value_clear(e->d);
2728 *e = p->arr[1];
2729 free(p);
2730 free_evalue_refs(&inc);
2731 r = 1;
2732 } else if (bounded && value_one_p(d) && p->size > 3) {
2733 evalue rem;
2734 evalue inc;
2735 evalue t;
2736 evalue f;
2737 evalue factor;
2738 value_init(rem.d);
2739 value_set_si(rem.d, 0);
2740 rem.x.p = new_enode(fractional, 3, -1);
2741 evalue_copy(&rem.x.p->arr[0], &p->arr[0]);
2742 value_clear(rem.x.p->arr[1].d);
2743 value_clear(rem.x.p->arr[2].d);
2744 rem.x.p->arr[1] = p->arr[1];
2745 rem.x.p->arr[2] = p->arr[2];
2746 for (i = 3; i < p->size; ++i)
2747 p->arr[i-2] = p->arr[i];
2748 p->size -= 2;
2749
2750 value_init(inc.d);
2751 value_init(inc.x.n);
2752 value_set_si(inc.d, 1);
2753 value_oppose(inc.x.n, min);
2754
2755 value_init(t.d);
2756 evalue_copy(&t, &p->arr[0]);
2757 eadd(&inc, &t);
2758
2759 value_init(f.d);
2760 value_set_si(f.d, 0);
2761 f.x.p = new_enode(fractional, 3, -1);
2762 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
2763 evalue_set_si(&f.x.p->arr[1], 1, 1);
2764 evalue_set_si(&f.x.p->arr[2], 2, 1);
2765
2766 value_init(factor.d);
2767 evalue_set_si(&factor, -1, 1);
2768 emul(&t, &factor);
2769
2770 eadd(&f, &factor);
2771 emul(&t, &factor);
2772
2773 value_clear(f.x.p->arr[1].x.n);
2774 value_clear(f.x.p->arr[2].x.n);
2775 evalue_set_si(&f.x.p->arr[1], 0, 1);
2776 evalue_set_si(&f.x.p->arr[2], -1, 1);
2777 eadd(&f, &factor);
2778
2779 emul(&factor, e);
2780 eadd(&rem, e);
2781
2782 free_evalue_refs(&inc);
2783 free_evalue_refs(&t);
2784 free_evalue_refs(&f);
2785 free_evalue_refs(&factor);
2786 free_evalue_refs(&rem);
2787
2788 evalue_range_reduction_in_domain(e, D);
2789
2790 r = 1;
2791 } else {
2792 _reduce_evalue(&p->arr[0], 0, 1);
2793 if (r == 1) {
2794 evalue f;
2795 value_init(f.d);
2796 value_set_si(f.d, 0);
2797 f.x.p = new_enode(fractional, 3, -1);
2798 value_clear(f.x.p->arr[0].d);
2799 f.x.p->arr[0] = p->arr[0];
2800 evalue_set_si(&f.x.p->arr[1], 0, 1);
2801 evalue_set_si(&f.x.p->arr[2], 1, 1);
2802 reorder_terms(p, &f);
2803 value_clear(e->d);
2804 *e = p->arr[1];
2805 free(p);
2806 }
2807 }
2808
2809 value_clear(d);
2810 value_clear(min);
2811 value_clear(max);
2812
2813 return r;
2814 }
2815
2816 void evalue_range_reduction(evalue *e)
2817 {
2818 int i;
2819 if (value_notzero_p(e->d) || e->x.p->type != partition)
2820 return;
2821
2822 for (i = 0; i < e->x.p->size/2; ++i)
2823 if (evalue_range_reduction_in_domain(&e->x.p->arr[2*i+1],
2824 EVALUE_DOMAIN(e->x.p->arr[2*i]))) {
2825 reduce_evalue(&e->x.p->arr[2*i+1]);
2826
2827 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2828 free_evalue_refs(&e->x.p->arr[2*i+1]);
2829 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
2830 value_clear(e->x.p->arr[2*i].d);
2831 e->x.p->size -= 2;
2832 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2833 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2834 --i;
2835 }
2836 }
2837 }
2838
2839 /* Frees EP
2840 */
2841 Enumeration* partition2enumeration(evalue *EP)
2842 {
2843 int i;
2844 Enumeration *en, *res = NULL;
2845
2846 if (EVALUE_IS_ZERO(*EP)) {
2847 free(EP);
2848 return res;
2849 }
2850
2851 for (i = 0; i < EP->x.p->size/2; ++i) {
2852 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[2*i])->Dimension);
2853 en = (Enumeration *)malloc(sizeof(Enumeration));
2854 en->next = res;
2855 res = en;
2856 res->ValidityDomain = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
2857 value_clear(EP->x.p->arr[2*i].d);
2858 res->EP = EP->x.p->arr[2*i+1];
2859 }
2860 free(EP->x.p);
2861 value_clear(EP->d);
2862 free(EP);
2863 return res;
2864 }
2865
2866 int evalue_frac2floor_in_domain(evalue *e, Polyhedron *D)
2867 {
2868 enode *p;
2869 int r = 0;
2870 int i;
2871 Polyhedron *I;
2872 Matrix *T;
2873 Value d, min;
2874 evalue fl;
2875
2876 if (value_notzero_p(e->d))
2877 return r;
2878
2879 p = e->x.p;
2880
2881 i = p->type == relation ? 1 :
2882 p->type == fractional ? 1 : 0;
2883 for (; i<p->size; i++)
2884 r |= evalue_frac2floor_in_domain(&p->arr[i], D);
2885
2886 if (p->type != fractional) {
2887 if (r && p->type == polynomial) {
2888 evalue f;
2889 value_init(f.d);
2890 value_set_si(f.d, 0);
2891 f.x.p = new_enode(polynomial, 2, p->pos);
2892 evalue_set_si(&f.x.p->arr[0], 0, 1);
2893 evalue_set_si(&f.x.p->arr[1], 1, 1);
2894 reorder_terms(p, &f);
2895 value_clear(e->d);
2896 *e = p->arr[0];
2897 free(p);
2898 }
2899 return r;
2900 }
2901
2902 value_init(d);
2903 I = polynomial_projection(p, D, &d, &T, 0);
2904
2905 /*
2906 Polyhedron_Print(stderr, P_VALUE_FMT, I);
2907 */
2908
2909 assert(I->NbEq == 0); /* Should have been reduced */
2910
2911 /* Find minimum */
2912 for (i = 0; i < I->NbConstraints; ++i)
2913 if (value_pos_p(I->Constraint[i][1]))
2914 break;
2915
2916 if (i < I->NbConstraints) {
2917 value_init(min);
2918 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
2919 mpz_cdiv_q(min, I->Constraint[i][2], I->Constraint[i][1]);
2920 if (value_neg_p(min)) {
2921 evalue offset;
2922 mpz_fdiv_q(min, min, d);
2923 value_init(offset.d);
2924 value_set_si(offset.d, 1);
2925 value_init(offset.x.n);
2926 value_oppose(offset.x.n, min);
2927 eadd(&offset, &p->arr[0]);
2928 free_evalue_refs(&offset);
2929 }
2930 value_clear(min);
2931 }
2932
2933 Polyhedron_Free(I);
2934 Matrix_Free(T);
2935 value_clear(d);
2936
2937 value_init(fl.d);
2938 value_set_si(fl.d, 0);
2939 fl.x.p = new_enode(flooring, 3, -1);
2940 evalue_set_si(&fl.x.p->arr[1], 0, 1);
2941 evalue_set_si(&fl.x.p->arr[2], -1, 1);
2942 evalue_copy(&fl.x.p->arr[0], &p->arr[0]);
2943
2944 eadd(&fl, &p->arr[0]);
2945 reorder_terms(p, &p->arr[0]);
2946 value_clear(e->d);
2947 *e = p->arr[1];
2948 free(p);
2949 free_evalue_refs(&fl);
2950
2951 return 1;
2952 }
2953
2954 void evalue_frac2floor(evalue *e)
2955 {
2956 int i;
2957 if (value_notzero_p(e->d) || e->x.p->type != partition)
2958 return;
2959
2960 for (i = 0; i < e->x.p->size/2; ++i)
2961 if (evalue_frac2floor_in_domain(&e->x.p->arr[2*i+1],
2962 EVALUE_DOMAIN(e->x.p->arr[2*i])))
2963 reduce_evalue(&e->x.p->arr[2*i+1]);
2964 }
2965
2966 static Matrix *esum_add_constraint(int nvar, Polyhedron *D, Matrix *C,
2967 Vector *row)
2968 {
2969 int nr, nc;
2970 int i;
2971 int nparam = D->Dimension - nvar;
2972
2973 if (C == 0) {
2974 nr = D->NbConstraints + 2;
2975 nc = D->Dimension + 2 + 1;
2976 C = Matrix_Alloc(nr, nc);
2977 for (i = 0; i < D->NbConstraints; ++i) {
2978 Vector_Copy(D->Constraint[i], C->p[i], 1 + nvar);
2979 Vector_Copy(D->Constraint[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
2980 D->Dimension + 1 - nvar);
2981 }
2982 } else {
2983 Matrix *oldC = C;
2984 nr = C->NbRows + 2;
2985 nc = C->NbColumns + 1;
2986 C = Matrix_Alloc(nr, nc);
2987 for (i = 0; i < oldC->NbRows; ++i) {
2988 Vector_Copy(oldC->p[i], C->p[i], 1 + nvar);
2989 Vector_Copy(oldC->p[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
2990 oldC->NbColumns - 1 - nvar);
2991 }
2992 }
2993 value_set_si(C->p[nr-2][0], 1);
2994 value_set_si(C->p[nr-2][1 + nvar], 1);
2995 value_set_si(C->p[nr-2][nc - 1], -1);
2996
2997 Vector_Copy(row->p, C->p[nr-1], 1 + nvar + 1);
2998 Vector_Copy(row->p + 1 + nvar + 1, C->p[nr-1] + C->NbColumns - 1 - nparam,
2999 1 + nparam);
3000
3001 return C;
3002 }
3003
3004 static void floor2frac_r(evalue *e, int nvar)
3005 {
3006 enode *p;
3007 int i;
3008 evalue f;
3009 evalue *pp;
3010
3011 if (value_notzero_p(e->d))
3012 return;
3013
3014 p = e->x.p;
3015
3016 assert(p->type == flooring);
3017 for (i = 1; i < p->size; i++)
3018 floor2frac_r(&p->arr[i], nvar);
3019
3020 for (pp = &p->arr[0]; value_zero_p(pp->d); pp = &pp->x.p->arr[0]) {
3021 assert(pp->x.p->type == polynomial);
3022 pp->x.p->pos -= nvar;
3023 }
3024
3025 value_init(f.d);
3026 value_set_si(f.d, 0);
3027 f.x.p = new_enode(fractional, 3, -1);
3028 evalue_set_si(&f.x.p->arr[1], 0, 1);
3029 evalue_set_si(&f.x.p->arr[2], -1, 1);
3030 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
3031
3032 eadd(&f, &p->arr[0]);
3033 reorder_terms(p, &p->arr[0]);
3034 value_clear(e->d);
3035 *e = p->arr[1];
3036 free(p);
3037 free_evalue_refs(&f);
3038 }
3039
3040 /* Convert flooring back to fractional and shift position
3041 * of the parameters by nvar
3042 */
3043 static void floor2frac(evalue *e, int nvar)
3044 {
3045 floor2frac_r(e, nvar);
3046 reduce_evalue(e);
3047 }
3048
3049 evalue *esum_over_domain_cst(int nvar, Polyhedron *D, Matrix *C)
3050 {
3051 evalue *t;
3052 int nparam = D->Dimension - nvar;
3053
3054 if (C != 0) {
3055 C = Matrix_Copy(C);
3056 D = Constraints2Polyhedron(C, 0);
3057 Matrix_Free(C);
3058 }
3059
3060 t = barvinok_enumerate_e(D, 0, nparam, 0);
3061
3062 /* Double check that D was not unbounded. */
3063 assert(!(value_pos_p(t->d) && value_neg_p(t->x.n)));
3064
3065 if (C != 0)
3066 Polyhedron_Free(D);
3067
3068 return t;
3069 }
3070
3071 evalue *esum_over_domain(evalue *e, int nvar, Polyhedron *D,
3072 Matrix *C)
3073 {
3074 Vector *row = NULL;
3075 int i;
3076 evalue *res;
3077 Matrix *origC;
3078 evalue *factor = NULL;
3079 evalue cum;
3080
3081 if (EVALUE_IS_ZERO(*e))
3082 return 0;
3083
3084 if (D->next) {
3085 Polyhedron *DD = Disjoint_Domain(D, 0, 0);
3086 Polyhedron *Q;
3087
3088 Q = DD;
3089 DD = Q->next;
3090 Q->next = 0;
3091
3092 res = esum_over_domain(e, nvar, Q, C);
3093 Polyhedron_Free(Q);
3094
3095 for (Q = DD; Q; Q = DD) {
3096 evalue *t;
3097
3098 DD = Q->next;
3099 Q->next = 0;
3100
3101 t = esum_over_domain(e, nvar, Q, C);
3102 Polyhedron_Free(Q);
3103
3104 if (!res)
3105 res = t;
3106 else if (t) {
3107 eadd(t, res);
3108 free_evalue_refs(t);
3109 free(t);
3110 }
3111 }
3112 return res;
3113 }
3114
3115 if (value_notzero_p(e->d)) {
3116 evalue *t;
3117
3118 t = esum_over_domain_cst(nvar, D, C);
3119
3120 if (!EVALUE_IS_ONE(*e))
3121 emul(e, t);
3122
3123 return t;
3124 }
3125
3126 switch (e->x.p->type) {
3127 case flooring: {
3128 evalue *pp = &e->x.p->arr[0];
3129
3130 if (pp->x.p->pos > nvar) {
3131 /* remainder is independent of the summated vars */
3132 evalue f;
3133 evalue *t;
3134
3135 value_init(f.d);
3136 evalue_copy(&f, e);
3137 floor2frac(&f, nvar);
3138
3139 t = esum_over_domain_cst(nvar, D, C);
3140
3141 emul(&f, t);
3142
3143 free_evalue_refs(&f);
3144
3145 return t;
3146 }
3147
3148 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3149 poly_denom(pp, &row->p[1 + nvar]);
3150 value_set_si(row->p[0], 1);
3151 for (pp = &e->x.p->arr[0]; value_zero_p(pp->d);
3152 pp = &pp->x.p->arr[0]) {
3153 int pos;
3154 assert(pp->x.p->type == polynomial);
3155 pos = pp->x.p->pos;
3156 if (pos >= 1 + nvar)
3157 ++pos;
3158 value_assign(row->p[pos], row->p[1+nvar]);
3159 value_division(row->p[pos], row->p[pos], pp->x.p->arr[1].d);
3160 value_multiply(row->p[pos], row->p[pos], pp->x.p->arr[1].x.n);
3161 }
3162 value_assign(row->p[1 + D->Dimension + 1], row->p[1+nvar]);
3163 value_division(row->p[1 + D->Dimension + 1],
3164 row->p[1 + D->Dimension + 1],
3165 pp->d);
3166 value_multiply(row->p[1 + D->Dimension + 1],
3167 row->p[1 + D->Dimension + 1],
3168 pp->x.n);
3169 value_oppose(row->p[1 + nvar], row->p[1 + nvar]);
3170 break;
3171 }
3172 case polynomial: {
3173 int pos = e->x.p->pos;
3174
3175 if (pos > nvar) {
3176 factor = ALLOC(evalue);
3177 value_init(factor->d);
3178 value_set_si(factor->d, 0);
3179 factor->x.p = new_enode(polynomial, 2, pos - nvar);
3180 evalue_set_si(&factor->x.p->arr[0], 0, 1);
3181 evalue_set_si(&factor->x.p->arr[1], 1, 1);
3182 break;
3183 }
3184
3185 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3186 for (i = 0; i < D->NbRays; ++i)
3187 if (value_notzero_p(D->Ray[i][pos]))
3188 break;
3189 assert(i < D->NbRays);
3190 if (value_neg_p(D->Ray[i][pos])) {
3191 factor = ALLOC(evalue);
3192 value_init(factor->d);
3193 evalue_set_si(factor, -1, 1);
3194 }
3195 value_set_si(row->p[0], 1);
3196 value_set_si(row->p[pos], 1);
3197 value_set_si(row->p[1 + nvar], -1);
3198 break;
3199 }
3200 default:
3201 assert(0);
3202 }
3203
3204 i = type_offset(e->x.p);
3205
3206 res = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3207 ++i;
3208
3209 if (factor) {
3210 value_init(cum.d);
3211 evalue_copy(&cum, factor);
3212 }
3213
3214 origC = C;
3215 for (; i < e->x.p->size; ++i) {
3216 evalue *t;
3217 if (row) {
3218 Matrix *prevC = C;
3219 C = esum_add_constraint(nvar, D, C, row);
3220 if (prevC != origC)
3221 Matrix_Free(prevC);
3222 }
3223 /*
3224 if (row)
3225 Vector_Print(stderr, P_VALUE_FMT, row);
3226 if (C)
3227 Matrix_Print(stderr, P_VALUE_FMT, C);
3228 */
3229 t = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3230
3231 if (t && factor)
3232 emul(&cum, t);
3233
3234 if (!res)
3235 res = t;
3236 else if (t) {
3237 eadd(t, res);
3238 free_evalue_refs(t);
3239 free(t);
3240 }
3241 if (factor && i+1 < e->x.p->size)
3242 emul(factor, &cum);
3243 }
3244 if (C != origC)
3245 Matrix_Free(C);
3246
3247 if (factor) {
3248 free_evalue_refs(factor);
3249 free_evalue_refs(&cum);
3250 free(factor);
3251 }
3252
3253 if (row)
3254 Vector_Free(row);
3255
3256 reduce_evalue(res);
3257
3258 return res;
3259 }
3260
3261 evalue *esum(evalue *e, int nvar)
3262 {
3263 int i;
3264 evalue *res = ALLOC(evalue);
3265 value_init(res->d);
3266
3267 assert(nvar >= 0);
3268 if (nvar == 0 || EVALUE_IS_ZERO(*e)) {
3269 evalue_copy(res, e);
3270 return res;
3271 }
3272
3273 evalue_set_si(res, 0, 1);
3274
3275 assert(value_zero_p(e->d));
3276 assert(e->x.p->type == partition);
3277
3278 for (i = 0; i < e->x.p->size/2; ++i) {
3279 evalue *t;
3280 t = esum_over_domain(&e->x.p->arr[2*i+1], nvar,
3281 EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
3282 eadd(t, res);
3283 free_evalue_refs(t);
3284 free(t);
3285 }
3286
3287 reduce_evalue(res);
3288
3289 return res;
3290 }
3291
3292 /* Initial silly implementation */
3293 void eor(evalue *e1, evalue *res)
3294 {
3295 evalue E;
3296 evalue mone;
3297 value_init(E.d);
3298 value_init(mone.d);
3299 evalue_set_si(&mone, -1, 1);
3300
3301 evalue_copy(&E, res);
3302 eadd(e1, &E);
3303 emul(e1, res);
3304 emul(&mone, res);
3305 eadd(&E, res);
3306
3307 free_evalue_refs(&E);
3308 free_evalue_refs(&mone);
3309 }
lattice point enumeration library