2 #include <barvinok/util.h>
4 #include "lattice_point.h"
10 struct OrthogonalException Orthogonal
;
12 void np_base::handle(const signed_cone
& sc
, barvinok_options
*options
)
14 assert(sc
.rays
.NumRows() == dim
);
16 handle(sc
.rays
, current_vertex
, factor
, sc
.det
, sc
.closed
, options
);
20 void np_base::start(Polyhedron
*P
, barvinok_options
*options
)
26 for (int i
= 0; i
< P
->NbRays
; ++i
) {
27 if (!value_pos_p(P
->Ray
[i
][dim
+1]))
30 Polyhedron
*C
= supporting_cone(P
, i
);
31 do_vertex_cone(factor
, C
, P
->Ray
[i
]+1, options
);
34 } catch (OrthogonalException
&e
) {
41 * f: the powers in the denominator for the remaining vars
42 * each row refers to a factor
43 * den_s: for each factor, the power of (s+1)
45 * num_s: powers in the numerator corresponding to the summed vars
46 * num_p: powers in the numerator corresponding to the remaining vars
47 * number of rays in cone: "dim" = "k"
48 * length of each ray: "dim" = "d"
49 * for now, it is assumed: k == d
51 * den_p: for each factor
52 * 0: independent of remaining vars
53 * 1: power corresponds to corresponding row in f
55 * all inputs are subject to change
57 void normalize(ZZ
& sign
, vec_ZZ
& num_s
, mat_ZZ
& num_p
, vec_ZZ
& den_s
, vec_ZZ
& den_p
,
60 unsigned dim
= f
.NumRows();
61 unsigned nparam
= num_p
.NumCols();
62 unsigned nvar
= dim
- nparam
;
66 for (int j
= 0; j
< den_s
.length(); ++j
) {
72 for (k
= 0; k
< nparam
; ++k
)
79 for (int i
= 0; i
< num_p
.NumRows(); ++i
)
87 den_s
[j
] = abs(den_s
[j
]);
88 for (int i
= 0; i
< num_p
.NumRows(); ++i
)
97 void reducer::base(const vec_QQ
& c
, const mat_ZZ
& num
, const mat_ZZ
& den_f
)
99 for (int i
= 0; i
< num
.NumRows(); ++i
)
100 base(c
[i
], num
[i
], den_f
);
103 struct dpoly_r_scanner
{
105 const dpoly
* const *num
;
108 dpoly_r_term_list::iterator
*iter
;
112 dpoly_r_scanner(const dpoly
* const *num
, int n
, const dpoly_r
*rc
, int dim
)
113 : num(num
), rc(rc
), n(n
), dim(dim
), powers(dim
, 0) {
115 iter
= new dpoly_r_term_list::iterator
[rc
->len
];
116 for (int i
= 0; i
< rc
->len
; ++i
) {
118 for (k
= 0; k
< n
; ++k
)
119 if (num
[k
]->coeff
[rc
->len
-1-i
] != 0)
122 iter
[i
] = rc
->c
[i
].begin();
124 iter
[i
] = rc
->c
[i
].end();
131 for (int i
= 0; i
< rc
->len
; ++i
) {
132 if (iter
[i
] == rc
->c
[i
].end())
137 if ((*iter
[i
])->powers
< (*iter
[pos
[0]])->powers
) {
140 } else if ((*iter
[i
])->powers
== (*iter
[pos
[0]])->powers
)
148 powers
= (*iter
[pos
[0]])->powers
;
149 for (int k
= 0; k
< n
; ++k
)
150 mul(coeff
[k
], (*iter
[pos
[0]])->coeff
, num
[k
]->coeff
[rc
->len
-1-pos
[0]]);
152 for (int i
= 1; i
< len
; ++i
) {
153 for (int k
= 0; k
< n
; ++k
) {
154 mul(tmp
, (*iter
[pos
[i
]])->coeff
, num
[k
]->coeff
[rc
->len
-1-pos
[i
]]);
155 add(coeff
[k
], coeff
[k
], tmp
);
169 void reducer::reduce(const vec_QQ
& c
, const mat_ZZ
& num
, const mat_ZZ
& den_f
)
171 assert(c
.length() == num
.NumRows());
172 unsigned len
= den_f
.NumRows(); // number of factors in den
175 if (num
.NumCols() == lower
) {
179 assert(num
.NumCols() > 1);
180 assert(num
.NumRows() > 0);
187 split(num
, num_s
, num_p
, den_f
, den_s
, den_r
);
190 den_p
.SetLength(len
);
192 ZZ
sign(INIT_VAL
, 1);
193 normalize(sign
, num_s
, num_p
, den_s
, den_p
, den_r
);
196 int only_param
= 0; // k-r-s from text
197 int no_param
= 0; // r from text
198 for (int k
= 0; k
< len
; ++k
) {
201 else if (den_s
[k
] == 0)
205 reduce(c2
, num_p
, den_r
);
209 pden
.SetDims(only_param
, den_r
.NumCols());
211 for (k
= 0, l
= 0; k
< len
; ++k
)
213 pden
[l
++] = den_r
[k
];
215 for (k
= 0; k
< len
; ++k
)
219 dpoly
*n
[num_s
.length()];
220 for (int i
= 0; i
< num_s
.length(); ++i
) {
221 n
[i
] = new dpoly(no_param
, num_s
[i
]);
222 /* Search for other numerator (j) with same num_p.
223 * If found, replace a[j]/b[j] * n[j] and a[i]/b[i] * n[i]
224 * by 1/(b[j]*b[i]/g) * (a[j]*b[i]/g * n[j] + a[i]*b[j]/g * n[i])
225 * where g = gcd(b[i], b[j].
227 for (int j
= 0; j
< i
; ++j
) {
228 if (num_p
[i
] != num_p
[j
])
230 ZZ g
= GCD(c2
[i
].d
, c2
[j
].d
);
231 *n
[j
] *= c2
[j
].n
* c2
[i
].d
/g
;
232 *n
[i
] *= c2
[i
].n
* c2
[j
].d
/g
;
235 c2
[j
].d
*= c2
[i
].d
/g
;
237 if (i
< num_s
.length()-1) {
238 num_s
[i
] = num_s
[num_s
.length()-1];
239 num_p
[i
] = num_p
[num_s
.length()-1];
240 c2
[i
] = c2
[num_s
.length()-1];
242 num_s
.SetLength(num_s
.length()-1);
243 c2
.SetLength(c2
.length()-1);
244 num_p
.SetDims(num_p
.NumRows()-1, num_p
.NumCols());
249 dpoly
D(no_param
, den_s
[k
], 1);
252 dpoly
fact(no_param
, den_s
[k
], 1);
256 if (no_param
+ only_param
== len
) {
258 q
.SetLength(num_s
.length());
259 for (int i
= 0; i
< num_s
.length(); ++i
) {
260 mpq_set_si(tcount
, 0, 1);
261 n
[i
]->div(D
, tcount
, one
);
263 value2zz(mpq_numref(tcount
), q
[i
].n
);
264 value2zz(mpq_denref(tcount
), q
[i
].d
);
267 for (int i
= q
.length()-1; i
>= 0; --i
) {
269 q
[i
] = q
[q
.length()-1];
270 num_p
[i
] = num_p
[q
.length()-1];
271 q
.SetLength(q
.length()-1);
272 num_p
.SetDims(num_p
.NumRows()-1, num_p
.NumCols());
277 reduce(q
, num_p
, pden
);
279 ZZ
zz_zero(INIT_VAL
, 0);
280 dpoly
one(no_param
, zz_zero
);
283 for (k
= 0; k
< len
; ++k
) {
284 if (den_s
[k
] == 0 || den_p
[k
] == 0)
287 dpoly
pd(no_param
-1, den_s
[k
], 1);
290 for (l
= 0; l
< k
; ++l
)
291 if (den_r
[l
] == den_r
[k
])
295 r
= new dpoly_r(one
, pd
, l
, len
);
297 dpoly_r
*nr
= new dpoly_r(r
, pd
, l
, len
);
304 factor
.SetLength(c2
.length());
305 int common
= pden
.NumRows();
306 dpoly_r
*rc
= r
->div(D
);
307 for (int i
= 0; i
< num_s
.length(); ++i
) {
308 factor
[i
].d
= c2
[i
].d
;
309 factor
[i
].d
*= rc
->denom
;
312 dpoly_r_scanner
scanner(n
, num_s
.length(), rc
, len
);
314 while (scanner
.next()) {
316 for (i
= 0; i
< num_s
.length(); ++i
)
317 if (scanner
.coeff
[i
] != 0)
319 if (i
== num_s
.length())
322 pden
.SetDims(rows
, pden
.NumCols());
323 for (int k
= 0; k
< rc
->dim
; ++k
) {
324 int n
= scanner
.powers
[k
];
327 pden
.SetDims(rows
+n
, pden
.NumCols());
328 for (int l
= 0; l
< n
; ++l
)
329 pden
[rows
+l
] = den_r
[k
];
332 for (int i
= 0; i
< num_s
.length(); ++i
) {
333 factor
[i
].n
= c2
[i
].n
;
334 factor
[i
].n
*= scanner
.coeff
[i
];
336 reduce(factor
, num_p
, pden
);
342 for (int i
= 0; i
< num_s
.length(); ++i
)
347 void reducer::handle(const mat_ZZ
& den
, Value
*V
, const QQ
& c
, unsigned long det
,
348 int *closed
, barvinok_options
*options
)
352 lattice_point(V
, den
, vertex
, det
, closed
);
354 vc
.SetLength(vertex
.NumRows());
355 for (int i
= 0; i
< vc
.length(); ++i
)
358 reduce(vc
, vertex
, den
);
361 void split_one(const mat_ZZ
& num
, vec_ZZ
& num_s
, mat_ZZ
& num_p
,
362 const mat_ZZ
& den_f
, vec_ZZ
& den_s
, mat_ZZ
& den_r
)
364 unsigned len
= den_f
.NumRows(); // number of factors in den
365 unsigned d
= num
.NumCols() - 1;
367 den_s
.SetLength(len
);
368 den_r
.SetDims(len
, d
);
370 for (int r
= 0; r
< len
; ++r
) {
371 den_s
[r
] = den_f
[r
][0];
372 for (int k
= 1; k
<= d
; ++k
)
373 den_r
[r
][k
-1] = den_f
[r
][k
];
376 num_s
.SetLength(num
.NumRows());
377 num_p
.SetDims(num
.NumRows(), d
);
378 for (int i
= 0; i
< num
.NumRows(); ++i
) {
379 num_s
[i
] = num
[i
][0];
380 for (int k
= 1 ; k
<= d
; ++k
)
381 num_p
[i
][k
-1] = num
[i
][k
];
385 void normalize(ZZ
& sign
, ZZ
& num
, vec_ZZ
& den
)
387 unsigned dim
= den
.length();
391 for (int j
= 0; j
< den
.length(); ++j
) {
395 den
[j
] = abs(den
[j
]);
403 void icounter::base(const QQ
& c
, const vec_ZZ
& num
, const mat_ZZ
& den_f
)
406 unsigned len
= den_f
.NumRows(); // number of factors in den
408 den_s
.SetLength(len
);
409 assert(num
.length() == 1);
411 for (r
= 0; r
< len
; ++r
)
412 den_s
[r
] = den_f
[r
][0];
413 ZZ sign
= ZZ(INIT_VAL
, 1);
414 normalize(sign
, num_s
, den_s
);
417 dpoly
D(len
, den_s
[0], 1);
418 for (int k
= 1; k
< len
; ++k
) {
419 dpoly
fact(len
, den_s
[k
], 1);
422 mpq_set_si(tcount
, 0, 1);
423 n
.div(D
, tcount
, one
);
426 value_oppose(tn
, tn
);
428 mpz_mul(mpq_numref(tcount
), mpq_numref(tcount
), tn
);
429 mpz_mul(mpq_denref(tcount
), mpq_denref(tcount
), td
);
430 mpq_canonicalize(tcount
);
431 mpq_add(count
, count
, tcount
);
434 void infinite_icounter::base(const QQ
& c
, const vec_ZZ
& num
, const mat_ZZ
& den_f
)
437 unsigned len
= den_f
.NumRows(); // number of factors in den
439 den_s
.SetLength(len
);
440 assert(num
.length() == 1);
443 for (r
= 0; r
< len
; ++r
)
444 den_s
[r
] = den_f
[r
][0];
445 ZZ sign
= ZZ(INIT_VAL
, 1);
446 normalize(sign
, num_s
, den_s
);
449 dpoly
D(len
, den_s
[0], 1);
450 for (int k
= 1; k
< len
; ++k
) {
451 dpoly
fact(len
, den_s
[k
], 1);
461 value_oppose(tmp
, tmp
);
462 value_assign(mpq_numref(factor
), tmp
);
464 value_assign(mpq_denref(factor
), tmp
);
466 n
.div(D
, count
, factor
);