2 #include <barvinok/barvinok.h>
3 #include <barvinok/util.h>
8 static Matrix
*extract_matrix(Polyhedron
*P
, unsigned dim
)
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]))
19 A
= Matrix_Alloc(2, n_col
+2);
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]))
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]);
29 value_set_si(A
->p
[0][n_col
], 1);
30 value_set_si(A
->p
[1][n_col
+1], 1);
35 static int lex_sign(Value
*v
, int len
)
39 first
= First_Non_Zero(v
, len
);
40 return first
== -1 ? 0 : value_sign(v
[first
]);
43 static void set_pos(int pos
[4], int actual
, int wanted
)
48 pos
[actual
] = pos
[wanted
];
52 static Matrix
*normalize_matrix(Matrix
*A
, int pos
[4], int *n
)
55 Value tmp
, tmp2
, factor
;
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]);
71 B
= Matrix_Alloc(2, A
->NbColumns
);
72 Matrix_Product(T
, A
, B
);
75 /* Make zero in first position negative */
76 if (lex_sign(B
->p
[1], B
->NbColumns
) > 0) {
77 value_set_si(tmp
, -1);
78 Vector_Scale(B
->p
[1], B
->p
[1], tmp
, B
->NbColumns
);
81 /* First determine whether the matrix is of sign pattern I or II
85 assert(value_neg_p(B
->p
[1][pos
[1]]));
86 assert(value_pos_p(B
->p
[1][pos
[2]]));
88 value_set_si(factor
, 0);
89 for (int i
= 1; i
<= 2; ++i
) {
90 value_pdivision(tmp
, B
->p
[0][pos
[i
]], B
->p
[1][pos
[i
]]);
91 value_increment(tmp
, tmp
);
92 if (value_gt(tmp
, factor
))
93 value_assign(factor
, tmp
);
95 value_oppose(factor
, factor
);
97 Vector_Combine(B
->p
[0], B
->p
[1], B
->p
[0],
98 tmp
, factor
, B
->NbColumns
);
99 Vector_Exchange(B
->p
[0], B
->p
[1], B
->NbColumns
);
100 /* problems with three constraints are considered
101 * to be of sign pattern II
106 for (i
= 1; i
<= 3; ++i
)
107 if (value_zero_p(B
->p
[1][pos
[i
]]))
110 /* put zero in position 3 */
113 /* put positive one in position 1 */
114 for (i
= 1; i
<= 3; ++i
)
115 if (value_pos_p(B
->p
[1][pos
[i
]]))
119 value_set_si(factor
, 0);
120 for (int i
= 1; i
<= 2; ++i
) {
121 value_pdivision(tmp
, B
->p
[0][pos
[i
]], B
->p
[1][pos
[i
]]);
122 value_increment(tmp
, tmp
);
123 if (value_gt(tmp
, factor
))
124 value_assign(factor
, tmp
);
126 value_oppose(factor
, factor
);
127 value_set_si(tmp
, 1);
128 Vector_Combine(B
->p
[0], B
->p
[1], B
->p
[0], tmp
, factor
, B
->NbColumns
);
130 assert(value_notzero_p(B
->p
[0][pos
[3]]));
131 type
= value_pos_p(B
->p
[0][pos
[3]]) ? 1 : 2;
134 int sign
= lex_sign(B
->p
[1], B
->NbColumns
);
136 for (int i
= 1; i
<= 3; ++i
)
137 if (value_neg_p(B
->p
[1][pos
[i
]]))
139 assert(neg
== 1 || neg
== 2);
142 /* put negative one in position 1 */
143 for (i
= 1; i
<= 3; ++i
)
144 if (value_neg_p(B
->p
[1][pos
[i
]]))
148 value_set_si(factor
, 0);
149 for (int i
= 1; i
<= 3; ++i
) {
150 value_pdivision(tmp
, B
->p
[0][pos
[i
]], B
->p
[1][pos
[i
]]);
151 value_increment(tmp
, tmp
);
152 if (value_gt(tmp
, factor
))
153 value_assign(factor
, tmp
);
155 value_oppose(factor
, factor
);
156 value_set_si(tmp
, 1);
157 Vector_Combine(B
->p
[0], B
->p
[1], B
->p
[0],
158 tmp
, factor
, B
->NbColumns
);
159 Vector_Exchange(B
->p
[0], B
->p
[1], B
->NbColumns
);
163 /* put positive one in position 1 */
164 for (i
= 1; i
<= 3; ++i
)
165 if (value_pos_p(B
->p
[1][pos
[i
]]))
169 value_set_si(factor
, 0);
170 for (int i
= 1; i
<= 3; ++i
) {
171 value_pdivision(tmp
, B
->p
[0][pos
[i
]], B
->p
[1][pos
[i
]]);
172 value_increment(tmp
, tmp
);
173 if (value_gt(tmp
, factor
))
174 value_assign(factor
, tmp
);
176 value_oppose(factor
, factor
);
177 value_set_si(tmp
, 1);
178 Vector_Combine(B
->p
[0], B
->p
[1], B
->p
[0],
179 tmp
, factor
, B
->NbColumns
);
189 value_oppose(tmp
, B
->p
[0][pos
[1]]);
190 value_pdivision(factor
, tmp
, B
->p
[1][pos
[1]]);
191 value_oppose(tmp
, B
->p
[1][pos
[2]]);
192 value_pdivision(tmp
, tmp
, B
->p
[0][pos
[2]]);
193 if (value_zero_p(factor
) && value_zero_p(tmp
))
195 assert(value_zero_p(factor
) || value_zero_p(tmp
));
196 if (value_pos_p(factor
)) {
197 value_set_si(tmp
, 1);
198 Vector_Combine(B
->p
[0], B
->p
[1], B
->p
[0], tmp
, factor
, B
->NbColumns
);
199 if (value_zero_p(B
->p
[0][pos
[1]])) {
200 /* We will deal with this later */
201 assert(lex_sign(B
->p
[0], B
->NbColumns
) < 0);
204 value_set_si(factor
, 1);
205 Vector_Combine(B
->p
[0], B
->p
[1], B
->p
[1], tmp
, factor
, B
->NbColumns
);
206 if (value_zero_p(B
->p
[1][pos
[2]])) {
207 /* We will deal with this later */
208 assert(lex_sign(B
->p
[1], B
->NbColumns
) < 0);
215 bool progress
= true;
218 for (int i
= 0; i
<= 1; ++i
) {
219 value_set_si(factor
, -1);
220 for (int j
= 1; j
<= 3; ++j
) {
221 if (value_zero_p(B
->p
[1-i
][pos
[j
]]))
223 value_oppose(tmp
, B
->p
[i
][pos
[j
]]);
224 value_pdivision(tmp
, tmp
, B
->p
[1-i
][pos
[j
]]);
225 if (value_neg_p(factor
) || value_lt(tmp
, factor
))
226 value_assign(factor
, tmp
);
228 if (value_pos_p(factor
)) {
229 value_set_si(tmp
, 1);
230 Vector_Combine(B
->p
[i
], B
->p
[1-i
], B
->p
[i
], tmp
, factor
,
232 sign
= lex_sign(B
->p
[i
], B
->NbColumns
);
233 for (int j
= 1; j
<= 3; ++j
) {
234 if (value_notzero_p(B
->p
[i
][pos
[j
]]))
236 /* a zero is interpreted to be of sign sign */
237 if ((sign
> 0 && value_pos_p(B
->p
[1-i
][pos
[j
]])) ||
238 (sign
< 0 && value_neg_p(B
->p
[1-i
][pos
[j
]]))) {
239 /* the zero is of the wrong sign => back-off one */
240 value_set_si(tmp2
, -1);
241 Vector_Combine(B
->p
[i
], B
->p
[1-i
], B
->p
[i
], tmp
, tmp2
,
243 value_decrement(factor
, factor
);
246 /* We may have backed-off, so we need to check again. */
247 if (value_pos_p(factor
))
253 for (int i
= 0; i
< B
->NbColumns
; ++i
) {
254 value_addto(tmp
, B
->p
[0][i
], B
->p
[1][i
]);
255 if (value_zero_p(tmp
))
257 sign
= value_neg_p(tmp
) ? -1 : 1;
261 for (int i
= 1; i
<= 3; ++i
) {
262 value_addto(tmp
, B
->p
[0][pos
[i
]], B
->p
[1][pos
[i
]]);
263 if (value_neg_p(tmp
) || (sign
< 0 && value_zero_p(tmp
)))
270 /* cases 4 and 5 in Theorem 11.1 */
271 value_set_si(tmp
, 1);
272 Vector_Combine(B
->p
[0], B
->p
[1], B
->p
[1], tmp
, tmp
, B
->NbColumns
);
274 /* put positive pair in position 3 */
275 for (i
= 1; i
<= 3; ++i
)
276 if (value_pos_p(B
->p
[0][pos
[i
]]) && value_pos_p(B
->p
[1][pos
[i
]]))
283 /* cases 1 and 2 in Theorem 11.1 */
284 value_set_si(tmp
, 1);
285 Vector_Combine(B
->p
[0], B
->p
[1], B
->p
[0], tmp
, tmp
, B
->NbColumns
);
287 /* put positive one in position 2 */
288 for (i
= 1; i
<= 3; ++i
)
289 if (value_pos_p(B
->p
[0][pos
[i
]]))
294 /* fourth constraint is redundant with respect to neighborhoods */
298 /* We will deal with these later */
311 Value last
; // last multiple of offset in link
313 Matrix
*M
; // rows: elements different from (0,0)
317 M
= Matrix_Alloc(d
, 2);
320 simplex(int d
, int mask
, Value last
) {
321 M
= Matrix_Alloc(d
, 2);
322 offset
= Vector_Alloc(2);
323 value_init(this->last
);
324 value_assign(this->last
, last
);
327 void transform(Matrix
*T
);
329 Polyhedron
*shrunk_polyhedron(Polyhedron
*P
, int dim
, Matrix
*A
,
331 void print(FILE *out
);
334 void simplex::print(FILE *out
)
337 Matrix_Print(out
, P_VALUE_FMT
, M
);
339 fprintf(out
, "%d %d\n", M
->NbRows
, M
->NbColumns
);
340 for (int j
= 0; j
< M
->NbRows
; ++j
) {
341 for (int k
= 0; k
< M
->NbColumns
; ++k
)
342 value_print(out
, P_VALUE_FMT
, M
->p
[j
][k
]);
343 if (mask
& (1 << j
)) {
344 fprintf(out
, " + k * ");
345 for (int k
= 0; k
< M
->NbColumns
; ++k
)
346 value_print(out
, P_VALUE_FMT
, offset
->p
[k
]);
350 fprintf(out
, "\t0 <= k <= ");
351 value_print(out
, P_VALUE_FMT
, last
);
356 static bool lex_smaller(Value
*v1
, Value
*v2
, int n
)
358 for (int i
= 0; i
< n
; ++i
)
359 if (value_lt(v1
[i
], v2
[i
]))
361 else if (value_gt(v1
[i
], v2
[i
]))
366 void simplex::transform(Matrix
*T
)
369 M
= Matrix_Alloc(M2
->NbRows
, M2
->NbColumns
);
370 Matrix_Product(M2
, T
, M
);
374 Vector
*offset2
= offset
;
375 offset
= Vector_Alloc(offset2
->Size
);
376 Vector_Matrix_Product(offset2
->p
, T
, offset
->p
);
377 Vector_Free(offset2
);
381 void simplex::normalize()
384 for (int i
= 1; i
< M
->NbRows
; ++i
)
385 if (lex_smaller(M
->p
[i
], M
->p
[lexmin
], 2))
387 if (lex_sign(M
->p
[lexmin
], 2) < 0) {
390 value_set_si(tmp
, -1);
391 Vector_Scale(M
->p
[lexmin
], M
->p
[lexmin
], tmp
, 2);
392 value_set_si(tmp
, 1);
393 for (int i
= 0; i
< M
->NbRows
; ++i
) {
396 Vector_Combine(M
->p
[lexmin
], M
->p
[i
], M
->p
[i
], tmp
, tmp
, 2);
398 if (offset
&& (mask
& (1 << lexmin
))) {
399 value_set_si(tmp
, -1);
400 Vector_Scale(offset
->p
, offset
->p
, tmp
, 2);
401 mask
^= (1 << M
->NbRows
) - 1 - (1 << lexmin
);
407 static Matrix
*InsertColumn(Matrix
*M
, int pos
)
412 R
= Matrix_Alloc(M
->NbRows
, M
->NbColumns
+1);
413 for (i
= 0; i
< M
->NbRows
; ++i
) {
414 Vector_Copy(M
->p
[i
], R
->p
[i
], pos
);
415 Vector_Copy(M
->p
[i
]+pos
, R
->p
[i
]+pos
+1, M
->NbColumns
-pos
);
420 Polyhedron
*simplex::shrunk_polyhedron(Polyhedron
*P
, int dim
, Matrix
*A
,
423 Matrix
*Constraints
, *b
;
424 Vector
*b_offset
= NULL
;
426 Value min
, min_var
, tmp
;
432 b
= Matrix_Alloc(M
->NbRows
, A
->NbColumns
);
433 Matrix_Product(M
, A
, b
);
436 b_offset
= Vector_Alloc(A
->NbColumns
);
437 Vector_Matrix_Product(offset
->p
, A
, b_offset
->p
);
441 Constraints
= Polyhedron2Constraints(P
);
443 Constraints
= Matrix_Alloc(P
->NbConstraints
+2, P
->Dimension
+2+1);
444 for (int i
= 0; i
< P
->NbConstraints
; ++i
) {
445 Vector_Copy(P
->Constraint
[i
], Constraints
->p
[i
], 1+dim
+2);
446 Vector_Copy(P
->Constraint
[i
]+1+dim
+2, Constraints
->p
[i
]+1+dim
+2+1,
447 (P
->Dimension
+2)-(1+dim
+2));
449 value_set_si(Constraints
->p
[P
->NbConstraints
][0], 1);
450 value_set_si(Constraints
->p
[P
->NbConstraints
][1+dim
+2], 1);
451 value_set_si(Constraints
->p
[P
->NbConstraints
+1][0], 1);
452 value_set_si(Constraints
->p
[P
->NbConstraints
+1][1+dim
+2], -1);
453 value_assign(Constraints
->p
[P
->NbConstraints
+1][Constraints
->NbColumns
-1],
456 constant
= Constraints
->NbColumns
- 1;
458 for (int i
= 0, j
= 0; i
< P
->NbConstraints
; ++i
) {
459 if (value_zero_p(Constraints
->p
[i
][1+dim
]) &&
460 value_zero_p(Constraints
->p
[i
][1+dim
+1]))
462 value_set_si(min
, 0);
463 for (int k
= 0; k
< b
->NbRows
; ++k
) {
464 if (offset
&& (mask
& (1 << k
)))
466 if (value_lt(b
->p
[k
][j
], min
))
467 value_assign(min
, b
->p
[k
][j
]);
469 value_set_si(min_var
, 0);
471 if (value_neg_p(b_offset
->p
[j
])) {
472 value_oppose(min_var
, b_offset
->p
[j
]);
473 value_multiply(min_var
, min_var
, last
);
474 value_increment(min_var
, min_var
);
476 for (int k
= 0; k
< b
->NbRows
; ++k
) {
477 if (!(mask
& (1 << k
)))
479 if (value_lt(b
->p
[k
][j
], min_var
))
480 value_assign(min_var
, b
->p
[k
][j
]);
483 if (!offset
|| value_pos_p(b_offset
->p
[j
])) {
484 if (value_le(min
, min_var
))
485 value_addto(Constraints
->p
[i
][constant
],
486 Constraints
->p
[i
][constant
], min
);
488 value_assign(tmp
, min_var
);
489 value_addmul(tmp
, last
, b_offset
->p
[j
]);
490 if (value_le(tmp
, min
)) {
491 value_addto(Constraints
->p
[i
][constant
],
492 Constraints
->p
[i
][constant
], min_var
);
493 value_addto(Constraints
->p
[i
][1+dim
+2],
494 Constraints
->p
[i
][1+dim
+2], b_offset
->p
[j
]);
496 int lastrow
= Constraints
->NbRows
;
497 int cols
= Constraints
->NbColumns
;
498 Matrix
*C
= Constraints
;
499 Constraints
= AddANullRow(Constraints
);
501 Vector_Copy(Constraints
->p
[i
], Constraints
->p
[lastrow
], cols
);
502 value_addto(Constraints
->p
[i
][constant
],
503 Constraints
->p
[i
][constant
], min_var
);
504 value_addto(Constraints
->p
[i
][1+dim
+2],
505 Constraints
->p
[i
][1+dim
+2], b_offset
->p
[j
]);
506 value_addto(Constraints
->p
[lastrow
][constant
],
507 Constraints
->p
[lastrow
][constant
], min
);
511 if (value_le(min_var
, min
)) {
512 value_addto(Constraints
->p
[i
][constant
],
513 Constraints
->p
[i
][constant
], min_var
);
514 value_addto(Constraints
->p
[i
][1+dim
+2],
515 Constraints
->p
[i
][1+dim
+2], b_offset
->p
[j
]);
517 value_assign(tmp
, min_var
);
518 value_addmul(tmp
, last
, b_offset
->p
[j
]);
519 if (value_le(min
, tmp
)) {
520 value_addto(Constraints
->p
[i
][constant
],
521 Constraints
->p
[i
][constant
], min
);
523 int lastrow
= Constraints
->NbRows
;
524 int cols
= Constraints
->NbColumns
;
525 Matrix
*C
= Constraints
;
526 Constraints
= AddANullRow(Constraints
);
528 Vector_Copy(Constraints
->p
[i
], Constraints
->p
[lastrow
], cols
);
529 value_addto(Constraints
->p
[i
][constant
],
530 Constraints
->p
[i
][constant
], min_var
);
531 value_addto(Constraints
->p
[i
][1+dim
+2],
532 Constraints
->p
[i
][1+dim
+2], b_offset
->p
[j
]);
533 value_addto(Constraints
->p
[lastrow
][constant
],
534 Constraints
->p
[lastrow
][constant
], min
);
540 Q
= Constraints2Polyhedron(Constraints
, MaxRays
);
543 Vector_Free(b_offset
);
545 Matrix_Free(Constraints
);
548 value_clear(min_var
);
553 struct scarf_complex
{
554 vector
<simplex
> simplices
;
555 void add(Matrix
*B
, int pos
[4], int n
);
556 void add(Matrix
*T
, simplex s
);
557 void print(FILE *out
);
559 for (int i
= 0; i
< simplices
.size(); ++i
) {
560 Matrix_Free(simplices
[i
].M
);
561 if (simplices
[i
].offset
) {
562 Vector_Free(simplices
[i
].offset
);
563 value_clear(simplices
[i
].last
);
569 void scarf_complex::add(Matrix
*T
, simplex s
)
573 if (s
.offset
&& lex_sign(s
.offset
->p
, 2) < 0) {
578 /* compute the smallest multiple (factor) of the offset that
579 * makes on of the vertices lexico-negative.
582 for (int i
= 0; i
< s
.M
->NbRows
; ++i
) {
583 if (!(s
.mask
& (1 << i
)))
585 if (lexmin
== -1 || lex_smaller(s
.M
->p
[i
], s
.M
->p
[lexmin
], 2))
588 if (value_zero_p(s
.offset
->p
[0])) {
589 if (value_pos_p(s
.M
->p
[lexmin
][0]))
590 value_increment(factor
, s
.last
);
592 value_oppose(factor
, s
.M
->p
[lexmin
][1]);
593 mpz_cdiv_q(factor
, factor
, s
.offset
->p
[1]);
596 value_oppose(factor
, s
.M
->p
[lexmin
][0]);
597 mpz_cdiv_q(factor
, factor
, s
.offset
->p
[0]);
598 if (mpz_divisible_p(s
.M
->p
[lexmin
][0], s
.offset
->p
[0])) {
599 value_assign(tmp
, s
.M
->p
[lexmin
][1]);
600 value_addmul(tmp
, factor
, s
.offset
->p
[1]);
601 if (value_pos_p(tmp
))
602 value_increment(factor
, factor
);
605 if (value_le(factor
, s
.last
)) {
606 simplex
part(s
.M
->NbRows
, s
.mask
, s
.last
);
607 Vector_Copy(s
.offset
->p
, part
.offset
->p
, 2);
608 value_set_si(tmp
, 1);
609 for (int i
= 0; i
< s
.M
->NbRows
; ++i
) {
610 if (s
.mask
& (1 << i
))
611 Vector_Combine(s
.M
->p
[i
], s
.offset
->p
, part
.M
->p
[i
],
614 Vector_Copy(s
.M
->p
[i
], part
.M
->p
[i
], 2);
616 value_subtract(part
.last
, part
.last
, factor
);
617 value_decrement(s
.last
, factor
);
619 simplices
.push_back(part
);
624 simplices
.push_back(s
);
627 void scarf_complex::add(Matrix
*B
, int pos
[4], int n
)
631 T
= Matrix_Alloc(2, 2);
632 Vector_Copy(B
->p
[0]+B
->NbColumns
-2, T
->p
[0], 2);
633 Vector_Copy(B
->p
[1]+B
->NbColumns
-2, T
->p
[1], 2);
635 if (n
== 3 || value_neg_p(B
->p
[0][pos
[3]])) {
636 assert(n
== 3 || value_neg_p(B
->p
[1][pos
[3]]));
639 value_set_si(s1
.M
->p
[0][0], 0);
640 value_set_si(s1
.M
->p
[0][1], 1);
644 value_set_si(s2
.M
->p
[0][0], 1);
645 value_set_si(s2
.M
->p
[0][1], 1);
649 value_set_si(s3
.M
->p
[0][0], 1);
650 value_set_si(s3
.M
->p
[0][1], 0);
654 value_set_si(s4
.M
->p
[0][0], 0);
655 value_set_si(s4
.M
->p
[0][1], 1);
656 value_set_si(s4
.M
->p
[1][0], 1);
657 value_set_si(s4
.M
->p
[1][1], 1);
661 value_set_si(s5
.M
->p
[0][0], 1);
662 value_set_si(s5
.M
->p
[0][1], 0);
663 value_set_si(s5
.M
->p
[1][0], 1);
664 value_set_si(s5
.M
->p
[1][1], 1);
670 bool progress
= true;
671 Value tmp
, tmp2
, factor
;
678 assert(value_pos_p(B
->p
[0][pos
[3]]));
679 assert(value_pos_p(B
->p
[1][pos
[3]]));
681 h
= Matrix_Alloc(3, 2);
682 value_set_si(h
->p
[0][0], 1);
683 value_set_si(h
->p
[0][1], 0);
684 value_set_si(h
->p
[1][0], 0);
685 value_set_si(h
->p
[1][1], 1);
686 value_set_si(h
->p
[2][0], 1);
687 value_set_si(h
->p
[2][1], 1);
689 offset
= Vector_Alloc(2);
693 for (int i
= 0; i
<= 1; ++i
) {
694 value_set_si(factor
, -1);
695 for (int j
= 1; j
<= 2; ++j
) {
696 if (value_zero_p(B
->p
[1-i
][pos
[j
]]))
698 value_oppose(tmp
, B
->p
[i
][pos
[j
]]);
699 value_pdivision(tmp
, tmp
, B
->p
[1-i
][pos
[j
]]);
700 if (value_neg_p(factor
) || value_lt(tmp
, factor
))
701 value_assign(factor
, tmp
);
703 if (value_pos_p(factor
)) {
704 value_set_si(tmp
, 1);
705 Vector_Combine(B
->p
[i
], B
->p
[1-i
], B
->p
[i
], tmp
, factor
,
707 sign
= lex_sign(B
->p
[i
], B
->NbColumns
);
708 for (int j
= 1; j
<= 2; ++j
) {
709 if (value_notzero_p(B
->p
[i
][pos
[j
]]))
711 /* a zero is interpreted to be of sign sign */
712 if ((sign
> 0 && value_pos_p(B
->p
[1-i
][pos
[j
]])) ||
713 (sign
< 0 && value_neg_p(B
->p
[1-i
][pos
[j
]]))) {
714 /* the zero is of the wrong sign => back-off one */
715 value_set_si(tmp2
, -1);
716 Vector_Combine(B
->p
[i
], B
->p
[1-i
], B
->p
[i
], tmp
, tmp2
,
718 value_decrement(factor
, factor
);
721 /* We may have backed-off, so we need to check again. */
722 if (value_pos_p(factor
)) {
724 value_set_si(tmp
, 1);
725 value_set_si(tmp2
, -1);
727 Vector_Combine(h
->p
[2], h
->p
[i
], offset
->p
, tmp
, tmp2
, 2);
730 /* the initial simplices not in any link */
732 Vector_Copy(h
->p
[0], l1
.M
->p
[0], 2);
736 Vector_Copy(h
->p
[1], l2
.M
->p
[0], 2);
740 Vector_Combine(h
->p
[0], h
->p
[1], l3
.M
->p
[0],
745 Vector_Copy(h
->p
[0], t1
.M
->p
[0], 2);
746 Vector_Copy(h
->p
[1], t1
.M
->p
[1], 2);
750 Vector_Combine(h
->p
[0], h
->p
[1], t2
.M
->p
[0],
752 Vector_Combine(h
->p
[2], h
->p
[1], t2
.M
->p
[1],
757 Vector_Combine(h
->p
[i
], offset
->p
, h
->p
[i
],
759 Vector_Combine(h
->p
[2], offset
->p
, h
->p
[2],
761 value_decrement(factor
, factor
);
764 simplex
q(3, 0x4 | (1 << i
), factor
);
765 Vector_Copy(h
->p
[0], q
.M
->p
[0], 2);
766 Vector_Copy(h
->p
[1], q
.M
->p
[1], 2);
767 Vector_Copy(h
->p
[2], q
.M
->p
[2], 2);
768 Vector_Copy(offset
->p
, q
.offset
->p
, 2);
771 simplex
t1(2, 0x3, factor
);
772 Vector_Copy(h
->p
[i
], t1
.M
->p
[0], 2);
773 Vector_Copy(h
->p
[2], t1
.M
->p
[1], 2);
774 Vector_Copy(offset
->p
, t1
.offset
->p
, 2);
777 simplex
t2(2, 0x2, factor
);
778 Vector_Copy(h
->p
[1-i
], t2
.M
->p
[0], 2);
779 Vector_Copy(h
->p
[2], t2
.M
->p
[1], 2);
780 Vector_Copy(offset
->p
, t2
.offset
->p
, 2);
783 simplex
l(1, 0x1, factor
);
784 Vector_Copy(h
->p
[2], l
.M
->p
[0], 2);
785 Vector_Copy(offset
->p
, l
.offset
->p
, 2);
789 Vector_Combine(h
->p
[i
], offset
->p
, h
->p
[i
], tmp
, factor
, 2);
790 Vector_Combine(h
->p
[2], offset
->p
, h
->p
[2], tmp
, factor
, 2);
798 /* the initial simplices not in any link */
800 Vector_Copy(h
->p
[0], l1
.M
->p
[0], 2);
804 Vector_Copy(h
->p
[1], l2
.M
->p
[0], 2);
808 Vector_Combine(h
->p
[0], h
->p
[1], l3
.M
->p
[0],
813 Vector_Copy(h
->p
[0], t1
.M
->p
[0], 2);
814 Vector_Copy(h
->p
[1], t1
.M
->p
[1], 2);
818 Vector_Combine(h
->p
[0], h
->p
[1], t2
.M
->p
[0],
820 Vector_Combine(h
->p
[2], h
->p
[1], t2
.M
->p
[1],
825 /* the simplices in a link, here of length 1 */
827 Vector_Copy(h
->p
[0], q
.M
->p
[0], 2);
828 Vector_Copy(h
->p
[1], q
.M
->p
[1], 2);
829 Vector_Copy(h
->p
[2], q
.M
->p
[2], 2);
833 Vector_Copy(h
->p
[0], t1
.M
->p
[0], 2);
834 Vector_Copy(h
->p
[2], t1
.M
->p
[1], 2);
838 Vector_Copy(h
->p
[1], t2
.M
->p
[0], 2);
839 Vector_Copy(h
->p
[2], t2
.M
->p
[1], 2);
843 Vector_Copy(h
->p
[2], l
.M
->p
[0], 2);
859 void scarf_complex::print(FILE *out
)
861 for (int i
= 0; i
< simplices
.size(); ++i
)
862 simplices
[i
].print(out
);
865 struct scarf_collector
{
866 virtual void add(Polyhedron
*P
, int sign
, Polyhedron
*C
, unsigned MaxRays
) = 0;
869 static void scarf(Polyhedron
*P
, unsigned exist
, unsigned nparam
, unsigned MaxRays
,
870 scarf_collector
& col
)
873 int dim
= P
->Dimension
- exist
- nparam
;
880 A
= extract_matrix(P
, dim
);
882 n
= A
->NbColumns
- 2;
883 assert(n
>= 3 && n
<= 4);
886 for (int i
= 0; i
< n
; ++i
) {
888 for (j
= 0; j
< l
; ++j
)
889 if (value_eq(A
->p
[0][pos
[j
]], A
->p
[0][i
]) &&
890 value_eq(A
->p
[1][pos
[j
]], A
->p
[1][i
]))
897 assert(l
>= 3 && l
<= 4);
898 B
= normalize_matrix(A
, pos
, &l
);
901 scarf
.add(B
, pos
, l
);
903 U
= Universe_Polyhedron(nparam
);
904 col
.add(P
, 0, U
, MaxRays
);
905 for (int i
= 0; i
< scarf
.simplices
.size(); ++i
) {
907 int sign
= (scarf
.simplices
[i
].M
->NbRows
% 2) ? -1 : 1;
908 Q
= scarf
.simplices
[i
].shrunk_polyhedron(P
, dim
, A
, MaxRays
);
909 col
.add(Q
, sign
, U
, MaxRays
);
919 struct scarf_collector_gf
: public scarf_collector
{
923 scarf_collector_gf() {
926 virtual void add(Polyhedron
*P
, int sign
, Polyhedron
*C
, unsigned MaxRays
);
929 void scarf_collector_gf::add(Polyhedron
*P
, int sign
, Polyhedron
*C
,
933 gf
= barvinok_series(P
, C
, MaxRays
);
937 gf2
= barvinok_series(P
, C
, MaxRays
);
943 gen_fun
*barvinok_enumerate_scarf_series(Polyhedron
*P
,
944 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
946 scarf_collector_gf scgf
;
947 scarf(P
, exist
, nparam
, MaxRays
, scgf
);
951 struct scarf_collector_ev
: public scarf_collector
{
955 scarf_collector_ev() {
957 evalue_set_si(&mone
, -1, 1);
959 ~scarf_collector_ev() {
960 free_evalue_refs(&mone
);
962 virtual void add(Polyhedron
*P
, int sign
, Polyhedron
*C
, unsigned MaxRays
);
965 void scarf_collector_ev::add(Polyhedron
*P
, int sign
, Polyhedron
*C
,
969 EP
= barvinok_enumerate_ev(P
, C
, MaxRays
);
972 E2
= barvinok_enumerate_ev(P
, C
, MaxRays
);
976 free_evalue_refs(E2
);
981 evalue
*barvinok_enumerate_scarf(Polyhedron
*P
,
982 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
984 scarf_collector_ev scev
;
985 scarf(P
, exist
, nparam
, MaxRays
, scev
);