1 #include <assert.h>
2 #include <NTL/ZZ.h>
3 #include <NTL/vec_ZZ.h>
4 #include <NTL/mat_ZZ.h>
5 #include <barvinok/polylib.h>
6 #include <barvinok/util.h>
7 #include "counter.h"
8 #include "lattice_point.h"
9
10 /* Computes the integer points in the fundamental parallelepiped and
11 * passes them along (in num) to the counter specific (i.e., specialization
12 * specific) add_lattice_points.
13 */
handle(const mat_ZZ & rays,Value * V,const QQ & c,unsigned long det,barvinok_options * options)14 void counter_base::handle(const mat_ZZ& rays, Value *V, const QQ& c,
15 unsigned long det, barvinok_options *options)
16 {
17 Matrix* Rays = zz2matrix(rays);
18
19 assert(c.d == 1);
20 assert(c.n == 1 || c.n == -1);
21 int sign = to_int(c.n);
22
23 Matrix_Vector_Product(Rays, lambda->p, den->p_Init);
24 for (int k = 0; k < dim; ++k)
25 if (value_zero_p(den->p_Init[k])) {
26 Matrix_Free(Rays);
27 throw Orthogonal;
28 }
29 Inner_Product(lambda->p, V, dim, &tmp);
30 lattice_points_fixed(V, &tmp, Rays, den, num, det);
31 num->NbRows = det;
32 Matrix_Free(Rays);
33
34 if (dim % 2)
35 sign = -sign;
36
37 add_lattice_points(sign);
38 }
39
add_falling_powers(dpoly & n,Value degree)40 static void add_falling_powers(dpoly& n, Value degree)
41 {
42 value_increment(n.coeff->p[0], n.coeff->p[0]);
43 if (n.coeff->Size == 1)
44 return;
45
46 int min = n.coeff->Size-1;
47 if (value_posz_p(degree) && value_cmp_si(degree, min) < 0)
48 min = VALUE_TO_INT(degree);
49
50 Value tmp;
51 value_init(tmp);
52 value_assign(tmp, degree);
53 value_addto(n.coeff->p[1], n.coeff->p[1], tmp);
54 for (int i = 2; i <= min; ++i) {
55 value_decrement(degree, degree);
56 value_multiply(tmp, tmp, degree);
57 mpz_divexact_ui(tmp, tmp, i);
58 value_addto(n.coeff->p[i], n.coeff->p[i], tmp);
59 }
60 value_clear(tmp);
61 }
62
add_lattice_points(int sign)63 void counter::add_lattice_points(int sign)
64 {
65 dpoly d(dim);
66 for (int k = 0; k < num->NbRows; ++k)
67 add_falling_powers(d, num->p_Init[k]);
68 dpoly n(dim, den->p_Init[0], 1);
69 for (int k = 1; k < dim; ++k) {
70 dpoly fact(dim, den->p_Init[k], 1);
71 n *= fact;
72 }
73 d.div(n, count, sign);
74 }
75
76
77
78
setup_todd(unsigned dim)79 void tcounter::setup_todd(unsigned dim)
80 {
81 value_set_si(todd.coeff->p[0], 1);
82
83 dpoly d(dim);
84 value_set_si(d.coeff->p[dim], 1);
85 for (int i = dim-1; i >= 0; --i)
86 mpz_mul_ui(d.coeff->p[i], d.coeff->p[i+1], i+2);
87
88 todd_denom = todd.div(d);
89 /* shift denominators up -> divide by (dim+1)! */
90 for (int i = todd.coeff->Size-1; i >= 1; --i)
91 value_assign(todd_denom->p[i], todd_denom->p[i-1]);
92 value_set_si(todd_denom->p[0], 1);
93 }
94
adapt_todd(dpoly & t,const Value c)95 void tcounter::adapt_todd(dpoly& t, const Value c)
96 {
97 if (t.coeff->Size <= 1)
98 return;
99 value_assign(tmp, c);
100 value_multiply(t.coeff->p[1], t.coeff->p[1], tmp);
101 for (int i = 2; i < t.coeff->Size; ++i) {
102 value_multiply(tmp, tmp, c);
103 value_multiply(t.coeff->p[i], t.coeff->p[i], tmp);
104 }
105 }
106
add_powers(dpoly & n,const Value c)107 void tcounter::add_powers(dpoly& n, const Value c)
108 {
109 value_increment(n.coeff->p[0], n.coeff->p[0]);
110 if (n.coeff->Size == 1)
111 return;
112 value_assign(tmp, c);
113 value_addto(n.coeff->p[1], n.coeff->p[1], tmp);
114 for (int i = 2; i < n.coeff->Size; ++i) {
115 value_multiply(tmp, tmp, c);
116 value_addto(n.coeff->p[i], n.coeff->p[i], tmp);
117 }
118 }
119
add_lattice_points(int sign)120 void tcounter::add_lattice_points(int sign)
121 {
122 dpoly t(todd);
123 value_assign(denom, den->p_Init[0]);
124 adapt_todd(t, den->p_Init[0]);
125 for (int k = 1; k < dim; ++k) {
126 dpoly fact(todd);
127 value_multiply(denom, denom, den->p_Init[k]);
128 adapt_todd(fact, den->p_Init[k]);
129 t *= fact;
130 }
131
132 dpoly n(dim);
133 for (int k = 0; k < num->NbRows; ++k)
134 add_powers(n, num->p_Init[k]);
135
136 for (int i = 0; i < n.coeff->Size; ++i)
137 value_multiply(n.coeff->p[i], n.coeff->p[i], todd_denom->p[i]);
138 value_multiply(denom, denom, todd_denom->p[todd_denom->Size-1]);
139
140 value_set_si(tmp, 1);
141 for (int i = 2; i < n.coeff->Size; ++i) {
142 mpz_mul_ui(tmp, tmp, i);
143 mpz_divexact(n.coeff->p[i], n.coeff->p[i], tmp);
144 }
145
146 value_multiply(tmp, t.coeff->p[0], n.coeff->p[n.coeff->Size-1]);
147 for (int i = 1; i < n.coeff->Size; ++i)
148 value_addmul(tmp, t.coeff->p[i], n.coeff->p[n.coeff->Size-1-i]);
149
150 value_assign(mpq_numref(tcount), tmp);
151 value_assign(mpq_denref(tcount), denom);
152 mpq_canonicalize(tcount);
153 if (sign == -1)
154 mpq_sub(count, count, tcount);
155 else
156 mpq_add(count, count, tcount);
157 }
158
159
160 /*
161 * Set lambda to a random vector that has a positive inner product
162 * with all the rays of the context { x | A x + b >= 0 }.
163 *
164 * To do so, we take d rows A' from the constraint matrix A.
165 * For every ray, we have
166 * A' r >= 0
167 * We compute a random positive row vector x' and set x = x' A'.
168 * We then have, for each ray r,
169 * x r = x' A' r >= 0
170 * Although we can take any d rows from A, we choose linearly
171 * independent rows from A to avoid the elements of the transformed
172 * random vector to have simple linear relations, which would
173 * increase the risk of the vector being orthogonal to one of
174 * powers in the denominator of one of the terms in the generating
175 * function.
176 */
init(Polyhedron * context,int n_try)177 void infinite_counter::init(Polyhedron *context, int n_try)
178 {
179 Matrix *M, *H, *Q, *U;
180 mat_ZZ A;
181
182 randomvector(context, lambda, context->Dimension, n_try);
183
184 M = Matrix_Alloc(context->NbConstraints, context->Dimension);
185 for (int i = 0; i < context->NbConstraints; ++i)
186 Vector_Copy(context->Constraint[i]+1, M->p[i], context->Dimension);
187 left_hermite(M, &H, &Q, &U);
188 Matrix_Free(Q);
189 Matrix_Free(U);
190
191 for (int col = 0, row = 0; col < H->NbColumns; ++col, ++row) {
192 for (; row < H->NbRows; ++row)
193 if (value_notzero_p(H->p[row][col]))
194 break;
195 assert(row < H->NbRows);
196 Vector_Copy(M->p[row], M->p[col], M->NbColumns);
197 }
198 matrix2zz(M, A, context->Dimension, context->Dimension);
199 Matrix_Free(H);
200 Matrix_Free(M);
201
202 for (int i = 0; i < lambda.length(); ++i)
203 lambda[i] = abs(lambda[i]);
204 lambda = lambda * A;
205 }
206
LCM(const ZZ & a,const ZZ & b)207 static ZZ LCM(const ZZ& a, const ZZ& b)
208 {
209 return a * b / GCD(a, b);
210 }
211
212 /* Normalize the powers in the denominator to be positive
213 * and return -1 is the sign has to be changed.
214 */
normalized_sign(vec_ZZ & num,vec_ZZ & den)215 static int normalized_sign(vec_ZZ& num, vec_ZZ& den)
216 {
217 int change = 0;
218
219 for (int j = 0; j < den.length(); ++j) {
220 if (den[j] > 0)
221 change ^= 1;
222 else {
223 den[j] = abs(den[j]);
224 for (int k = 0; k < num.length(); ++k)
225 num[k] += den[j];
226 }
227 }
228 return change ? -1 : 1;
229 }
230
reduce(const vec_QQ & c,const mat_ZZ & num,const mat_ZZ & den_f)231 void infinite_counter::reduce(const vec_QQ& c, const mat_ZZ& num,
232 const mat_ZZ& den_f)
233 {
234 mpq_t factor;
235 mpq_init(factor);
236 unsigned len = den_f.NumRows();
237
238 ZZ lcm = c[0].d;
239 for (int i = 1; i < c.length(); ++i)
240 lcm = LCM(lcm, c[i].d);
241
242 vec_ZZ coeff;
243 coeff.SetLength(c.length());
244 for (int i = 0; i < c.length(); ++i)
245 coeff[i] = c[i].n * lcm/c[i].d;
246
247 if (len == 0) {
248 for (int i = 0; i < c.length(); ++i) {
249 zz2value(coeff[i], tz);
250 value_addto(mpq_numref(factor), mpq_numref(factor), tz);
251 }
252 zz2value(lcm, tz);
253 value_assign(mpq_denref(factor), tz);
254 mpq_add(count[0], count[0], factor);
255 mpq_clear(factor);
256 return;
257 }
258
259 vec_ZZ num_s = num * lambda;
260 vec_ZZ den_s = den_f * lambda;
261 for (int i = 0; i < den_s.length(); ++i)
262 assert(den_s[i] != 0);
263 int sign = normalized_sign(num_s, den_s);
264 if (sign < 0)
265 coeff = -coeff;
266
267 dpoly n(len);
268 zz2value(num_s[0], tz);
269 add_falling_powers(n, tz);
270 zz2value(coeff[0], tz);
271 n *= tz;
272 for (int i = 1; i < c.length(); ++i) {
273 dpoly t(len);
274 zz2value(num_s[i], tz);
275 add_falling_powers(t, tz);
276 zz2value(coeff[i], tz);
277 t *= tz;
278 n += t;
279 }
280 zz2value(den_s[0], tz);
281 dpoly d(len, tz, 1);
282 for (int i = 1; i < len; ++i) {
283 zz2value(den_s[i], tz);
284 dpoly fact(len, tz, 1);
285 d *= fact;
286 }
287 value_set_si(mpq_numref(factor), 1);
288 zz2value(lcm, tz);
289 value_assign(mpq_denref(factor), tz);
290 n.div(d, count, factor);
291
292 mpq_clear(factor);
293 }
294