4 #include <barvinok/genfun.h>
5 #include <barvinok/barvinok.h>
10 static int lex_cmp(mat_ZZ
& a
, mat_ZZ
& b
)
12 assert(a
.NumCols() == b
.NumCols());
13 int alen
= a
.NumRows();
14 int blen
= b
.NumRows();
15 int len
= alen
< blen
? alen
: blen
;
17 for (int i
= 0; i
< len
; ++i
) {
18 int s
= lex_cmp(a
[i
], b
[i
]);
25 static void lex_order_terms(struct short_rat
* rat
)
27 for (int i
= 0; i
< rat
->n
.power
.NumRows(); ++i
) {
29 for (int j
= i
+1; j
< rat
->n
.power
.NumRows(); ++j
)
30 if (lex_cmp(rat
->n
.power
[j
], rat
->n
.power
[m
]) < 0)
33 vec_ZZ tmp
= rat
->n
.power
[m
];
34 rat
->n
.power
[m
] = rat
->n
.power
[i
];
35 rat
->n
.power
[i
] = tmp
;
36 tmp
= rat
->n
.coeff
[m
];
37 rat
->n
.coeff
[m
] = rat
->n
.coeff
[i
];
38 rat
->n
.coeff
[i
] = tmp
;
43 void short_rat::add(short_rat
*r
)
45 for (int i
= 0; i
< r
->n
.power
.NumRows(); ++i
) {
46 int len
= n
.coeff
.NumRows();
48 for (j
= 0; j
< len
; ++j
)
49 if (r
->n
.power
[i
] == n
.power
[j
])
52 ZZ g
= GCD(r
->n
.coeff
[i
][1], n
.coeff
[j
][1]);
53 ZZ num
= n
.coeff
[j
][0] * (r
->n
.coeff
[i
][1] / g
) +
54 (n
.coeff
[j
][1] / g
) * r
->n
.coeff
[i
][0];
55 ZZ d
= n
.coeff
[j
][1] / g
* r
->n
.coeff
[i
][1];
58 n
.coeff
[j
][0] = num
/g
;
62 n
.power
[j
] = n
.power
[len
-1];
63 n
.coeff
[j
] = n
.coeff
[len
-1];
65 int dim
= n
.power
.NumCols();
66 n
.coeff
.SetDims(len
-1, 2);
67 n
.power
.SetDims(len
-1, dim
);
70 int dim
= n
.power
.NumCols();
71 n
.coeff
.SetDims(len
+1, 2);
72 n
.power
.SetDims(len
+1, dim
);
73 n
.coeff
[len
] = r
->n
.coeff
[i
];
74 n
.power
[len
] = r
->n
.power
[i
];
79 bool short_rat::reduced()
81 int dim
= n
.power
.NumCols();
82 lex_order_terms(this);
83 if (n
.power
.NumRows() % 2 == 0) {
84 if (n
.coeff
[0][0] == -n
.coeff
[1][0] &&
85 n
.coeff
[0][1] == n
.coeff
[1][1]) {
86 vec_ZZ step
= n
.power
[1] - n
.power
[0];
88 for (k
= 1; k
< n
.power
.NumRows()/2; ++k
) {
89 if (n
.coeff
[2*k
][0] != -n
.coeff
[2*k
+1][0] ||
90 n
.coeff
[2*k
][1] != n
.coeff
[2*k
+1][1])
92 if (step
!= n
.power
[2*k
+1] - n
.power
[2*k
])
95 if (k
== n
.power
.NumRows()/2) {
96 for (k
= 0; k
< d
.power
.NumRows(); ++k
)
97 if (d
.power
[k
] == step
)
99 if (k
< d
.power
.NumRows()) {
100 for (++k
; k
< d
.power
.NumRows(); ++k
)
101 d
.power
[k
-1] = d
.power
[k
];
102 d
.power
.SetDims(k
-1, dim
);
103 for (k
= 1; k
< n
.power
.NumRows()/2; ++k
) {
104 n
.coeff
[k
] = n
.coeff
[2*k
];
105 n
.power
[k
] = n
.power
[2*k
];
107 n
.coeff
.SetDims(k
, 2);
108 n
.power
.SetDims(k
, dim
);
117 void gen_fun::add(const ZZ
& cn
, const ZZ
& cd
, const vec_ZZ
& num
,
123 short_rat
* r
= new short_rat
;
124 r
->n
.coeff
.SetDims(1, 2);
126 r
->n
.coeff
[0][0] = cn
/g
;
127 r
->n
.coeff
[0][1] = cd
/g
;
128 r
->n
.power
.SetDims(1, num
.length());
132 /* Make all powers in denominator lexico-positive */
133 for (int i
= 0; i
< r
->d
.power
.NumRows(); ++i
) {
135 for (j
= 0; j
< r
->d
.power
.NumCols(); ++j
)
136 if (r
->d
.power
[i
][j
] != 0)
138 if (r
->d
.power
[i
][j
] < 0) {
139 r
->d
.power
[i
] = -r
->d
.power
[i
];
140 r
->n
.coeff
[0][0] = -r
->n
.coeff
[0][0];
141 r
->n
.power
[0] += r
->d
.power
[i
];
145 /* Order powers in denominator */
146 lex_order_rows(r
->d
.power
);
148 for (int i
= 0; i
< term
.size(); ++i
)
149 if (lex_cmp(term
[i
]->d
.power
, r
->d
.power
) == 0) {
151 if (term
[i
]->n
.coeff
.NumRows() == 0) {
153 if (i
!= term
.size()-1)
154 term
[i
] = term
[term
.size()-1];
156 } else if (term
[i
]->reduced()) {
158 /* we've modified term[i], so removed it
159 * and add it back again
162 if (i
!= term
.size()-1)
163 term
[i
] = term
[term
.size()-1];
175 void gen_fun::add(const ZZ
& cn
, const ZZ
& cd
, gen_fun
*gf
)
178 for (int i
= 0; i
< gf
->term
.size(); ++i
) {
179 for (int j
= 0; j
< gf
->term
[i
]->n
.power
.NumRows(); ++j
) {
180 n
= cn
* gf
->term
[i
]->n
.coeff
[j
][0];
181 d
= cd
* gf
->term
[i
]->n
.coeff
[j
][1];
182 add(n
, d
, gf
->term
[i
]->n
.power
[j
], gf
->term
[i
]->d
.power
);
188 * Perform the substitution specified by CP and (map, offset)
190 * CP is a homogeneous matrix that maps a set of "compressed parameters"
191 * to the original set of parameters.
193 * This function is applied to a gen_fun computed with the compressed parameters
194 * and adapts it to refer to the original parameters.
196 * That is, if y are the compressed parameters and x = A y + b are the original
197 * parameters, then we want the coefficient of the monomial t^y in the original
198 * generating function to be the coefficient of the monomial u^x in the resulting
199 * generating function.
200 * The original generating function has the form
202 * a t^m/(1-t^n) = a t^m + a t^{m+n} + a t^{m+2n} + ...
204 * Since each term t^y should correspond to a term u^x, with x = A y + b, we want
206 * a u^{A m + b} + a u^{A (m+n) + b} + a u^{A (m+2n) +b} + ... =
208 * = a u^{A m + b}/(1-u^{A n})
210 * Therefore, we multiply the powers m and n in both numerator and denominator by A
211 * and add b to the power in the numerator.
212 * Since the above powers are stored as row vectors m^T and n^T,
213 * we compute, say, m'^T = m^T A^T to obtain m' = A m.
215 * The pair (map, offset) contains the same information as CP.
216 * map is the transpose of the linear part of CP, while offset is the constant part.
218 void gen_fun::substitute(Matrix
*CP
, const mat_ZZ
& map
, const vec_ZZ
& offset
)
220 Polyhedron
*C
= Polyhedron_Image(context
, CP
, 0);
221 Polyhedron_Free(context
);
223 for (int i
= 0; i
< term
.size(); ++i
) {
224 term
[i
]->d
.power
*= map
;
225 term
[i
]->n
.power
*= map
;
226 for (int j
= 0; j
< term
[i
]->n
.power
.NumRows(); ++j
)
227 term
[i
]->n
.power
[j
] += offset
;
231 gen_fun
*gen_fun::Hadamard_product(gen_fun
*gf
, unsigned MaxRays
)
233 Polyhedron
*C
= DomainIntersection(context
, gf
->context
, MaxRays
);
234 Polyhedron
*U
= Universe_Polyhedron(C
->Dimension
);
235 gen_fun
*sum
= new gen_fun(C
);
236 for (int i
= 0; i
< term
.size(); ++i
) {
237 for (int i2
= 0; i2
< gf
->term
.size(); ++i2
) {
238 int d
= term
[i
]->d
.power
.NumCols();
239 int k1
= term
[i
]->d
.power
.NumRows();
240 int k2
= gf
->term
[i2
]->d
.power
.NumRows();
241 assert(term
[i
]->d
.power
.NumCols() == gf
->term
[i2
]->d
.power
.NumCols());
242 for (int j
= 0; j
< term
[i
]->n
.power
.NumRows(); ++j
) {
243 for (int j2
= 0; j2
< gf
->term
[i2
]->n
.power
.NumRows(); ++j2
) {
244 Matrix
*M
= Matrix_Alloc(k1
+k2
+d
+d
, 1+k1
+k2
+d
+1);
245 for (int k
= 0; k
< k1
+k2
; ++k
) {
246 value_set_si(M
->p
[k
][0], 1);
247 value_set_si(M
->p
[k
][1+k
], 1);
249 for (int k
= 0; k
< d
; ++k
) {
250 value_set_si(M
->p
[k1
+k2
+k
][1+k1
+k2
+k
], -1);
251 zz2value(term
[i
]->n
.power
[j
][k
], M
->p
[k1
+k2
+k
][1+k1
+k2
+d
]);
252 for (int l
= 0; l
< k1
; ++l
)
253 zz2value(term
[i
]->d
.power
[l
][k
], M
->p
[k1
+k2
+k
][1+l
]);
255 for (int k
= 0; k
< d
; ++k
) {
256 value_set_si(M
->p
[k1
+k2
+d
+k
][1+k1
+k2
+k
], -1);
257 zz2value(gf
->term
[i2
]->n
.power
[j2
][k
],
258 M
->p
[k1
+k2
+d
+k
][1+k1
+k2
+d
]);
259 for (int l
= 0; l
< k2
; ++l
)
260 zz2value(gf
->term
[i2
]->d
.power
[l
][k
],
261 M
->p
[k1
+k2
+d
+k
][1+k1
+l
]);
263 Polyhedron
*P
= Constraints2Polyhedron(M
, MaxRays
);
266 gen_fun
*t
= barvinok_series(P
, U
, MaxRays
);
268 ZZ cn
= term
[i
]->n
.coeff
[j
][0] * gf
->term
[i2
]->n
.coeff
[j2
][0];
269 ZZ cd
= term
[i
]->n
.coeff
[j
][1] * gf
->term
[i2
]->n
.coeff
[j2
][1];
282 void gen_fun::add_union(gen_fun
*gf
, unsigned MaxRays
)
288 gen_fun
*hp
= Hadamard_product(gf
, MaxRays
);
294 static void print_power(vec_ZZ
& c
, vec_ZZ
& p
,
295 unsigned int nparam
, char **param_name
)
299 for (int i
= 0; i
< p
.length(); ++i
) {
303 if (c
[0] == -1 && c
[1] == 1)
305 else if (c
[0] != 1 || c
[1] != 1) {
308 cout
<< " / " << c
[1];
315 cout
<< param_name
[i
];
321 cout
<< "^(" << p
[i
] << ")";
328 cout
<< " / " << c
[1];
332 void gen_fun::print(unsigned int nparam
, char **param_name
) const
338 for (int i
= 0; i
< term
.size(); ++i
) {
342 for (int j
= 0; j
< term
[i
]->n
.coeff
.NumRows(); ++j
) {
343 if (j
!= 0 && term
[i
]->n
.coeff
[j
][0] > 0)
345 print_power(term
[i
]->n
.coeff
[j
], term
[i
]->n
.power
[j
],
349 for (int j
= 0; j
< term
[i
]->d
.power
.NumRows(); ++j
) {
353 print_power(mone
, term
[i
]->d
.power
[j
],
361 gen_fun::operator evalue
*() const
365 value_init(factor
.d
);
366 value_init(factor
.x
.n
);
367 for (int i
= 0; i
< term
.size(); ++i
) {
368 unsigned nvar
= term
[i
]->d
.power
.NumRows();
369 unsigned nparam
= term
[i
]->d
.power
.NumCols();
370 Matrix
*C
= Matrix_Alloc(nparam
+ nvar
, 1 + nvar
+ nparam
+ 1);
371 mat_ZZ
& d
= term
[i
]->d
.power
;
372 Polyhedron
*U
= context
? context
: Universe_Polyhedron(nparam
);
374 for (int j
= 0; j
< term
[i
]->n
.coeff
.NumRows(); ++j
) {
375 for (int r
= 0; r
< nparam
; ++r
) {
376 value_set_si(C
->p
[r
][0], 0);
377 for (int c
= 0; c
< nvar
; ++c
) {
378 zz2value(d
[c
][r
], C
->p
[r
][1+c
]);
380 Vector_Set(&C
->p
[r
][1+nvar
], 0, nparam
);
381 value_set_si(C
->p
[r
][1+nvar
+r
], -1);
382 zz2value(term
[i
]->n
.power
[j
][r
], C
->p
[r
][1+nvar
+nparam
]);
384 for (int r
= 0; r
< nvar
; ++r
) {
385 value_set_si(C
->p
[nparam
+r
][0], 1);
386 Vector_Set(&C
->p
[nparam
+r
][1], 0, nvar
+ nparam
+ 1);
387 value_set_si(C
->p
[nparam
+r
][1+r
], 1);
389 Polyhedron
*P
= Constraints2Polyhedron(C
, 0);
390 evalue
*E
= barvinok_enumerate_ev(P
, U
, 0);
392 if (EVALUE_IS_ZERO(*E
)) {
397 zz2value(term
[i
]->n
.coeff
[j
][0], factor
.x
.n
);
398 zz2value(term
[i
]->n
.coeff
[j
][1], factor
.d
);
401 Matrix_Print(stdout, P_VALUE_FMT, C);
402 char *test[] = { "A", "B", "C", "D", "E", "F", "G" };
403 print_evalue(stdout, E, test);
417 value_clear(factor
.d
);
418 value_clear(factor
.x
.n
);
422 void gen_fun::coefficient(Value
* params
, Value
* c
) const
424 if (context
&& !in_domain(context
, params
)) {
431 value_init(part
.x
.n
);
434 evalue_set_si(&sum
, 0, 1);
438 for (int i
= 0; i
< term
.size(); ++i
) {
439 unsigned nvar
= term
[i
]->d
.power
.NumRows();
440 unsigned nparam
= term
[i
]->d
.power
.NumCols();
441 Matrix
*C
= Matrix_Alloc(nparam
+ nvar
, 1 + nvar
+ 1);
442 mat_ZZ
& d
= term
[i
]->d
.power
;
444 for (int j
= 0; j
< term
[i
]->n
.coeff
.NumRows(); ++j
) {
445 for (int r
= 0; r
< nparam
; ++r
) {
446 value_set_si(C
->p
[r
][0], 0);
447 for (int c
= 0; c
< nvar
; ++c
) {
448 zz2value(d
[c
][r
], C
->p
[r
][1+c
]);
450 zz2value(term
[i
]->n
.power
[j
][r
], C
->p
[r
][1+nvar
]);
451 value_subtract(C
->p
[r
][1+nvar
], C
->p
[r
][1+nvar
], params
[r
]);
453 for (int r
= 0; r
< nvar
; ++r
) {
454 value_set_si(C
->p
[nparam
+r
][0], 1);
455 Vector_Set(&C
->p
[nparam
+r
][1], 0, nvar
+ 1);
456 value_set_si(C
->p
[nparam
+r
][1+r
], 1);
458 Polyhedron
*P
= Constraints2Polyhedron(C
, 0);
463 barvinok_count(P
, &tmp
, 0);
465 if (value_zero_p(tmp
))
467 zz2value(term
[i
]->n
.coeff
[j
][0], part
.x
.n
);
468 zz2value(term
[i
]->n
.coeff
[j
][1], part
.d
);
469 value_multiply(part
.x
.n
, part
.x
.n
, tmp
);
475 assert(value_one_p(sum
.d
));
476 value_assign(*c
, sum
.x
.n
);
480 value_clear(part
.x
.n
);
482 value_clear(sum
.x
.n
);