2 #include <barvinok/util.h>
4 #include "lattice_point.h"
8 struct OrthogonalException Orthogonal
;
10 void np_base::handle(const signed_cone
& sc
, barvinok_options
*options
)
12 assert(sc
.rays
.NumRows() == dim
);
14 handle(sc
.rays
, current_vertex
, factor
, sc
.closed
, options
);
18 void np_base::start(Polyhedron
*P
, barvinok_options
*options
)
24 for (int i
= 0; i
< P
->NbRays
; ++i
) {
25 if (!value_pos_p(P
->Ray
[i
][dim
+1]))
28 Polyhedron
*C
= supporting_cone(P
, i
);
29 do_vertex_cone(factor
, C
, P
->Ray
[i
]+1, options
);
32 } catch (OrthogonalException
&e
) {
39 * f: the powers in the denominator for the remaining vars
40 * each row refers to a factor
41 * den_s: for each factor, the power of (s+1)
43 * num_s: powers in the numerator corresponding to the summed vars
44 * num_p: powers in the numerator corresponding to the remaining vars
45 * number of rays in cone: "dim" = "k"
46 * length of each ray: "dim" = "d"
47 * for now, it is assumed: k == d
49 * den_p: for each factor
50 * 0: independent of remaining vars
51 * 1: power corresponds to corresponding row in f
53 * all inputs are subject to change
55 void normalize(ZZ
& sign
, ZZ
& num_s
, vec_ZZ
& num_p
, vec_ZZ
& den_s
, vec_ZZ
& den_p
,
58 unsigned dim
= f
.NumRows();
59 unsigned nparam
= num_p
.length();
60 unsigned nvar
= dim
- nparam
;
64 for (int j
= 0; j
< den_s
.length(); ++j
) {
70 for (k
= 0; k
< nparam
; ++k
)
84 den_s
[j
] = abs(den_s
[j
]);
93 void reducer::reduce(QQ c
, vec_ZZ
& num
, const mat_ZZ
& den_f
)
95 unsigned len
= den_f
.NumRows(); // number of factors in den
97 if (num
.length() == lower
) {
101 assert(num
.length() > 1);
108 split(num
, num_s
, num_p
, den_f
, den_s
, den_r
);
111 den_p
.SetLength(len
);
113 normalize(c
.n
, num_s
, num_p
, den_s
, den_p
, den_r
);
115 int only_param
= 0; // k-r-s from text
116 int no_param
= 0; // r from text
117 for (int k
= 0; k
< len
; ++k
) {
120 else if (den_s
[k
] == 0)
124 reduce(c
, num_p
, den_r
);
128 pden
.SetDims(only_param
, den_r
.NumCols());
130 for (k
= 0, l
= 0; k
< len
; ++k
)
132 pden
[l
++] = den_r
[k
];
134 for (k
= 0; k
< len
; ++k
)
138 dpoly
n(no_param
, num_s
);
139 dpoly
D(no_param
, den_s
[k
], 1);
142 dpoly
fact(no_param
, den_s
[k
], 1);
146 if (no_param
+ only_param
== len
) {
147 mpq_set_si(tcount
, 0, 1);
148 n
.div(D
, tcount
, one
);
151 value2zz(mpq_numref(tcount
), q
.n
);
152 value2zz(mpq_denref(tcount
), q
.d
);
157 reduce(q
, num_p
, pden
);
161 for (k
= 0; k
< len
; ++k
) {
162 if (den_s
[k
] == 0 || den_p
[k
] == 0)
165 dpoly
pd(no_param
-1, den_s
[k
], 1);
168 for (l
= 0; l
< k
; ++l
)
169 if (den_r
[l
] == den_r
[k
])
173 r
= new dpoly_r(n
, pd
, l
, len
);
175 dpoly_r
*nr
= new dpoly_r(r
, pd
, l
, len
);
181 dpoly_r
*rc
= r
->div(D
);
184 factor
.d
= rc
->denom
* c
.d
;
186 int common
= pden
.NumRows();
187 dpoly_r_term_list
& final
= rc
->c
[rc
->len
-1];
189 dpoly_r_term_list::iterator j
;
190 for (j
= final
.begin(); j
!= final
.end(); ++j
) {
191 if ((*j
)->coeff
== 0)
194 pden
.SetDims(rows
, pden
.NumCols());
195 for (int k
= 0; k
< rc
->dim
; ++k
) {
196 int n
= (*j
)->powers
[k
];
199 pden
.SetDims(rows
+n
, pden
.NumCols());
200 for (int l
= 0; l
< n
; ++l
)
201 pden
[rows
+l
] = den_r
[k
];
204 factor
.n
= (*j
)->coeff
*= c
.n
;
205 reduce(factor
, num_p
, pden
);
214 void reducer::handle(const mat_ZZ
& den
, Value
*V
, QQ c
, int *closed
,
215 barvinok_options
*options
)
217 lattice_point(V
, den
, vertex
, closed
);
219 reduce(c
, vertex
, den
);
222 void ireducer::split(vec_ZZ
& num
, ZZ
& num_s
, vec_ZZ
& num_p
,
223 const mat_ZZ
& den_f
, vec_ZZ
& den_s
, mat_ZZ
& den_r
)
225 unsigned len
= den_f
.NumRows(); // number of factors in den
226 unsigned d
= num
.length() - 1;
228 den_s
.SetLength(len
);
229 den_r
.SetDims(len
, d
);
231 for (int r
= 0; r
< len
; ++r
) {
232 den_s
[r
] = den_f
[r
][0];
233 for (int k
= 1; k
<= d
; ++k
)
234 den_r
[r
][k
-1] = den_f
[r
][k
];
239 for (int k
= 1 ; k
<= d
; ++k
)
243 void normalize(ZZ
& sign
, ZZ
& num
, vec_ZZ
& den
)
245 unsigned dim
= den
.length();
249 for (int j
= 0; j
< den
.length(); ++j
) {
253 den
[j
] = abs(den
[j
]);
261 void icounter::base(QQ
& c
, const vec_ZZ
& num
, const mat_ZZ
& den_f
)
264 unsigned len
= den_f
.NumRows(); // number of factors in den
266 den_s
.SetLength(len
);
268 for (r
= 0; r
< len
; ++r
)
269 den_s
[r
] = den_f
[r
][0];
270 normalize(c
.n
, num_s
, den_s
);
273 dpoly
D(len
, den_s
[0], 1);
274 for (int k
= 1; k
< len
; ++k
) {
275 dpoly
fact(len
, den_s
[k
], 1);
278 mpq_set_si(tcount
, 0, 1);
279 n
.div(D
, tcount
, one
);
282 mpz_mul(mpq_numref(tcount
), mpq_numref(tcount
), tn
);
283 mpz_mul(mpq_denref(tcount
), mpq_denref(tcount
), td
);
284 mpq_canonicalize(tcount
);
285 mpq_add(count
, count
, tcount
);
288 void infinite_icounter::base(QQ
& c
, const vec_ZZ
& num
, const mat_ZZ
& den_f
)
291 unsigned len
= den_f
.NumRows(); // number of factors in den
293 den_s
.SetLength(len
);
296 for (r
= 0; r
< len
; ++r
)
297 den_s
[r
] = den_f
[r
][0];
298 normalize(c
.n
, num_s
, den_s
);
301 dpoly
D(len
, den_s
[0], 1);
302 for (int k
= 1; k
< len
; ++k
) {
303 dpoly
fact(len
, den_s
[k
], 1);
312 value_assign(mpq_numref(factor
), tmp
);
314 value_assign(mpq_denref(factor
), tmp
);
316 n
.div(D
, count
, factor
);