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