4 #include <barvinok/genfun.h>
5 #include <barvinok/barvinok.h>
6 #include "conversion.h"
11 static int lex_cmp(mat_ZZ
& a
, mat_ZZ
& b
)
13 assert(a
.NumCols() == b
.NumCols());
14 int alen
= a
.NumRows();
15 int blen
= b
.NumRows();
16 int len
= alen
< blen
? alen
: blen
;
18 for (int i
= 0; i
< len
; ++i
) {
19 int s
= lex_cmp(a
[i
], b
[i
]);
26 static void lex_order_terms(struct short_rat
* rat
)
28 for (int i
= 0; i
< rat
->n
.power
.NumRows(); ++i
) {
30 for (int j
= i
+1; j
< rat
->n
.power
.NumRows(); ++j
)
31 if (lex_cmp(rat
->n
.power
[j
], rat
->n
.power
[m
]) < 0)
34 vec_ZZ tmp
= rat
->n
.power
[m
];
35 rat
->n
.power
[m
] = rat
->n
.power
[i
];
36 rat
->n
.power
[i
] = tmp
;
37 QQ tmp_coeff
= rat
->n
.coeff
[m
];
38 rat
->n
.coeff
[m
] = rat
->n
.coeff
[i
];
39 rat
->n
.coeff
[i
] = tmp_coeff
;
44 void short_rat::add(short_rat
*r
)
46 for (int i
= 0; i
< r
->n
.power
.NumRows(); ++i
) {
47 int len
= n
.coeff
.length();
49 for (j
= 0; j
< len
; ++j
)
50 if (r
->n
.power
[i
] == n
.power
[j
])
53 n
.coeff
[j
] += r
->n
.coeff
[i
];
54 if (n
.coeff
[j
].n
== 0) {
56 n
.power
[j
] = n
.power
[len
-1];
57 n
.coeff
[j
] = n
.coeff
[len
-1];
59 int dim
= n
.power
.NumCols();
60 n
.coeff
.SetLength(len
-1);
61 n
.power
.SetDims(len
-1, dim
);
64 int dim
= n
.power
.NumCols();
65 n
.coeff
.SetLength(len
+1);
66 n
.power
.SetDims(len
+1, dim
);
67 n
.coeff
[len
] = r
->n
.coeff
[i
];
68 n
.power
[len
] = r
->n
.power
[i
];
73 bool short_rat::reduced()
75 int dim
= n
.power
.NumCols();
76 lex_order_terms(this);
77 if (n
.power
.NumRows() % 2 == 0) {
78 if (n
.coeff
[0].n
== -n
.coeff
[1].n
&&
79 n
.coeff
[0].d
== n
.coeff
[1].d
) {
80 vec_ZZ step
= n
.power
[1] - n
.power
[0];
82 for (k
= 1; k
< n
.power
.NumRows()/2; ++k
) {
83 if (n
.coeff
[2*k
].n
!= -n
.coeff
[2*k
+1].n
||
84 n
.coeff
[2*k
].d
!= n
.coeff
[2*k
+1].d
)
86 if (step
!= n
.power
[2*k
+1] - n
.power
[2*k
])
89 if (k
== n
.power
.NumRows()/2) {
90 for (k
= 0; k
< d
.power
.NumRows(); ++k
)
91 if (d
.power
[k
] == step
)
93 if (k
< d
.power
.NumRows()) {
94 for (++k
; k
< d
.power
.NumRows(); ++k
)
95 d
.power
[k
-1] = d
.power
[k
];
96 d
.power
.SetDims(k
-1, dim
);
97 for (k
= 1; k
< n
.power
.NumRows()/2; ++k
) {
98 n
.coeff
[k
] = n
.coeff
[2*k
];
99 n
.power
[k
] = n
.power
[2*k
];
101 n
.coeff
.SetLength(k
);
102 n
.power
.SetDims(k
, dim
);
111 void gen_fun::add(const QQ
& c
, const vec_ZZ
& num
, const mat_ZZ
& den
)
116 short_rat
* r
= new short_rat
;
117 r
->n
.coeff
.SetLength(1);
118 ZZ g
= GCD(c
.n
, c
.d
);
119 r
->n
.coeff
[0].n
= c
.n
/g
;
120 r
->n
.coeff
[0].d
= c
.d
/g
;
121 r
->n
.power
.SetDims(1, num
.length());
125 /* Make all powers in denominator lexico-positive */
126 for (int i
= 0; i
< r
->d
.power
.NumRows(); ++i
) {
128 for (j
= 0; j
< r
->d
.power
.NumCols(); ++j
)
129 if (r
->d
.power
[i
][j
] != 0)
131 if (r
->d
.power
[i
][j
] < 0) {
132 r
->d
.power
[i
] = -r
->d
.power
[i
];
133 r
->n
.coeff
[0].n
= -r
->n
.coeff
[0].n
;
134 r
->n
.power
[0] += r
->d
.power
[i
];
138 /* Order powers in denominator */
139 lex_order_rows(r
->d
.power
);
141 for (int i
= 0; i
< term
.size(); ++i
)
142 if (lex_cmp(term
[i
]->d
.power
, r
->d
.power
) == 0) {
144 if (term
[i
]->n
.coeff
.length() == 0) {
146 if (i
!= term
.size()-1)
147 term
[i
] = term
[term
.size()-1];
149 } else if (term
[i
]->reduced()) {
151 /* we've modified term[i], so removed it
152 * and add it back again
155 if (i
!= term
.size()-1)
156 term
[i
] = term
[term
.size()-1];
168 void gen_fun::add(const QQ
& c
, const gen_fun
*gf
)
171 for (int i
= 0; i
< gf
->term
.size(); ++i
) {
172 for (int j
= 0; j
< gf
->term
[i
]->n
.power
.NumRows(); ++j
) {
174 p
*= gf
->term
[i
]->n
.coeff
[j
];
175 add(p
, gf
->term
[i
]->n
.power
[j
], gf
->term
[i
]->d
.power
);
181 * Perform the substitution specified by CP and (map, offset)
183 * CP is a homogeneous matrix that maps a set of "compressed parameters"
184 * to the original set of parameters.
186 * This function is applied to a gen_fun computed with the compressed parameters
187 * and adapts it to refer to the original parameters.
189 * That is, if y are the compressed parameters and x = A y + b are the original
190 * parameters, then we want the coefficient of the monomial t^y in the original
191 * generating function to be the coefficient of the monomial u^x in the resulting
192 * generating function.
193 * The original generating function has the form
195 * a t^m/(1-t^n) = a t^m + a t^{m+n} + a t^{m+2n} + ...
197 * Since each term t^y should correspond to a term u^x, with x = A y + b, we want
199 * a u^{A m + b} + a u^{A (m+n) + b} + a u^{A (m+2n) +b} + ... =
201 * = a u^{A m + b}/(1-u^{A n})
203 * Therefore, we multiply the powers m and n in both numerator and denominator by A
204 * and add b to the power in the numerator.
205 * Since the above powers are stored as row vectors m^T and n^T,
206 * we compute, say, m'^T = m^T A^T to obtain m' = A m.
208 * The pair (map, offset) contains the same information as CP.
209 * map is the transpose of the linear part of CP, while offset is the constant part.
211 void gen_fun::substitute(Matrix
*CP
, const mat_ZZ
& map
, const vec_ZZ
& offset
)
213 Polyhedron
*C
= Polyhedron_Image(context
, CP
, 0);
214 Polyhedron_Free(context
);
216 for (int i
= 0; i
< term
.size(); ++i
) {
217 term
[i
]->d
.power
*= map
;
218 term
[i
]->n
.power
*= map
;
219 for (int j
= 0; j
< term
[i
]->n
.power
.NumRows(); ++j
)
220 term
[i
]->n
.power
[j
] += offset
;
224 gen_fun
*gen_fun::Hadamard_product(gen_fun
*gf
, unsigned MaxRays
)
226 Polyhedron
*C
= DomainIntersection(context
, gf
->context
, MaxRays
);
227 Polyhedron
*U
= Universe_Polyhedron(C
->Dimension
);
228 gen_fun
*sum
= new gen_fun(C
);
229 for (int i
= 0; i
< term
.size(); ++i
) {
230 for (int i2
= 0; i2
< gf
->term
.size(); ++i2
) {
231 int d
= term
[i
]->d
.power
.NumCols();
232 int k1
= term
[i
]->d
.power
.NumRows();
233 int k2
= gf
->term
[i2
]->d
.power
.NumRows();
234 assert(term
[i
]->d
.power
.NumCols() == gf
->term
[i2
]->d
.power
.NumCols());
235 for (int j
= 0; j
< term
[i
]->n
.power
.NumRows(); ++j
) {
236 for (int j2
= 0; j2
< gf
->term
[i2
]->n
.power
.NumRows(); ++j2
) {
237 Matrix
*M
= Matrix_Alloc(k1
+k2
+d
+d
, 1+k1
+k2
+d
+1);
238 for (int k
= 0; k
< k1
+k2
; ++k
) {
239 value_set_si(M
->p
[k
][0], 1);
240 value_set_si(M
->p
[k
][1+k
], 1);
242 for (int k
= 0; k
< d
; ++k
) {
243 value_set_si(M
->p
[k1
+k2
+k
][1+k1
+k2
+k
], -1);
244 zz2value(term
[i
]->n
.power
[j
][k
], M
->p
[k1
+k2
+k
][1+k1
+k2
+d
]);
245 for (int l
= 0; l
< k1
; ++l
)
246 zz2value(term
[i
]->d
.power
[l
][k
], M
->p
[k1
+k2
+k
][1+l
]);
248 for (int k
= 0; k
< d
; ++k
) {
249 value_set_si(M
->p
[k1
+k2
+d
+k
][1+k1
+k2
+k
], -1);
250 zz2value(gf
->term
[i2
]->n
.power
[j2
][k
],
251 M
->p
[k1
+k2
+d
+k
][1+k1
+k2
+d
]);
252 for (int l
= 0; l
< k2
; ++l
)
253 zz2value(gf
->term
[i2
]->d
.power
[l
][k
],
254 M
->p
[k1
+k2
+d
+k
][1+k1
+l
]);
256 Polyhedron
*P
= Constraints2Polyhedron(M
, MaxRays
);
259 gen_fun
*t
= barvinok_series(P
, U
, MaxRays
);
261 QQ c
= term
[i
]->n
.coeff
[j
];
262 c
*= gf
->term
[i2
]->n
.coeff
[j2
];
275 void gen_fun::add_union(gen_fun
*gf
, unsigned MaxRays
)
277 QQ
one(1, 1), mone(-1, 1);
279 gen_fun
*hp
= Hadamard_product(gf
, MaxRays
);
285 static void Polyhedron_Shift(Polyhedron
*P
, Vector
*offset
)
289 for (int i
= 0; i
< P
->NbConstraints
; ++i
) {
290 Inner_Product(P
->Constraint
[i
]+1, offset
->p
, P
->Dimension
, &tmp
);
291 value_subtract(P
->Constraint
[i
][1+P
->Dimension
],
292 P
->Constraint
[i
][1+P
->Dimension
], tmp
);
294 for (int i
= 0; i
< P
->NbRays
; ++i
) {
295 if (value_notone_p(P
->Ray
[i
][0]))
297 if (value_zero_p(P
->Ray
[i
][1+P
->Dimension
]))
299 Vector_Combine(P
->Ray
[i
]+1, offset
->p
, P
->Ray
[i
]+1,
300 P
->Ray
[i
][0], P
->Ray
[i
][1+P
->Dimension
], P
->Dimension
);
305 void gen_fun::shift(const vec_ZZ
& offset
)
307 for (int i
= 0; i
< term
.size(); ++i
)
308 for (int j
= 0; j
< term
[i
]->n
.power
.NumRows(); ++j
)
309 term
[i
]->n
.power
[j
] += offset
;
311 Vector
*v
= Vector_Alloc(offset
.length());
312 zz2values(offset
, v
->p
);
313 Polyhedron_Shift(context
, v
);
317 static void print_power(std::ostream
& os
, QQ
& c
, vec_ZZ
& p
,
318 unsigned int nparam
, char **param_name
)
322 for (int i
= 0; i
< p
.length(); ++i
) {
326 if (c
.n
== -1 && c
.d
== 1)
328 else if (c
.n
!= 1 || c
.d
!= 1) {
344 os
<< "^(" << p
[i
] << ")";
355 void gen_fun::print(std::ostream
& os
, unsigned int nparam
, char **param_name
) const
358 for (int i
= 0; i
< term
.size(); ++i
) {
362 for (int j
= 0; j
< term
[i
]->n
.coeff
.length(); ++j
) {
363 if (j
!= 0 && term
[i
]->n
.coeff
[j
].n
> 0)
365 print_power(os
, term
[i
]->n
.coeff
[j
], term
[i
]->n
.power
[j
],
369 for (int j
= 0; j
< term
[i
]->d
.power
.NumRows(); ++j
) {
373 print_power(os
, mone
, term
[i
]->d
.power
[j
], nparam
, param_name
);
380 gen_fun::operator evalue
*() const
384 value_init(factor
.d
);
385 value_init(factor
.x
.n
);
386 for (int i
= 0; i
< term
.size(); ++i
) {
387 unsigned nvar
= term
[i
]->d
.power
.NumRows();
388 unsigned nparam
= term
[i
]->d
.power
.NumCols();
389 Matrix
*C
= Matrix_Alloc(nparam
+ nvar
, 1 + nvar
+ nparam
+ 1);
390 mat_ZZ
& d
= term
[i
]->d
.power
;
391 Polyhedron
*U
= context
? context
: Universe_Polyhedron(nparam
);
393 for (int j
= 0; j
< term
[i
]->n
.coeff
.length(); ++j
) {
394 for (int r
= 0; r
< nparam
; ++r
) {
395 value_set_si(C
->p
[r
][0], 0);
396 for (int c
= 0; c
< nvar
; ++c
) {
397 zz2value(d
[c
][r
], C
->p
[r
][1+c
]);
399 Vector_Set(&C
->p
[r
][1+nvar
], 0, nparam
);
400 value_set_si(C
->p
[r
][1+nvar
+r
], -1);
401 zz2value(term
[i
]->n
.power
[j
][r
], C
->p
[r
][1+nvar
+nparam
]);
403 for (int r
= 0; r
< nvar
; ++r
) {
404 value_set_si(C
->p
[nparam
+r
][0], 1);
405 Vector_Set(&C
->p
[nparam
+r
][1], 0, nvar
+ nparam
+ 1);
406 value_set_si(C
->p
[nparam
+r
][1+r
], 1);
408 Polyhedron
*P
= Constraints2Polyhedron(C
, 0);
409 evalue
*E
= barvinok_enumerate_ev(P
, U
, 0);
411 if (EVALUE_IS_ZERO(*E
)) {
416 zz2value(term
[i
]->n
.coeff
[j
].n
, factor
.x
.n
);
417 zz2value(term
[i
]->n
.coeff
[j
].d
, factor
.d
);
420 Matrix_Print(stdout, P_VALUE_FMT, C);
421 char *test[] = { "A", "B", "C", "D", "E", "F", "G" };
422 print_evalue(stdout, E, test);
436 value_clear(factor
.d
);
437 value_clear(factor
.x
.n
);
441 void gen_fun::coefficient(Value
* params
, Value
* c
) const
443 if (context
&& !in_domain(context
, params
)) {
450 value_init(part
.x
.n
);
453 evalue_set_si(&sum
, 0, 1);
457 for (int i
= 0; i
< term
.size(); ++i
) {
458 unsigned nvar
= term
[i
]->d
.power
.NumRows();
459 unsigned nparam
= term
[i
]->d
.power
.NumCols();
460 Matrix
*C
= Matrix_Alloc(nparam
+ nvar
, 1 + nvar
+ 1);
461 mat_ZZ
& d
= term
[i
]->d
.power
;
463 for (int j
= 0; j
< term
[i
]->n
.coeff
.length(); ++j
) {
464 for (int r
= 0; r
< nparam
; ++r
) {
465 value_set_si(C
->p
[r
][0], 0);
466 for (int c
= 0; c
< nvar
; ++c
) {
467 zz2value(d
[c
][r
], C
->p
[r
][1+c
]);
469 zz2value(term
[i
]->n
.power
[j
][r
], C
->p
[r
][1+nvar
]);
470 value_subtract(C
->p
[r
][1+nvar
], C
->p
[r
][1+nvar
], params
[r
]);
472 for (int r
= 0; r
< nvar
; ++r
) {
473 value_set_si(C
->p
[nparam
+r
][0], 1);
474 Vector_Set(&C
->p
[nparam
+r
][1], 0, nvar
+ 1);
475 value_set_si(C
->p
[nparam
+r
][1+r
], 1);
477 Polyhedron
*P
= Constraints2Polyhedron(C
, 0);
482 barvinok_count(P
, &tmp
, 0);
484 if (value_zero_p(tmp
))
486 zz2value(term
[i
]->n
.coeff
[j
].n
, part
.x
.n
);
487 zz2value(term
[i
]->n
.coeff
[j
].d
, part
.d
);
488 value_multiply(part
.x
.n
, part
.x
.n
, tmp
);
494 assert(value_one_p(sum
.d
));
495 value_assign(*c
, sum
.x
.n
);
499 value_clear(part
.x
.n
);
501 value_clear(sum
.x
.n
);