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