?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
omega/convert.cc: add conversion from PolyLib to Omega
[barvinok.git] / genfun.cc
blobe4bfff9a346f9f34f7f75a4432e200913d33ecf6
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 r->n.coeff[0][0] = cn;
43 r->n.coeff[0][1] = cd;
44 r->n.power.SetDims(1, num.length());
45 r->n.power[0] = num;
46 r->d.power = den;
47
48 /* Make all powers in denominator lexico-positive */
49 for (int i = 0; i < r->d.power.NumRows(); ++i) {
50 int j;
51 for (j = 0; j < r->d.power.NumCols(); ++j)
52 if (r->d.power[i][j] != 0)
53 break;
54 if (r->d.power[i][j] < 0) {
55 r->d.power[i] = -r->d.power[i];
56 r->n.coeff[0][0] = -r->n.coeff[0][0];
57 r->n.power[0] += r->d.power[i];
58 }
59 }
60
61 /* Order powers in denominator */
62 for (int i = 0; i < r->d.power.NumRows(); ++i) {
63 int m = i;
64 for (int j = i+1; j < r->d.power.NumRows(); ++j)
65 if (lex_cmp(r->d.power[j], r->d.power[m]) < 0)
66 m = j;
67 if (m != i) {
68 vec_ZZ tmp = r->d.power[m];
69 r->d.power[m] = r->d.power[i];
70 r->d.power[i] = tmp;
71 }
72 }
73
74 for (int i = 0; i < term.size(); ++i)
75 if (lex_cmp(term[i]->d.power, r->d.power) == 0) {
76 int len = term[i]->n.coeff.NumRows();
77 int j;
78 for (j = 0; j < len; ++j)
79 if (r->n.power[0] == term[i]->n.power[j])
80 break;
81 if (j < len) {
82 ZZ g = GCD(r->n.coeff[0][1], term[i]->n.coeff[j][1]);
83 ZZ n = term[i]->n.coeff[j][0] * (r->n.coeff[0][1] / g) +
84 (term[i]->n.coeff[j][1] / g) * r->n.coeff[0][0];
85 ZZ d = term[i]->n.coeff[j][1] / g * r->n.coeff[0][1];
86 if (n != 0) {
87 term[i]->n.coeff[j][0] = n;
88 term[i]->n.coeff[j][1] = d;
89 } else {
90 if (len > 1) {
91 if (j < len-1) {
92 term[i]->n.power[j] = term[i]->n.power[len-1];
93 term[i]->n.coeff[j] = term[i]->n.coeff[len-1];
94 }
95 int dim = term[i]->n.power.NumCols();
96 term[i]->n.coeff.SetDims(len-1, 2);
97 term[i]->n.power.SetDims(len-1, dim);
98 } else {
99 delete term[i];
100 if (i != term.size()-1)
101 term[i] = term[term.size()-1];
102 term.pop_back();
103 }
104 }
105 } else {
106 int dim = term[i]->n.power.NumCols();
107 term[i]->n.coeff.SetDims(len+1, 2);
108 term[i]->n.power.SetDims(len+1, dim);
109 term[i]->n.coeff[len] = r->n.coeff[0];
110 term[i]->n.power[len] = r->n.power[0];
111 }
112 delete r;
113 return;
114 }
115
116 term.push_back(r);
117 }
118
119 static void print_power(vec_ZZ& c, vec_ZZ& p,
120 unsigned int nparam, char **param_name)
121 {
122 bool first = true;
123
124 for (int i = 0; i < p.length(); ++i) {
125 if (p[i] == 0)
126 continue;
127 if (first) {
128 if (c[0] == -1 && c[1] == 1)
129 cout << "-";
130 else if (c[0] != 1 || c[1] != 1) {
131 cout << c[0];
132 if (c[1] != 1)
133 cout << " / " << c[1];
134 cout << " * ";
135 }
136 first = false;
137 } else
138 cout << " * ";
139 if (i < nparam)
140 cout << param_name[i];
141 else
142 cout << "x" << i;
143 if (p[i] == 1)
144 continue;
145 if (p[i] < 0)
146 cout << "^(" << p[i] << ")";
147 else
148 cout << "^" << p[i];
149 }
150 if (first) {
151 cout << c[0];
152 if (c[1] != 1)
153 cout << " / " << c[1];
154 }
155 }
156
157 void gen_fun::print(unsigned int nparam, char **param_name) const
158 {
159 vec_ZZ mone;
160 mone.SetLength(2);
161 mone[0] = -1;
162 mone[1] = 1;
163 for (int i = 0; i < term.size(); ++i) {
164 if (i != 0)
165 cout << " + ";
166 cout << "(";
167 for (int j = 0; j < term[i]->n.coeff.NumRows(); ++j) {
168 if (j != 0 && term[i]->n.coeff[j][0] > 0)
169 cout << "+";
170 print_power(term[i]->n.coeff[j], term[i]->n.power[j],
171 nparam, param_name);
172 }
173 cout << ")/(";
174 for (int j = 0; j < term[i]->d.power.NumRows(); ++j) {
175 if (j != 0)
176 cout << " * ";
177 cout << "(1";
178 print_power(mone, term[i]->d.power[j],
179 nparam, param_name);
180 cout << ")";
181 }
182 cout << ")";
183 }
184 }
185
186 gen_fun::operator evalue *() const
187 {
188 evalue *EP = NULL;
189 evalue factor;
190 value_init(factor.d);
191 value_init(factor.x.n);
192 for (int i = 0; i < term.size(); ++i) {
193 unsigned nvar = term[i]->d.power.NumRows();
194 unsigned nparam = term[i]->d.power.NumCols();
195 Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + nparam + 1);
196 mat_ZZ& d = term[i]->d.power;
197 Polyhedron *U = context ? context : Universe_Polyhedron(nparam);
198
199 for (int j = 0; j < term[i]->n.coeff.NumRows(); ++j) {
200 for (int r = 0; r < nparam; ++r) {
201 value_set_si(C->p[r][0], 0);
202 for (int c = 0; c < nvar; ++c) {
203 zz2value(d[c][r], C->p[r][1+c]);
204 }
205 Vector_Set(&C->p[r][1+nvar], 0, nparam);
206 value_set_si(C->p[r][1+nvar+r], -1);
207 zz2value(term[i]->n.power[j][r], C->p[r][1+nvar+nparam]);
208 }
209 for (int r = 0; r < nvar; ++r) {
210 value_set_si(C->p[nparam+r][0], 1);
211 Vector_Set(&C->p[nparam+r][1], 0, nvar + nparam + 1);
212 value_set_si(C->p[nparam+r][1+r], 1);
213 }
214 Polyhedron *P = Constraints2Polyhedron(C, 0);
215 evalue *E = barvinok_enumerate_ev(P, U, 0);
216 Polyhedron_Free(P);
217 if (EVALUE_IS_ZERO(*E))
218 continue;
219 zz2value(term[i]->n.coeff[j][0], factor.x.n);
220 zz2value(term[i]->n.coeff[j][1], factor.d);
221 emul(&factor, E);
222 /*
223 Matrix_Print(stdout, P_VALUE_FMT, C);
224 char *test[] = { "A", "B", "C", "D", "E", "F", "G" };
225 print_evalue(stdout, E, test);
226 */
227 if (!EP)
228 EP = E;
229 else {
230 eadd(E, EP);
231 free_evalue_refs(E);
232 }
233 }
234 Matrix_Free(C);
235 if (!context)
236 Polyhedron_Free(U);
237 }
238 value_clear(factor.d);
239 value_clear(factor.x.n);
240 return EP;
241 }
242
243 void gen_fun::coefficient(Value* params, Value* c) const
244 {
245 if (context && !in_domain(context, params)) {
246 value_set_si(*c, 0);
247 return;
248 }
249
250 evalue part;
251 value_init(part.d);
252 value_init(part.x.n);
253 evalue sum;
254 value_init(sum.d);
255 evalue_set_si(&sum, 0, 1);
256 Value tmp;
257 value_init(tmp);
258
259 for (int i = 0; i < term.size(); ++i) {
260 unsigned nvar = term[i]->d.power.NumRows();
261 unsigned nparam = term[i]->d.power.NumCols();
262 Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + 1);
263 mat_ZZ& d = term[i]->d.power;
264
265 for (int j = 0; j < term[i]->n.coeff.NumRows(); ++j) {
266 for (int r = 0; r < nparam; ++r) {
267 value_set_si(C->p[r][0], 0);
268 for (int c = 0; c < nvar; ++c) {
269 zz2value(d[c][r], C->p[r][1+c]);
270 }
271 zz2value(term[i]->n.power[j][r], C->p[r][1+nvar]);
272 value_subtract(C->p[r][1+nvar], C->p[r][1+nvar], params[r]);
273 }
274 for (int r = 0; r < nvar; ++r) {
275 value_set_si(C->p[nparam+r][0], 1);
276 Vector_Set(&C->p[nparam+r][1], 0, nvar + 1);
277 value_set_si(C->p[nparam+r][1+r], 1);
278 }
279 Polyhedron *P = Constraints2Polyhedron(C, 0);
280 if (emptyQ(P)) {
281 Polyhedron_Free(P);
282 continue;
283 }
284 barvinok_count(P, &tmp, 0);
285 Polyhedron_Free(P);
286 if (value_zero_p(tmp))
287 continue;
288 zz2value(term[i]->n.coeff[j][0], part.x.n);
289 zz2value(term[i]->n.coeff[j][1], part.d);
290 value_multiply(part.x.n, part.x.n, tmp);
291 eadd(&part, &sum);
292 }
293 Matrix_Free(C);
294 }
295
296 assert(value_one_p(sum.d));
297 value_assign(*c, sum.x.n);
298
299 value_clear(tmp);
300 value_clear(part.d);
301 value_clear(part.x.n);
302 value_clear(sum.d);
303 value_clear(sum.x.n);
304 }
lattice point enumeration library