?  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] / scarf.cc
blob7524f356e034bd88d4590e96d3980de8b66c1093
1 #include <vector>
2 #include <barvinok/barvinok.h>
3 #include <barvinok/util.h>
4 #include "scarf.h"
5
6 using std::vector;
7
8 static Matrix *extract_matrix(Polyhedron *P, unsigned dim)
9 {
10 Matrix *A;
11 int n_col;
12
13 n_col = 0;
14 for (int i = 0; i < P->NbConstraints; ++i)
15 if (value_notzero_p(P->Constraint[i][1+dim]) ||
16 value_notzero_p(P->Constraint[i][1+dim+1]))
17 ++n_col;
18
19 A = Matrix_Alloc(2, n_col+2);
20 n_col = 0;
21 for (int i = 0; i < P->NbConstraints; ++i) {
22 if (value_zero_p(P->Constraint[i][1+dim]) &&
23 value_zero_p(P->Constraint[i][1+dim+1]))
24 continue;
25 value_assign(A->p[0][n_col], P->Constraint[i][1+dim]);
26 value_assign(A->p[1][n_col], P->Constraint[i][1+dim+1]);
27 ++n_col;
28 }
29 value_set_si(A->p[0][n_col], 1);
30 value_set_si(A->p[1][n_col+1], 1);
31
32 return A;
33 }
34
35 static int lex_sign(Value *v, int len)
36 {
37 int first;
38
39 first = First_Non_Zero(v, len);
40 return first == -1 ? 0 : value_sign(v[first]);
41 }
42
43 static void set_pos(int pos[4], int actual, int wanted)
44 {
45 if (actual == wanted)
46 return;
47 int t = pos[actual];
48 pos[actual] = pos[wanted];
49 pos[wanted] = t;
50 }
51
52 static Matrix *normalize_matrix(Matrix *A, int pos[4], int n)
53 {
54 Matrix *T, *B;
55 Value tmp, tmp2, factor;
56 int type = -1;
57
58 value_init(tmp);
59 value_init(tmp2);
60 value_init(factor);
61
62 T = Matrix_Alloc(2, 2);
63 Extended_Euclid(A->p[0][pos[0]], A->p[1][pos[0]],
64 &T->p[0][0], &T->p[0][1], &tmp);
65 value_division(T->p[1][0], A->p[1][pos[0]], tmp);
66 value_division(T->p[1][1], A->p[0][pos[0]], tmp);
67 value_oppose(T->p[0][0], T->p[0][0]);
68 value_oppose(T->p[0][1], T->p[0][1]);
69 value_oppose(T->p[1][0], T->p[1][0]);
70
71 B = Matrix_Alloc(2, A->NbColumns);
72 Matrix_Product(T, A, B);
73 Matrix_Free(T);
74
75 if (lex_sign(B->p[1], B->NbColumns) > 0) {
76 value_set_si(tmp, -1);
77 Vector_Scale(B->p[1], B->p[1], tmp, B->NbColumns);
78 }
79
80 if (n == 3) {
81 assert(value_neg_p(B->p[1][pos[1]]));
82 assert(value_pos_p(B->p[1][pos[2]]));
83 assert(value_pos_p(B->p[0][pos[1]]));
84 assert(value_neg_p(B->p[0][pos[2]]));
85 Vector_Exchange(B->p[0], B->p[1], B->NbColumns);
86 type = 2;
87 } else {
88 int i;
89 for (i = 1; i <= 3; ++i)
90 if (value_zero_p(B->p[1][pos[i]]))
91 break;
92 if (i <= 3) {
93 /* put zero in position 3 */
94 set_pos(pos, i, 3);
95
96 /* put positive one in position 1 */
97 for (i = 1; i <= 3; ++i)
98 if (value_pos_p(B->p[1][pos[i]]))
99 break;
100 set_pos(pos, i, 1);
101
102 value_set_si(factor, 0);
103 for (int i = 1; i <= 2; ++i) {
104 value_pdivision(tmp, B->p[0][pos[i]], B->p[1][pos[i]]);
105 value_increment(tmp, tmp);
106 if (value_gt(tmp, factor))
107 value_assign(factor, tmp);
108 }
109 value_oppose(factor, factor);
110 value_set_si(tmp, 1);
111 Vector_Combine(B->p[0], B->p[1], B->p[0], tmp, factor, B->NbColumns);
112
113 assert(value_notzero_p(B->p[0][pos[3]]));
114 type = value_pos_p(B->p[0][pos[3]]) ? 1 : 2;
115 } else {
116 int neg = 0;
117 for (int i = 1; i <= 3; ++i)
118 if (value_neg_p(B->p[1][pos[i]]))
119 ++neg;
120 assert(neg == 1 || neg == 2);
121 if (neg == 1) {
122 /* We will deal with this later */
123 assert(0);
124 } else {
125 int i;
126 /* put positive one in position 1 */
127 for (i = 1; i <= 3; ++i)
128 if (value_pos_p(B->p[1][pos[i]]))
129 break;
130 set_pos(pos, i, 1);
131
132 value_set_si(factor, 0);
133 for (int i = 1; i <= 3; ++i) {
134 value_pdivision(tmp, B->p[0][pos[i]], B->p[1][pos[i]]);
135 value_increment(tmp, tmp);
136 if (value_gt(tmp, factor))
137 value_assign(factor, tmp);
138 }
139 value_oppose(factor, factor);
140 value_set_si(tmp, 1);
141 Vector_Combine(B->p[0], B->p[1], B->p[0],
142 tmp, factor, B->NbColumns);
143 type = 1;
144 }
145 }
146 }
147
148 assert(type != -1);
149
150 if (type == 2) {
151 for (;;) {
152 value_oppose(tmp, B->p[0][pos[1]]);
153 value_pdivision(factor, tmp, B->p[1][pos[1]]);
154 value_oppose(tmp, B->p[1][pos[2]]);
155 value_pdivision(tmp, tmp, B->p[0][pos[2]]);
156 if (value_zero_p(factor) && value_zero_p(tmp))
157 break;
158 assert(value_zero_p(factor) || value_zero_p(tmp));
159 if (value_pos_p(factor)) {
160 value_set_si(tmp, 1);
161 Vector_Combine(B->p[0], B->p[1], B->p[0], tmp, factor, B->NbColumns);
162 if (value_zero_p(B->p[0][pos[1]])) {
163 /* We will deal with this later */
164 assert(lex_sign(B->p[0], B->NbColumns) < 0);
165 }
166 } else {
167 value_set_si(factor, 1);
168 Vector_Combine(B->p[0], B->p[1], B->p[1], tmp, factor, B->NbColumns);
169 if (value_zero_p(B->p[1][pos[2]])) {
170 /* We will deal with this later */
171 assert(lex_sign(B->p[1], B->NbColumns) < 0);
172 }
173 }
174 }
175 } else {
176 int neg;
177 int sign;
178 bool progress = true;
179 while (progress) {
180 progress = false;
181 for (int i = 0; i <= 1; ++i) {
182 value_set_si(factor, -1);
183 for (int j = 1; j <= 3; ++j) {
184 if (value_zero_p(B->p[1-i][pos[j]]))
185 continue;
186 value_oppose(tmp, B->p[i][pos[j]]);
187 value_pdivision(tmp, tmp, B->p[1-i][pos[j]]);
188 if (value_neg_p(factor) || value_lt(tmp, factor))
189 value_assign(factor, tmp);
190 }
191 if (value_pos_p(factor)) {
192 value_set_si(tmp, 1);
193 Vector_Combine(B->p[i], B->p[1-i], B->p[i], tmp, factor,
194 B->NbColumns);
195 sign = lex_sign(B->p[i], B->NbColumns);
196 for (int j = 1; j <= 3; ++j) {
197 if (value_notzero_p(B->p[i][pos[j]]))
198 continue;
199 /* a zero is interpreted to be of sign sign */
200 if ((sign > 0 && value_pos_p(B->p[1-i][pos[j]])) ||
201 (sign < 0 && value_neg_p(B->p[1-i][pos[j]]))) {
202 /* the zero is of the wrong sign => back-off one */
203 value_set_si(tmp2, -1);
204 Vector_Combine(B->p[i], B->p[1-i], B->p[i], tmp, tmp2,
205 B->NbColumns);
206 value_decrement(factor, factor);
207 }
208 }
209 /* We may have backed-off, so we need to check again. */
210 if (value_pos_p(factor))
211 progress = true;
212 }
213 }
214 }
215 sign = 0;
216 for (int i = 0; i < B->NbColumns; ++i) {
217 value_addto(tmp, B->p[0][i], B->p[1][i]);
218 if (value_zero_p(tmp))
219 continue;
220 sign = value_neg_p(tmp) ? -1 : 1;
221 break;
222 }
223 neg = 0;
224 for (int i = 1; i <= 3; ++i) {
225 value_addto(tmp, B->p[0][pos[i]], B->p[1][pos[i]]);
226 if (value_neg_p(tmp) || (sign < 0 && value_zero_p(tmp)))
227 ++neg;
228 }
229 assert(neg <= 2);
230 switch(neg) {
231 int i;
232 case 1:
233 value_set_si(tmp, 1);
234 Vector_Combine(B->p[0], B->p[1], B->p[1], tmp, tmp, B->NbColumns);
235
236 /* put positive pair in position 3 */
237 for (i = 1; i <= 3; ++i)
238 if (value_pos_p(B->p[0][pos[i]]) && value_pos_p(B->p[1][pos[i]]))
239 break;
240 assert(i <= 3);
241 if (i != 3) {
242 int t = pos[i];
243 pos[i] = pos[3];
244 pos[3] = t;
245 }
246
247 break;
248 case 0:
249 case 2:
250 /* We will deal with these later */
251 assert(0);
252 }
253 }
254
255 value_clear(tmp);
256 value_clear(tmp2);
257 value_clear(factor);
258
259 return B;
260 }
261
262 struct simplex {
263 Value last; // last multiple of offset in link
264 Vector *offset;
265 Matrix *M; // rows: elements different from (0,0)
266 int mask;
267
268 simplex(int d) {
269 M = Matrix_Alloc(d, 2);
270 offset = NULL;
271 }
272 simplex(int d, int mask, Value last) {
273 M = Matrix_Alloc(d, 2);
274 offset = Vector_Alloc(2);
275 value_init(this->last);
276 value_assign(this->last, last);
277 this->mask = mask;
278 }
279 void transform(Matrix *T);
280 void normalize();
281 Polyhedron *shrunk_polyhedron(Polyhedron *P, int dim, Matrix *A,
282 unsigned MaxRays);
283 void print(FILE *out);
284 };
285
286 void simplex::print(FILE *out)
287 {
288 if (!offset)
289 Matrix_Print(out, P_VALUE_FMT, M);
290 else {
291 fprintf(out, "%d %d\n", M->NbRows, M->NbColumns);
292 for (int j = 0; j < M->NbRows; ++j) {
293 for (int k = 0; k < M->NbColumns; ++k)
294 value_print(out, P_VALUE_FMT, M->p[j][k]);
295 if (mask & (1 << j)) {
296 fprintf(out, " + k * ");
297 for (int k = 0; k < M->NbColumns; ++k)
298 value_print(out, P_VALUE_FMT, offset->p[k]);
299 }
300 fprintf(out, "\n");
301 }
302 fprintf(out, "\t0 <= k <= ");
303 value_print(out, P_VALUE_FMT, last);
304 fprintf(out, "\n");
305 }
306 }
307
308 static bool lex_smaller(Value *v1, Value *v2, int n)
309 {
310 for (int i = 0; i < n; ++i)
311 if (value_lt(v1[i], v2[i]))
312 return true;
313 else if (value_gt(v1[i], v2[i]))
314 return false;
315 return false;
316 }
317
318 void simplex::transform(Matrix *T)
319 {
320 Matrix *M2 = M;
321 M = Matrix_Alloc(M2->NbRows, M2->NbColumns);
322 Matrix_Product(M2, T, M);
323 Matrix_Free(M2);
324
325 if (offset) {
326 Vector *offset2 = offset;
327 offset = Vector_Alloc(offset2->Size);
328 Vector_Matrix_Product(offset2->p, T, offset->p);
329 Vector_Free(offset2);
330 }
331 }
332
333 void simplex::normalize()
334 {
335 int lexmin = 0;
336 for (int i = 1; i < M->NbRows; ++i)
337 if (lex_smaller(M->p[i], M->p[lexmin], 2))
338 lexmin = i;
339 if (lex_sign(M->p[lexmin], 2) < 0) {
340 Value tmp;
341 value_init(tmp);
342 value_set_si(tmp, -1);
343 Vector_Scale(M->p[lexmin], M->p[lexmin], tmp, 2);
344 value_set_si(tmp, 1);
345 for (int i = 0; i < M->NbRows; ++i) {
346 if (i == lexmin)
347 continue;
348 Vector_Combine(M->p[lexmin], M->p[i], M->p[i], tmp, tmp, 2);
349 }
350 if (offset && (mask & (1 << lexmin))) {
351 value_set_si(tmp, -1);
352 Vector_Scale(offset->p, offset->p, tmp, 2);
353 mask ^= (1 << M->NbRows) - 1 - (1 << lexmin);
354 }
355 value_clear(tmp);
356 }
357 }
358
359 static Matrix *InsertColumn(Matrix *M, int pos)
360 {
361 Matrix *R;
362 int i;
363
364 R = Matrix_Alloc(M->NbRows, M->NbColumns+1);
365 for (i = 0; i < M->NbRows; ++i) {
366 Vector_Copy(M->p[i], R->p[i], pos);
367 Vector_Copy(M->p[i]+pos, R->p[i]+pos+1, M->NbColumns-pos);
368 }
369 return R;
370 }
371
372 Polyhedron *simplex::shrunk_polyhedron(Polyhedron *P, int dim, Matrix *A,
373 unsigned MaxRays)
374 {
375 Matrix *Constraints, *b;
376 Vector *b_offset = NULL;
377 Polyhedron *Q;
378 Value min, min_var, tmp;
379 value_init(tmp);
380 value_init(min);
381 value_init(min_var);
382 int constant;
383
384 b = Matrix_Alloc(M->NbRows, A->NbColumns);
385 Matrix_Product(M, A, b);
386
387 if (offset) {
388 b_offset = Vector_Alloc(A->NbColumns);
389 Vector_Matrix_Product(offset->p, A, b_offset->p);
390 }
391
392 if (!offset)
393 Constraints = Polyhedron2Constraints(P);
394 else {
395 Constraints = Matrix_Alloc(P->NbConstraints+2, P->Dimension+2+1);
396 for (int i = 0; i < P->NbConstraints; ++i) {
397 Vector_Copy(P->Constraint[i], Constraints->p[i], 1+dim+2);
398 Vector_Copy(P->Constraint[i]+1+dim+2, Constraints->p[i]+1+dim+2+1,
399 (P->Dimension+2)-(1+dim+2));
400 }
401 value_set_si(Constraints->p[P->NbConstraints][0], 1);
402 value_set_si(Constraints->p[P->NbConstraints][1+dim+2], 1);
403 value_set_si(Constraints->p[P->NbConstraints+1][0], 1);
404 value_set_si(Constraints->p[P->NbConstraints+1][1+dim+2], -1);
405 value_assign(Constraints->p[P->NbConstraints+1][Constraints->NbColumns-1],
406 last);
407 }
408 constant = Constraints->NbColumns - 1;
409
410 for (int i = 0, j = 0; i < P->NbConstraints; ++i) {
411 if (value_zero_p(Constraints->p[i][1+dim]) &&
412 value_zero_p(Constraints->p[i][1+dim+1]))
413 continue;
414 value_set_si(min, 0);
415 for (int k = 0; k < b->NbRows; ++k) {
416 if (offset && (mask & (1 << k)))
417 continue;
418 if (value_lt(b->p[k][j], min))
419 value_assign(min, b->p[k][j]);
420 }
421 value_set_si(min_var, 0);
422 if (offset) {
423 if (value_neg_p(b_offset->p[j])) {
424 value_oppose(min_var, b_offset->p[j]);
425 value_multiply(min_var, min_var, last);
426 value_increment(min_var, min_var);
427 }
428 for (int k = 0; k < b->NbRows; ++k) {
429 if (!(mask & (1 << k)))
430 continue;
431 if (value_lt(b->p[k][j], min_var))
432 value_assign(min_var, b->p[k][j]);
433 }
434 }
435 if (!offset || value_pos_p(b_offset->p[j])) {
436 if (value_le(min, min_var))
437 value_addto(Constraints->p[i][constant],
438 Constraints->p[i][constant], min);
439 else {
440 value_assign(tmp, min_var);
441 value_addmul(tmp, last, b_offset->p[j]);
442 if (value_le(tmp, min)) {
443 value_addto(Constraints->p[i][constant],
444 Constraints->p[i][constant], min_var);
445 value_addto(Constraints->p[i][1+dim+2],
446 Constraints->p[i][1+dim+2], b_offset->p[j]);
447 } else {
448 int lastrow = Constraints->NbRows;
449 int cols = Constraints->NbColumns;
450 Matrix *C = Constraints;
451 Constraints = AddANullRow(Constraints);
452 Matrix_Free(C);
453 Vector_Copy(Constraints->p[i], Constraints->p[lastrow], cols);
454 value_addto(Constraints->p[i][constant],
455 Constraints->p[i][constant], min_var);
456 value_addto(Constraints->p[i][1+dim+2],
457 Constraints->p[i][1+dim+2], b_offset->p[j]);
458 value_addto(Constraints->p[lastrow][constant],
459 Constraints->p[lastrow][constant], min);
460 }
461 }
462 } else {
463 if (value_le(min_var, min)) {
464 value_addto(Constraints->p[i][constant],
465 Constraints->p[i][constant], min_var);
466 value_addto(Constraints->p[i][1+dim+2],
467 Constraints->p[i][1+dim+2], b_offset->p[j]);
468 } else {
469 value_assign(tmp, min_var);
470 value_addmul(tmp, last, b_offset->p[j]);
471 if (value_le(min, tmp)) {
472 value_addto(Constraints->p[i][constant],
473 Constraints->p[i][constant], min);
474 } else {
475 int lastrow = Constraints->NbRows;
476 int cols = Constraints->NbColumns;
477 Matrix *C = Constraints;
478 Constraints = AddANullRow(Constraints);
479 Matrix_Free(C);
480 Vector_Copy(Constraints->p[i], Constraints->p[lastrow], cols);
481 value_addto(Constraints->p[i][constant],
482 Constraints->p[i][constant], min_var);
483 value_addto(Constraints->p[i][1+dim+2],
484 Constraints->p[i][1+dim+2], b_offset->p[j]);
485 value_addto(Constraints->p[lastrow][constant],
486 Constraints->p[lastrow][constant], min);
487 }
488 }
489 }
490 ++j;
491 }
492 Q = Constraints2Polyhedron(Constraints, MaxRays);
493
494 if (b_offset)
495 Vector_Free(b_offset);
496 Matrix_Free(b);
497 Matrix_Free(Constraints);
498 value_clear(tmp);
499 value_clear(min);
500 value_clear(min_var);
501
502 return Q;
503 }
504
505 struct scarf_complex {
506 vector<simplex> simplices;
507 void add(Matrix *B, int pos[4], int n);
508 void add(Matrix *T, simplex s);
509 void print(FILE *out);
510 ~scarf_complex() {
511 for (int i = 0; i < simplices.size(); ++i) {
512 Matrix_Free(simplices[i].M);
513 if (simplices[i].offset) {
514 Vector_Free(simplices[i].offset);
515 value_clear(simplices[i].last);
516 }
517 }
518 }
519 };
520
521 void scarf_complex::add(Matrix *T, simplex s)
522 {
523 s.transform(T);
524 s.normalize();
525 if (s.offset && lex_sign(s.offset->p, 2) < 0) {
526 Value factor;
527 Value tmp;
528 value_init(factor);
529 value_init(tmp);
530 /* compute the smallest multiple (factor) of the offset that
531 * makes on of the vertices lexico-negative.
532 */
533 int lexmin = -1;
534 for (int i = 0; i < s.M->NbRows; ++i) {
535 if (!(s.mask & (1 << i)))
536 continue;
537 if (lexmin == -1 || lex_smaller(s.M->p[i], s.M->p[lexmin], 2))
538 lexmin = i;
539 }
540 if (value_zero_p(s.offset->p[0])) {
541 if (value_pos_p(s.M->p[lexmin][0]))
542 value_increment(factor, s.last);
543 else {
544 value_oppose(factor, s.M->p[lexmin][1]);
545 mpz_cdiv_q(factor, factor, s.offset->p[1]);
546 }
547 } else {
548 value_oppose(factor, s.M->p[lexmin][0]);
549 mpz_cdiv_q(factor, factor, s.offset->p[0]);
550 if (mpz_divisible_p(s.M->p[lexmin][0], s.offset->p[0])) {
551 value_assign(tmp, s.M->p[lexmin][1]);
552 value_addmul(tmp, factor, s.offset->p[1]);
553 if (value_pos_p(tmp))
554 value_increment(factor, factor);
555 }
556 }
557 if (value_le(factor, s.last)) {
558 simplex part(s.M->NbRows, s.mask, s.last);
559 Vector_Copy(s.offset->p, part.offset->p, 2);
560 value_set_si(tmp, 1);
561 for (int i = 0; i < s.M->NbRows; ++i) {
562 if (s.mask & (1 << i))
563 Vector_Combine(s.M->p[i], s.offset->p, part.M->p[i],
564 tmp, factor, 2);
565 else
566 Vector_Copy(s.M->p[i], part.M->p[i], 2);
567 }
568 value_subtract(part.last, part.last, factor);
569 value_decrement(s.last, factor);
570 part.normalize();
571 simplices.push_back(part);
572 }
573 value_clear(tmp);
574 value_clear(factor);
575 }
576 simplices.push_back(s);
577 }
578
579 void scarf_complex::add(Matrix *B, int pos[4], int n)
580 {
581 Matrix *T;
582
583 T = Matrix_Alloc(2, 2);
584 Vector_Copy(B->p[0]+B->NbColumns-2, T->p[0], 2);
585 Vector_Copy(B->p[1]+B->NbColumns-2, T->p[1], 2);
586
587 if (n == 3 || value_neg_p(B->p[0][pos[3]])) {
588 assert(n == 3 || value_neg_p(B->p[1][pos[3]]));
589
590 simplex s1(1);
591 value_set_si(s1.M->p[0][0], 0);
592 value_set_si(s1.M->p[0][1], 1);
593 add(T, s1);
594
595 simplex s2(1);
596 value_set_si(s2.M->p[0][0], 1);
597 value_set_si(s2.M->p[0][1], 1);
598 add(T, s2);
599
600 simplex s3(1);
601 value_set_si(s3.M->p[0][0], 1);
602 value_set_si(s3.M->p[0][1], 0);
603 add(T, s3);
604
605 simplex s4(2);
606 value_set_si(s4.M->p[0][0], 0);
607 value_set_si(s4.M->p[0][1], 1);
608 value_set_si(s4.M->p[1][0], 1);
609 value_set_si(s4.M->p[1][1], 1);
610 add(T, s4);
611
612 simplex s5(2);
613 value_set_si(s5.M->p[0][0], 1);
614 value_set_si(s5.M->p[0][1], 0);
615 value_set_si(s5.M->p[1][0], 1);
616 value_set_si(s5.M->p[1][1], 1);
617 add(T, s5);
618 } else {
619 Matrix *h;
620 Vector *offset;
621 bool initial = true;
622 bool progress = true;
623 Value tmp, tmp2, factor;
624 int sign;
625
626 value_init(tmp);
627 value_init(tmp2);
628 value_init(factor);
629
630 assert(value_pos_p(B->p[0][pos[3]]));
631 assert(value_pos_p(B->p[1][pos[3]]));
632
633 h = Matrix_Alloc(3, 2);
634 value_set_si(h->p[0][0], 1);
635 value_set_si(h->p[0][1], 0);
636 value_set_si(h->p[1][0], 0);
637 value_set_si(h->p[1][1], 1);
638 value_set_si(h->p[2][0], 1);
639 value_set_si(h->p[2][1], 1);
640
641 offset = Vector_Alloc(2);
642
643 while (progress) {
644 progress = false;
645 for (int i = 0; i <= 1; ++i) {
646 value_set_si(factor, -1);
647 for (int j = 1; j <= 2; ++j) {
648 if (value_zero_p(B->p[1-i][pos[j]]))
649 continue;
650 value_oppose(tmp, B->p[i][pos[j]]);
651 value_pdivision(tmp, tmp, B->p[1-i][pos[j]]);
652 if (value_neg_p(factor) || value_lt(tmp, factor))
653 value_assign(factor, tmp);
654 }
655 if (value_pos_p(factor)) {
656 value_set_si(tmp, 1);
657 Vector_Combine(B->p[i], B->p[1-i], B->p[i], tmp, factor,
658 B->NbColumns);
659 sign = lex_sign(B->p[i], B->NbColumns);
660 for (int j = 1; j <= 2; ++j) {
661 if (value_notzero_p(B->p[i][pos[j]]))
662 continue;
663 /* a zero is interpreted to be of sign sign */
664 if ((sign > 0 && value_pos_p(B->p[1-i][pos[j]])) ||
665 (sign < 0 && value_neg_p(B->p[1-i][pos[j]]))) {
666 /* the zero is of the wrong sign => back-off one */
667 value_set_si(tmp2, -1);
668 Vector_Combine(B->p[i], B->p[1-i], B->p[i], tmp, tmp2,
669 B->NbColumns);
670 value_decrement(factor, factor);
671 }
672 }
673 /* We may have backed-off, so we need to check again. */
674 if (value_pos_p(factor)) {
675 progress = true;
676 value_set_si(tmp, 1);
677 value_set_si(tmp2, -1);
678
679 Vector_Combine(h->p[2], h->p[i], offset->p, tmp, tmp2, 2);
680
681 if (initial) {
682 simplex l1(1);
683 Vector_Copy(h->p[0], l1.M->p[0], 2);
684 add(T, l1);
685
686 simplex l2(1);
687 Vector_Copy(h->p[1], l2.M->p[0], 2);
688 add(T, l2);
689
690 simplex l3(1);
691 Vector_Combine(h->p[0], h->p[1], l3.M->p[0],
692 tmp, tmp2, 2);
693 add(T, l3);
694
695 simplex t1(2);
696 Vector_Copy(h->p[0], t1.M->p[0], 2);
697 Vector_Copy(h->p[1], t1.M->p[1], 2);
698 add(T, t1);
699
700 simplex t2(2);
701 Vector_Combine(h->p[0], h->p[1], t2.M->p[0],
702 tmp, tmp2, 2);
703 Vector_Combine(h->p[2], h->p[1], t2.M->p[1],
704 tmp, tmp2, 2);
705 add(T, t2);
706 } else {
707 /* update h */
708 Vector_Combine(h->p[i], offset->p, h->p[i],
709 tmp, tmp, 2);
710 Vector_Combine(h->p[2], offset->p, h->p[2],
711 tmp, tmp, 2);
712 value_decrement(factor, factor);
713 }
714
715 simplex q(3, 0x4 | (1 << i), factor);
716 Vector_Copy(h->p[0], q.M->p[0], 2);
717 Vector_Copy(h->p[1], q.M->p[1], 2);
718 Vector_Copy(h->p[2], q.M->p[2], 2);
719 Vector_Copy(offset->p, q.offset->p, 2);
720 add(T, q);
721
722 simplex t1(2, 0x3, factor);
723 Vector_Copy(h->p[i], t1.M->p[0], 2);
724 Vector_Copy(h->p[2], t1.M->p[1], 2);
725 Vector_Copy(offset->p, t1.offset->p, 2);
726 add(T, t1);
727
728 simplex t2(2, 0x2, factor);
729 Vector_Copy(h->p[1-i], t2.M->p[0], 2);
730 Vector_Copy(h->p[2], t2.M->p[1], 2);
731 Vector_Copy(offset->p, t2.offset->p, 2);
732 add(T, t2);
733
734 simplex l(1, 0x1, factor);
735 Vector_Copy(h->p[2], l.M->p[0], 2);
736 Vector_Copy(offset->p, l.offset->p, 2);
737 add(T, l);
738
739 /* update h */
740 Vector_Combine(h->p[i], offset->p, h->p[i], tmp, factor, 2);
741 Vector_Combine(h->p[2], offset->p, h->p[2], tmp, factor, 2);
742
743 initial = false;
744 }
745 }
746 }
747 }
748 assert(!initial);
749
750 Vector_Free(offset);
751 Matrix_Free(h);
752
753 value_clear(tmp);
754 value_clear(tmp2);
755 value_clear(factor);
756 }
757
758 Matrix_Free(T);
759 }
760
761 void scarf_complex::print(FILE *out)
762 {
763 for (int i = 0; i < simplices.size(); ++i)
764 simplices[i].print(out);
765 }
766
767 struct scarf_collector {
768 virtual void add(Polyhedron *P, int sign, Polyhedron *C, unsigned MaxRays) = 0;
769 };
770
771 static void scarf(Polyhedron *P, unsigned exist, unsigned nparam, unsigned MaxRays,
772 scarf_collector& col)
773 {
774 Matrix *A, *B;
775 int dim = P->Dimension - exist - nparam;
776 assert(exist == 2);
777 int pos[4];
778 Polyhedron *U;
779 gen_fun *gf;
780 int n;
781
782 A = extract_matrix(P, dim);
783
784 n = A->NbColumns - 2;
785 assert(n >= 3 && n <= 4);
786
787 for (int i = 0; i < n; ++i)
788 pos[i] = i;
789 B = normalize_matrix(A, pos, n);
790
791 scarf_complex scarf;
792 scarf.add(B, pos, n);
793
794 U = Universe_Polyhedron(nparam);
795 col.add(P, 0, U, MaxRays);
796 for (int i = 0; i < scarf.simplices.size(); ++i) {
797 Polyhedron *Q;
798 int sign = (scarf.simplices[i].M->NbRows % 2) ? -1 : 1;
799 Q = scarf.simplices[i].shrunk_polyhedron(P, dim, A, MaxRays);
800 col.add(Q, sign, U, MaxRays);
801 Polyhedron_Free(Q);
802 }
803 Polyhedron_Free(U);
804
805 Matrix_Free(B);
806
807 Matrix_Free(A);
808 }
809
810 struct scarf_collector_gf : public scarf_collector {
811 QQ c;
812 gen_fun *gf;
813
814 scarf_collector_gf() {
815 c.d = 1;
816 }
817 virtual void add(Polyhedron *P, int sign, Polyhedron *C, unsigned MaxRays);
818 };
819
820 void scarf_collector_gf::add(Polyhedron *P, int sign, Polyhedron *C,
821 unsigned MaxRays)
822 {
823 if (!sign)
824 gf = barvinok_series(P, C, MaxRays);
825 else {
826 gen_fun *gf2;
827 c.n = sign;
828 gf2 = barvinok_series(P, C, MaxRays);
829 gf->add(c, gf2);
830 delete gf2;
831 }
832 }
833
834 gen_fun *barvinok_enumerate_scarf_series(Polyhedron *P,
835 unsigned exist, unsigned nparam, unsigned MaxRays)
836 {
837 scarf_collector_gf scgf;
838 scarf(P, exist, nparam, MaxRays, scgf);
839 return scgf.gf;
840 }
841
842 struct scarf_collector_ev : public scarf_collector {
843 evalue mone;
844 evalue *EP;
845
846 scarf_collector_ev() {
847 value_init(mone.d);
848 evalue_set_si(&mone, -1, 1);
849 }
850 ~scarf_collector_ev() {
851 free_evalue_refs(&mone);
852 }
853 virtual void add(Polyhedron *P, int sign, Polyhedron *C, unsigned MaxRays);
854 };
855
856 void scarf_collector_ev::add(Polyhedron *P, int sign, Polyhedron *C,
857 unsigned MaxRays)
858 {
859 if (!sign)
860 EP = barvinok_enumerate_ev(P, C, MaxRays);
861 else {
862 evalue *E2;
863 E2 = barvinok_enumerate_ev(P, C, MaxRays);
864 if (sign < 0)
865 emul(&mone, E2);
866 eadd(E2, EP);
867 free_evalue_refs(E2);
868 free(E2);
869 }
870 }
871
872 evalue *barvinok_enumerate_scarf(Polyhedron *P,
873 unsigned exist, unsigned nparam, unsigned MaxRays)
874 {
875 scarf_collector_ev scev;
876 scarf(P, exist, nparam, MaxRays, scev);
877 return scev.EP;
878 }
lattice point enumeration library