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 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
);
81 assert(value_neg_p(B
->p
[1][pos
[1]]));
82 assert(value_pos_p(B
->p
[1][pos
[2]]));
83 assert(value_pos_p(B
->p
[0][pos
[1]]));
84 assert(value_neg_p(B
->p
[0][pos
[2]]));
85 Vector_Exchange(B
->p
[0], B
->p
[1], B
->NbColumns
);
89 for (i
= 1; i
<= 3; ++i
)
90 if (value_zero_p(B
->p
[1][pos
[i
]]))
93 /* put zero in position 3 */
96 /* put positive one in position 1 */
97 for (i
= 1; i
<= 3; ++i
)
98 if (value_pos_p(B
->p
[1][pos
[i
]]))
102 value_set_si(factor
, 0);
103 for (int i
= 1; i
<= 2; ++i
) {
104 value_pdivision(tmp
, B
->p
[0][pos
[i
]], B
->p
[1][pos
[i
]]);
105 value_increment(tmp
, tmp
);
106 if (value_gt(tmp
, factor
))
107 value_assign(factor
, tmp
);
109 value_oppose(factor
, factor
);
110 value_set_si(tmp
, 1);
111 Vector_Combine(B
->p
[0], B
->p
[1], B
->p
[0], tmp
, factor
, B
->NbColumns
);
113 assert(value_notzero_p(B
->p
[0][pos
[3]]));
114 type
= value_pos_p(B
->p
[0][pos
[3]]) ? 1 : 2;
117 for (int i
= 1; i
<= 3; ++i
)
118 if (value_neg_p(B
->p
[1][pos
[i
]]))
120 assert(neg
== 1 || neg
== 2);
122 /* We will deal with this later */
126 /* put positive one in position 1 */
127 for (i
= 1; i
<= 3; ++i
)
128 if (value_pos_p(B
->p
[1][pos
[i
]]))
132 value_set_si(factor
, 0);
133 for (int i
= 1; i
<= 3; ++i
) {
134 value_pdivision(tmp
, B
->p
[0][pos
[i
]], B
->p
[1][pos
[i
]]);
135 value_increment(tmp
, tmp
);
136 if (value_gt(tmp
, factor
))
137 value_assign(factor
, tmp
);
139 value_oppose(factor
, factor
);
140 value_set_si(tmp
, 1);
141 Vector_Combine(B
->p
[0], B
->p
[1], B
->p
[0],
142 tmp
, factor
, B
->NbColumns
);
152 value_oppose(tmp
, B
->p
[0][pos
[1]]);
153 value_pdivision(factor
, tmp
, B
->p
[1][pos
[1]]);
154 value_oppose(tmp
, B
->p
[1][pos
[2]]);
155 value_pdivision(tmp
, tmp
, B
->p
[0][pos
[2]]);
156 if (value_zero_p(factor
) && value_zero_p(tmp
))
158 assert(value_zero_p(factor
) || value_zero_p(tmp
));
159 if (value_pos_p(factor
)) {
160 value_set_si(tmp
, 1);
161 Vector_Combine(B
->p
[0], B
->p
[1], B
->p
[0], tmp
, factor
, B
->NbColumns
);
162 if (value_zero_p(B
->p
[0][pos
[1]])) {
163 /* We will deal with this later */
164 assert(lex_sign(B
->p
[0], B
->NbColumns
) < 0);
167 value_set_si(factor
, 1);
168 Vector_Combine(B
->p
[0], B
->p
[1], B
->p
[1], tmp
, factor
, B
->NbColumns
);
169 if (value_zero_p(B
->p
[1][pos
[2]])) {
170 /* We will deal with this later */
171 assert(lex_sign(B
->p
[1], B
->NbColumns
) < 0);
178 bool progress
= true;
181 for (int i
= 0; i
<= 1; ++i
) {
182 value_set_si(factor
, -1);
183 for (int j
= 1; j
<= 3; ++j
) {
184 if (value_zero_p(B
->p
[1-i
][pos
[j
]]))
186 value_oppose(tmp
, B
->p
[i
][pos
[j
]]);
187 value_pdivision(tmp
, tmp
, B
->p
[1-i
][pos
[j
]]);
188 if (value_neg_p(factor
) || value_lt(tmp
, factor
))
189 value_assign(factor
, tmp
);
191 if (value_pos_p(factor
)) {
192 value_set_si(tmp
, 1);
193 Vector_Combine(B
->p
[i
], B
->p
[1-i
], B
->p
[i
], tmp
, factor
,
195 sign
= lex_sign(B
->p
[i
], B
->NbColumns
);
196 for (int j
= 1; j
<= 3; ++j
) {
197 if (value_notzero_p(B
->p
[i
][pos
[j
]]))
199 /* a zero is interpreted to be of sign sign */
200 if ((sign
> 0 && value_pos_p(B
->p
[1-i
][pos
[j
]])) ||
201 (sign
< 0 && value_neg_p(B
->p
[1-i
][pos
[j
]]))) {
202 /* the zero is of the wrong sign => back-off one */
203 value_set_si(tmp2
, -1);
204 Vector_Combine(B
->p
[i
], B
->p
[1-i
], B
->p
[i
], tmp
, tmp2
,
206 value_decrement(factor
, factor
);
209 /* We may have backed-off, so we need to check again. */
210 if (value_pos_p(factor
))
216 for (int i
= 0; i
< B
->NbColumns
; ++i
) {
217 value_addto(tmp
, B
->p
[0][i
], B
->p
[1][i
]);
218 if (value_zero_p(tmp
))
220 sign
= value_neg_p(tmp
) ? -1 : 1;
224 for (int i
= 1; i
<= 3; ++i
) {
225 value_addto(tmp
, B
->p
[0][pos
[i
]], B
->p
[1][pos
[i
]]);
226 if (value_neg_p(tmp
) || (sign
< 0 && value_zero_p(tmp
)))
233 value_set_si(tmp
, 1);
234 Vector_Combine(B
->p
[0], B
->p
[1], B
->p
[1], tmp
, tmp
, B
->NbColumns
);
236 /* put positive pair in position 3 */
237 for (i
= 1; i
<= 3; ++i
)
238 if (value_pos_p(B
->p
[0][pos
[i
]]) && value_pos_p(B
->p
[1][pos
[i
]]))
250 /* We will deal with these later */
263 Value last
; // last multiple of offset in link
265 Matrix
*M
; // rows: elements different from (0,0)
269 M
= Matrix_Alloc(d
, 2);
272 simplex(int d
, int mask
, Value last
) {
273 M
= Matrix_Alloc(d
, 2);
274 offset
= Vector_Alloc(2);
275 value_init(this->last
);
276 value_assign(this->last
, last
);
279 void transform(Matrix
*T
);
281 Polyhedron
*shrunk_polyhedron(Polyhedron
*P
, int dim
, Matrix
*A
,
283 void print(FILE *out
);
286 void simplex::print(FILE *out
)
289 Matrix_Print(out
, P_VALUE_FMT
, M
);
291 fprintf(out
, "%d %d\n", M
->NbRows
, M
->NbColumns
);
292 for (int j
= 0; j
< M
->NbRows
; ++j
) {
293 for (int k
= 0; k
< M
->NbColumns
; ++k
)
294 value_print(out
, P_VALUE_FMT
, M
->p
[j
][k
]);
295 if (mask
& (1 << j
)) {
296 fprintf(out
, " + k * ");
297 for (int k
= 0; k
< M
->NbColumns
; ++k
)
298 value_print(out
, P_VALUE_FMT
, offset
->p
[k
]);
302 fprintf(out
, "\t0 <= k <= ");
303 value_print(out
, P_VALUE_FMT
, last
);
308 static bool lex_smaller(Value
*v1
, Value
*v2
, int n
)
310 for (int i
= 0; i
< n
; ++i
)
311 if (value_lt(v1
[i
], v2
[i
]))
313 else if (value_gt(v1
[i
], v2
[i
]))
318 void simplex::transform(Matrix
*T
)
321 M
= Matrix_Alloc(M2
->NbRows
, M2
->NbColumns
);
322 Matrix_Product(M2
, T
, M
);
326 Vector
*offset2
= offset
;
327 offset
= Vector_Alloc(offset2
->Size
);
328 Vector_Matrix_Product(offset2
->p
, T
, offset
->p
);
329 Vector_Free(offset2
);
333 void simplex::normalize()
336 for (int i
= 1; i
< M
->NbRows
; ++i
)
337 if (lex_smaller(M
->p
[i
], M
->p
[lexmin
], 2))
339 if (lex_sign(M
->p
[lexmin
], 2) < 0) {
342 value_set_si(tmp
, -1);
343 Vector_Scale(M
->p
[lexmin
], M
->p
[lexmin
], tmp
, 2);
344 value_set_si(tmp
, 1);
345 for (int i
= 0; i
< M
->NbRows
; ++i
) {
348 Vector_Combine(M
->p
[lexmin
], M
->p
[i
], M
->p
[i
], tmp
, tmp
, 2);
350 if (offset
&& (mask
& (1 << lexmin
))) {
351 value_set_si(tmp
, -1);
352 Vector_Scale(offset
->p
, offset
->p
, tmp
, 2);
353 mask
^= (1 << M
->NbRows
) - 1 - (1 << lexmin
);
359 static Matrix
*InsertColumn(Matrix
*M
, int pos
)
364 R
= Matrix_Alloc(M
->NbRows
, M
->NbColumns
+1);
365 for (i
= 0; i
< M
->NbRows
; ++i
) {
366 Vector_Copy(M
->p
[i
], R
->p
[i
], pos
);
367 Vector_Copy(M
->p
[i
]+pos
, R
->p
[i
]+pos
+1, M
->NbColumns
-pos
);
372 Polyhedron
*simplex::shrunk_polyhedron(Polyhedron
*P
, int dim
, Matrix
*A
,
375 Matrix
*Constraints
, *b
;
376 Vector
*b_offset
= NULL
;
378 Value min
, min_var
, tmp
;
384 b
= Matrix_Alloc(M
->NbRows
, A
->NbColumns
);
385 Matrix_Product(M
, A
, b
);
388 b_offset
= Vector_Alloc(A
->NbColumns
);
389 Vector_Matrix_Product(offset
->p
, A
, b_offset
->p
);
393 Constraints
= Polyhedron2Constraints(P
);
395 Constraints
= Matrix_Alloc(P
->NbConstraints
+2, P
->Dimension
+2+1);
396 for (int i
= 0; i
< P
->NbConstraints
; ++i
) {
397 Vector_Copy(P
->Constraint
[i
], Constraints
->p
[i
], 1+dim
+2);
398 Vector_Copy(P
->Constraint
[i
]+1+dim
+2, Constraints
->p
[i
]+1+dim
+2+1,
399 (P
->Dimension
+2)-(1+dim
+2));
401 value_set_si(Constraints
->p
[P
->NbConstraints
][0], 1);
402 value_set_si(Constraints
->p
[P
->NbConstraints
][1+dim
+2], 1);
403 value_set_si(Constraints
->p
[P
->NbConstraints
+1][0], 1);
404 value_set_si(Constraints
->p
[P
->NbConstraints
+1][1+dim
+2], -1);
405 value_assign(Constraints
->p
[P
->NbConstraints
+1][Constraints
->NbColumns
-1],
408 constant
= Constraints
->NbColumns
- 1;
410 for (int i
= 0, j
= 0; i
< P
->NbConstraints
; ++i
) {
411 if (value_zero_p(Constraints
->p
[i
][1+dim
]) &&
412 value_zero_p(Constraints
->p
[i
][1+dim
+1]))
414 value_set_si(min
, 0);
415 for (int k
= 0; k
< b
->NbRows
; ++k
) {
416 if (offset
&& (mask
& (1 << k
)))
418 if (value_lt(b
->p
[k
][j
], min
))
419 value_assign(min
, b
->p
[k
][j
]);
421 value_set_si(min_var
, 0);
423 if (value_neg_p(b_offset
->p
[j
])) {
424 value_oppose(min_var
, b_offset
->p
[j
]);
425 value_multiply(min_var
, min_var
, last
);
426 value_increment(min_var
, min_var
);
428 for (int k
= 0; k
< b
->NbRows
; ++k
) {
429 if (!(mask
& (1 << k
)))
431 if (value_lt(b
->p
[k
][j
], min_var
))
432 value_assign(min_var
, b
->p
[k
][j
]);
435 if (!offset
|| value_pos_p(b_offset
->p
[j
])) {
436 if (value_le(min
, min_var
))
437 value_addto(Constraints
->p
[i
][constant
],
438 Constraints
->p
[i
][constant
], min
);
440 value_assign(tmp
, min_var
);
441 value_addmul(tmp
, last
, b_offset
->p
[j
]);
442 if (value_le(tmp
, min
)) {
443 value_addto(Constraints
->p
[i
][constant
],
444 Constraints
->p
[i
][constant
], min_var
);
445 value_addto(Constraints
->p
[i
][1+dim
+2],
446 Constraints
->p
[i
][1+dim
+2], b_offset
->p
[j
]);
448 int lastrow
= Constraints
->NbRows
;
449 int cols
= Constraints
->NbColumns
;
450 Matrix
*C
= Constraints
;
451 Constraints
= AddANullRow(Constraints
);
453 Vector_Copy(Constraints
->p
[i
], Constraints
->p
[lastrow
], cols
);
454 value_addto(Constraints
->p
[i
][constant
],
455 Constraints
->p
[i
][constant
], min_var
);
456 value_addto(Constraints
->p
[i
][1+dim
+2],
457 Constraints
->p
[i
][1+dim
+2], b_offset
->p
[j
]);
458 value_addto(Constraints
->p
[lastrow
][constant
],
459 Constraints
->p
[lastrow
][constant
], min
);
463 if (value_le(min_var
, min
)) {
464 value_addto(Constraints
->p
[i
][constant
],
465 Constraints
->p
[i
][constant
], min_var
);
466 value_addto(Constraints
->p
[i
][1+dim
+2],
467 Constraints
->p
[i
][1+dim
+2], b_offset
->p
[j
]);
469 value_assign(tmp
, min_var
);
470 value_addmul(tmp
, last
, b_offset
->p
[j
]);
471 if (value_le(min
, tmp
)) {
472 value_addto(Constraints
->p
[i
][constant
],
473 Constraints
->p
[i
][constant
], min
);
475 int lastrow
= Constraints
->NbRows
;
476 int cols
= Constraints
->NbColumns
;
477 Matrix
*C
= Constraints
;
478 Constraints
= AddANullRow(Constraints
);
480 Vector_Copy(Constraints
->p
[i
], Constraints
->p
[lastrow
], cols
);
481 value_addto(Constraints
->p
[i
][constant
],
482 Constraints
->p
[i
][constant
], min_var
);
483 value_addto(Constraints
->p
[i
][1+dim
+2],
484 Constraints
->p
[i
][1+dim
+2], b_offset
->p
[j
]);
485 value_addto(Constraints
->p
[lastrow
][constant
],
486 Constraints
->p
[lastrow
][constant
], min
);
492 Q
= Constraints2Polyhedron(Constraints
, MaxRays
);
495 Vector_Free(b_offset
);
497 Matrix_Free(Constraints
);
500 value_clear(min_var
);
505 struct scarf_complex
{
506 vector
<simplex
> simplices
;
507 void add(Matrix
*B
, int pos
[4], int n
);
508 void add(Matrix
*T
, simplex s
);
509 void print(FILE *out
);
511 for (int i
= 0; i
< simplices
.size(); ++i
) {
512 Matrix_Free(simplices
[i
].M
);
513 if (simplices
[i
].offset
) {
514 Vector_Free(simplices
[i
].offset
);
515 value_clear(simplices
[i
].last
);
521 void scarf_complex::add(Matrix
*T
, simplex s
)
525 if (s
.offset
&& lex_sign(s
.offset
->p
, 2) < 0) {
530 /* compute the smallest multiple (factor) of the offset that
531 * makes on of the vertices lexico-negative.
534 for (int i
= 0; i
< s
.M
->NbRows
; ++i
) {
535 if (!(s
.mask
& (1 << i
)))
537 if (lexmin
== -1 || lex_smaller(s
.M
->p
[i
], s
.M
->p
[lexmin
], 2))
540 if (value_zero_p(s
.offset
->p
[0])) {
541 if (value_pos_p(s
.M
->p
[lexmin
][0]))
542 value_increment(factor
, s
.last
);
544 value_oppose(factor
, s
.M
->p
[lexmin
][1]);
545 mpz_cdiv_q(factor
, factor
, s
.offset
->p
[1]);
548 value_oppose(factor
, s
.M
->p
[lexmin
][0]);
549 mpz_cdiv_q(factor
, factor
, s
.offset
->p
[0]);
550 if (mpz_divisible_p(s
.M
->p
[lexmin
][0], s
.offset
->p
[0])) {
551 value_assign(tmp
, s
.M
->p
[lexmin
][1]);
552 value_addmul(tmp
, factor
, s
.offset
->p
[1]);
553 if (value_pos_p(tmp
))
554 value_increment(factor
, factor
);
557 if (value_le(factor
, s
.last
)) {
558 simplex
part(s
.M
->NbRows
, s
.mask
, s
.last
);
559 Vector_Copy(s
.offset
->p
, part
.offset
->p
, 2);
560 value_set_si(tmp
, 1);
561 for (int i
= 0; i
< s
.M
->NbRows
; ++i
) {
562 if (s
.mask
& (1 << i
))
563 Vector_Combine(s
.M
->p
[i
], s
.offset
->p
, part
.M
->p
[i
],
566 Vector_Copy(s
.M
->p
[i
], part
.M
->p
[i
], 2);
568 value_subtract(part
.last
, part
.last
, factor
);
569 value_decrement(s
.last
, factor
);
571 simplices
.push_back(part
);
576 simplices
.push_back(s
);
579 void scarf_complex::add(Matrix
*B
, int pos
[4], int n
)
583 T
= Matrix_Alloc(2, 2);
584 Vector_Copy(B
->p
[0]+B
->NbColumns
-2, T
->p
[0], 2);
585 Vector_Copy(B
->p
[1]+B
->NbColumns
-2, T
->p
[1], 2);
587 if (n
== 3 || value_neg_p(B
->p
[0][pos
[3]])) {
588 assert(n
== 3 || value_neg_p(B
->p
[1][pos
[3]]));
591 value_set_si(s1
.M
->p
[0][0], 0);
592 value_set_si(s1
.M
->p
[0][1], 1);
596 value_set_si(s2
.M
->p
[0][0], 1);
597 value_set_si(s2
.M
->p
[0][1], 1);
601 value_set_si(s3
.M
->p
[0][0], 1);
602 value_set_si(s3
.M
->p
[0][1], 0);
606 value_set_si(s4
.M
->p
[0][0], 0);
607 value_set_si(s4
.M
->p
[0][1], 1);
608 value_set_si(s4
.M
->p
[1][0], 1);
609 value_set_si(s4
.M
->p
[1][1], 1);
613 value_set_si(s5
.M
->p
[0][0], 1);
614 value_set_si(s5
.M
->p
[0][1], 0);
615 value_set_si(s5
.M
->p
[1][0], 1);
616 value_set_si(s5
.M
->p
[1][1], 1);
622 bool progress
= true;
623 Value tmp
, tmp2
, factor
;
630 assert(value_pos_p(B
->p
[0][pos
[3]]));
631 assert(value_pos_p(B
->p
[1][pos
[3]]));
633 h
= Matrix_Alloc(3, 2);
634 value_set_si(h
->p
[0][0], 1);
635 value_set_si(h
->p
[0][1], 0);
636 value_set_si(h
->p
[1][0], 0);
637 value_set_si(h
->p
[1][1], 1);
638 value_set_si(h
->p
[2][0], 1);
639 value_set_si(h
->p
[2][1], 1);
641 offset
= Vector_Alloc(2);
645 for (int i
= 0; i
<= 1; ++i
) {
646 value_set_si(factor
, -1);
647 for (int j
= 1; j
<= 2; ++j
) {
648 if (value_zero_p(B
->p
[1-i
][pos
[j
]]))
650 value_oppose(tmp
, B
->p
[i
][pos
[j
]]);
651 value_pdivision(tmp
, tmp
, B
->p
[1-i
][pos
[j
]]);
652 if (value_neg_p(factor
) || value_lt(tmp
, factor
))
653 value_assign(factor
, tmp
);
655 if (value_pos_p(factor
)) {
656 value_set_si(tmp
, 1);
657 Vector_Combine(B
->p
[i
], B
->p
[1-i
], B
->p
[i
], tmp
, factor
,
659 sign
= lex_sign(B
->p
[i
], B
->NbColumns
);
660 for (int j
= 1; j
<= 2; ++j
) {
661 if (value_notzero_p(B
->p
[i
][pos
[j
]]))
663 /* a zero is interpreted to be of sign sign */
664 if ((sign
> 0 && value_pos_p(B
->p
[1-i
][pos
[j
]])) ||
665 (sign
< 0 && value_neg_p(B
->p
[1-i
][pos
[j
]]))) {
666 /* the zero is of the wrong sign => back-off one */
667 value_set_si(tmp2
, -1);
668 Vector_Combine(B
->p
[i
], B
->p
[1-i
], B
->p
[i
], tmp
, tmp2
,
670 value_decrement(factor
, factor
);
673 /* We may have backed-off, so we need to check again. */
674 if (value_pos_p(factor
)) {
676 value_set_si(tmp
, 1);
677 value_set_si(tmp2
, -1);
679 Vector_Combine(h
->p
[2], h
->p
[i
], offset
->p
, tmp
, tmp2
, 2);
683 Vector_Copy(h
->p
[0], l1
.M
->p
[0], 2);
687 Vector_Copy(h
->p
[1], l2
.M
->p
[0], 2);
691 Vector_Combine(h
->p
[0], h
->p
[1], l3
.M
->p
[0],
696 Vector_Copy(h
->p
[0], t1
.M
->p
[0], 2);
697 Vector_Copy(h
->p
[1], t1
.M
->p
[1], 2);
701 Vector_Combine(h
->p
[0], h
->p
[1], t2
.M
->p
[0],
703 Vector_Combine(h
->p
[2], h
->p
[1], t2
.M
->p
[1],
708 Vector_Combine(h
->p
[i
], offset
->p
, h
->p
[i
],
710 Vector_Combine(h
->p
[2], offset
->p
, h
->p
[2],
712 value_decrement(factor
, factor
);
715 simplex
q(3, 0x4 | (1 << i
), factor
);
716 Vector_Copy(h
->p
[0], q
.M
->p
[0], 2);
717 Vector_Copy(h
->p
[1], q
.M
->p
[1], 2);
718 Vector_Copy(h
->p
[2], q
.M
->p
[2], 2);
719 Vector_Copy(offset
->p
, q
.offset
->p
, 2);
722 simplex
t1(2, 0x3, factor
);
723 Vector_Copy(h
->p
[i
], t1
.M
->p
[0], 2);
724 Vector_Copy(h
->p
[2], t1
.M
->p
[1], 2);
725 Vector_Copy(offset
->p
, t1
.offset
->p
, 2);
728 simplex
t2(2, 0x2, factor
);
729 Vector_Copy(h
->p
[1-i
], t2
.M
->p
[0], 2);
730 Vector_Copy(h
->p
[2], t2
.M
->p
[1], 2);
731 Vector_Copy(offset
->p
, t2
.offset
->p
, 2);
734 simplex
l(1, 0x1, factor
);
735 Vector_Copy(h
->p
[2], l
.M
->p
[0], 2);
736 Vector_Copy(offset
->p
, l
.offset
->p
, 2);
740 Vector_Combine(h
->p
[i
], offset
->p
, h
->p
[i
], tmp
, factor
, 2);
741 Vector_Combine(h
->p
[2], offset
->p
, h
->p
[2], tmp
, factor
, 2);
761 void scarf_complex::print(FILE *out
)
763 for (int i
= 0; i
< simplices
.size(); ++i
)
764 simplices
[i
].print(out
);
767 struct scarf_collector
{
768 virtual void add(Polyhedron
*P
, int sign
, Polyhedron
*C
, unsigned MaxRays
) = 0;
771 static void scarf(Polyhedron
*P
, unsigned exist
, unsigned nparam
, unsigned MaxRays
,
772 scarf_collector
& col
)
775 int dim
= P
->Dimension
- exist
- nparam
;
782 A
= extract_matrix(P
, dim
);
784 n
= A
->NbColumns
- 2;
785 assert(n
>= 3 && n
<= 4);
787 for (int i
= 0; i
< n
; ++i
)
789 B
= normalize_matrix(A
, pos
, n
);
792 scarf
.add(B
, pos
, n
);
794 U
= Universe_Polyhedron(nparam
);
795 col
.add(P
, 0, U
, MaxRays
);
796 for (int i
= 0; i
< scarf
.simplices
.size(); ++i
) {
798 int sign
= (scarf
.simplices
[i
].M
->NbRows
% 2) ? -1 : 1;
799 Q
= scarf
.simplices
[i
].shrunk_polyhedron(P
, dim
, A
, MaxRays
);
800 col
.add(Q
, sign
, U
, MaxRays
);
810 struct scarf_collector_gf
: public scarf_collector
{
814 scarf_collector_gf() {
817 virtual void add(Polyhedron
*P
, int sign
, Polyhedron
*C
, unsigned MaxRays
);
820 void scarf_collector_gf::add(Polyhedron
*P
, int sign
, Polyhedron
*C
,
824 gf
= barvinok_series(P
, C
, MaxRays
);
828 gf2
= barvinok_series(P
, C
, MaxRays
);
834 gen_fun
*barvinok_enumerate_scarf_series(Polyhedron
*P
,
835 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
837 scarf_collector_gf scgf
;
838 scarf(P
, exist
, nparam
, MaxRays
, scgf
);
842 struct scarf_collector_ev
: public scarf_collector
{
846 scarf_collector_ev() {
848 evalue_set_si(&mone
, -1, 1);
850 ~scarf_collector_ev() {
851 free_evalue_refs(&mone
);
853 virtual void add(Polyhedron
*P
, int sign
, Polyhedron
*C
, unsigned MaxRays
);
856 void scarf_collector_ev::add(Polyhedron
*P
, int sign
, Polyhedron
*C
,
860 EP
= barvinok_enumerate_ev(P
, C
, MaxRays
);
863 E2
= barvinok_enumerate_ev(P
, C
, MaxRays
);
867 free_evalue_refs(E2
);
872 evalue
*barvinok_enumerate_scarf(Polyhedron
*P
,
873 unsigned exist
, unsigned nparam
, unsigned MaxRays
)
875 scarf_collector_ev scev
;
876 scarf(P
, exist
, nparam
, MaxRays
, scev
);