?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
add tests/ehrhart.README explaining origin of some of the ehrhart inputs
[barvinok.git] / scarf.cc
blob0585a98e192be54e722806a9460b87bcb680c8e7
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 Polyhedron *simplex::shrunk_polyhedron(Polyhedron *P, int dim, Matrix *A,
409 unsigned MaxRays)
410 {
411 Matrix *Constraints, *b;
412 Vector *b_offset = NULL;
413 Polyhedron *Q;
414 Value min, min_var, tmp;
415 value_init(tmp);
416 value_init(min);
417 value_init(min_var);
418 int constant;
419
420 b = Matrix_Alloc(M->NbRows, A->NbColumns);
421 Matrix_Product(M, A, b);
422
423 if (offset) {
424 b_offset = Vector_Alloc(A->NbColumns);
425 Vector_Matrix_Product(offset->p, A, b_offset->p);
426 }
427
428 if (!offset)
429 Constraints = Polyhedron2Constraints(P);
430 else {
431 Constraints = Matrix_Alloc(P->NbConstraints+2, P->Dimension+2+1);
432 for (int i = 0; i < P->NbConstraints; ++i) {
433 Vector_Copy(P->Constraint[i], Constraints->p[i], 1+dim+2);
434 Vector_Copy(P->Constraint[i]+1+dim+2, Constraints->p[i]+1+dim+2+1,
435 (P->Dimension+2)-(1+dim+2));
436 }
437 value_set_si(Constraints->p[P->NbConstraints][0], 1);
438 value_set_si(Constraints->p[P->NbConstraints][1+dim+2], 1);
439 value_set_si(Constraints->p[P->NbConstraints+1][0], 1);
440 value_set_si(Constraints->p[P->NbConstraints+1][1+dim+2], -1);
441 value_assign(Constraints->p[P->NbConstraints+1][Constraints->NbColumns-1],
442 last);
443 }
444 constant = Constraints->NbColumns - 1;
445
446 for (int i = 0, j = 0; i < P->NbConstraints; ++i) {
447 if (value_zero_p(Constraints->p[i][1+dim]) &&
448 value_zero_p(Constraints->p[i][1+dim+1]))
449 continue;
450 value_set_si(min, 0);
451 for (int k = 0; k < b->NbRows; ++k) {
452 if (offset && (mask & (1 << k)))
453 continue;
454 if (value_lt(b->p[k][j], min))
455 value_assign(min, b->p[k][j]);
456 }
457 value_set_si(min_var, 0);
458 if (offset) {
459 if (value_neg_p(b_offset->p[j])) {
460 value_oppose(min_var, b_offset->p[j]);
461 value_multiply(min_var, min_var, last);
462 value_increment(min_var, min_var);
463 }
464 for (int k = 0; k < b->NbRows; ++k) {
465 if (!(mask & (1 << k)))
466 continue;
467 if (value_lt(b->p[k][j], min_var))
468 value_assign(min_var, b->p[k][j]);
469 }
470 }
471 if (!offset || value_pos_p(b_offset->p[j])) {
472 if (value_le(min, min_var))
473 value_addto(Constraints->p[i][constant],
474 Constraints->p[i][constant], min);
475 else {
476 value_assign(tmp, min_var);
477 value_addmul(tmp, last, b_offset->p[j]);
478 if (value_le(tmp, min)) {
479 value_addto(Constraints->p[i][constant],
480 Constraints->p[i][constant], min_var);
481 value_addto(Constraints->p[i][1+dim+2],
482 Constraints->p[i][1+dim+2], b_offset->p[j]);
483 } else {
484 int lastrow = Constraints->NbRows;
485 int cols = Constraints->NbColumns;
486 Matrix *C = Constraints;
487 Constraints = AddANullRow(Constraints);
488 Matrix_Free(C);
489 Vector_Copy(Constraints->p[i], Constraints->p[lastrow], cols);
490 value_addto(Constraints->p[i][constant],
491 Constraints->p[i][constant], min_var);
492 value_addto(Constraints->p[i][1+dim+2],
493 Constraints->p[i][1+dim+2], b_offset->p[j]);
494 value_addto(Constraints->p[lastrow][constant],
495 Constraints->p[lastrow][constant], min);
496 }
497 }
498 } else {
499 if (value_le(min_var, min)) {
500 value_addto(Constraints->p[i][constant],
501 Constraints->p[i][constant], min_var);
502 value_addto(Constraints->p[i][1+dim+2],
503 Constraints->p[i][1+dim+2], b_offset->p[j]);
504 } else {
505 value_assign(tmp, min_var);
506 value_addmul(tmp, last, b_offset->p[j]);
507 if (value_le(min, tmp)) {
508 value_addto(Constraints->p[i][constant],
509 Constraints->p[i][constant], min);
510 } else {
511 int lastrow = Constraints->NbRows;
512 int cols = Constraints->NbColumns;
513 Matrix *C = Constraints;
514 Constraints = AddANullRow(Constraints);
515 Matrix_Free(C);
516 Vector_Copy(Constraints->p[i], Constraints->p[lastrow], cols);
517 value_addto(Constraints->p[i][constant],
518 Constraints->p[i][constant], min_var);
519 value_addto(Constraints->p[i][1+dim+2],
520 Constraints->p[i][1+dim+2], b_offset->p[j]);
521 value_addto(Constraints->p[lastrow][constant],
522 Constraints->p[lastrow][constant], min);
523 }
524 }
525 }
526 ++j;
527 }
528 Q = Constraints2Polyhedron(Constraints, MaxRays);
529
530 if (b_offset)
531 Vector_Free(b_offset);
532 Matrix_Free(b);
533 Matrix_Free(Constraints);
534 value_clear(tmp);
535 value_clear(min);
536 value_clear(min_var);
537
538 return Q;
539 }
540
541 struct scarf_complex {
542 vector<simplex> simplices;
543 void add(Matrix *B, int pos[4], int n);
544 void add(Matrix *T, simplex s);
545 void print(FILE *out);
546 ~scarf_complex() {
547 for (int i = 0; i < simplices.size(); ++i) {
548 Matrix_Free(simplices[i].M);
549 if (simplices[i].offset) {
550 Vector_Free(simplices[i].offset);
551 value_clear(simplices[i].last);
552 }
553 }
554 }
555 };
556
557 void scarf_complex::add(Matrix *T, simplex s)
558 {
559 s.transform(T);
560 s.normalize();
561 if (s.offset && lex_sign(s.offset->p, 2) < 0) {
562 Value factor;
563 Value tmp;
564 value_init(factor);
565 value_init(tmp);
566 /* compute the smallest multiple (factor) of the offset that
567 * makes on of the vertices lexico-negative.
568 */
569 int lexmin = -1;
570 for (int i = 0; i < s.M->NbRows; ++i) {
571 if (!(s.mask & (1 << i)))
572 continue;
573 if (lexmin == -1 || lex_smaller(s.M->p[i], s.M->p[lexmin], 2))
574 lexmin = i;
575 }
576 if (value_zero_p(s.offset->p[0])) {
577 if (value_pos_p(s.M->p[lexmin][0]))
578 value_increment(factor, s.last);
579 else {
580 value_oppose(factor, s.M->p[lexmin][1]);
581 mpz_cdiv_q(factor, factor, s.offset->p[1]);
582 }
583 } else {
584 value_oppose(factor, s.M->p[lexmin][0]);
585 mpz_cdiv_q(factor, factor, s.offset->p[0]);
586 if (mpz_divisible_p(s.M->p[lexmin][0], s.offset->p[0])) {
587 value_assign(tmp, s.M->p[lexmin][1]);
588 value_addmul(tmp, factor, s.offset->p[1]);
589 if (value_pos_p(tmp))
590 value_increment(factor, factor);
591 }
592 }
593 if (value_le(factor, s.last)) {
594 simplex part(s.M->NbRows, s.mask, s.last);
595 Vector_Copy(s.offset->p, part.offset->p, 2);
596 value_set_si(tmp, 1);
597 for (int i = 0; i < s.M->NbRows; ++i) {
598 if (s.mask & (1 << i))
599 Vector_Combine(s.M->p[i], s.offset->p, part.M->p[i],
600 tmp, factor, 2);
601 else
602 Vector_Copy(s.M->p[i], part.M->p[i], 2);
603 }
604 value_subtract(part.last, part.last, factor);
605 value_decrement(s.last, factor);
606 part.normalize();
607 simplices.push_back(part);
608 }
609 value_clear(tmp);
610 value_clear(factor);
611 }
612 simplices.push_back(s);
613 }
614
615 void scarf_complex::add(Matrix *B, int pos[4], int n)
616 {
617 Matrix *T;
618
619 T = Matrix_Alloc(2, 2);
620 Vector_Copy(B->p[0]+B->NbColumns-2, T->p[0], 2);
621 Vector_Copy(B->p[1]+B->NbColumns-2, T->p[1], 2);
622
623 if (n == 3 || value_neg_p(B->p[0][pos[3]])) {
624 assert(n == 3 || value_neg_p(B->p[1][pos[3]]));
625
626 simplex s1(1);
627 value_set_si(s1.M->p[0][0], 0);
628 value_set_si(s1.M->p[0][1], 1);
629 add(T, s1);
630
631 simplex s2(1);
632 value_set_si(s2.M->p[0][0], 1);
633 value_set_si(s2.M->p[0][1], 1);
634 add(T, s2);
635
636 simplex s3(1);
637 value_set_si(s3.M->p[0][0], 1);
638 value_set_si(s3.M->p[0][1], 0);
639 add(T, s3);
640
641 simplex s4(2);
642 value_set_si(s4.M->p[0][0], 0);
643 value_set_si(s4.M->p[0][1], 1);
644 value_set_si(s4.M->p[1][0], 1);
645 value_set_si(s4.M->p[1][1], 1);
646 add(T, s4);
647
648 simplex s5(2);
649 value_set_si(s5.M->p[0][0], 1);
650 value_set_si(s5.M->p[0][1], 0);
651 value_set_si(s5.M->p[1][0], 1);
652 value_set_si(s5.M->p[1][1], 1);
653 add(T, s5);
654 } else {
655 Matrix *h;
656 Vector *offset;
657 bool initial = true;
658 bool progress = true;
659 Value tmp, tmp2, factor;
660 int sign;
661
662 value_init(tmp);
663 value_init(tmp2);
664 value_init(factor);
665
666 assert(value_pos_p(B->p[0][pos[3]]));
667 assert(value_pos_p(B->p[1][pos[3]]));
668
669 h = Matrix_Alloc(3, 2);
670 value_set_si(h->p[0][0], 1);
671 value_set_si(h->p[0][1], 0);
672 value_set_si(h->p[1][0], 0);
673 value_set_si(h->p[1][1], 1);
674 value_set_si(h->p[2][0], 1);
675 value_set_si(h->p[2][1], 1);
676
677 offset = Vector_Alloc(2);
678
679 while (progress) {
680 progress = false;
681 for (int i = 0; i <= 1; ++i) {
682 value_set_si(factor, -1);
683 for (int j = 1; j <= 2; ++j) {
684 if (value_zero_p(B->p[1-i][pos[j]]))
685 continue;
686 value_oppose(tmp, B->p[i][pos[j]]);
687 value_pdivision(tmp, tmp, B->p[1-i][pos[j]]);
688 if (value_neg_p(factor) || value_lt(tmp, factor))
689 value_assign(factor, tmp);
690 }
691 if (value_pos_p(factor)) {
692 value_set_si(tmp, 1);
693 Vector_Combine(B->p[i], B->p[1-i], B->p[i], tmp, factor,
694 B->NbColumns);
695 sign = lex_sign(B->p[i], B->NbColumns);
696 for (int j = 1; j <= 2; ++j) {
697 if (value_notzero_p(B->p[i][pos[j]]))
698 continue;
699 /* a zero is interpreted to be of sign sign */
700 if ((sign > 0 && value_pos_p(B->p[1-i][pos[j]])) ||
701 (sign < 0 && value_neg_p(B->p[1-i][pos[j]]))) {
702 /* the zero is of the wrong sign => back-off one */
703 value_set_si(tmp2, -1);
704 Vector_Combine(B->p[i], B->p[1-i], B->p[i], tmp, tmp2,
705 B->NbColumns);
706 value_decrement(factor, factor);
707 }
708 }
709 /* We may have backed-off, so we need to check again. */
710 if (value_pos_p(factor)) {
711 progress = true;
712 value_set_si(tmp, 1);
713 value_set_si(tmp2, -1);
714
715 Vector_Combine(h->p[2], h->p[i], offset->p, tmp, tmp2, 2);
716
717 if (initial) {
718 /* the initial simplices not in any link */
719 simplex l1(1);
720 Vector_Copy(h->p[0], l1.M->p[0], 2);
721 add(T, l1);
722
723 simplex l2(1);
724 Vector_Copy(h->p[1], l2.M->p[0], 2);
725 add(T, l2);
726
727 simplex l3(1);
728 Vector_Combine(h->p[0], h->p[1], l3.M->p[0],
729 tmp, tmp2, 2);
730 add(T, l3);
731
732 simplex t1(2);
733 Vector_Copy(h->p[0], t1.M->p[0], 2);
734 Vector_Copy(h->p[1], t1.M->p[1], 2);
735 add(T, t1);
736
737 simplex t2(2);
738 Vector_Combine(h->p[0], h->p[1], t2.M->p[0],
739 tmp, tmp2, 2);
740 Vector_Combine(h->p[2], h->p[1], t2.M->p[1],
741 tmp, tmp2, 2);
742 add(T, t2);
743 } else {
744 /* update h */
745 Vector_Combine(h->p[i], offset->p, h->p[i],
746 tmp, tmp, 2);
747 Vector_Combine(h->p[2], offset->p, h->p[2],
748 tmp, tmp, 2);
749 value_decrement(factor, factor);
750 }
751
752 simplex q(3, 0x4 | (1 << i), factor);
753 Vector_Copy(h->p[0], q.M->p[0], 2);
754 Vector_Copy(h->p[1], q.M->p[1], 2);
755 Vector_Copy(h->p[2], q.M->p[2], 2);
756 Vector_Copy(offset->p, q.offset->p, 2);
757 add(T, q);
758
759 simplex t1(2, 0x3, factor);
760 Vector_Copy(h->p[i], t1.M->p[0], 2);
761 Vector_Copy(h->p[2], t1.M->p[1], 2);
762 Vector_Copy(offset->p, t1.offset->p, 2);
763 add(T, t1);
764
765 simplex t2(2, 0x2, factor);
766 Vector_Copy(h->p[1-i], t2.M->p[0], 2);
767 Vector_Copy(h->p[2], t2.M->p[1], 2);
768 Vector_Copy(offset->p, t2.offset->p, 2);
769 add(T, t2);
770
771 simplex l(1, 0x1, factor);
772 Vector_Copy(h->p[2], l.M->p[0], 2);
773 Vector_Copy(offset->p, l.offset->p, 2);
774 add(T, l);
775
776 /* update h */
777 Vector_Combine(h->p[i], offset->p, h->p[i], tmp, factor, 2);
778 Vector_Combine(h->p[2], offset->p, h->p[2], tmp, factor, 2);
779
780 initial = false;
781 }
782 }
783 }
784 }
785 if (initial) {
786 /* the initial simplices not in any link */
787 simplex l1(1);
788 Vector_Copy(h->p[0], l1.M->p[0], 2);
789 add(T, l1);
790
791 simplex l2(1);
792 Vector_Copy(h->p[1], l2.M->p[0], 2);
793 add(T, l2);
794
795 simplex l3(1);
796 Vector_Combine(h->p[0], h->p[1], l3.M->p[0],
797 tmp, tmp2, 2);
798 add(T, l3);
799
800 simplex t1(2);
801 Vector_Copy(h->p[0], t1.M->p[0], 2);
802 Vector_Copy(h->p[1], t1.M->p[1], 2);
803 add(T, t1);
804
805 simplex t2(2);
806 Vector_Combine(h->p[0], h->p[1], t2.M->p[0],
807 tmp, tmp2, 2);
808 Vector_Combine(h->p[2], h->p[1], t2.M->p[1],
809 tmp, tmp2, 2);
810 add(T, t2);
811
812 {
813 /* the simplices in a link, here of length 1 */
814 simplex q(3);
815 Vector_Copy(h->p[0], q.M->p[0], 2);
816 Vector_Copy(h->p[1], q.M->p[1], 2);
817 Vector_Copy(h->p[2], q.M->p[2], 2);
818 add(T, q);
819
820 simplex t1(2);
821 Vector_Copy(h->p[0], t1.M->p[0], 2);
822 Vector_Copy(h->p[2], t1.M->p[1], 2);
823 add(T, t1);
824
825 simplex t2(2);
826 Vector_Copy(h->p[1], t2.M->p[0], 2);
827 Vector_Copy(h->p[2], t2.M->p[1], 2);
828 add(T, t2);
829
830 simplex l(1);
831 Vector_Copy(h->p[2], l.M->p[0], 2);
832 add(T, l);
833 }
834 }
835
836 Vector_Free(offset);
837 Matrix_Free(h);
838
839 value_clear(tmp);
840 value_clear(tmp2);
841 value_clear(factor);
842 }
843
844 Matrix_Free(T);
845 }
846
847 void scarf_complex::print(FILE *out)
848 {
849 for (int i = 0; i < simplices.size(); ++i)
850 simplices[i].print(out);
851 }
852
853 struct scarf_collector {
854 virtual void add(Polyhedron *P, int sign, Polyhedron *C,
855 barvinok_options *options) = 0;
856 virtual ~scarf_collector() {}
857 };
858
859 static void scarf(Polyhedron *P, unsigned exist, unsigned nparam,
860 barvinok_options *options, scarf_collector& col)
861 {
862 Matrix *A, *B;
863 int dim = P->Dimension - exist - nparam;
864 assert(exist == 2);
865 int pos[4];
866 Polyhedron *U;
867 int n;
868
869 A = extract_matrix(P, dim);
870
871 n = A->NbColumns - 2;
872 assert(n >= 3 && n <= 4);
873
874 int l = 0;
875 for (int i = 0; i < n; ++i) {
876 int j;
877 for (j = 0; j < l; ++j)
878 if (value_eq(A->p[0][pos[j]], A->p[0][i]) &&
879 value_eq(A->p[1][pos[j]], A->p[1][i]))
880 break;
881 if (j < l)
882 continue;
883 pos[l++] = i;
884 }
885
886 assert(l >= 3 && l <= 4);
887 B = normalize_matrix(A, pos, &l);
888
889 scarf_complex scarf;
890 scarf.add(B, pos, l);
891
892 U = Universe_Polyhedron(nparam);
893 col.add(P, 0, U, options);
894 for (int i = 0; i < scarf.simplices.size(); ++i) {
895 Polyhedron *Q;
896 int sign = (scarf.simplices[i].M->NbRows % 2) ? -1 : 1;
897 Q = scarf.simplices[i].shrunk_polyhedron(P, dim, A, options->MaxRays);
898 col.add(Q, sign, U, options);
899 Polyhedron_Free(Q);
900 }
901 Polyhedron_Free(U);
902
903 Matrix_Free(B);
904
905 Matrix_Free(A);
906 }
907
908 struct scarf_collector_gf : public scarf_collector {
909 QQ c;
910 gen_fun *gf;
911
912 scarf_collector_gf() {
913 c.d = 1;
914 }
915 virtual void add(Polyhedron *P, int sign, Polyhedron *C,
916 barvinok_options *options);
917 };
918
919 void scarf_collector_gf::add(Polyhedron *P, int sign, Polyhedron *C,
920 barvinok_options *options)
921 {
922 if (!sign)
923 gf = barvinok_series_with_options(P, C, options);
924 else {
925 gen_fun *gf2;
926 c.n = sign;
927 gf2 = barvinok_series_with_options(P, C, options);
928 gf->add(c, gf2, options);
929 delete gf2;
930 }
931 }
932
933 gen_fun *barvinok_enumerate_scarf_series(Polyhedron *P,
934 unsigned exist, unsigned nparam, barvinok_options *options)
935 {
936 scarf_collector_gf scgf;
937 scarf(P, exist, nparam, options, scgf);
938 return scgf.gf;
939 }
940
941 struct scarf_collector_ev : public scarf_collector {
942 evalue mone;
943 evalue *EP;
944
945 scarf_collector_ev() {
946 value_init(mone.d);
947 evalue_set_si(&mone, -1, 1);
948 }
949 ~scarf_collector_ev() {
950 free_evalue_refs(&mone);
951 }
952 virtual void add(Polyhedron *P, int sign, Polyhedron *C,
953 barvinok_options *options);
954 };
955
956 void scarf_collector_ev::add(Polyhedron *P, int sign, Polyhedron *C,
957 barvinok_options *options)
958 {
959 if (!sign)
960 EP = barvinok_enumerate_with_options(P, C, options);
961 else {
962 evalue *E2;
963 E2 = barvinok_enumerate_with_options(P, C, options);
964 if (sign < 0)
965 emul(&mone, E2);
966 eadd(E2, EP);
967 evalue_free(E2);
968 }
969 }
970
971 evalue *barvinok_enumerate_scarf(Polyhedron *P,
972 unsigned exist, unsigned nparam, barvinok_options *options)
973 {
974 scarf_collector_ev scev;
975 scarf(P, exist, nparam, options, scev);
976 return scev.EP;
977 }
lattice point enumeration library