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