?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
gen_fun::substitute: add more detailed explanation
[barvinok.git] / genfun.cc
bloba8437fa1ed45df3a6edc2d6918af6d71cceba214
1 #include <iostream>
2 #include <assert.h>
3 #include "config.h"
4 #include <barvinok/genfun.h>
5 #include <barvinok/barvinok.h>
6
7 using std::cout;
8
9 static int lex_cmp(vec_ZZ& a, vec_ZZ& b)
10 {
11 assert(a.length() == b.length());
12
13 for (int j = 0; j < a.length(); ++j)
14 if (a[j] != b[j])
15 return a[j] < b[j] ? -1 : 1;
16 return 0;
17 }
18
19 static int lex_cmp(mat_ZZ& a, mat_ZZ& b)
20 {
21 assert(a.NumCols() == b.NumCols());
22 int alen = a.NumRows();
23 int blen = b.NumRows();
24 int len = alen < blen ? alen : blen;
25
26 for (int i = 0; i < len; ++i) {
27 int s = lex_cmp(a[i], b[i]);
28 if (s)
29 return s;
30 }
31 return alen-blen;
32 }
33
34 void gen_fun::add(const ZZ& cn, const ZZ& cd, const vec_ZZ& num,
35 const mat_ZZ& den)
36 {
37 if (cn == 0)
38 return;
39
40 short_rat * r = new short_rat;
41 r->n.coeff.SetDims(1, 2);
42 ZZ g = GCD(cn, cd);
43 r->n.coeff[0][0] = cn/g;
44 r->n.coeff[0][1] = cd/g;
45 r->n.power.SetDims(1, num.length());
46 r->n.power[0] = num;
47 r->d.power = den;
48
49 /* Make all powers in denominator lexico-positive */
50 for (int i = 0; i < r->d.power.NumRows(); ++i) {
51 int j;
52 for (j = 0; j < r->d.power.NumCols(); ++j)
53 if (r->d.power[i][j] != 0)
54 break;
55 if (r->d.power[i][j] < 0) {
56 r->d.power[i] = -r->d.power[i];
57 r->n.coeff[0][0] = -r->n.coeff[0][0];
58 r->n.power[0] += r->d.power[i];
59 }
60 }
61
62 /* Order powers in denominator */
63 for (int i = 0; i < r->d.power.NumRows(); ++i) {
64 int m = i;
65 for (int j = i+1; j < r->d.power.NumRows(); ++j)
66 if (lex_cmp(r->d.power[j], r->d.power[m]) < 0)
67 m = j;
68 if (m != i) {
69 vec_ZZ tmp = r->d.power[m];
70 r->d.power[m] = r->d.power[i];
71 r->d.power[i] = tmp;
72 }
73 }
74
75 for (int i = 0; i < term.size(); ++i)
76 if (lex_cmp(term[i]->d.power, r->d.power) == 0) {
77 int len = term[i]->n.coeff.NumRows();
78 int j;
79 for (j = 0; j < len; ++j)
80 if (r->n.power[0] == term[i]->n.power[j])
81 break;
82 if (j < len) {
83 ZZ g = GCD(r->n.coeff[0][1], term[i]->n.coeff[j][1]);
84 ZZ n = term[i]->n.coeff[j][0] * (r->n.coeff[0][1] / g) +
85 (term[i]->n.coeff[j][1] / g) * r->n.coeff[0][0];
86 ZZ d = term[i]->n.coeff[j][1] / g * r->n.coeff[0][1];
87 if (n != 0) {
88 g = GCD(n,d);
89 term[i]->n.coeff[j][0] = n/g;
90 term[i]->n.coeff[j][1] = d/g;
91 } else {
92 if (len > 1) {
93 if (j < len-1) {
94 term[i]->n.power[j] = term[i]->n.power[len-1];
95 term[i]->n.coeff[j] = term[i]->n.coeff[len-1];
96 }
97 int dim = term[i]->n.power.NumCols();
98 term[i]->n.coeff.SetDims(len-1, 2);
99 term[i]->n.power.SetDims(len-1, dim);
100 } else {
101 delete term[i];
102 if (i != term.size()-1)
103 term[i] = term[term.size()-1];
104 term.pop_back();
105 }
106 }
107 } else {
108 int dim = term[i]->n.power.NumCols();
109 term[i]->n.coeff.SetDims(len+1, 2);
110 term[i]->n.power.SetDims(len+1, dim);
111 term[i]->n.coeff[len] = r->n.coeff[0];
112 term[i]->n.power[len] = r->n.power[0];
113 }
114 delete r;
115 return;
116 }
117
118 term.push_back(r);
119 }
120
121 void gen_fun::add(const ZZ& cn, const ZZ& cd, gen_fun *gf)
122 {
123 ZZ n, d;
124 for (int i = 0; i < gf->term.size(); ++i) {
125 for (int j = 0; j < gf->term[i]->n.power.NumRows(); ++j) {
126 n = cn * gf->term[i]->n.coeff[j][0];
127 d = cd * gf->term[i]->n.coeff[j][1];
128 add(n, d, gf->term[i]->n.power[j], gf->term[i]->d.power);
129 }
130 }
131 }
132
133 /*
134 * Perform the substitution specified by CP and (map, offset)
135 *
136 * CP is a homogeneous matrix that maps a set of "compressed parameters"
137 * to the original set of parameters.
138 *
139 * This function is applied to a gen_fun computed with the compressed parameters
140 * and adapts it to refer to the original parameters.
141 *
142 * That is, if y are the compressed parameters and x = A y + b are the original
143 * parameters, then we want the coefficient of the monomial t^y in the original
144 * generating function to be the coefficient of the monomial u^x in the resulting
145 * generating function.
146 * The original generating function has the form
147 *
148 * a t^m/(1-t^n) = a t^m + a t^{m+n} + a t^{m+2n} + ...
149 *
150 * Since each term t^y should correspond to a term u^x, with x = A y + b, we want
151 *
152 * a u^{A m + b} + a u^{A (m+n) + b} + a u^{A (m+2n) +b} + ... =
153 *
154 * = a u^{A m + b}/(1-u^{A n})
155 *
156 * Therefore, we multiply the powers m and n in both numerator and denominator by A
157 * and add b to the power in the numerator.
158 * Since the above powers are stored as row vectors m^T and n^T,
159 * we compute, say, m'^T = m^T A^T to obtain m' = A m.
160 *
161 * The pair (map, offset) contains the same information as CP.
162 * map is the transpose of the linear part of CP, while offset is the constant part.
163 */
164 void gen_fun::substitute(Matrix *CP, const mat_ZZ& map, const vec_ZZ& offset)
165 {
166 Polyhedron *C = Polyhedron_Image(context, CP, 0);
167 Polyhedron_Free(context);
168 context = C;
169 for (int i = 0; i < term.size(); ++i) {
170 term[i]->d.power *= map;
171 term[i]->n.power *= map;
172 for (int j = 0; j < term[i]->n.power.NumRows(); ++j)
173 term[i]->n.power[j] += offset;
174 }
175 }
176
177 gen_fun *gen_fun::Hadamard_product(gen_fun *gf, unsigned MaxRays)
178 {
179 Polyhedron *C = DomainIntersection(context, gf->context, MaxRays);
180 Polyhedron *U = Universe_Polyhedron(C->Dimension);
181 gen_fun *sum = new gen_fun(C);
182 for (int i = 0; i < term.size(); ++i) {
183 for (int i2 = 0; i2 < gf->term.size(); ++i2) {
184 int d = term[i]->d.power.NumCols();
185 int k1 = term[i]->d.power.NumRows();
186 int k2 = gf->term[i2]->d.power.NumRows();
187 assert(term[i]->d.power.NumCols() == gf->term[i2]->d.power.NumCols());
188 for (int j = 0; j < term[i]->n.power.NumRows(); ++j) {
189 for (int j2 = 0; j2 < gf->term[i2]->n.power.NumRows(); ++j2) {
190 Matrix *M = Matrix_Alloc(k1+k2+d+d, 1+k1+k2+d+1);
191 for (int k = 0; k < k1+k2; ++k) {
192 value_set_si(M->p[k][0], 1);
193 value_set_si(M->p[k][1+k], 1);
194 }
195 for (int k = 0; k < d; ++k) {
196 value_set_si(M->p[k1+k2+k][1+k1+k2+k], -1);
197 zz2value(term[i]->n.power[j][k], M->p[k1+k2+k][1+k1+k2+d]);
198 for (int l = 0; l < k1; ++l)
199 zz2value(term[i]->d.power[l][k], M->p[k1+k2+k][1+l]);
200 }
201 for (int k = 0; k < d; ++k) {
202 value_set_si(M->p[k1+k2+d+k][1+k1+k2+k], -1);
203 zz2value(gf->term[i2]->n.power[j2][k],
204 M->p[k1+k2+d+k][1+k1+k2+d]);
205 for (int l = 0; l < k2; ++l)
206 zz2value(gf->term[i2]->d.power[l][k],
207 M->p[k1+k2+d+k][1+k1+l]);
208 }
209 Polyhedron *P = Constraints2Polyhedron(M, MaxRays);
210 Matrix_Free(M);
211
212 gen_fun *t = barvinok_series(P, U, MaxRays);
213
214 ZZ cn = term[i]->n.coeff[j][0] * gf->term[i2]->n.coeff[j2][0];
215 ZZ cd = term[i]->n.coeff[j][1] * gf->term[i2]->n.coeff[j2][1];
216 sum->add(cn, cd, t);
217 delete t;
218
219 Polyhedron_Free(P);
220 }
221 }
222 }
223 }
224 Polyhedron_Free(U);
225 return sum;
226 }
227
228 void gen_fun::add_union(gen_fun *gf, unsigned MaxRays)
229 {
230 ZZ one, mone;
231 one = 1;
232 mone = -1;
233
234 gen_fun *hp = gf->Hadamard_product(gf, MaxRays);
235 add(one, one, gf);
236 add(mone, one, hp);
237 delete hp;
238 }
239
240 static void print_power(vec_ZZ& c, vec_ZZ& p,
241 unsigned int nparam, char **param_name)
242 {
243 bool first = true;
244
245 for (int i = 0; i < p.length(); ++i) {
246 if (p[i] == 0)
247 continue;
248 if (first) {
249 if (c[0] == -1 && c[1] == 1)
250 cout << "-";
251 else if (c[0] != 1 || c[1] != 1) {
252 cout << c[0];
253 if (c[1] != 1)
254 cout << " / " << c[1];
255 cout << "*";
256 }
257 first = false;
258 } else
259 cout << "*";
260 if (i < nparam)
261 cout << param_name[i];
262 else
263 cout << "x" << i;
264 if (p[i] == 1)
265 continue;
266 if (p[i] < 0)
267 cout << "^(" << p[i] << ")";
268 else
269 cout << "^" << p[i];
270 }
271 if (first) {
272 cout << c[0];
273 if (c[1] != 1)
274 cout << " / " << c[1];
275 }
276 }
277
278 void gen_fun::print(unsigned int nparam, char **param_name) const
279 {
280 vec_ZZ mone;
281 mone.SetLength(2);
282 mone[0] = -1;
283 mone[1] = 1;
284 for (int i = 0; i < term.size(); ++i) {
285 if (i != 0)
286 cout << " + ";
287 cout << "(";
288 for (int j = 0; j < term[i]->n.coeff.NumRows(); ++j) {
289 if (j != 0 && term[i]->n.coeff[j][0] > 0)
290 cout << "+";
291 print_power(term[i]->n.coeff[j], term[i]->n.power[j],
292 nparam, param_name);
293 }
294 cout << ")/(";
295 for (int j = 0; j < term[i]->d.power.NumRows(); ++j) {
296 if (j != 0)
297 cout << " * ";
298 cout << "(1";
299 print_power(mone, term[i]->d.power[j],
300 nparam, param_name);
301 cout << ")";
302 }
303 cout << ")";
304 }
305 }
306
307 gen_fun::operator evalue *() const
308 {
309 evalue *EP = NULL;
310 evalue factor;
311 value_init(factor.d);
312 value_init(factor.x.n);
313 for (int i = 0; i < term.size(); ++i) {
314 unsigned nvar = term[i]->d.power.NumRows();
315 unsigned nparam = term[i]->d.power.NumCols();
316 Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + nparam + 1);
317 mat_ZZ& d = term[i]->d.power;
318 Polyhedron *U = context ? context : Universe_Polyhedron(nparam);
319
320 for (int j = 0; j < term[i]->n.coeff.NumRows(); ++j) {
321 for (int r = 0; r < nparam; ++r) {
322 value_set_si(C->p[r][0], 0);
323 for (int c = 0; c < nvar; ++c) {
324 zz2value(d[c][r], C->p[r][1+c]);
325 }
326 Vector_Set(&C->p[r][1+nvar], 0, nparam);
327 value_set_si(C->p[r][1+nvar+r], -1);
328 zz2value(term[i]->n.power[j][r], C->p[r][1+nvar+nparam]);
329 }
330 for (int r = 0; r < nvar; ++r) {
331 value_set_si(C->p[nparam+r][0], 1);
332 Vector_Set(&C->p[nparam+r][1], 0, nvar + nparam + 1);
333 value_set_si(C->p[nparam+r][1+r], 1);
334 }
335 Polyhedron *P = Constraints2Polyhedron(C, 0);
336 evalue *E = barvinok_enumerate_ev(P, U, 0);
337 Polyhedron_Free(P);
338 if (EVALUE_IS_ZERO(*E)) {
339 free_evalue_refs(E);
340 free(E);
341 continue;
342 }
343 zz2value(term[i]->n.coeff[j][0], factor.x.n);
344 zz2value(term[i]->n.coeff[j][1], factor.d);
345 emul(&factor, E);
346 /*
347 Matrix_Print(stdout, P_VALUE_FMT, C);
348 char *test[] = { "A", "B", "C", "D", "E", "F", "G" };
349 print_evalue(stdout, E, test);
350 */
351 if (!EP)
352 EP = E;
353 else {
354 eadd(E, EP);
355 free_evalue_refs(E);
356 free(E);
357 }
358 }
359 Matrix_Free(C);
360 if (!context)
361 Polyhedron_Free(U);
362 }
363 value_clear(factor.d);
364 value_clear(factor.x.n);
365 return EP;
366 }
367
368 void gen_fun::coefficient(Value* params, Value* c) const
369 {
370 if (context && !in_domain(context, params)) {
371 value_set_si(*c, 0);
372 return;
373 }
374
375 evalue part;
376 value_init(part.d);
377 value_init(part.x.n);
378 evalue sum;
379 value_init(sum.d);
380 evalue_set_si(&sum, 0, 1);
381 Value tmp;
382 value_init(tmp);
383
384 for (int i = 0; i < term.size(); ++i) {
385 unsigned nvar = term[i]->d.power.NumRows();
386 unsigned nparam = term[i]->d.power.NumCols();
387 Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + 1);
388 mat_ZZ& d = term[i]->d.power;
389
390 for (int j = 0; j < term[i]->n.coeff.NumRows(); ++j) {
391 for (int r = 0; r < nparam; ++r) {
392 value_set_si(C->p[r][0], 0);
393 for (int c = 0; c < nvar; ++c) {
394 zz2value(d[c][r], C->p[r][1+c]);
395 }
396 zz2value(term[i]->n.power[j][r], C->p[r][1+nvar]);
397 value_subtract(C->p[r][1+nvar], C->p[r][1+nvar], params[r]);
398 }
399 for (int r = 0; r < nvar; ++r) {
400 value_set_si(C->p[nparam+r][0], 1);
401 Vector_Set(&C->p[nparam+r][1], 0, nvar + 1);
402 value_set_si(C->p[nparam+r][1+r], 1);
403 }
404 Polyhedron *P = Constraints2Polyhedron(C, 0);
405 if (emptyQ(P)) {
406 Polyhedron_Free(P);
407 continue;
408 }
409 barvinok_count(P, &tmp, 0);
410 Polyhedron_Free(P);
411 if (value_zero_p(tmp))
412 continue;
413 zz2value(term[i]->n.coeff[j][0], part.x.n);
414 zz2value(term[i]->n.coeff[j][1], part.d);
415 value_multiply(part.x.n, part.x.n, tmp);
416 eadd(&part, &sum);
417 }
418 Matrix_Free(C);
419 }
420
421 assert(value_one_p(sum.d));
422 value_assign(*c, sum.x.n);
423
424 value_clear(tmp);
425 value_clear(part.d);
426 value_clear(part.x.n);
427 value_clear(sum.d);
428 value_clear(sum.x.n);
429 }
lattice point enumeration library