?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
barvinok.cc: remove ancient debugging code
[barvinok.git] / evalue.c
blobd1db85a41f8ed0c44e8309827110d7f13904bbec
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(evalue *e1, 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 void reduce_evalue (evalue *e) {
770 struct subst s = { NULL, 0, 0 };
771
772 if (value_notzero_p(e->d))
773 return; /* a rational number, its already reduced */
774
775 if (e->x.p->type == partition) {
776 int i;
777 unsigned dim = -1;
778 for (i = 0; i < e->x.p->size/2; ++i) {
779 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
780 s.n = 0;
781 /* This shouldn't really happen;
782 * Empty domains should not be added.
783 */
784 POL_ENSURE_VERTICES(D);
785 if (emptyQ(D))
786 goto discard;
787
788 dim = D->Dimension;
789 if (D->next)
790 D = DomainConvex(D, 0);
791 if (!D->next && D->NbEq) {
792 int j, k;
793 if (s.max < dim) {
794 if (s.max != 0)
795 realloc_substitution(&s, dim);
796 else {
797 int d = relations_depth(&e->x.p->arr[2*i+1]);
798 s.max = dim+d;
799 NALLOC(s.fixed, s.max);
800 }
801 }
802 for (j = 0; j < D->NbEq; ++j)
803 add_substitution(&s, D->Constraint[j], dim);
804 }
805 if (D != EVALUE_DOMAIN(e->x.p->arr[2*i]))
806 Domain_Free(D);
807 _reduce_evalue(&e->x.p->arr[2*i+1], &s, 0);
808 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
809 discard:
810 free_evalue_refs(&e->x.p->arr[2*i+1]);
811 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
812 value_clear(e->x.p->arr[2*i].d);
813 e->x.p->size -= 2;
814 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
815 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
816 --i;
817 }
818 if (s.n != 0) {
819 int j;
820 for (j = 0; j < s.n; ++j) {
821 value_clear(s.fixed[j].d);
822 value_clear(s.fixed[j].m);
823 free_evalue_refs(&s.fixed[j].s);
824 }
825 }
826 }
827 if (e->x.p->size == 0) {
828 free(e->x.p);
829 evalue_set_si(e, 0, 1);
830 }
831 } else
832 _reduce_evalue(e, &s, 0);
833 if (s.max != 0)
834 free(s.fixed);
835 }
836
837 void print_evalue(FILE *DST,evalue *e,char **pname) {
838
839 if(value_notzero_p(e->d)) {
840 if(value_notone_p(e->d)) {
841 value_print(DST,VALUE_FMT,e->x.n);
842 fprintf(DST,"/");
843 value_print(DST,VALUE_FMT,e->d);
844 }
845 else {
846 value_print(DST,VALUE_FMT,e->x.n);
847 }
848 }
849 else
850 print_enode(DST,e->x.p,pname);
851 return;
852 } /* print_evalue */
853
854 void print_enode(FILE *DST,enode *p,char **pname) {
855
856 int i;
857
858 if (!p) {
859 fprintf(DST, "NULL");
860 return;
861 }
862 switch (p->type) {
863 case evector:
864 fprintf(DST, "{ ");
865 for (i=0; i<p->size; i++) {
866 print_evalue(DST, &p->arr[i], pname);
867 if (i!=(p->size-1))
868 fprintf(DST, ", ");
869 }
870 fprintf(DST, " }\n");
871 break;
872 case polynomial:
873 fprintf(DST, "( ");
874 for (i=p->size-1; i>=0; i--) {
875 print_evalue(DST, &p->arr[i], pname);
876 if (i==1) fprintf(DST, " * %s + ", pname[p->pos-1]);
877 else if (i>1)
878 fprintf(DST, " * %s^%d + ", pname[p->pos-1], i);
879 }
880 fprintf(DST, " )\n");
881 break;
882 case periodic:
883 fprintf(DST, "[ ");
884 for (i=0; i<p->size; i++) {
885 print_evalue(DST, &p->arr[i], pname);
886 if (i!=(p->size-1)) fprintf(DST, ", ");
887 }
888 fprintf(DST," ]_%s", pname[p->pos-1]);
889 break;
890 case flooring:
891 case fractional:
892 fprintf(DST, "( ");
893 for (i=p->size-1; i>=1; i--) {
894 print_evalue(DST, &p->arr[i], pname);
895 if (i >= 2) {
896 fprintf(DST, " * ");
897 fprintf(DST, p->type == flooring ? "[" : "{");
898 print_evalue(DST, &p->arr[0], pname);
899 fprintf(DST, p->type == flooring ? "]" : "}");
900 if (i>2)
901 fprintf(DST, "^%d + ", i-1);
902 else
903 fprintf(DST, " + ");
904 }
905 }
906 fprintf(DST, " )\n");
907 break;
908 case relation:
909 fprintf(DST, "[ ");
910 print_evalue(DST, &p->arr[0], pname);
911 fprintf(DST, "= 0 ] * \n");
912 print_evalue(DST, &p->arr[1], pname);
913 if (p->size > 2) {
914 fprintf(DST, " +\n [ ");
915 print_evalue(DST, &p->arr[0], pname);
916 fprintf(DST, "!= 0 ] * \n");
917 print_evalue(DST, &p->arr[2], pname);
918 }
919 break;
920 case partition: {
921 char **names = pname;
922 int maxdim = EVALUE_DOMAIN(p->arr[0])->Dimension;
923 if (!pname || p->pos < maxdim) {
924 NALLOC(names, maxdim);
925 for (i = 0; i < p->pos; ++i) {
926 if (pname)
927 names[i] = pname[i];
928 else {
929 NALLOC(names[i], 10);
930 snprintf(names[i], 10, "%c", 'P'+i);
931 }
932 }
933 for ( ; i < maxdim; ++i) {
934 NALLOC(names[i], 10);
935 snprintf(names[i], 10, "_p%d", i);
936 }
937 }
938
939 for (i=0; i<p->size/2; i++) {
940 Print_Domain(DST, EVALUE_DOMAIN(p->arr[2*i]), names);
941 print_evalue(DST, &p->arr[2*i+1], names);
942 }
943
944 if (!pname || p->pos < maxdim) {
945 for (i = pname ? p->pos : 0; i < maxdim; ++i)
946 free(names[i]);
947 free(names);
948 }
949
950 break;
951 }
952 default:
953 assert(0);
954 }
955 return;
956 } /* print_enode */
957
958 static void eadd_rev(evalue *e1, evalue *res)
959 {
960 evalue ev;
961 value_init(ev.d);
962 evalue_copy(&ev, e1);
963 eadd(res, &ev);
964 free_evalue_refs(res);
965 *res = ev;
966 }
967
968 static void eadd_rev_cst (evalue *e1, evalue *res)
969 {
970 evalue ev;
971 value_init(ev.d);
972 evalue_copy(&ev, e1);
973 eadd(res, &ev.x.p->arr[type_offset(ev.x.p)]);
974 free_evalue_refs(res);
975 *res = ev;
976 }
977
978 static int is_zero_on(evalue *e, Polyhedron *D)
979 {
980 int is_zero;
981 evalue tmp;
982 value_init(tmp.d);
983 tmp.x.p = new_enode(partition, 2, D->Dimension);
984 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Domain_Copy(D));
985 evalue_copy(&tmp.x.p->arr[1], e);
986 reduce_evalue(&tmp);
987 is_zero = EVALUE_IS_ZERO(tmp);
988 free_evalue_refs(&tmp);
989 return is_zero;
990 }
991
992 struct section { Polyhedron * D; evalue E; };
993
994 void eadd_partitions (evalue *e1,evalue *res)
995 {
996 int n, i, j;
997 Polyhedron *d, *fd;
998 struct section *s;
999 s = (struct section *)
1000 malloc((e1->x.p->size/2+1) * (res->x.p->size/2+1) *
1001 sizeof(struct section));
1002 assert(s);
1003 assert(e1->x.p->pos == res->x.p->pos);
1004 assert(e1->x.p->pos == EVALUE_DOMAIN(e1->x.p->arr[0])->Dimension);
1005 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1006
1007 n = 0;
1008 for (j = 0; j < e1->x.p->size/2; ++j) {
1009 assert(res->x.p->size >= 2);
1010 fd = DomainDifference(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1011 EVALUE_DOMAIN(res->x.p->arr[0]), 0);
1012 if (!emptyQ(fd))
1013 for (i = 1; i < res->x.p->size/2; ++i) {
1014 Polyhedron *t = fd;
1015 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1016 Domain_Free(t);
1017 if (emptyQ(fd))
1018 break;
1019 }
1020 if (emptyQ(fd)) {
1021 Domain_Free(fd);
1022 continue;
1023 }
1024 /* See if we can extend one of the domains in res to cover fd */
1025 for (i = 0; i < res->x.p->size/2; ++i)
1026 if (is_zero_on(&res->x.p->arr[2*i+1], fd))
1027 break;
1028 if (i < res->x.p->size/2) {
1029 EVALUE_SET_DOMAIN(res->x.p->arr[2*i],
1030 DomainConcat(fd, EVALUE_DOMAIN(res->x.p->arr[2*i])));
1031 continue;
1032 }
1033 value_init(s[n].E.d);
1034 evalue_copy(&s[n].E, &e1->x.p->arr[2*j+1]);
1035 s[n].D = fd;
1036 ++n;
1037 }
1038 for (i = 0; i < res->x.p->size/2; ++i) {
1039 fd = EVALUE_DOMAIN(res->x.p->arr[2*i]);
1040 for (j = 0; j < e1->x.p->size/2; ++j) {
1041 Polyhedron *t;
1042 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1043 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1044 if (emptyQ(d)) {
1045 Domain_Free(d);
1046 continue;
1047 }
1048 t = fd;
1049 fd = DomainDifference(fd, EVALUE_DOMAIN(e1->x.p->arr[2*j]), 0);
1050 if (t != EVALUE_DOMAIN(res->x.p->arr[2*i]))
1051 Domain_Free(t);
1052 value_init(s[n].E.d);
1053 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
1054 eadd(&e1->x.p->arr[2*j+1], &s[n].E);
1055 if (!emptyQ(fd) && is_zero_on(&e1->x.p->arr[2*j+1], fd)) {
1056 d = DomainConcat(fd, d);
1057 fd = Empty_Polyhedron(fd->Dimension);
1058 }
1059 s[n].D = d;
1060 ++n;
1061 }
1062 if (!emptyQ(fd)) {
1063 s[n].E = res->x.p->arr[2*i+1];
1064 s[n].D = fd;
1065 ++n;
1066 } else {
1067 free_evalue_refs(&res->x.p->arr[2*i+1]);
1068 Domain_Free(fd);
1069 }
1070 if (fd != EVALUE_DOMAIN(res->x.p->arr[2*i]))
1071 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
1072 value_clear(res->x.p->arr[2*i].d);
1073 }
1074
1075 free(res->x.p);
1076 assert(n > 0);
1077 res->x.p = new_enode(partition, 2*n, e1->x.p->pos);
1078 for (j = 0; j < n; ++j) {
1079 s[j].D = DomainConstraintSimplify(s[j].D, 0);
1080 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1081 value_clear(res->x.p->arr[2*j+1].d);
1082 res->x.p->arr[2*j+1] = s[j].E;
1083 }
1084
1085 free(s);
1086 }
1087
1088 static void explicit_complement(evalue *res)
1089 {
1090 enode *rel = new_enode(relation, 3, 0);
1091 assert(rel);
1092 value_clear(rel->arr[0].d);
1093 rel->arr[0] = res->x.p->arr[0];
1094 value_clear(rel->arr[1].d);
1095 rel->arr[1] = res->x.p->arr[1];
1096 value_set_si(rel->arr[2].d, 1);
1097 value_init(rel->arr[2].x.n);
1098 value_set_si(rel->arr[2].x.n, 0);
1099 free(res->x.p);
1100 res->x.p = rel;
1101 }
1102
1103 void eadd(evalue *e1,evalue *res) {
1104
1105 int i;
1106 if (value_notzero_p(e1->d) && value_notzero_p(res->d)) {
1107 /* Add two rational numbers */
1108 Value g,m1,m2;
1109 value_init(g);
1110 value_init(m1);
1111 value_init(m2);
1112
1113 value_multiply(m1,e1->x.n,res->d);
1114 value_multiply(m2,res->x.n,e1->d);
1115 value_addto(res->x.n,m1,m2);
1116 value_multiply(res->d,e1->d,res->d);
1117 Gcd(res->x.n,res->d,&g);
1118 if (value_notone_p(g)) {
1119 value_division(res->d,res->d,g);
1120 value_division(res->x.n,res->x.n,g);
1121 }
1122 value_clear(g); value_clear(m1); value_clear(m2);
1123 return ;
1124 }
1125 else if (value_notzero_p(e1->d) && value_zero_p(res->d)) {
1126 switch (res->x.p->type) {
1127 case polynomial:
1128 /* Add the constant to the constant term of a polynomial*/
1129 eadd(e1, &res->x.p->arr[0]);
1130 return ;
1131 case periodic:
1132 /* Add the constant to all elements of a periodic number */
1133 for (i=0; i<res->x.p->size; i++) {
1134 eadd(e1, &res->x.p->arr[i]);
1135 }
1136 return ;
1137 case evector:
1138 fprintf(stderr, "eadd: cannot add const with vector\n");
1139 return;
1140 case flooring:
1141 case fractional:
1142 eadd(e1, &res->x.p->arr[1]);
1143 return ;
1144 case partition:
1145 assert(EVALUE_IS_ZERO(*e1));
1146 break; /* Do nothing */
1147 case relation:
1148 /* Create (zero) complement if needed */
1149 if (res->x.p->size < 3 && !EVALUE_IS_ZERO(*e1))
1150 explicit_complement(res);
1151 for (i = 1; i < res->x.p->size; ++i)
1152 eadd(e1, &res->x.p->arr[i]);
1153 break;
1154 default:
1155 assert(0);
1156 }
1157 }
1158 /* add polynomial or periodic to constant
1159 * you have to exchange e1 and res, before doing addition */
1160
1161 else if (value_zero_p(e1->d) && value_notzero_p(res->d)) {
1162 eadd_rev(e1, res);
1163 return;
1164 }
1165 else { // ((e1->d==0) && (res->d==0))
1166 assert(!((e1->x.p->type == partition) ^
1167 (res->x.p->type == partition)));
1168 if (e1->x.p->type == partition) {
1169 eadd_partitions(e1, res);
1170 return;
1171 }
1172 if (e1->x.p->type == relation &&
1173 (res->x.p->type != relation ||
1174 mod_term_smaller(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1175 eadd_rev(e1, res);
1176 return;
1177 }
1178 if (res->x.p->type == relation) {
1179 if (e1->x.p->type == relation &&
1180 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1181 if (res->x.p->size < 3 && e1->x.p->size == 3)
1182 explicit_complement(res);
1183 if (e1->x.p->size < 3 && res->x.p->size == 3)
1184 explicit_complement(e1);
1185 for (i = 1; i < res->x.p->size; ++i)
1186 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1187 return;
1188 }
1189 if (res->x.p->size < 3)
1190 explicit_complement(res);
1191 for (i = 1; i < res->x.p->size; ++i)
1192 eadd(e1, &res->x.p->arr[i]);
1193 return;
1194 }
1195 if ((e1->x.p->type != res->x.p->type) ) {
1196 /* adding to evalues of different type. two cases are possible
1197 * res is periodic and e1 is polynomial, you have to exchange
1198 * e1 and res then to add e1 to the constant term of res */
1199 if (e1->x.p->type == polynomial) {
1200 eadd_rev_cst(e1, res);
1201 }
1202 else if (res->x.p->type == polynomial) {
1203 /* res is polynomial and e1 is periodic,
1204 add e1 to the constant term of res */
1205
1206 eadd(e1,&res->x.p->arr[0]);
1207 } else
1208 assert(0);
1209
1210 return;
1211 }
1212 else if (e1->x.p->pos != res->x.p->pos ||
1213 ((res->x.p->type == fractional ||
1214 res->x.p->type == flooring) &&
1215 !eequal(&e1->x.p->arr[0], &res->x.p->arr[0]))) {
1216 /* adding evalues of different position (i.e function of different unknowns
1217 * to case are possible */
1218
1219 switch (res->x.p->type) {
1220 case flooring:
1221 case fractional:
1222 if(mod_term_smaller(res, e1))
1223 eadd(e1,&res->x.p->arr[1]);
1224 else
1225 eadd_rev_cst(e1, res);
1226 return;
1227 case polynomial: // res and e1 are polynomials
1228 // add e1 to the constant term of res
1229
1230 if(res->x.p->pos < e1->x.p->pos)
1231 eadd(e1,&res->x.p->arr[0]);
1232 else
1233 eadd_rev_cst(e1, res);
1234 // value_clear(g); value_clear(m1); value_clear(m2);
1235 return;
1236 case periodic: // res and e1 are pointers to periodic numbers
1237 //add e1 to all elements of res
1238
1239 if(res->x.p->pos < e1->x.p->pos)
1240 for (i=0;i<res->x.p->size;i++) {
1241 eadd(e1,&res->x.p->arr[i]);
1242 }
1243 else
1244 eadd_rev(e1, res);
1245 return;
1246 default:
1247 assert(0);
1248 }
1249 }
1250
1251
1252 //same type , same pos and same size
1253 if (e1->x.p->size == res->x.p->size) {
1254 // add any element in e1 to the corresponding element in res
1255 i = type_offset(res->x.p);
1256 if (i == 1)
1257 assert(eequal(&e1->x.p->arr[0], &res->x.p->arr[0]));
1258 for (; i<res->x.p->size; i++) {
1259 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1260 }
1261 return ;
1262 }
1263
1264 /* Sizes are different */
1265 switch(res->x.p->type) {
1266 case polynomial:
1267 case flooring:
1268 case fractional:
1269 /* VIN100: if e1-size > res-size you have to copy e1 in a */
1270 /* new enode and add res to that new node. If you do not do */
1271 /* that, you lose the the upper weight part of e1 ! */
1272
1273 if(e1->x.p->size > res->x.p->size)
1274 eadd_rev(e1, res);
1275 else {
1276 i = type_offset(res->x.p);
1277 if (i == 1)
1278 assert(eequal(&e1->x.p->arr[0],
1279 &res->x.p->arr[0]));
1280 for (; i<e1->x.p->size ; i++) {
1281 eadd(&e1->x.p->arr[i], &res->x.p->arr[i]);
1282 }
1283
1284 return ;
1285 }
1286 break;
1287
1288 /* add two periodics of the same pos (unknown) but whith different sizes (periods) */
1289 case periodic:
1290 {
1291 /* you have to create a new evalue 'ne' in whitch size equals to the lcm
1292 of the sizes of e1 and res, then to copy res periodicaly in ne, after
1293 to add periodicaly elements of e1 to elements of ne, and finaly to
1294 return ne. */
1295 int x,y,p;
1296 Value ex, ey ,ep;
1297 evalue *ne;
1298 value_init(ex); value_init(ey);value_init(ep);
1299 x=e1->x.p->size;
1300 y= res->x.p->size;
1301 value_set_si(ex,e1->x.p->size);
1302 value_set_si(ey,res->x.p->size);
1303 value_assign (ep,*Lcm(ex,ey));
1304 p=(int)mpz_get_si(ep);
1305 ne= (evalue *) malloc (sizeof(evalue));
1306 value_init(ne->d);
1307 value_set_si( ne->d,0);
1308
1309 ne->x.p=new_enode(res->x.p->type,p, res->x.p->pos);
1310 for(i=0;i<p;i++) {
1311 evalue_copy(&ne->x.p->arr[i], &res->x.p->arr[i%y]);
1312 }
1313 for(i=0;i<p;i++) {
1314 eadd(&e1->x.p->arr[i%x], &ne->x.p->arr[i]);
1315 }
1316
1317 value_assign(res->d, ne->d);
1318 res->x.p=ne->x.p;
1319
1320 return ;
1321 }
1322 case evector:
1323 fprintf(stderr, "eadd: ?cannot add vectors of different length\n");
1324 return ;
1325 default:
1326 assert(0);
1327 }
1328 }
1329 return ;
1330 }/* eadd */
1331
1332 static void emul_rev (evalue *e1, evalue *res)
1333 {
1334 evalue ev;
1335 value_init(ev.d);
1336 evalue_copy(&ev, e1);
1337 emul(res, &ev);
1338 free_evalue_refs(res);
1339 *res = ev;
1340 }
1341
1342 static void emul_poly (evalue *e1, evalue *res)
1343 {
1344 int i, j, o = type_offset(res->x.p);
1345 evalue tmp;
1346 int size=(e1->x.p->size + res->x.p->size - o - 1);
1347 value_init(tmp.d);
1348 value_set_si(tmp.d,0);
1349 tmp.x.p=new_enode(res->x.p->type, size, res->x.p->pos);
1350 if (o)
1351 evalue_copy(&tmp.x.p->arr[0], &e1->x.p->arr[0]);
1352 for (i=o; i < e1->x.p->size; i++) {
1353 evalue_copy(&tmp.x.p->arr[i], &e1->x.p->arr[i]);
1354 emul(&res->x.p->arr[o], &tmp.x.p->arr[i]);
1355 }
1356 for (; i<size; i++)
1357 evalue_set_si(&tmp.x.p->arr[i], 0, 1);
1358 for (i=o+1; i<res->x.p->size; i++)
1359 for (j=o; j<e1->x.p->size; j++) {
1360 evalue ev;
1361 value_init(ev.d);
1362 evalue_copy(&ev, &e1->x.p->arr[j]);
1363 emul(&res->x.p->arr[i], &ev);
1364 eadd(&ev, &tmp.x.p->arr[i+j-o]);
1365 free_evalue_refs(&ev);
1366 }
1367 free_evalue_refs(res);
1368 *res = tmp;
1369 }
1370
1371 void emul_partitions (evalue *e1,evalue *res)
1372 {
1373 int n, i, j, k;
1374 Polyhedron *d;
1375 struct section *s;
1376 s = (struct section *)
1377 malloc((e1->x.p->size/2) * (res->x.p->size/2) *
1378 sizeof(struct section));
1379 assert(s);
1380 assert(e1->x.p->pos == res->x.p->pos);
1381 assert(e1->x.p->pos == EVALUE_DOMAIN(e1->x.p->arr[0])->Dimension);
1382 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1383
1384 n = 0;
1385 for (i = 0; i < res->x.p->size/2; ++i) {
1386 for (j = 0; j < e1->x.p->size/2; ++j) {
1387 d = DomainIntersection(EVALUE_DOMAIN(e1->x.p->arr[2*j]),
1388 EVALUE_DOMAIN(res->x.p->arr[2*i]), 0);
1389 if (emptyQ(d)) {
1390 Domain_Free(d);
1391 continue;
1392 }
1393
1394 /* This code is only needed because the partitions
1395 are not true partitions.
1396 */
1397 for (k = 0; k < n; ++k) {
1398 if (DomainIncludes(s[k].D, d))
1399 break;
1400 if (DomainIncludes(d, s[k].D)) {
1401 Domain_Free(s[k].D);
1402 free_evalue_refs(&s[k].E);
1403 if (n > k)
1404 s[k] = s[--n];
1405 --k;
1406 }
1407 }
1408 if (k < n) {
1409 Domain_Free(d);
1410 continue;
1411 }
1412
1413 value_init(s[n].E.d);
1414 evalue_copy(&s[n].E, &res->x.p->arr[2*i+1]);
1415 emul(&e1->x.p->arr[2*j+1], &s[n].E);
1416 s[n].D = d;
1417 ++n;
1418 }
1419 Domain_Free(EVALUE_DOMAIN(res->x.p->arr[2*i]));
1420 value_clear(res->x.p->arr[2*i].d);
1421 free_evalue_refs(&res->x.p->arr[2*i+1]);
1422 }
1423
1424 free(res->x.p);
1425 if (n == 0)
1426 evalue_set_si(res, 0, 1);
1427 else {
1428 res->x.p = new_enode(partition, 2*n, e1->x.p->pos);
1429 for (j = 0; j < n; ++j) {
1430 s[j].D = DomainConstraintSimplify(s[j].D, 0);
1431 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1432 value_clear(res->x.p->arr[2*j+1].d);
1433 res->x.p->arr[2*j+1] = s[j].E;
1434 }
1435 }
1436
1437 free(s);
1438 }
1439
1440 #define value_two_p(val) (mpz_cmp_si(val,2) == 0)
1441
1442 /* Computes the product of two evalues "e1" and "res" and puts the result in "res". you must
1443 * do a copy of "res" befor calling this function if you nead it after. The vector type of
1444 * evalues is not treated here */
1445
1446 void emul (evalue *e1, evalue *res ){
1447 int i,j;
1448
1449 if((value_zero_p(e1->d)&&e1->x.p->type==evector)||(value_zero_p(res->d)&&(res->x.p->type==evector))) {
1450 fprintf(stderr, "emul: do not proced on evector type !\n");
1451 return;
1452 }
1453
1454 if (EVALUE_IS_ZERO(*res))
1455 return;
1456
1457 if (value_zero_p(e1->d) && e1->x.p->type == partition) {
1458 if (value_zero_p(res->d) && res->x.p->type == partition)
1459 emul_partitions(e1, res);
1460 else
1461 emul_rev(e1, res);
1462 } else if (value_zero_p(res->d) && res->x.p->type == partition) {
1463 for (i = 0; i < res->x.p->size/2; ++i)
1464 emul(e1, &res->x.p->arr[2*i+1]);
1465 } else
1466 if (value_zero_p(res->d) && res->x.p->type == relation) {
1467 if (value_zero_p(e1->d) && e1->x.p->type == relation &&
1468 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1469 if (res->x.p->size < 3 && e1->x.p->size == 3)
1470 explicit_complement(res);
1471 if (e1->x.p->size < 3 && res->x.p->size == 3)
1472 explicit_complement(e1);
1473 for (i = 1; i < res->x.p->size; ++i)
1474 emul(&e1->x.p->arr[i], &res->x.p->arr[i]);
1475 return;
1476 }
1477 for (i = 1; i < res->x.p->size; ++i)
1478 emul(e1, &res->x.p->arr[i]);
1479 } else
1480 if(value_zero_p(e1->d)&& value_zero_p(res->d)) {
1481 switch(e1->x.p->type) {
1482 case polynomial:
1483 switch(res->x.p->type) {
1484 case polynomial:
1485 if(e1->x.p->pos == res->x.p->pos) {
1486 /* Product of two polynomials of the same variable */
1487 emul_poly(e1, res);
1488 return;
1489 }
1490 else {
1491 /* Product of two polynomials of different variables */
1492
1493 if(res->x.p->pos < e1->x.p->pos)
1494 for( i=0; i<res->x.p->size ; i++)
1495 emul(e1, &res->x.p->arr[i]);
1496 else
1497 emul_rev(e1, res);
1498
1499 return ;
1500 }
1501 case periodic:
1502 case flooring:
1503 case fractional:
1504 /* Product of a polynomial and a periodic or fractional */
1505 emul_rev(e1, res);
1506 return;
1507 default:
1508 assert(0);
1509 }
1510 case periodic:
1511 switch(res->x.p->type) {
1512 case periodic:
1513 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size==res->x.p->size) {
1514 /* Product of two periodics of the same parameter and period */
1515
1516 for(i=0; i<res->x.p->size;i++)
1517 emul(&(e1->x.p->arr[i]), &(res->x.p->arr[i]));
1518
1519 return;
1520 }
1521 else{
1522 if(e1->x.p->pos==res->x.p->pos && e1->x.p->size!=res->x.p->size) {
1523 /* Product of two periodics of the same parameter and different periods */
1524 evalue *newp;
1525 Value x,y,z;
1526 int ix,iy,lcm;
1527 value_init(x); value_init(y);value_init(z);
1528 ix=e1->x.p->size;
1529 iy=res->x.p->size;
1530 value_set_si(x,e1->x.p->size);
1531 value_set_si(y,res->x.p->size);
1532 value_assign (z,*Lcm(x,y));
1533 lcm=(int)mpz_get_si(z);
1534 newp= (evalue *) malloc (sizeof(evalue));
1535 value_init(newp->d);
1536 value_set_si( newp->d,0);
1537 newp->x.p=new_enode(periodic,lcm, e1->x.p->pos);
1538 for(i=0;i<lcm;i++) {
1539 evalue_copy(&newp->x.p->arr[i],
1540 &res->x.p->arr[i%iy]);
1541 }
1542 for(i=0;i<lcm;i++)
1543 emul(&e1->x.p->arr[i%ix], &newp->x.p->arr[i]);
1544
1545 value_assign(res->d,newp->d);
1546 res->x.p=newp->x.p;
1547
1548 value_clear(x); value_clear(y);value_clear(z);
1549 return ;
1550 }
1551 else {
1552 /* Product of two periodics of different parameters */
1553
1554 if(res->x.p->pos < e1->x.p->pos)
1555 for(i=0; i<res->x.p->size; i++)
1556 emul(e1, &(res->x.p->arr[i]));
1557 else
1558 emul_rev(e1, res);
1559
1560 return;
1561 }
1562 }
1563 case polynomial:
1564 /* Product of a periodic and a polynomial */
1565
1566 for(i=0; i<res->x.p->size ; i++)
1567 emul(e1, &(res->x.p->arr[i]));
1568
1569 return;
1570
1571 }
1572 case flooring:
1573 case fractional:
1574 switch(res->x.p->type) {
1575 case polynomial:
1576 for(i=0; i<res->x.p->size ; i++)
1577 emul(e1, &(res->x.p->arr[i]));
1578 return;
1579 default:
1580 case periodic:
1581 assert(0);
1582 case flooring:
1583 case fractional:
1584 assert(e1->x.p->type == res->x.p->type);
1585 if (e1->x.p->pos == res->x.p->pos &&
1586 eequal(&e1->x.p->arr[0], &res->x.p->arr[0])) {
1587 evalue d;
1588 value_init(d.d);
1589 poly_denom(&e1->x.p->arr[0], &d.d);
1590 if (e1->x.p->type != fractional || !value_two_p(d.d))
1591 emul_poly(e1, res);
1592 else {
1593 evalue tmp;
1594 value_init(d.x.n);
1595 value_set_si(d.x.n, 1);
1596 /* { x }^2 == { x }/2 */
1597 /* a0 b0 + (a0 b1 + a1 b0 + a1 b1/2) { x } */
1598 assert(e1->x.p->size == 3);
1599 assert(res->x.p->size == 3);
1600 value_init(tmp.d);
1601 evalue_copy(&tmp, &res->x.p->arr[2]);
1602 emul(&d, &tmp);
1603 eadd(&res->x.p->arr[1], &tmp);
1604 emul(&e1->x.p->arr[2], &tmp);
1605 emul(&e1->x.p->arr[1], res);
1606 eadd(&tmp, &res->x.p->arr[2]);
1607 free_evalue_refs(&tmp);
1608 value_clear(d.x.n);
1609 }
1610 value_clear(d.d);
1611 } else {
1612 if(mod_term_smaller(res, e1))
1613 for(i=1; i<res->x.p->size ; i++)
1614 emul(e1, &(res->x.p->arr[i]));
1615 else
1616 emul_rev(e1, res);
1617 return;
1618 }
1619 }
1620 break;
1621 case relation:
1622 emul_rev(e1, res);
1623 break;
1624 default:
1625 assert(0);
1626 }
1627 }
1628 else {
1629 if (value_notzero_p(e1->d)&& value_notzero_p(res->d)) {
1630 /* Product of two rational numbers */
1631
1632 Value g;
1633 value_init(g);
1634 value_multiply(res->d,e1->d,res->d);
1635 value_multiply(res->x.n,e1->x.n,res->x.n );
1636 Gcd(res->x.n, res->d,&g);
1637 if (value_notone_p(g)) {
1638 value_division(res->d,res->d,g);
1639 value_division(res->x.n,res->x.n,g);
1640 }
1641 value_clear(g);
1642 return ;
1643 }
1644 else {
1645 if(value_zero_p(e1->d)&& value_notzero_p(res->d)) {
1646 /* Product of an expression (polynomial or peririodic) and a rational number */
1647
1648 emul_rev(e1, res);
1649 return ;
1650 }
1651 else {
1652 /* Product of a rationel number and an expression (polynomial or peririodic) */
1653
1654 i = type_offset(res->x.p);
1655 for (; i<res->x.p->size; i++)
1656 emul(e1, &res->x.p->arr[i]);
1657
1658 return ;
1659 }
1660 }
1661 }
1662
1663 return ;
1664 }
1665
1666 /* Frees mask content ! */
1667 void emask(evalue *mask, evalue *res) {
1668 int n, i, j;
1669 Polyhedron *d, *fd;
1670 struct section *s;
1671 evalue mone;
1672 int pos;
1673
1674 if (EVALUE_IS_ZERO(*res)) {
1675 free_evalue_refs(mask);
1676 return;
1677 }
1678
1679 assert(value_zero_p(mask->d));
1680 assert(mask->x.p->type == partition);
1681 assert(value_zero_p(res->d));
1682 assert(res->x.p->type == partition);
1683 assert(mask->x.p->pos == res->x.p->pos);
1684 assert(res->x.p->pos == EVALUE_DOMAIN(res->x.p->arr[0])->Dimension);
1685 assert(mask->x.p->pos == EVALUE_DOMAIN(mask->x.p->arr[0])->Dimension);
1686 pos = res->x.p->pos;
1687
1688 s = (struct section *)
1689 malloc((mask->x.p->size/2+1) * (res->x.p->size/2) *
1690 sizeof(struct section));
1691 assert(s);
1692
1693 value_init(mone.d);
1694 evalue_set_si(&mone, -1, 1);
1695
1696 n = 0;
1697 for (j = 0; j < res->x.p->size/2; ++j) {
1698 assert(mask->x.p->size >= 2);
1699 fd = DomainDifference(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1700 EVALUE_DOMAIN(mask->x.p->arr[0]), 0);
1701 if (!emptyQ(fd))
1702 for (i = 1; i < mask->x.p->size/2; ++i) {
1703 Polyhedron *t = fd;
1704 fd = DomainDifference(fd, EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1705 Domain_Free(t);
1706 if (emptyQ(fd))
1707 break;
1708 }
1709 if (emptyQ(fd)) {
1710 Domain_Free(fd);
1711 continue;
1712 }
1713 value_init(s[n].E.d);
1714 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1715 s[n].D = fd;
1716 ++n;
1717 }
1718 for (i = 0; i < mask->x.p->size/2; ++i) {
1719 if (EVALUE_IS_ONE(mask->x.p->arr[2*i+1]))
1720 continue;
1721
1722 fd = EVALUE_DOMAIN(mask->x.p->arr[2*i]);
1723 eadd(&mone, &mask->x.p->arr[2*i+1]);
1724 emul(&mone, &mask->x.p->arr[2*i+1]);
1725 for (j = 0; j < res->x.p->size/2; ++j) {
1726 Polyhedron *t;
1727 d = DomainIntersection(EVALUE_DOMAIN(res->x.p->arr[2*j]),
1728 EVALUE_DOMAIN(mask->x.p->arr[2*i]), 0);
1729 if (emptyQ(d)) {
1730 Domain_Free(d);
1731 continue;
1732 }
1733 t = fd;
1734 fd = DomainDifference(fd, EVALUE_DOMAIN(res->x.p->arr[2*j]), 0);
1735 if (t != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1736 Domain_Free(t);
1737 value_init(s[n].E.d);
1738 evalue_copy(&s[n].E, &res->x.p->arr[2*j+1]);
1739 emul(&mask->x.p->arr[2*i+1], &s[n].E);
1740 s[n].D = d;
1741 ++n;
1742 }
1743
1744 if (!emptyQ(fd)) {
1745 /* Just ignore; this may have been previously masked off */
1746 }
1747 if (fd != EVALUE_DOMAIN(mask->x.p->arr[2*i]))
1748 Domain_Free(fd);
1749 }
1750
1751 free_evalue_refs(&mone);
1752 free_evalue_refs(mask);
1753 free_evalue_refs(res);
1754 value_init(res->d);
1755 if (n == 0)
1756 evalue_set_si(res, 0, 1);
1757 else {
1758 res->x.p = new_enode(partition, 2*n, pos);
1759 for (j = 0; j < n; ++j) {
1760 EVALUE_SET_DOMAIN(res->x.p->arr[2*j], s[j].D);
1761 value_clear(res->x.p->arr[2*j+1].d);
1762 res->x.p->arr[2*j+1] = s[j].E;
1763 }
1764 }
1765
1766 free(s);
1767 }
1768
1769 void evalue_copy(evalue *dst, evalue *src)
1770 {
1771 value_assign(dst->d, src->d);
1772 if(value_notzero_p(src->d)) {
1773 value_init(dst->x.n);
1774 value_assign(dst->x.n, src->x.n);
1775 } else
1776 dst->x.p = ecopy(src->x.p);
1777 }
1778
1779 enode *new_enode(enode_type type,int size,int pos) {
1780
1781 enode *res;
1782 int i;
1783
1784 if(size == 0) {
1785 fprintf(stderr, "Allocating enode of size 0 !\n" );
1786 return NULL;
1787 }
1788 res = (enode *) malloc(sizeof(enode) + (size-1)*sizeof(evalue));
1789 res->type = type;
1790 res->size = size;
1791 res->pos = pos;
1792 for(i=0; i<size; i++) {
1793 value_init(res->arr[i].d);
1794 value_set_si(res->arr[i].d,0);
1795 res->arr[i].x.p = 0;
1796 }
1797 return res;
1798 } /* new_enode */
1799
1800 enode *ecopy(enode *e) {
1801
1802 enode *res;
1803 int i;
1804
1805 res = new_enode(e->type,e->size,e->pos);
1806 for(i=0;i<e->size;++i) {
1807 value_assign(res->arr[i].d,e->arr[i].d);
1808 if(value_zero_p(res->arr[i].d))
1809 res->arr[i].x.p = ecopy(e->arr[i].x.p);
1810 else if (EVALUE_IS_DOMAIN(res->arr[i]))
1811 EVALUE_SET_DOMAIN(res->arr[i], Domain_Copy(EVALUE_DOMAIN(e->arr[i])));
1812 else {
1813 value_init(res->arr[i].x.n);
1814 value_assign(res->arr[i].x.n,e->arr[i].x.n);
1815 }
1816 }
1817 return(res);
1818 } /* ecopy */
1819
1820 int ecmp(const evalue *e1, const evalue *e2)
1821 {
1822 enode *p1, *p2;
1823 int i;
1824 int r;
1825
1826 if (value_notzero_p(e1->d) && value_notzero_p(e2->d)) {
1827 Value m, m2;
1828 value_init(m);
1829 value_init(m2);
1830 value_multiply(m, e1->x.n, e2->d);
1831 value_multiply(m2, e2->x.n, e1->d);
1832
1833 if (value_lt(m, m2))
1834 r = -1;
1835 else if (value_gt(m, m2))
1836 r = 1;
1837 else
1838 r = 0;
1839
1840 value_clear(m);
1841 value_clear(m2);
1842
1843 return r;
1844 }
1845 if (value_notzero_p(e1->d))
1846 return -1;
1847 if (value_notzero_p(e2->d))
1848 return 1;
1849
1850 p1 = e1->x.p;
1851 p2 = e2->x.p;
1852
1853 if (p1->type != p2->type)
1854 return p1->type - p2->type;
1855 if (p1->pos != p2->pos)
1856 return p1->pos - p2->pos;
1857 if (p1->size != p2->size)
1858 return p1->size - p2->size;
1859
1860 for (i = p1->size-1; i >= 0; --i)
1861 if ((r = ecmp(&p1->arr[i], &p2->arr[i])) != 0)
1862 return r;
1863
1864 return 0;
1865 }
1866
1867 int eequal(evalue *e1,evalue *e2) {
1868
1869 int i;
1870 enode *p1, *p2;
1871
1872 if (value_ne(e1->d,e2->d))
1873 return 0;
1874
1875 /* e1->d == e2->d */
1876 if (value_notzero_p(e1->d)) {
1877 if (value_ne(e1->x.n,e2->x.n))
1878 return 0;
1879
1880 /* e1->d == e2->d != 0 AND e1->n == e2->n */
1881 return 1;
1882 }
1883
1884 /* e1->d == e2->d == 0 */
1885 p1 = e1->x.p;
1886 p2 = e2->x.p;
1887 if (p1->type != p2->type) return 0;
1888 if (p1->size != p2->size) return 0;
1889 if (p1->pos != p2->pos) return 0;
1890 for (i=0; i<p1->size; i++)
1891 if (!eequal(&p1->arr[i], &p2->arr[i]) )
1892 return 0;
1893 return 1;
1894 } /* eequal */
1895
1896 void free_evalue_refs(evalue *e) {
1897
1898 enode *p;
1899 int i;
1900
1901 if (EVALUE_IS_DOMAIN(*e)) {
1902 Domain_Free(EVALUE_DOMAIN(*e));
1903 value_clear(e->d);
1904 return;
1905 } else if (value_pos_p(e->d)) {
1906
1907 /* 'e' stores a constant */
1908 value_clear(e->d);
1909 value_clear(e->x.n);
1910 return;
1911 }
1912 assert(value_zero_p(e->d));
1913 value_clear(e->d);
1914 p = e->x.p;
1915 if (!p) return; /* null pointer */
1916 for (i=0; i<p->size; i++) {
1917 free_evalue_refs(&(p->arr[i]));
1918 }
1919 free(p);
1920 return;
1921 } /* free_evalue_refs */
1922
1923 static void mod2table_r(evalue *e, Vector *periods, Value m, int p,
1924 Vector * val, evalue *res)
1925 {
1926 unsigned nparam = periods->Size;
1927
1928 if (p == nparam) {
1929 double d = compute_evalue(e, val->p);
1930 d *= VALUE_TO_DOUBLE(m);
1931 if (d > 0)
1932 d += .25;
1933 else
1934 d -= .25;
1935 value_assign(res->d, m);
1936 value_init(res->x.n);
1937 value_set_double(res->x.n, d);
1938 mpz_fdiv_r(res->x.n, res->x.n, m);
1939 return;
1940 }
1941 if (value_one_p(periods->p[p]))
1942 mod2table_r(e, periods, m, p+1, val, res);
1943 else {
1944 Value tmp;
1945 value_init(tmp);
1946
1947 value_assign(tmp, periods->p[p]);
1948 value_set_si(res->d, 0);
1949 res->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
1950 do {
1951 value_decrement(tmp, tmp);
1952 value_assign(val->p[p], tmp);
1953 mod2table_r(e, periods, m, p+1, val,
1954 &res->x.p->arr[VALUE_TO_INT(tmp)]);
1955 } while (value_pos_p(tmp));
1956
1957 value_clear(tmp);
1958 }
1959 }
1960
1961 static void rel2table(evalue *e, int zero)
1962 {
1963 if (value_pos_p(e->d)) {
1964 if (value_zero_p(e->x.n) == zero)
1965 value_set_si(e->x.n, 1);
1966 else
1967 value_set_si(e->x.n, 0);
1968 value_set_si(e->d, 1);
1969 } else {
1970 int i;
1971 for (i = 0; i < e->x.p->size; ++i)
1972 rel2table(&e->x.p->arr[i], zero);
1973 }
1974 }
1975
1976 void evalue_mod2table(evalue *e, int nparam)
1977 {
1978 enode *p;
1979 int i;
1980
1981 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
1982 return;
1983 p = e->x.p;
1984 for (i=0; i<p->size; i++) {
1985 evalue_mod2table(&(p->arr[i]), nparam);
1986 }
1987 if (p->type == relation) {
1988 evalue copy;
1989
1990 if (p->size > 2) {
1991 value_init(copy.d);
1992 evalue_copy(&copy, &p->arr[0]);
1993 }
1994 rel2table(&p->arr[0], 1);
1995 emul(&p->arr[0], &p->arr[1]);
1996 if (p->size > 2) {
1997 rel2table(&copy, 0);
1998 emul(&copy, &p->arr[2]);
1999 eadd(&p->arr[2], &p->arr[1]);
2000 free_evalue_refs(&p->arr[2]);
2001 free_evalue_refs(&copy);
2002 }
2003 free_evalue_refs(&p->arr[0]);
2004 value_clear(e->d);
2005 *e = p->arr[1];
2006 free(p);
2007 } else if (p->type == fractional) {
2008 Vector *periods = Vector_Alloc(nparam);
2009 Vector *val = Vector_Alloc(nparam);
2010 Value tmp;
2011 evalue *ev;
2012 evalue EP, res;
2013
2014 value_init(tmp);
2015 value_set_si(tmp, 1);
2016 Vector_Set(periods->p, 1, nparam);
2017 Vector_Set(val->p, 0, nparam);
2018 for (ev = &p->arr[0]; value_zero_p(ev->d); ev = &ev->x.p->arr[0]) {
2019 enode *p = ev->x.p;
2020
2021 assert(p->type == polynomial);
2022 assert(p->size == 2);
2023 value_assign(periods->p[p->pos-1], p->arr[1].d);
2024 value_lcm(tmp, p->arr[1].d, &tmp);
2025 }
2026 value_lcm(tmp, ev->d, &tmp);
2027 value_init(EP.d);
2028 mod2table_r(&p->arr[0], periods, tmp, 0, val, &EP);
2029
2030 value_init(res.d);
2031 evalue_set_si(&res, 0, 1);
2032 /* Compute the polynomial using Horner's rule */
2033 for (i=p->size-1;i>1;i--) {
2034 eadd(&p->arr[i], &res);
2035 emul(&EP, &res);
2036 }
2037 eadd(&p->arr[1], &res);
2038
2039 free_evalue_refs(e);
2040 free_evalue_refs(&EP);
2041 *e = res;
2042
2043 value_clear(tmp);
2044 Vector_Free(val);
2045 Vector_Free(periods);
2046 }
2047 } /* evalue_mod2table */
2048
2049 /********************************************************/
2050 /* function in domain */
2051 /* check if the parameters in list_args */
2052 /* verifies the constraints of Domain P */
2053 /********************************************************/
2054 int in_domain(Polyhedron *P, Value *list_args) {
2055
2056 int col,row;
2057 Value v; /* value of the constraint of a row when
2058 parameters are instanciated*/
2059 Value tmp;
2060
2061 value_init(v);
2062 value_init(tmp);
2063
2064 /*P->Constraint constraint matrice of polyhedron P*/
2065 for(row=0;row<P->NbConstraints;row++) {
2066 value_assign(v,P->Constraint[row][P->Dimension+1]); /*constant part*/
2067 for(col=1;col<P->Dimension+1;col++) {
2068 value_multiply(tmp,P->Constraint[row][col],list_args[col-1]);
2069 value_addto(v,v,tmp);
2070 }
2071 if (value_notzero_p(P->Constraint[row][0])) {
2072
2073 /*if v is not >=0 then this constraint is not respected */
2074 if (value_neg_p(v)) {
2075 next:
2076 value_clear(v);
2077 value_clear(tmp);
2078 return P->next ? in_domain(P->next, list_args) : 0;
2079 }
2080 }
2081 else {
2082
2083 /*if v is not = 0 then this constraint is not respected */
2084 if (value_notzero_p(v))
2085 goto next;
2086 }
2087 }
2088
2089 /*if not return before this point => all
2090 the constraints are respected */
2091 value_clear(v);
2092 value_clear(tmp);
2093 return 1;
2094 } /* in_domain */
2095
2096 /****************************************************/
2097 /* function compute enode */
2098 /* compute the value of enode p with parameters */
2099 /* list "list_args */
2100 /* compute the polynomial or the periodic */
2101 /****************************************************/
2102
2103 static double compute_enode(enode *p, Value *list_args) {
2104
2105 int i;
2106 Value m, param;
2107 double res=0.0;
2108
2109 if (!p)
2110 return(0.);
2111
2112 value_init(m);
2113 value_init(param);
2114
2115 if (p->type == polynomial) {
2116 if (p->size > 1)
2117 value_assign(param,list_args[p->pos-1]);
2118
2119 /* Compute the polynomial using Horner's rule */
2120 for (i=p->size-1;i>0;i--) {
2121 res +=compute_evalue(&p->arr[i],list_args);
2122 res *=VALUE_TO_DOUBLE(param);
2123 }
2124 res +=compute_evalue(&p->arr[0],list_args);
2125 }
2126 else if (p->type == fractional) {
2127 double d = compute_evalue(&p->arr[0], list_args);
2128 d -= floor(d+1e-10);
2129
2130 /* Compute the polynomial using Horner's rule */
2131 for (i=p->size-1;i>1;i--) {
2132 res +=compute_evalue(&p->arr[i],list_args);
2133 res *=d;
2134 }
2135 res +=compute_evalue(&p->arr[1],list_args);
2136 }
2137 else if (p->type == flooring) {
2138 double d = compute_evalue(&p->arr[0], list_args);
2139 d = floor(d+1e-10);
2140
2141 /* Compute the polynomial using Horner's rule */
2142 for (i=p->size-1;i>1;i--) {
2143 res +=compute_evalue(&p->arr[i],list_args);
2144 res *=d;
2145 }
2146 res +=compute_evalue(&p->arr[1],list_args);
2147 }
2148 else if (p->type == periodic) {
2149 value_assign(m,list_args[p->pos-1]);
2150
2151 /* Choose the right element of the periodic */
2152 value_set_si(param,p->size);
2153 value_pmodulus(m,m,param);
2154 res = compute_evalue(&p->arr[VALUE_TO_INT(m)],list_args);
2155 }
2156 else if (p->type == relation) {
2157 if (fabs(compute_evalue(&p->arr[0], list_args)) < 1e-10)
2158 res = compute_evalue(&p->arr[1], list_args);
2159 else if (p->size > 2)
2160 res = compute_evalue(&p->arr[2], list_args);
2161 }
2162 else if (p->type == partition) {
2163 int dim = EVALUE_DOMAIN(p->arr[0])->Dimension;
2164 Value *vals = list_args;
2165 if (p->pos < dim) {
2166 NALLOC(vals, dim);
2167 for (i = 0; i < dim; ++i) {
2168 value_init(vals[i]);
2169 if (i < p->pos)
2170 value_assign(vals[i], list_args[i]);
2171 }
2172 }
2173 for (i = 0; i < p->size/2; ++i)
2174 if (DomainContains(EVALUE_DOMAIN(p->arr[2*i]), vals, p->pos, 0, 1)) {
2175 res = compute_evalue(&p->arr[2*i+1], vals);
2176 break;
2177 }
2178 if (p->pos < dim) {
2179 for (i = 0; i < dim; ++i)
2180 value_clear(vals[i]);
2181 free(vals);
2182 }
2183 }
2184 else
2185 assert(0);
2186 value_clear(m);
2187 value_clear(param);
2188 return res;
2189 } /* compute_enode */
2190
2191 /*************************************************/
2192 /* return the value of Ehrhart Polynomial */
2193 /* It returns a double, because since it is */
2194 /* a recursive function, some intermediate value */
2195 /* might not be integral */
2196 /*************************************************/
2197
2198 double compute_evalue(evalue *e,Value *list_args) {
2199
2200 double res;
2201
2202 if (value_notzero_p(e->d)) {
2203 if (value_notone_p(e->d))
2204 res = VALUE_TO_DOUBLE(e->x.n) / VALUE_TO_DOUBLE(e->d);
2205 else
2206 res = VALUE_TO_DOUBLE(e->x.n);
2207 }
2208 else
2209 res = compute_enode(e->x.p,list_args);
2210 return res;
2211 } /* compute_evalue */
2212
2213
2214 /****************************************************/
2215 /* function compute_poly : */
2216 /* Check for the good validity domain */
2217 /* return the number of point in the Polyhedron */
2218 /* in allocated memory */
2219 /* Using the Ehrhart pseudo-polynomial */
2220 /****************************************************/
2221 Value *compute_poly(Enumeration *en,Value *list_args) {
2222
2223 Value *tmp;
2224 /* double d; int i; */
2225
2226 tmp = (Value *) malloc (sizeof(Value));
2227 assert(tmp != NULL);
2228 value_init(*tmp);
2229 value_set_si(*tmp,0);
2230
2231 if(!en)
2232 return(tmp); /* no ehrhart polynomial */
2233 if(en->ValidityDomain) {
2234 if(!en->ValidityDomain->Dimension) { /* no parameters */
2235 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2236 return(tmp);
2237 }
2238 }
2239 else
2240 return(tmp); /* no Validity Domain */
2241 while(en) {
2242 if(in_domain(en->ValidityDomain,list_args)) {
2243
2244 #ifdef EVAL_EHRHART_DEBUG
2245 Print_Domain(stdout,en->ValidityDomain);
2246 print_evalue(stdout,&en->EP);
2247 #endif
2248
2249 /* d = compute_evalue(&en->EP,list_args);
2250 i = d;
2251 printf("(double)%lf = %d\n", d, i ); */
2252 value_set_double(*tmp,compute_evalue(&en->EP,list_args)+.25);
2253 return(tmp);
2254 }
2255 else
2256 en=en->next;
2257 }
2258 value_set_si(*tmp,0);
2259 return(tmp); /* no compatible domain with the arguments */
2260 } /* compute_poly */
2261
2262 size_t value_size(Value v) {
2263 return (v[0]._mp_size > 0 ? v[0]._mp_size : -v[0]._mp_size)
2264 * sizeof(v[0]._mp_d[0]);
2265 }
2266
2267 size_t domain_size(Polyhedron *D)
2268 {
2269 int i, j;
2270 size_t s = sizeof(*D);
2271
2272 for (i = 0; i < D->NbConstraints; ++i)
2273 for (j = 0; j < D->Dimension+2; ++j)
2274 s += value_size(D->Constraint[i][j]);
2275
2276 /*
2277 for (i = 0; i < D->NbRays; ++i)
2278 for (j = 0; j < D->Dimension+2; ++j)
2279 s += value_size(D->Ray[i][j]);
2280 */
2281
2282 return D->next ? s+domain_size(D->next) : s;
2283 }
2284
2285 size_t enode_size(enode *p) {
2286 size_t s = sizeof(*p) - sizeof(p->arr[0]);
2287 int i;
2288
2289 if (p->type == partition)
2290 for (i = 0; i < p->size/2; ++i) {
2291 s += domain_size(EVALUE_DOMAIN(p->arr[2*i]));
2292 s += evalue_size(&p->arr[2*i+1]);
2293 }
2294 else
2295 for (i = 0; i < p->size; ++i) {
2296 s += evalue_size(&p->arr[i]);
2297 }
2298 return s;
2299 }
2300
2301 size_t evalue_size(evalue *e)
2302 {
2303 size_t s = sizeof(*e);
2304 s += value_size(e->d);
2305 if (value_notzero_p(e->d))
2306 s += value_size(e->x.n);
2307 else
2308 s += enode_size(e->x.p);
2309 return s;
2310 }
2311
2312 static evalue *find_second(evalue *base, evalue *cst, evalue *e, Value m)
2313 {
2314 evalue *found = NULL;
2315 evalue offset;
2316 evalue copy;
2317 int i;
2318
2319 if (value_pos_p(e->d) || e->x.p->type != fractional)
2320 return NULL;
2321
2322 value_init(offset.d);
2323 value_init(offset.x.n);
2324 poly_denom(&e->x.p->arr[0], &offset.d);
2325 value_lcm(m, offset.d, &offset.d);
2326 value_set_si(offset.x.n, 1);
2327
2328 value_init(copy.d);
2329 evalue_copy(&copy, cst);
2330
2331 eadd(&offset, cst);
2332 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2333
2334 if (eequal(base, &e->x.p->arr[0]))
2335 found = &e->x.p->arr[0];
2336 else {
2337 value_set_si(offset.x.n, -2);
2338
2339 eadd(&offset, cst);
2340 mpz_fdiv_r(cst->x.n, cst->x.n, cst->d);
2341
2342 if (eequal(base, &e->x.p->arr[0]))
2343 found = base;
2344 }
2345 free_evalue_refs(cst);
2346 free_evalue_refs(&offset);
2347 *cst = copy;
2348
2349 for (i = 1; !found && i < e->x.p->size; ++i)
2350 found = find_second(base, cst, &e->x.p->arr[i], m);
2351
2352 return found;
2353 }
2354
2355 static evalue *find_relation_pair(evalue *e)
2356 {
2357 int i;
2358 evalue *found = NULL;
2359
2360 if (EVALUE_IS_DOMAIN(*e) || value_pos_p(e->d))
2361 return NULL;
2362
2363 if (e->x.p->type == fractional) {
2364 Value m;
2365 evalue *cst;
2366
2367 value_init(m);
2368 poly_denom(&e->x.p->arr[0], &m);
2369
2370 for (cst = &e->x.p->arr[0]; value_zero_p(cst->d);
2371 cst = &cst->x.p->arr[0])
2372 ;
2373
2374 for (i = 1; !found && i < e->x.p->size; ++i)
2375 found = find_second(&e->x.p->arr[0], cst, &e->x.p->arr[i], m);
2376
2377 value_clear(m);
2378 }
2379
2380 i = e->x.p->type == relation;
2381 for (; !found && i < e->x.p->size; ++i)
2382 found = find_relation_pair(&e->x.p->arr[i]);
2383
2384 return found;
2385 }
2386
2387 void evalue_mod2relation(evalue *e) {
2388 evalue *d;
2389
2390 if (value_zero_p(e->d) && e->x.p->type == partition) {
2391 int i;
2392
2393 for (i = 0; i < e->x.p->size/2; ++i) {
2394 evalue_mod2relation(&e->x.p->arr[2*i+1]);
2395 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2396 value_clear(e->x.p->arr[2*i].d);
2397 free_evalue_refs(&e->x.p->arr[2*i+1]);
2398 e->x.p->size -= 2;
2399 if (2*i < e->x.p->size) {
2400 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2401 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2402 }
2403 --i;
2404 }
2405 }
2406 if (e->x.p->size == 0) {
2407 free(e->x.p);
2408 evalue_set_si(e, 0, 1);
2409 }
2410
2411 return;
2412 }
2413
2414 while ((d = find_relation_pair(e)) != NULL) {
2415 evalue split;
2416 evalue *ev;
2417
2418 value_init(split.d);
2419 value_set_si(split.d, 0);
2420 split.x.p = new_enode(relation, 3, 0);
2421 evalue_set_si(&split.x.p->arr[1], 1, 1);
2422 evalue_set_si(&split.x.p->arr[2], 1, 1);
2423
2424 ev = &split.x.p->arr[0];
2425 value_set_si(ev->d, 0);
2426 ev->x.p = new_enode(fractional, 3, -1);
2427 evalue_set_si(&ev->x.p->arr[1], 0, 1);
2428 evalue_set_si(&ev->x.p->arr[2], 1, 1);
2429 evalue_copy(&ev->x.p->arr[0], d);
2430
2431 emul(&split, e);
2432
2433 reduce_evalue(e);
2434
2435 free_evalue_refs(&split);
2436 }
2437 }
2438
2439 static int evalue_comp(const void * a, const void * b)
2440 {
2441 const evalue *e1 = *(const evalue **)a;
2442 const evalue *e2 = *(const evalue **)b;
2443 return ecmp(e1, e2);
2444 }
2445
2446 void evalue_combine(evalue *e)
2447 {
2448 evalue **evs;
2449 int i, k;
2450 enode *p;
2451 evalue tmp;
2452
2453 if (value_notzero_p(e->d) || e->x.p->type != partition)
2454 return;
2455
2456 NALLOC(evs, e->x.p->size/2);
2457 for (i = 0; i < e->x.p->size/2; ++i)
2458 evs[i] = &e->x.p->arr[2*i+1];
2459 qsort(evs, e->x.p->size/2, sizeof(evs[0]), evalue_comp);
2460 p = new_enode(partition, e->x.p->size, e->x.p->pos);
2461 for (i = 0, k = 0; i < e->x.p->size/2; ++i) {
2462 if (k == 0 || ecmp(&p->arr[2*k-1], evs[i]) != 0) {
2463 value_clear(p->arr[2*k].d);
2464 value_clear(p->arr[2*k+1].d);
2465 p->arr[2*k] = *(evs[i]-1);
2466 p->arr[2*k+1] = *(evs[i]);
2467 ++k;
2468 } else {
2469 Polyhedron *D = EVALUE_DOMAIN(*(evs[i]-1));
2470 Polyhedron *L = D;
2471
2472 value_clear((evs[i]-1)->d);
2473
2474 while (L->next)
2475 L = L->next;
2476 L->next = EVALUE_DOMAIN(p->arr[2*k-2]);
2477 EVALUE_SET_DOMAIN(p->arr[2*k-2], D);
2478 free_evalue_refs(evs[i]);
2479 }
2480 }
2481
2482 for (i = 2*k ; i < p->size; ++i)
2483 value_clear(p->arr[i].d);
2484
2485 free(evs);
2486 free(e->x.p);
2487 p->size = 2*k;
2488 e->x.p = p;
2489
2490 for (i = 0; i < e->x.p->size/2; ++i) {
2491 Polyhedron *H;
2492 if (value_notzero_p(e->x.p->arr[2*i+1].d))
2493 continue;
2494 H = DomainConvex(EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
2495 if (H == NULL)
2496 continue;
2497 for (k = 0; k < e->x.p->size/2; ++k) {
2498 Polyhedron *D, *N, **P;
2499 if (i == k)
2500 continue;
2501 P = &EVALUE_DOMAIN(e->x.p->arr[2*k]);
2502 D = *P;
2503 if (D == NULL)
2504 continue;
2505 for (; D; D = N) {
2506 *P = D;
2507 N = D->next;
2508 if (D->NbEq <= H->NbEq) {
2509 P = &D->next;
2510 continue;
2511 }
2512
2513 value_init(tmp.d);
2514 tmp.x.p = new_enode(partition, 2, e->x.p->pos);
2515 EVALUE_SET_DOMAIN(tmp.x.p->arr[0], Polyhedron_Copy(D));
2516 evalue_copy(&tmp.x.p->arr[1], &e->x.p->arr[2*i+1]);
2517 reduce_evalue(&tmp);
2518 if (value_notzero_p(tmp.d) ||
2519 ecmp(&tmp.x.p->arr[1], &e->x.p->arr[2*k+1]) != 0)
2520 P = &D->next;
2521 else {
2522 D->next = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2523 EVALUE_DOMAIN(e->x.p->arr[2*i]) = D;
2524 *P = NULL;
2525 }
2526 free_evalue_refs(&tmp);
2527 }
2528 }
2529 Polyhedron_Free(H);
2530 }
2531
2532 for (i = 0; i < e->x.p->size/2; ++i) {
2533 Polyhedron *H, *E;
2534 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
2535 if (!D) {
2536 value_clear(e->x.p->arr[2*i].d);
2537 free_evalue_refs(&e->x.p->arr[2*i+1]);
2538 e->x.p->size -= 2;
2539 if (2*i < e->x.p->size) {
2540 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2541 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2542 }
2543 --i;
2544 continue;
2545 }
2546 if (!D->next)
2547 continue;
2548 H = DomainConvex(D, 0);
2549 E = DomainDifference(H, D, 0);
2550 Domain_Free(D);
2551 D = DomainDifference(H, E, 0);
2552 Domain_Free(H);
2553 Domain_Free(E);
2554 EVALUE_SET_DOMAIN(p->arr[2*i], D);
2555 }
2556 }
2557
2558 /* May change coefficients to become non-standard if fiddle is set
2559 * => reduce p afterwards to correct
2560 */
2561 static Polyhedron *polynomial_projection(enode *p, Polyhedron *D, Value *d,
2562 Matrix **R, int fiddle)
2563 {
2564 Polyhedron *I, *H;
2565 evalue *pp;
2566 unsigned dim = D->Dimension;
2567 Matrix *T = Matrix_Alloc(2, dim+1);
2568 Value twice;
2569 value_init(twice);
2570 assert(T);
2571
2572 assert(p->type == fractional);
2573 pp = &p->arr[0];
2574 value_set_si(T->p[1][dim], 1);
2575 poly_denom(pp, d);
2576 while (value_zero_p(pp->d)) {
2577 assert(pp->x.p->type == polynomial);
2578 assert(pp->x.p->size == 2);
2579 assert(value_notzero_p(pp->x.p->arr[1].d));
2580 value_division(T->p[0][pp->x.p->pos-1], *d, pp->x.p->arr[1].d);
2581 mpz_mul_ui(twice, pp->x.p->arr[1].x.n, 2);
2582 if (fiddle && value_gt(twice, pp->x.p->arr[1].d))
2583 value_subtract(pp->x.p->arr[1].x.n,
2584 pp->x.p->arr[1].x.n, pp->x.p->arr[1].d);
2585 value_multiply(T->p[0][pp->x.p->pos-1],
2586 T->p[0][pp->x.p->pos-1], pp->x.p->arr[1].x.n);
2587 pp = &pp->x.p->arr[0];
2588 }
2589 value_division(T->p[0][dim], *d, pp->d);
2590 value_multiply(T->p[0][dim], T->p[0][dim], pp->x.n);
2591 I = DomainImage(D, T, 0);
2592 H = DomainConvex(I, 0);
2593 Domain_Free(I);
2594 *R = T;
2595
2596 value_clear(twice);
2597
2598 return H;
2599 }
2600
2601 int reduce_in_domain(evalue *e, Polyhedron *D)
2602 {
2603 int i;
2604 enode *p;
2605 Value d, min, max;
2606 int r = 0;
2607 Polyhedron *I;
2608 Matrix *T;
2609 int bounded;
2610
2611 if (value_notzero_p(e->d))
2612 return r;
2613
2614 p = e->x.p;
2615
2616 if (p->type == relation) {
2617 int equal;
2618 value_init(d);
2619 value_init(min);
2620 value_init(max);
2621
2622 I = polynomial_projection(p->arr[0].x.p, D, &d, &T, 1);
2623 bounded = line_minmax(I, &min, &max); /* frees I */
2624 equal = value_eq(min, max);
2625 mpz_cdiv_q(min, min, d);
2626 mpz_fdiv_q(max, max, d);
2627
2628 if (bounded && value_gt(min, max)) {
2629 /* Never zero */
2630 if (p->size == 3) {
2631 value_clear(e->d);
2632 *e = p->arr[2];
2633 } else {
2634 evalue_set_si(e, 0, 1);
2635 r = 1;
2636 }
2637 free_evalue_refs(&(p->arr[1]));
2638 free_evalue_refs(&(p->arr[0]));
2639 free(p);
2640 value_clear(d);
2641 value_clear(min);
2642 value_clear(max);
2643 Matrix_Free(T);
2644 return r ? r : reduce_in_domain(e, D);
2645 } else if (bounded && equal) {
2646 /* Always zero */
2647 if (p->size == 3)
2648 free_evalue_refs(&(p->arr[2]));
2649 value_clear(e->d);
2650 *e = p->arr[1];
2651 free_evalue_refs(&(p->arr[0]));
2652 free(p);
2653 value_clear(d);
2654 value_clear(min);
2655 value_clear(max);
2656 Matrix_Free(T);
2657 return reduce_in_domain(e, D);
2658 } else if (bounded && value_eq(min, max)) {
2659 /* zero for a single value */
2660 Polyhedron *E;
2661 Matrix *M = Matrix_Alloc(1, D->Dimension+2);
2662 Vector_Copy(T->p[0], M->p[0]+1, D->Dimension+1);
2663 value_multiply(min, min, d);
2664 value_subtract(M->p[0][D->Dimension+1],
2665 M->p[0][D->Dimension+1], min);
2666 E = DomainAddConstraints(D, M, 0);
2667 value_clear(d);
2668 value_clear(min);
2669 value_clear(max);
2670 Matrix_Free(T);
2671 Matrix_Free(M);
2672 r = reduce_in_domain(&p->arr[1], E);
2673 if (p->size == 3)
2674 r |= reduce_in_domain(&p->arr[2], D);
2675 Domain_Free(E);
2676 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2677 return r;
2678 }
2679
2680 value_clear(d);
2681 value_clear(min);
2682 value_clear(max);
2683 Matrix_Free(T);
2684 _reduce_evalue(&p->arr[0].x.p->arr[0], 0, 1);
2685 }
2686
2687 i = p->type == relation ? 1 :
2688 p->type == fractional ? 1 : 0;
2689 for (; i<p->size; i++)
2690 r |= reduce_in_domain(&p->arr[i], D);
2691
2692 if (p->type != fractional) {
2693 if (r && p->type == polynomial) {
2694 evalue f;
2695 value_init(f.d);
2696 value_set_si(f.d, 0);
2697 f.x.p = new_enode(polynomial, 2, p->pos);
2698 evalue_set_si(&f.x.p->arr[0], 0, 1);
2699 evalue_set_si(&f.x.p->arr[1], 1, 1);
2700 reorder_terms(p, &f);
2701 value_clear(e->d);
2702 *e = p->arr[0];
2703 free(p);
2704 }
2705 return r;
2706 }
2707
2708 value_init(d);
2709 value_init(min);
2710 value_init(max);
2711 I = polynomial_projection(p, D, &d, &T, 1);
2712 Matrix_Free(T);
2713 bounded = line_minmax(I, &min, &max); /* frees I */
2714 mpz_fdiv_q(min, min, d);
2715 mpz_fdiv_q(max, max, d);
2716 value_subtract(d, max, min);
2717
2718 if (bounded && value_eq(min, max)) {
2719 evalue inc;
2720 value_init(inc.d);
2721 value_init(inc.x.n);
2722 value_set_si(inc.d, 1);
2723 value_oppose(inc.x.n, min);
2724 eadd(&inc, &p->arr[0]);
2725 reorder_terms(p, &p->arr[0]); /* frees arr[0] */
2726 value_clear(e->d);
2727 *e = p->arr[1];
2728 free(p);
2729 free_evalue_refs(&inc);
2730 r = 1;
2731 } else if (bounded && value_one_p(d) && p->size > 3) {
2732 evalue rem;
2733 evalue inc;
2734 evalue t;
2735 evalue f;
2736 evalue factor;
2737 value_init(rem.d);
2738 value_set_si(rem.d, 0);
2739 rem.x.p = new_enode(fractional, 3, -1);
2740 evalue_copy(&rem.x.p->arr[0], &p->arr[0]);
2741 value_clear(rem.x.p->arr[1].d);
2742 value_clear(rem.x.p->arr[2].d);
2743 rem.x.p->arr[1] = p->arr[1];
2744 rem.x.p->arr[2] = p->arr[2];
2745 for (i = 3; i < p->size; ++i)
2746 p->arr[i-2] = p->arr[i];
2747 p->size -= 2;
2748
2749 value_init(inc.d);
2750 value_init(inc.x.n);
2751 value_set_si(inc.d, 1);
2752 value_oppose(inc.x.n, min);
2753
2754 value_init(t.d);
2755 evalue_copy(&t, &p->arr[0]);
2756 eadd(&inc, &t);
2757
2758 value_init(f.d);
2759 value_set_si(f.d, 0);
2760 f.x.p = new_enode(fractional, 3, -1);
2761 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
2762 evalue_set_si(&f.x.p->arr[1], 1, 1);
2763 evalue_set_si(&f.x.p->arr[2], 2, 1);
2764
2765 value_init(factor.d);
2766 evalue_set_si(&factor, -1, 1);
2767 emul(&t, &factor);
2768
2769 eadd(&f, &factor);
2770 emul(&t, &factor);
2771
2772 value_clear(f.x.p->arr[1].x.n);
2773 value_clear(f.x.p->arr[2].x.n);
2774 evalue_set_si(&f.x.p->arr[1], 0, 1);
2775 evalue_set_si(&f.x.p->arr[2], -1, 1);
2776 eadd(&f, &factor);
2777
2778 emul(&factor, e);
2779 eadd(&rem, e);
2780
2781 free_evalue_refs(&inc);
2782 free_evalue_refs(&t);
2783 free_evalue_refs(&f);
2784 free_evalue_refs(&factor);
2785 free_evalue_refs(&rem);
2786
2787 reduce_in_domain(e, D);
2788
2789 r = 1;
2790 } else {
2791 _reduce_evalue(&p->arr[0], 0, 1);
2792 if (r == 1) {
2793 evalue f;
2794 value_init(f.d);
2795 value_set_si(f.d, 0);
2796 f.x.p = new_enode(fractional, 3, -1);
2797 value_clear(f.x.p->arr[0].d);
2798 f.x.p->arr[0] = p->arr[0];
2799 evalue_set_si(&f.x.p->arr[1], 0, 1);
2800 evalue_set_si(&f.x.p->arr[2], 1, 1);
2801 reorder_terms(p, &f);
2802 value_clear(e->d);
2803 *e = p->arr[1];
2804 free(p);
2805 }
2806 }
2807
2808 value_clear(d);
2809 value_clear(min);
2810 value_clear(max);
2811
2812 return r;
2813 }
2814
2815 void evalue_range_reduction(evalue *e)
2816 {
2817 int i;
2818 if (value_notzero_p(e->d) || e->x.p->type != partition)
2819 return;
2820
2821 for (i = 0; i < e->x.p->size/2; ++i)
2822 if (reduce_in_domain(&e->x.p->arr[2*i+1],
2823 EVALUE_DOMAIN(e->x.p->arr[2*i]))) {
2824 reduce_evalue(&e->x.p->arr[2*i+1]);
2825
2826 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2827 free_evalue_refs(&e->x.p->arr[2*i+1]);
2828 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
2829 value_clear(e->x.p->arr[2*i].d);
2830 e->x.p->size -= 2;
2831 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2832 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2833 --i;
2834 }
2835 }
2836 }
2837
2838 /* Frees EP
2839 */
2840 Enumeration* partition2enumeration(evalue *EP)
2841 {
2842 int i;
2843 Enumeration *en, *res = NULL;
2844
2845 if (EVALUE_IS_ZERO(*EP)) {
2846 free(EP);
2847 return res;
2848 }
2849
2850 for (i = 0; i < EP->x.p->size/2; ++i) {
2851 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[2*i])->Dimension);
2852 en = (Enumeration *)malloc(sizeof(Enumeration));
2853 en->next = res;
2854 res = en;
2855 res->ValidityDomain = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
2856 value_clear(EP->x.p->arr[2*i].d);
2857 res->EP = EP->x.p->arr[2*i+1];
2858 }
2859 free(EP->x.p);
2860 value_clear(EP->d);
2861 free(EP);
2862 return res;
2863 }
2864
2865 int evalue_frac2floor_in_domain(evalue *e, Polyhedron *D)
2866 {
2867 enode *p;
2868 int r = 0;
2869 int i;
2870 Polyhedron *I;
2871 Matrix *T;
2872 Value d, min;
2873 evalue fl;
2874
2875 if (value_notzero_p(e->d))
2876 return r;
2877
2878 p = e->x.p;
2879
2880 i = p->type == relation ? 1 :
2881 p->type == fractional ? 1 : 0;
2882 for (; i<p->size; i++)
2883 r |= evalue_frac2floor_in_domain(&p->arr[i], D);
2884
2885 if (p->type != fractional) {
2886 if (r && p->type == polynomial) {
2887 evalue f;
2888 value_init(f.d);
2889 value_set_si(f.d, 0);
2890 f.x.p = new_enode(polynomial, 2, p->pos);
2891 evalue_set_si(&f.x.p->arr[0], 0, 1);
2892 evalue_set_si(&f.x.p->arr[1], 1, 1);
2893 reorder_terms(p, &f);
2894 value_clear(e->d);
2895 *e = p->arr[0];
2896 free(p);
2897 }
2898 return r;
2899 }
2900
2901 value_init(d);
2902 I = polynomial_projection(p, D, &d, &T, 0);
2903
2904 /*
2905 Polyhedron_Print(stderr, P_VALUE_FMT, I);
2906 */
2907
2908 assert(I->NbEq == 0); /* Should have been reduced */
2909
2910 /* Find minimum */
2911 for (i = 0; i < I->NbConstraints; ++i)
2912 if (value_pos_p(I->Constraint[i][1]))
2913 break;
2914
2915 assert(i < I->NbConstraints);
2916 value_init(min);
2917 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
2918 mpz_cdiv_q(min, I->Constraint[i][2], I->Constraint[i][1]);
2919 if (value_neg_p(min)) {
2920 evalue offset;
2921 mpz_fdiv_q(min, min, d);
2922 value_init(offset.d);
2923 value_set_si(offset.d, 1);
2924 value_init(offset.x.n);
2925 value_oppose(offset.x.n, min);
2926 eadd(&offset, &p->arr[0]);
2927 free_evalue_refs(&offset);
2928 }
2929
2930 Polyhedron_Free(I);
2931 Matrix_Free(T);
2932 value_clear(min);
2933 value_clear(d);
2934
2935 value_init(fl.d);
2936 value_set_si(fl.d, 0);
2937 fl.x.p = new_enode(flooring, 3, -1);
2938 evalue_set_si(&fl.x.p->arr[1], 0, 1);
2939 evalue_set_si(&fl.x.p->arr[2], -1, 1);
2940 evalue_copy(&fl.x.p->arr[0], &p->arr[0]);
2941
2942 eadd(&fl, &p->arr[0]);
2943 reorder_terms(p, &p->arr[0]);
2944 *e = p->arr[1];
2945 free(p);
2946 free_evalue_refs(&fl);
2947
2948 return 1;
2949 }
2950
2951 void evalue_frac2floor(evalue *e)
2952 {
2953 int i;
2954 if (value_notzero_p(e->d) || e->x.p->type != partition)
2955 return;
2956
2957 for (i = 0; i < e->x.p->size/2; ++i)
2958 if (evalue_frac2floor_in_domain(&e->x.p->arr[2*i+1],
2959 EVALUE_DOMAIN(e->x.p->arr[2*i])))
2960 reduce_evalue(&e->x.p->arr[2*i+1]);
2961 }
2962
2963 static Matrix *esum_add_constraint(int nvar, Polyhedron *D, Matrix *C,
2964 Vector *row)
2965 {
2966 int nr, nc;
2967 int i;
2968 int nparam = D->Dimension - nvar;
2969
2970 if (C == 0) {
2971 nr = D->NbConstraints + 2;
2972 nc = D->Dimension + 2 + 1;
2973 C = Matrix_Alloc(nr, nc);
2974 for (i = 0; i < D->NbConstraints; ++i) {
2975 Vector_Copy(D->Constraint[i], C->p[i], 1 + nvar);
2976 Vector_Copy(D->Constraint[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
2977 D->Dimension + 1 - nvar);
2978 }
2979 } else {
2980 Matrix *oldC = C;
2981 nr = C->NbRows + 2;
2982 nc = C->NbColumns + 1;
2983 C = Matrix_Alloc(nr, nc);
2984 for (i = 0; i < oldC->NbRows; ++i) {
2985 Vector_Copy(oldC->p[i], C->p[i], 1 + nvar);
2986 Vector_Copy(oldC->p[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
2987 oldC->NbColumns - 1 - nvar);
2988 }
2989 }
2990 value_set_si(C->p[nr-2][0], 1);
2991 value_set_si(C->p[nr-2][1 + nvar], 1);
2992 value_set_si(C->p[nr-2][nc - 1], -1);
2993
2994 Vector_Copy(row->p, C->p[nr-1], 1 + nvar + 1);
2995 Vector_Copy(row->p + 1 + nvar + 1, C->p[nr-1] + C->NbColumns - 1 - nparam,
2996 1 + nparam);
2997
2998 return C;
2999 }
3000
3001 static void floor2frac_r(evalue *e, int nvar)
3002 {
3003 enode *p;
3004 int i;
3005 evalue f;
3006 evalue *pp;
3007
3008 if (value_notzero_p(e->d))
3009 return;
3010
3011 p = e->x.p;
3012
3013 assert(p->type == flooring);
3014 for (i = 1; i < p->size; i++)
3015 floor2frac_r(&p->arr[i], nvar);
3016
3017 for (pp = &p->arr[0]; value_zero_p(pp->d); pp = &pp->x.p->arr[0]) {
3018 assert(pp->x.p->type == polynomial);
3019 pp->x.p->pos -= nvar;
3020 }
3021
3022 value_init(f.d);
3023 value_set_si(f.d, 0);
3024 f.x.p = new_enode(fractional, 3, -1);
3025 evalue_set_si(&f.x.p->arr[1], 0, 1);
3026 evalue_set_si(&f.x.p->arr[2], -1, 1);
3027 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
3028
3029 eadd(&f, &p->arr[0]);
3030 reorder_terms(p, &p->arr[0]);
3031 *e = p->arr[1];
3032 free(p);
3033 free_evalue_refs(&f);
3034 }
3035
3036 /* Convert flooring back to fractional and shift position
3037 * of the parameters by nvar
3038 */
3039 static void floor2frac(evalue *e, int nvar)
3040 {
3041 floor2frac_r(e, nvar);
3042 reduce_evalue(e);
3043 }
3044
3045 evalue *esum_over_domain_cst(int nvar, Polyhedron *D, Matrix *C)
3046 {
3047 evalue *t;
3048 int nparam = D->Dimension - nvar;
3049
3050 if (C != 0) {
3051 C = Matrix_Copy(C);
3052 D = Constraints2Polyhedron(C, 0);
3053 Matrix_Free(C);
3054 }
3055
3056 t = barvinok_enumerate_e(D, 0, nparam, 0);
3057
3058 /* Double check that D was not unbounded. */
3059 assert(!(value_pos_p(t->d) && value_neg_p(t->x.n)));
3060
3061 if (C != 0)
3062 Polyhedron_Free(D);
3063
3064 return t;
3065 }
3066
3067 evalue *esum_over_domain(evalue *e, int nvar, Polyhedron *D,
3068 Matrix *C)
3069 {
3070 Vector *row = NULL;
3071 int i;
3072 evalue *res;
3073 Matrix *origC;
3074 evalue *factor = NULL;
3075 evalue cum;
3076
3077 if (EVALUE_IS_ZERO(*e))
3078 return 0;
3079
3080 if (D->next) {
3081 Polyhedron *DD = Disjoint_Domain(D, 0, 0);
3082 Polyhedron *Q;
3083
3084 Q = DD;
3085 DD = Q->next;
3086 Q->next = 0;
3087
3088 res = esum_over_domain(e, nvar, Q, C);
3089 Polyhedron_Free(Q);
3090
3091 for (Q = DD; Q; Q = DD) {
3092 evalue *t;
3093
3094 DD = Q->next;
3095 Q->next = 0;
3096
3097 t = esum_over_domain(e, nvar, Q, C);
3098 Polyhedron_Free(Q);
3099
3100 if (!res)
3101 res = t;
3102 else if (t) {
3103 eadd(t, res);
3104 free_evalue_refs(t);
3105 free(t);
3106 }
3107 }
3108 return res;
3109 }
3110
3111 if (value_notzero_p(e->d)) {
3112 evalue *t;
3113
3114 t = esum_over_domain_cst(nvar, D, C);
3115
3116 if (!EVALUE_IS_ONE(*e))
3117 emul(e, t);
3118
3119 return t;
3120 }
3121
3122 switch (e->x.p->type) {
3123 case flooring: {
3124 evalue *pp = &e->x.p->arr[0];
3125
3126 if (pp->x.p->pos > nvar) {
3127 /* remainder is independent of the summated vars */
3128 evalue f;
3129 evalue *t;
3130
3131 value_init(f.d);
3132 evalue_copy(&f, e);
3133 floor2frac(&f, nvar);
3134
3135 t = esum_over_domain_cst(nvar, D, C);
3136
3137 emul(&f, t);
3138
3139 free_evalue_refs(&f);
3140
3141 return t;
3142 }
3143
3144 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3145 poly_denom(pp, &row->p[1 + nvar]);
3146 value_set_si(row->p[0], 1);
3147 for (pp = &e->x.p->arr[0]; value_zero_p(pp->d);
3148 pp = &pp->x.p->arr[0]) {
3149 int pos;
3150 assert(pp->x.p->type == polynomial);
3151 pos = pp->x.p->pos;
3152 if (pos >= 1 + nvar)
3153 ++pos;
3154 value_assign(row->p[pos], row->p[1+nvar]);
3155 value_division(row->p[pos], row->p[pos], pp->x.p->arr[1].d);
3156 value_multiply(row->p[pos], row->p[pos], pp->x.p->arr[1].x.n);
3157 }
3158 value_assign(row->p[1 + D->Dimension + 1], row->p[1+nvar]);
3159 value_division(row->p[1 + D->Dimension + 1],
3160 row->p[1 + D->Dimension + 1],
3161 pp->d);
3162 value_multiply(row->p[1 + D->Dimension + 1],
3163 row->p[1 + D->Dimension + 1],
3164 pp->x.n);
3165 value_oppose(row->p[1 + nvar], row->p[1 + nvar]);
3166 break;
3167 }
3168 case polynomial: {
3169 int pos = e->x.p->pos;
3170
3171 if (pos > nvar) {
3172 factor = ALLOC(evalue);
3173 value_init(factor->d);
3174 value_set_si(factor->d, 0);
3175 factor->x.p = new_enode(polynomial, 2, pos - nvar);
3176 evalue_set_si(&factor->x.p->arr[0], 0, 1);
3177 evalue_set_si(&factor->x.p->arr[1], 1, 1);
3178 break;
3179 }
3180
3181 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3182 for (i = 0; i < D->NbRays; ++i)
3183 if (value_notzero_p(D->Ray[i][pos]))
3184 break;
3185 assert(i < D->NbRays);
3186 if (value_neg_p(D->Ray[i][pos])) {
3187 factor = ALLOC(evalue);
3188 value_init(factor->d);
3189 evalue_set_si(factor, -1, 1);
3190 }
3191 value_set_si(row->p[0], 1);
3192 value_set_si(row->p[pos], 1);
3193 value_set_si(row->p[1 + nvar], -1);
3194 break;
3195 }
3196 default:
3197 assert(0);
3198 }
3199
3200 i = type_offset(e->x.p);
3201
3202 res = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3203 ++i;
3204
3205 if (factor) {
3206 value_init(cum.d);
3207 evalue_copy(&cum, factor);
3208 }
3209
3210 origC = C;
3211 for (; i < e->x.p->size; ++i) {
3212 evalue *t;
3213 if (row) {
3214 Matrix *prevC = C;
3215 C = esum_add_constraint(nvar, D, C, row);
3216 if (prevC != origC)
3217 Matrix_Free(prevC);
3218 }
3219 /*
3220 if (row)
3221 Vector_Print(stderr, P_VALUE_FMT, row);
3222 if (C)
3223 Matrix_Print(stderr, P_VALUE_FMT, C);
3224 */
3225 t = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3226
3227 if (t && factor)
3228 emul(&cum, t);
3229
3230 if (!res)
3231 res = t;
3232 else if (t) {
3233 eadd(t, res);
3234 free_evalue_refs(t);
3235 free(t);
3236 }
3237 if (factor && i+1 < e->x.p->size)
3238 emul(factor, &cum);
3239 }
3240 if (C != origC)
3241 Matrix_Free(C);
3242
3243 if (factor) {
3244 free_evalue_refs(factor);
3245 free_evalue_refs(&cum);
3246 free(factor);
3247 }
3248
3249 if (row)
3250 Vector_Free(row);
3251
3252 reduce_evalue(res);
3253
3254 return res;
3255 }
3256
3257 evalue *esum(evalue *e, int nvar)
3258 {
3259 int i;
3260 evalue *res = ALLOC(evalue);
3261 value_init(res->d);
3262
3263 assert(nvar >= 0);
3264 if (nvar == 0 || EVALUE_IS_ZERO(*e)) {
3265 evalue_copy(res, e);
3266 return res;
3267 }
3268
3269 evalue_set_si(res, 0, 1);
3270
3271 assert(value_zero_p(e->d));
3272 assert(e->x.p->type == partition);
3273
3274 for (i = 0; i < e->x.p->size/2; ++i) {
3275 evalue *t;
3276 t = esum_over_domain(&e->x.p->arr[2*i+1], nvar,
3277 EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
3278 eadd(t, res);
3279 free_evalue_refs(t);
3280 free(t);
3281 }
3282
3283 reduce_evalue(res);
3284
3285 return res;
3286 }
3287
3288 /* Initial silly implementation */
3289 void eor(evalue *e1, evalue *res)
3290 {
3291 evalue E;
3292 evalue mone;
3293 value_init(E.d);
3294 value_init(mone.d);
3295 evalue_set_si(&mone, -1, 1);
3296
3297 evalue_copy(&E, res);
3298 eadd(e1, &E);
3299 emul(e1, res);
3300 emul(&mone, res);
3301 eadd(&E, res);
3302
3303 free_evalue_refs(&E);
3304 free_evalue_refs(&mone);
3305 }
lattice point enumeration library