?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
omega/convert.cc: add conversion from PolyLib to Omega
[barvinok.git] / evalue.c
blobbf21fdeedc88cbb4db4fde0bb4837ebdccc41dc6
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 value_clear(rem.x.p->arr[1].d);
2734 value_clear(rem.x.p->arr[2].d);
2735 rem.x.p->arr[1] = p->arr[1];
2736 rem.x.p->arr[2] = p->arr[2];
2737 for (i = 3; i < p->size; ++i)
2738 p->arr[i-2] = p->arr[i];
2739 p->size -= 2;
2740
2741 value_init(inc.d);
2742 value_init(inc.x.n);
2743 value_set_si(inc.d, 1);
2744 value_oppose(inc.x.n, min);
2745
2746 value_init(t.d);
2747 evalue_copy(&t, &p->arr[0]);
2748 eadd(&inc, &t);
2749
2750 value_init(f.d);
2751 value_set_si(f.d, 0);
2752 f.x.p = new_enode(fractional, 3, -1);
2753 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
2754 evalue_set_si(&f.x.p->arr[1], 1, 1);
2755 evalue_set_si(&f.x.p->arr[2], 2, 1);
2756
2757 value_init(factor.d);
2758 evalue_set_si(&factor, -1, 1);
2759 emul(&t, &factor);
2760
2761 eadd(&f, &factor);
2762 emul(&t, &factor);
2763
2764 value_clear(f.x.p->arr[1].x.n);
2765 value_clear(f.x.p->arr[2].x.n);
2766 evalue_set_si(&f.x.p->arr[1], 0, 1);
2767 evalue_set_si(&f.x.p->arr[2], -1, 1);
2768 eadd(&f, &factor);
2769
2770 emul(&factor, e);
2771 eadd(&rem, e);
2772
2773 free_evalue_refs(&inc);
2774 free_evalue_refs(&t);
2775 free_evalue_refs(&f);
2776 free_evalue_refs(&factor);
2777 free_evalue_refs(&rem);
2778
2779 reduce_in_domain(e, D);
2780
2781 r = 1;
2782 } else {
2783 _reduce_evalue(&p->arr[0], 0, 1);
2784 if (r == 1) {
2785 evalue f;
2786 value_init(f.d);
2787 value_set_si(f.d, 0);
2788 f.x.p = new_enode(fractional, 3, -1);
2789 value_clear(f.x.p->arr[0].d);
2790 f.x.p->arr[0] = p->arr[0];
2791 evalue_set_si(&f.x.p->arr[1], 0, 1);
2792 evalue_set_si(&f.x.p->arr[2], 1, 1);
2793 reorder_terms(p, &f);
2794 value_clear(e->d);
2795 *e = p->arr[1];
2796 free(p);
2797 }
2798 }
2799
2800 value_clear(d);
2801 value_clear(min);
2802 value_clear(max);
2803
2804 return r;
2805 }
2806
2807 void evalue_range_reduction(evalue *e)
2808 {
2809 int i;
2810 if (value_notzero_p(e->d) || e->x.p->type != partition)
2811 return;
2812
2813 for (i = 0; i < e->x.p->size/2; ++i)
2814 if (reduce_in_domain(&e->x.p->arr[2*i+1],
2815 EVALUE_DOMAIN(e->x.p->arr[2*i]))) {
2816 reduce_evalue(&e->x.p->arr[2*i+1]);
2817
2818 if (EVALUE_IS_ZERO(e->x.p->arr[2*i+1])) {
2819 free_evalue_refs(&e->x.p->arr[2*i+1]);
2820 Domain_Free(EVALUE_DOMAIN(e->x.p->arr[2*i]));
2821 value_clear(e->x.p->arr[2*i].d);
2822 e->x.p->size -= 2;
2823 e->x.p->arr[2*i] = e->x.p->arr[e->x.p->size];
2824 e->x.p->arr[2*i+1] = e->x.p->arr[e->x.p->size+1];
2825 --i;
2826 }
2827 }
2828 }
2829
2830 /* Frees EP
2831 */
2832 Enumeration* partition2enumeration(evalue *EP)
2833 {
2834 int i;
2835 Enumeration *en, *res = NULL;
2836
2837 if (EVALUE_IS_ZERO(*EP)) {
2838 free(EP);
2839 return res;
2840 }
2841
2842 for (i = 0; i < EP->x.p->size/2; ++i) {
2843 assert(EP->x.p->pos == EVALUE_DOMAIN(EP->x.p->arr[2*i])->Dimension);
2844 en = (Enumeration *)malloc(sizeof(Enumeration));
2845 en->next = res;
2846 res = en;
2847 res->ValidityDomain = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
2848 value_clear(EP->x.p->arr[2*i].d);
2849 res->EP = EP->x.p->arr[2*i+1];
2850 }
2851 free(EP->x.p);
2852 value_clear(EP->d);
2853 free(EP);
2854 return res;
2855 }
2856
2857 static int frac2floor_in_domain(evalue *e, Polyhedron *D)
2858 {
2859 enode *p;
2860 int r = 0;
2861 int i;
2862 Polyhedron *I;
2863 Matrix *T;
2864 Value d, min;
2865 evalue fl;
2866
2867 if (value_notzero_p(e->d))
2868 return r;
2869
2870 p = e->x.p;
2871
2872 i = p->type == relation ? 1 :
2873 p->type == fractional ? 1 : 0;
2874 for (; i<p->size; i++)
2875 r |= frac2floor_in_domain(&p->arr[i], D);
2876
2877 if (p->type != fractional) {
2878 if (r && p->type == polynomial) {
2879 evalue f;
2880 value_init(f.d);
2881 value_set_si(f.d, 0);
2882 f.x.p = new_enode(polynomial, 2, p->pos);
2883 evalue_set_si(&f.x.p->arr[0], 0, 1);
2884 evalue_set_si(&f.x.p->arr[1], 1, 1);
2885 reorder_terms(p, &f);
2886 value_clear(e->d);
2887 *e = p->arr[0];
2888 free(p);
2889 }
2890 return r;
2891 }
2892
2893 value_init(d);
2894 I = polynomial_projection(p, D, &d, &T, 0);
2895
2896 /*
2897 Polyhedron_Print(stderr, P_VALUE_FMT, I);
2898 */
2899
2900 assert(I->NbEq == 0); /* Should have been reduced */
2901
2902 /* Find minimum */
2903 for (i = 0; i < I->NbConstraints; ++i)
2904 if (value_pos_p(I->Constraint[i][1]))
2905 break;
2906
2907 assert(i < I->NbConstraints);
2908 value_init(min);
2909 value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
2910 mpz_cdiv_q(min, I->Constraint[i][2], I->Constraint[i][1]);
2911 if (value_neg_p(min)) {
2912 evalue offset;
2913 mpz_fdiv_q(min, min, d);
2914 value_init(offset.d);
2915 value_set_si(offset.d, 1);
2916 value_init(offset.x.n);
2917 value_oppose(offset.x.n, min);
2918 eadd(&offset, &p->arr[0]);
2919 free_evalue_refs(&offset);
2920 }
2921
2922 Polyhedron_Free(I);
2923 Matrix_Free(T);
2924 value_clear(min);
2925 value_clear(d);
2926
2927 value_init(fl.d);
2928 value_set_si(fl.d, 0);
2929 fl.x.p = new_enode(flooring, 3, -1);
2930 evalue_set_si(&fl.x.p->arr[1], 0, 1);
2931 evalue_set_si(&fl.x.p->arr[2], -1, 1);
2932 evalue_copy(&fl.x.p->arr[0], &p->arr[0]);
2933
2934 eadd(&fl, &p->arr[0]);
2935 reorder_terms(p, &p->arr[0]);
2936 *e = p->arr[1];
2937 free(p);
2938 free_evalue_refs(&fl);
2939
2940 return 1;
2941 }
2942
2943 void evalue_frac2floor(evalue *e)
2944 {
2945 int i;
2946 if (value_notzero_p(e->d) || e->x.p->type != partition)
2947 return;
2948
2949 for (i = 0; i < e->x.p->size/2; ++i)
2950 if (frac2floor_in_domain(&e->x.p->arr[2*i+1],
2951 EVALUE_DOMAIN(e->x.p->arr[2*i])))
2952 reduce_evalue(&e->x.p->arr[2*i+1]);
2953 }
2954
2955 static Matrix *esum_add_constraint(int nvar, Polyhedron *D, Matrix *C,
2956 Vector *row)
2957 {
2958 int nr, nc;
2959 int i;
2960 int nparam = D->Dimension - nvar;
2961
2962 if (C == 0) {
2963 nr = D->NbConstraints + 2;
2964 nc = D->Dimension + 2 + 1;
2965 C = Matrix_Alloc(nr, nc);
2966 for (i = 0; i < D->NbConstraints; ++i) {
2967 Vector_Copy(D->Constraint[i], C->p[i], 1 + nvar);
2968 Vector_Copy(D->Constraint[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
2969 D->Dimension + 1 - nvar);
2970 }
2971 } else {
2972 Matrix *oldC = C;
2973 nr = C->NbRows + 2;
2974 nc = C->NbColumns + 1;
2975 C = Matrix_Alloc(nr, nc);
2976 for (i = 0; i < oldC->NbRows; ++i) {
2977 Vector_Copy(oldC->p[i], C->p[i], 1 + nvar);
2978 Vector_Copy(oldC->p[i] + 1 + nvar, C->p[i] + 1 + nvar + 1,
2979 oldC->NbColumns - 1 - nvar);
2980 }
2981 }
2982 value_set_si(C->p[nr-2][0], 1);
2983 value_set_si(C->p[nr-2][1 + nvar], 1);
2984 value_set_si(C->p[nr-2][nc - 1], -1);
2985
2986 Vector_Copy(row->p, C->p[nr-1], 1 + nvar + 1);
2987 Vector_Copy(row->p + 1 + nvar + 1, C->p[nr-1] + C->NbColumns - 1 - nparam,
2988 1 + nparam);
2989
2990 return C;
2991 }
2992
2993 static void floor2frac_r(evalue *e, int nvar)
2994 {
2995 enode *p;
2996 int i;
2997 evalue f;
2998 evalue *pp;
2999
3000 if (value_notzero_p(e->d))
3001 return;
3002
3003 p = e->x.p;
3004
3005 assert(p->type == flooring);
3006 for (i = 1; i < p->size; i++)
3007 floor2frac_r(&p->arr[i], nvar);
3008
3009 for (pp = &p->arr[0]; value_zero_p(pp->d); pp = &pp->x.p->arr[0]) {
3010 assert(pp->x.p->type == polynomial);
3011 pp->x.p->pos -= nvar;
3012 }
3013
3014 value_init(f.d);
3015 value_set_si(f.d, 0);
3016 f.x.p = new_enode(fractional, 3, -1);
3017 evalue_set_si(&f.x.p->arr[1], 0, 1);
3018 evalue_set_si(&f.x.p->arr[2], -1, 1);
3019 evalue_copy(&f.x.p->arr[0], &p->arr[0]);
3020
3021 eadd(&f, &p->arr[0]);
3022 reorder_terms(p, &p->arr[0]);
3023 *e = p->arr[1];
3024 free(p);
3025 free_evalue_refs(&f);
3026 }
3027
3028 /* Convert flooring back to fractional and shift position
3029 * of the parameters by nvar
3030 */
3031 static void floor2frac(evalue *e, int nvar)
3032 {
3033 floor2frac_r(e, nvar);
3034 reduce_evalue(e);
3035 }
3036
3037 evalue *esum_over_domain_cst(int nvar, Polyhedron *D, Matrix *C)
3038 {
3039 evalue *t;
3040 int nparam = D->Dimension - nvar;
3041
3042 if (C != 0) {
3043 C = Matrix_Copy(C);
3044 D = Constraints2Polyhedron(C, 0);
3045 Matrix_Free(C);
3046 }
3047
3048 t = barvinok_enumerate_e(D, 0, nparam, 0);
3049
3050 /* Double check that D was not unbounded. */
3051 assert(!(value_pos_p(t->d) && value_neg_p(t->x.n)));
3052
3053 if (C != 0)
3054 Polyhedron_Free(D);
3055
3056 return t;
3057 }
3058
3059 evalue *esum_over_domain(evalue *e, int nvar, Polyhedron *D,
3060 Matrix *C)
3061 {
3062 Vector *row = NULL;
3063 int i;
3064 evalue *res;
3065 Matrix *origC;
3066 evalue *factor = NULL;
3067 evalue cum;
3068
3069 if (EVALUE_IS_ZERO(*e))
3070 return 0;
3071
3072 if (D->next) {
3073 Polyhedron *DD = Disjoint_Domain(D, 0, 0);
3074 Polyhedron *Q;
3075
3076 Q = DD;
3077 DD = Q->next;
3078 Q->next = 0;
3079
3080 res = esum_over_domain(e, nvar, Q, C);
3081 Polyhedron_Free(Q);
3082
3083 for (Q = DD; Q; Q = DD) {
3084 evalue *t;
3085
3086 DD = Q->next;
3087 Q->next = 0;
3088
3089 t = esum_over_domain(e, nvar, Q, C);
3090 Polyhedron_Free(Q);
3091
3092 if (!res)
3093 res = t;
3094 else if (t) {
3095 eadd(t, res);
3096 free_evalue_refs(t);
3097 free(t);
3098 }
3099 }
3100 return res;
3101 }
3102
3103 if (value_notzero_p(e->d)) {
3104 evalue *t;
3105
3106 t = esum_over_domain_cst(nvar, D, C);
3107
3108 if (!EVALUE_IS_ONE(*e))
3109 emul(e, t);
3110
3111 return t;
3112 }
3113
3114 switch (e->x.p->type) {
3115 case flooring: {
3116 evalue *pp = &e->x.p->arr[0];
3117
3118 if (pp->x.p->pos > nvar) {
3119 /* remainder is independent of the summated vars */
3120 evalue f;
3121 evalue *t;
3122
3123 value_init(f.d);
3124 evalue_copy(&f, e);
3125 floor2frac(&f, nvar);
3126
3127 t = esum_over_domain_cst(nvar, D, C);
3128
3129 emul(&f, t);
3130
3131 free_evalue_refs(&f);
3132
3133 return t;
3134 }
3135
3136 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3137 poly_denom(pp, &row->p[1 + nvar]);
3138 value_set_si(row->p[0], 1);
3139 for (pp = &e->x.p->arr[0]; value_zero_p(pp->d);
3140 pp = &pp->x.p->arr[0]) {
3141 int pos;
3142 assert(pp->x.p->type == polynomial);
3143 pos = pp->x.p->pos;
3144 if (pos >= 1 + nvar)
3145 ++pos;
3146 value_assign(row->p[pos], row->p[1+nvar]);
3147 value_division(row->p[pos], row->p[pos], pp->x.p->arr[1].d);
3148 value_multiply(row->p[pos], row->p[pos], pp->x.p->arr[1].x.n);
3149 }
3150 value_assign(row->p[1 + D->Dimension + 1], row->p[1+nvar]);
3151 value_division(row->p[1 + D->Dimension + 1],
3152 row->p[1 + D->Dimension + 1],
3153 pp->d);
3154 value_multiply(row->p[1 + D->Dimension + 1],
3155 row->p[1 + D->Dimension + 1],
3156 pp->x.n);
3157 value_oppose(row->p[1 + nvar], row->p[1 + nvar]);
3158 break;
3159 }
3160 case polynomial: {
3161 int pos = e->x.p->pos;
3162
3163 if (pos > nvar) {
3164 ALLOC(factor);
3165 value_init(factor->d);
3166 value_set_si(factor->d, 0);
3167 factor->x.p = new_enode(polynomial, 2, pos - nvar);
3168 evalue_set_si(&factor->x.p->arr[0], 0, 1);
3169 evalue_set_si(&factor->x.p->arr[1], 1, 1);
3170 break;
3171 }
3172
3173 row = Vector_Alloc(1 + D->Dimension + 1 + 1);
3174 for (i = 0; i < D->NbRays; ++i)
3175 if (value_notzero_p(D->Ray[i][pos]))
3176 break;
3177 assert(i < D->NbRays);
3178 if (value_neg_p(D->Ray[i][pos])) {
3179 ALLOC(factor);
3180 value_init(factor->d);
3181 evalue_set_si(factor, -1, 1);
3182 }
3183 value_set_si(row->p[0], 1);
3184 value_set_si(row->p[pos], 1);
3185 value_set_si(row->p[1 + nvar], -1);
3186 break;
3187 }
3188 default:
3189 assert(0);
3190 }
3191
3192 i = type_offset(e->x.p);
3193
3194 res = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3195 ++i;
3196
3197 if (factor) {
3198 value_init(cum.d);
3199 evalue_copy(&cum, factor);
3200 }
3201
3202 origC = C;
3203 for (; i < e->x.p->size; ++i) {
3204 evalue *t;
3205 if (row) {
3206 Matrix *prevC = C;
3207 C = esum_add_constraint(nvar, D, C, row);
3208 if (prevC != origC)
3209 Matrix_Free(prevC);
3210 }
3211 /*
3212 if (row)
3213 Vector_Print(stderr, P_VALUE_FMT, row);
3214 if (C)
3215 Matrix_Print(stderr, P_VALUE_FMT, C);
3216 */
3217 t = esum_over_domain(&e->x.p->arr[i], nvar, D, C);
3218
3219 if (t && factor)
3220 emul(&cum, t);
3221
3222 if (!res)
3223 res = t;
3224 else if (t) {
3225 eadd(t, res);
3226 free_evalue_refs(t);
3227 free(t);
3228 }
3229 if (factor && i+1 < e->x.p->size)
3230 emul(factor, &cum);
3231 }
3232 if (C != origC)
3233 Matrix_Free(C);
3234
3235 if (factor) {
3236 free_evalue_refs(factor);
3237 free_evalue_refs(&cum);
3238 free(factor);
3239 }
3240
3241 if (row)
3242 Vector_Free(row);
3243
3244 reduce_evalue(res);
3245
3246 return res;
3247 }
3248
3249 evalue *esum(evalue *e, int nvar)
3250 {
3251 int i;
3252 evalue *res;
3253 ALLOC(res);
3254 value_init(res->d);
3255
3256 assert(nvar >= 0);
3257 if (nvar == 0 || EVALUE_IS_ZERO(*e)) {
3258 evalue_copy(res, e);
3259 return res;
3260 }
3261
3262 evalue_set_si(res, 0, 1);
3263
3264 assert(value_zero_p(e->d));
3265 assert(e->x.p->type == partition);
3266
3267 for (i = 0; i < e->x.p->size/2; ++i) {
3268 evalue *t;
3269 t = esum_over_domain(&e->x.p->arr[2*i+1], nvar,
3270 EVALUE_DOMAIN(e->x.p->arr[2*i]), 0);
3271 eadd(t, res);
3272 free_evalue_refs(t);
3273 free(t);
3274 }
3275
3276 reduce_evalue(res);
3277
3278 return res;
3279 }
3280
3281 /* Initial silly implementation */
3282 void eor(evalue *e1, evalue *res)
3283 {
3284 evalue E;
3285 evalue mone;
3286 value_init(E.d);
3287 value_init(mone.d);
3288 evalue_set_si(&mone, -1, 1);
3289
3290 evalue_copy(&E, res);
3291 eadd(e1, &E);
3292 emul(e1, res);
3293 emul(&mone, res);
3294 eadd(&E, res);
3295
3296 free_evalue_refs(&E);
3297 free_evalue_refs(&mone);
3298 }
lattice point enumeration library