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