?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
barvinok_count: postpone computation of dual if allowed by PolyLib.
[barvinok.git] / genfun.cc
blobfa0c0523fb4614a906c58681c075b402b91de132
1 #include <iostream>
2 #include <assert.h>
3 #include "config.h"
4 #include <genfun.h>
5 #include <barvinok.h>
6 #include <barvinok2.h>
7
8 using std::cout;
9
10 static int lex_cmp(vec_ZZ& a, vec_ZZ& b)
11 {
12 assert(a.length() == b.length());
13
14 for (int j = 0; j < a.length(); ++j)
15 if (a[j] != b[j])
16 return a[j] < b[j] ? -1 : 1;
17 return 0;
18 }
19
20 static int lex_cmp(mat_ZZ& a, mat_ZZ& b)
21 {
22 assert(a.NumCols() == b.NumCols());
23 int alen = a.NumRows();
24 int blen = b.NumRows();
25 int len = alen < blen ? alen : blen;
26
27 for (int i = 0; i < len; ++i) {
28 int s = lex_cmp(a[i], b[i]);
29 if (s)
30 return s;
31 }
32 return alen-blen;
33 }
34
35 void gen_fun::add(const ZZ& cn, const ZZ& cd, const vec_ZZ& num,
36 const mat_ZZ& den)
37 {
38 if (cn == 0)
39 return;
40
41 short_rat * r = new short_rat;
42 r->n.coeff.SetDims(1, 2);
43 r->n.coeff[0][0] = cn;
44 r->n.coeff[0][1] = cd;
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 term[i]->n.coeff[j][0] = n;
89 term[i]->n.coeff[j][1] = d;
90 } else {
91 if (len > 1) {
92 if (j < len-1) {
93 term[i]->n.power[j] = term[i]->n.power[len-1];
94 term[i]->n.coeff[j] = term[i]->n.coeff[len-1];
95 }
96 int dim = term[i]->n.power.NumCols();
97 term[i]->n.coeff.SetDims(len-1, 2);
98 term[i]->n.power.SetDims(len-1, dim);
99 } else {
100 delete term[i];
101 if (i != term.size()-1)
102 term[i] = term[term.size()-1];
103 term.pop_back();
104 }
105 }
106 } else {
107 int dim = term[i]->n.power.NumCols();
108 term[i]->n.coeff.SetDims(len+1, 2);
109 term[i]->n.power.SetDims(len+1, dim);
110 term[i]->n.coeff[len] = r->n.coeff[0];
111 term[i]->n.power[len] = r->n.power[0];
112 }
113 delete r;
114 return;
115 }
116
117 term.push_back(r);
118 }
119
120 static void print_power(vec_ZZ& c, vec_ZZ& p,
121 unsigned int nparam, char **param_name)
122 {
123 bool first = true;
124
125 for (int i = 0; i < p.length(); ++i) {
126 if (p[i] == 0)
127 continue;
128 if (first) {
129 if (c[0] == -1 && c[1] == 1)
130 cout << "-";
131 else if (c[0] != 1 || c[1] != 1) {
132 cout << c[0];
133 if (c[1] != 1)
134 cout << " / " << c[1];
135 cout << " * ";
136 }
137 first = false;
138 } else
139 cout << " * ";
140 if (i < nparam)
141 cout << param_name[i];
142 else
143 cout << "x" << i;
144 if (p[i] == 1)
145 continue;
146 if (p[i] < 0)
147 cout << "^(" << p[i] << ")";
148 else
149 cout << "^" << p[i];
150 }
151 if (first) {
152 cout << c[0];
153 if (c[1] != 1)
154 cout << " / " << c[1];
155 }
156 }
157
158 void gen_fun::print(unsigned int nparam, char **param_name) const
159 {
160 vec_ZZ mone;
161 mone.SetLength(2);
162 mone[0] = -1;
163 mone[1] = 1;
164 for (int i = 0; i < term.size(); ++i) {
165 if (i != 0)
166 cout << " + ";
167 cout << "(";
168 for (int j = 0; j < term[i]->n.coeff.NumRows(); ++j) {
169 if (j != 0 && term[i]->n.coeff[j][0] > 0)
170 cout << "+";
171 print_power(term[i]->n.coeff[j], term[i]->n.power[j],
172 nparam, param_name);
173 }
174 cout << ")/(";
175 for (int j = 0; j < term[i]->d.power.NumRows(); ++j) {
176 if (j != 0)
177 cout << " * ";
178 cout << "(1";
179 print_power(mone, term[i]->d.power[j],
180 nparam, param_name);
181 cout << ")";
182 }
183 cout << ")";
184 }
185 }
186
187 gen_fun::operator evalue *() const
188 {
189 evalue *EP = NULL;
190 evalue factor;
191 value_init(factor.d);
192 value_init(factor.x.n);
193 for (int i = 0; i < term.size(); ++i) {
194 unsigned nvar = term[i]->d.power.NumRows();
195 unsigned nparam = term[i]->d.power.NumCols();
196 Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + nparam + 1);
197 mat_ZZ& d = term[i]->d.power;
198 Polyhedron *U = context ? context : Universe_Polyhedron(nparam);
199
200 for (int j = 0; j < term[i]->n.coeff.NumRows(); ++j) {
201 for (int r = 0; r < nparam; ++r) {
202 value_set_si(C->p[r][0], 0);
203 for (int c = 0; c < nvar; ++c) {
204 zz2value(d[c][r], C->p[r][1+c]);
205 }
206 Vector_Set(&C->p[r][1+nvar], 0, nparam);
207 value_set_si(C->p[r][1+nvar+r], -1);
208 zz2value(term[i]->n.power[j][r], C->p[r][1+nvar+nparam]);
209 }
210 for (int r = 0; r < nvar; ++r) {
211 value_set_si(C->p[nparam+r][0], 1);
212 Vector_Set(&C->p[nparam+r][1], 0, nvar + nparam + 1);
213 value_set_si(C->p[nparam+r][1+r], 1);
214 }
215 Polyhedron *P = Constraints2Polyhedron(C, 0);
216 evalue *E = barvinok_enumerate_ev(P, U, 0);
217 Polyhedron_Free(P);
218 if (EVALUE_IS_ZERO(*E))
219 continue;
220 zz2value(term[i]->n.coeff[j][0], factor.x.n);
221 zz2value(term[i]->n.coeff[j][1], factor.d);
222 emul(&factor, E);
223 /*
224 Matrix_Print(stdout, P_VALUE_FMT, C);
225 char *test[] = { "A", "B", "C", "D", "E", "F", "G" };
226 print_evalue(stdout, E, test);
227 */
228 if (!EP)
229 EP = E;
230 else {
231 eadd(E, EP);
232 free_evalue_refs(E);
233 }
234 }
235 Matrix_Free(C);
236 if (!context)
237 Polyhedron_Free(U);
238 }
239 value_clear(factor.d);
240 value_clear(factor.x.n);
241 return EP;
242 }
243
244 void gen_fun::coefficient(Value* params, Value* c) const
245 {
246 if (context && !in_domain(context, params)) {
247 value_set_si(*c, 0);
248 return;
249 }
250
251 evalue part;
252 value_init(part.d);
253 value_init(part.x.n);
254 evalue sum;
255 value_init(sum.d);
256 evalue_set_si(&sum, 0, 1);
257 Value tmp;
258 value_init(tmp);
259
260 for (int i = 0; i < term.size(); ++i) {
261 unsigned nvar = term[i]->d.power.NumRows();
262 unsigned nparam = term[i]->d.power.NumCols();
263 Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + 1);
264 mat_ZZ& d = term[i]->d.power;
265
266 for (int j = 0; j < term[i]->n.coeff.NumRows(); ++j) {
267 for (int r = 0; r < nparam; ++r) {
268 value_set_si(C->p[r][0], 0);
269 for (int c = 0; c < nvar; ++c) {
270 zz2value(d[c][r], C->p[r][1+c]);
271 }
272 zz2value(term[i]->n.power[j][r], C->p[r][1+nvar]);
273 value_subtract(C->p[r][1+nvar], C->p[r][1+nvar], params[r]);
274 }
275 for (int r = 0; r < nvar; ++r) {
276 value_set_si(C->p[nparam+r][0], 1);
277 Vector_Set(&C->p[nparam+r][1], 0, nvar + 1);
278 value_set_si(C->p[nparam+r][1+r], 1);
279 }
280 Polyhedron *P = Constraints2Polyhedron(C, 0);
281 if (emptyQ(P)) {
282 Polyhedron_Free(P);
283 continue;
284 }
285 barvinok_count(P, &tmp, 0);
286 Polyhedron_Free(P);
287 if (value_zero_p(tmp))
288 continue;
289 zz2value(term[i]->n.coeff[j][0], part.x.n);
290 zz2value(term[i]->n.coeff[j][1], part.d);
291 value_multiply(part.x.n, part.x.n, tmp);
292 eadd(&part, &sum);
293 }
294 Matrix_Free(C);
295 }
296
297 assert(value_one_p(sum.d));
298 value_assign(*c, sum.x.n);
299
300 value_clear(tmp);
301 value_clear(part.d);
302 value_clear(part.x.n);
303 value_clear(sum.d);
304 value_clear(sum.x.n);
305 }
lattice point enumeration library