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