?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
gen_fun::Hadamard_product: optimize computation of terms with same denominator
[barvinok.git] / genfun.cc
blob674b6bb39c7c567dd06d14121f9c7976a78b3a6f
1 #include <iostream>
2 #include <vector>
3 #include <assert.h>
4 #include "config.h"
5 #include <barvinok/genfun.h>
6 #include <barvinok/barvinok.h>
7 #include "conversion.h"
8 #include "genfun_constructor.h"
9 #include "mat_util.h"
10
11 using std::cout;
12 using std::cerr;
13 using std::endl;
14 using std::pair;
15 using std::vector;
16
17 static int lex_cmp(mat_ZZ& a, mat_ZZ& b)
18 {
19 assert(a.NumCols() == b.NumCols());
20 int alen = a.NumRows();
21 int blen = b.NumRows();
22 int len = alen < blen ? alen : blen;
23
24 for (int i = 0; i < len; ++i) {
25 int s = lex_cmp(a[i], b[i]);
26 if (s)
27 return s;
28 }
29 return alen-blen;
30 }
31
32 static void lex_order_terms(struct short_rat* rat)
33 {
34 for (int i = 0; i < rat->n.power.NumRows(); ++i) {
35 int m = i;
36 for (int j = i+1; j < rat->n.power.NumRows(); ++j)
37 if (lex_cmp(rat->n.power[j], rat->n.power[m]) < 0)
38 m = j;
39 if (m != i) {
40 vec_ZZ tmp = rat->n.power[m];
41 rat->n.power[m] = rat->n.power[i];
42 rat->n.power[i] = tmp;
43 QQ tmp_coeff = rat->n.coeff[m];
44 rat->n.coeff[m] = rat->n.coeff[i];
45 rat->n.coeff[i] = tmp_coeff;
46 }
47 }
48 }
49
50 void short_rat::add(short_rat *r)
51 {
52 for (int i = 0; i < r->n.power.NumRows(); ++i) {
53 int len = n.coeff.length();
54 int j;
55 for (j = 0; j < len; ++j)
56 if (r->n.power[i] == n.power[j])
57 break;
58 if (j < len) {
59 n.coeff[j] += r->n.coeff[i];
60 if (n.coeff[j].n == 0) {
61 if (j < len-1) {
62 n.power[j] = n.power[len-1];
63 n.coeff[j] = n.coeff[len-1];
64 }
65 int dim = n.power.NumCols();
66 n.coeff.SetLength(len-1);
67 n.power.SetDims(len-1, dim);
68 }
69 } else {
70 int dim = n.power.NumCols();
71 n.coeff.SetLength(len+1);
72 n.power.SetDims(len+1, dim);
73 n.coeff[len] = r->n.coeff[i];
74 n.power[len] = r->n.power[i];
75 }
76 }
77 }
78
79 bool short_rat::reduced()
80 {
81 int dim = n.power.NumCols();
82 lex_order_terms(this);
83 if (n.power.NumRows() % 2 == 0) {
84 if (n.coeff[0].n == -n.coeff[1].n &&
85 n.coeff[0].d == n.coeff[1].d) {
86 vec_ZZ step = n.power[1] - n.power[0];
87 int k;
88 for (k = 1; k < n.power.NumRows()/2; ++k) {
89 if (n.coeff[2*k].n != -n.coeff[2*k+1].n ||
90 n.coeff[2*k].d != n.coeff[2*k+1].d)
91 break;
92 if (step != n.power[2*k+1] - n.power[2*k])
93 break;
94 }
95 if (k == n.power.NumRows()/2) {
96 for (k = 0; k < d.power.NumRows(); ++k)
97 if (d.power[k] == step)
98 break;
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];
106 }
107 n.coeff.SetLength(k);
108 n.power.SetDims(k, dim);
109 return true;
110 }
111 }
112 }
113 }
114 return false;
115 }
116
117 void gen_fun::add(const QQ& c, const vec_ZZ& num, const mat_ZZ& den)
118 {
119 if (c.n == 0)
120 return;
121
122 short_rat * r = new short_rat;
123 r->n.coeff.SetLength(1);
124 ZZ g = GCD(c.n, c.d);
125 r->n.coeff[0].n = c.n/g;
126 r->n.coeff[0].d = c.d/g;
127 r->n.power.SetDims(1, num.length());
128 r->n.power[0] = num;
129 r->d.power = den;
130
131 /* Make all powers in denominator lexico-positive */
132 for (int i = 0; i < r->d.power.NumRows(); ++i) {
133 int j;
134 for (j = 0; j < r->d.power.NumCols(); ++j)
135 if (r->d.power[i][j] != 0)
136 break;
137 if (r->d.power[i][j] < 0) {
138 r->d.power[i] = -r->d.power[i];
139 r->n.coeff[0].n = -r->n.coeff[0].n;
140 r->n.power[0] += r->d.power[i];
141 }
142 }
143
144 /* Order powers in denominator */
145 lex_order_rows(r->d.power);
146
147 for (int i = 0; i < term.size(); ++i)
148 if (lex_cmp(term[i]->d.power, r->d.power) == 0) {
149 term[i]->add(r);
150 if (term[i]->n.coeff.length() == 0) {
151 delete term[i];
152 if (i != term.size()-1)
153 term[i] = term[term.size()-1];
154 term.pop_back();
155 } else if (term[i]->reduced()) {
156 delete r;
157 /* we've modified term[i], so removed it
158 * and add it back again
159 */
160 r = term[i];
161 if (i != term.size()-1)
162 term[i] = term[term.size()-1];
163 term.pop_back();
164 i = -1;
165 continue;
166 }
167 delete r;
168 return;
169 }
170
171 term.push_back(r);
172 }
173
174 void gen_fun::add(const QQ& c, const gen_fun *gf)
175 {
176 QQ p;
177 for (int i = 0; i < gf->term.size(); ++i) {
178 for (int j = 0; j < gf->term[i]->n.power.NumRows(); ++j) {
179 p = c;
180 p *= gf->term[i]->n.coeff[j];
181 add(p, gf->term[i]->n.power[j], gf->term[i]->d.power);
182 }
183 }
184 }
185
186 /*
187 * Perform the substitution specified by CP and (map, offset)
188 *
189 * CP is a homogeneous matrix that maps a set of "compressed parameters"
190 * to the original set of parameters.
191 *
192 * This function is applied to a gen_fun computed with the compressed parameters
193 * and adapts it to refer to the original parameters.
194 *
195 * That is, if y are the compressed parameters and x = A y + b are the original
196 * parameters, then we want the coefficient of the monomial t^y in the original
197 * generating function to be the coefficient of the monomial u^x in the resulting
198 * generating function.
199 * The original generating function has the form
200 *
201 * a t^m/(1-t^n) = a t^m + a t^{m+n} + a t^{m+2n} + ...
202 *
203 * Since each term t^y should correspond to a term u^x, with x = A y + b, we want
204 *
205 * a u^{A m + b} + a u^{A (m+n) + b} + a u^{A (m+2n) +b} + ... =
206 *
207 * = a u^{A m + b}/(1-u^{A n})
208 *
209 * Therefore, we multiply the powers m and n in both numerator and denominator by A
210 * and add b to the power in the numerator.
211 * Since the above powers are stored as row vectors m^T and n^T,
212 * we compute, say, m'^T = m^T A^T to obtain m' = A m.
213 *
214 * The pair (map, offset) contains the same information as CP.
215 * map is the transpose of the linear part of CP, while offset is the constant part.
216 */
217 void gen_fun::substitute(Matrix *CP, const mat_ZZ& map, const vec_ZZ& offset)
218 {
219 Polyhedron *C = Polyhedron_Image(context, CP, 0);
220 Polyhedron_Free(context);
221 context = C;
222 for (int i = 0; i < term.size(); ++i) {
223 term[i]->d.power *= map;
224 term[i]->n.power *= map;
225 for (int j = 0; j < term[i]->n.power.NumRows(); ++j)
226 term[i]->n.power[j] += offset;
227 }
228 }
229
230 struct cone {
231 int *pos;
232 vector<pair<Vector *, QQ> > vertices;
233 cone(int *pos) : pos(pos) {}
234 };
235
236 struct parallel_polytopes {
237 gf_base *red;
238 Matrix *Constraints;
239 int dim;
240 int nparam;
241 vector<cone> cones;
242
243 parallel_polytopes(int n, Polyhedron *context, int dim, int nparam) :
244 dim(dim), nparam(nparam) {
245 red = gf_base::create(Polyhedron_Copy(context), dim, nparam);
246 Constraints = NULL;
247 }
248 void add(const QQ& c, Polyhedron *P, unsigned MaxRays) {
249 Polyhedron *Q = remove_equalities_p(Polyhedron_Copy(P), P->Dimension-nparam,
250 NULL);
251 assert(Q->Dimension == dim);
252
253 POL_ENSURE_VERTICES(Q);
254 if (emptyQ(Q)) {
255 Polyhedron_Free(Q);
256 return;
257 }
258
259 if (First_Non_Zero(Q->Constraint[Q->NbConstraints-1]+1, Q->Dimension) == -1)
260 Q->NbConstraints--;
261
262 if (!Constraints) {
263 red->base->init(Q);
264 Constraints = Matrix_Alloc(Q->NbConstraints, Q->Dimension);
265 for (int i = 0; i < Q->NbConstraints; ++i) {
266 Vector_Copy(Q->Constraint[i]+1, Constraints->p[i], Q->Dimension);
267 }
268 } else {
269 for (int i = 0; i < Q->NbConstraints; ++i) {
270 int j;
271 for (j = 0; j < Constraints->NbRows; ++j)
272 if (Vector_Equal(Q->Constraint[i]+1, Constraints->p[j],
273 Q->Dimension))
274 break;
275 assert(j < Constraints->NbRows);
276 }
277 }
278
279 for (int i = 0; i < Q->NbRays; ++i) {
280 if (!value_pos_p(Q->Ray[i][dim+1]))
281 continue;
282
283 Polyhedron *C = supporting_cone(Q, i);
284
285 if (First_Non_Zero(C->Constraint[C->NbConstraints-1]+1,
286 C->Dimension) == -1)
287 C->NbConstraints--;
288
289 int *pos = new int[1+C->NbConstraints];
290 pos[0] = C->NbConstraints;
291 int l = 0;
292 for (int k = 0; k < Constraints->NbRows; ++k) {
293 for (int j = 0; j < C->NbConstraints; ++j) {
294 if (Vector_Equal(C->Constraint[j]+1, Constraints->p[k],
295 C->Dimension)) {
296 pos[1+l++] = k;
297 break;
298 }
299 }
300 }
301 assert(l == C->NbConstraints);
302
303 int j;
304 for (j = 0; j < cones.size(); ++j)
305 if (!memcmp(pos, cones[j].pos, (1+C->NbConstraints)*sizeof(int)))
306 break;
307 if (j == cones.size())
308 cones.push_back(cone(pos));
309 else
310 delete [] pos;
311
312 Polyhedron_Free(C);
313
314 int k;
315 for (k = 0; k < cones[j].vertices.size(); ++k)
316 if (Vector_Equal(Q->Ray[i]+1, cones[j].vertices[k].first->p,
317 Q->Dimension+1))
318 break;
319
320 if (k == cones[j].vertices.size()) {
321 Vector *vertex = Vector_Alloc(Q->Dimension+1);
322 Vector_Copy(Q->Ray[i]+1, vertex->p, Q->Dimension+1);
323 cones[j].vertices.push_back(pair<Vector*,QQ>(vertex, c));
324 } else {
325 cones[j].vertices[k].second += c;
326 if (cones[j].vertices[k].second.n == 0) {
327 int size = cones[j].vertices.size();
328 Vector_Free(cones[j].vertices[k].first);
329 if (k < size-1)
330 cones[j].vertices[k] = cones[j].vertices[size-1];
331 cones[j].vertices.pop_back();
332 }
333 }
334 }
335
336 Polyhedron_Free(Q);
337 }
338 gen_fun *compute(unsigned MaxRays) {
339 for (int i = 0; i < cones.size(); ++i) {
340 Matrix *M = Matrix_Alloc(cones[i].pos[0], 1+Constraints->NbColumns+1);
341 Polyhedron *Cone;
342 for (int j = 0; j <cones[i].pos[0]; ++j) {
343 value_set_si(M->p[j][0], 1);
344 Vector_Copy(Constraints->p[cones[i].pos[1+j]], M->p[j]+1,
345 Constraints->NbColumns);
346 }
347 Cone = Constraints2Polyhedron(M, MaxRays);
348 Matrix_Free(M);
349 for (int j = 0; j < cones[i].vertices.size(); ++j) {
350 red->base->do_vertex_cone(cones[i].vertices[j].second,
351 Polyhedron_Copy(Cone),
352 cones[i].vertices[j].first->p,
353 MaxRays);
354 }
355 Polyhedron_Free(Cone);
356 }
357 return red->gf;
358 }
359 void print(std::ostream& os) const {
360 for (int i = 0; i < cones.size(); ++i) {
361 os << "[";
362 for (int j = 0; j < cones[i].pos[0]; ++j) {
363 if (j)
364 os << ", ";
365 os << cones[i].pos[1+j];
366 }
367 os << "]" << endl;
368 for (int j = 0; j < cones[i].vertices.size(); ++j) {
369 Vector_Print(stderr, P_VALUE_FMT, cones[i].vertices[j].first);
370 os << cones[i].vertices[j].second << endl;
371 }
372 }
373 }
374 ~parallel_polytopes() {
375 for (int i = 0; i < cones.size(); ++i) {
376 delete [] cones[i].pos;
377 for (int j = 0; j < cones[i].vertices.size(); ++j)
378 Vector_Free(cones[i].vertices[j].first);
379 }
380 if (Constraints)
381 Matrix_Free(Constraints);
382 delete red;
383 }
384 };
385
386 gen_fun *gen_fun::Hadamard_product(const gen_fun *gf, unsigned MaxRays)
387 {
388 QQ one(1, 1);
389 Polyhedron *C = DomainIntersection(context, gf->context, MaxRays);
390 Polyhedron *U = Universe_Polyhedron(C->Dimension);
391 gen_fun *sum = new gen_fun(C);
392 for (int i = 0; i < term.size(); ++i) {
393 for (int i2 = 0; i2 < gf->term.size(); ++i2) {
394 int d = term[i]->d.power.NumCols();
395 int k1 = term[i]->d.power.NumRows();
396 int k2 = gf->term[i2]->d.power.NumRows();
397 assert(term[i]->d.power.NumCols() == gf->term[i2]->d.power.NumCols());
398
399 parallel_polytopes pp(term[i]->n.power.NumRows() *
400 gf->term[i2]->n.power.NumRows(),
401 sum->context, k1+k2-d, d);
402
403 for (int j = 0; j < term[i]->n.power.NumRows(); ++j) {
404 for (int j2 = 0; j2 < gf->term[i2]->n.power.NumRows(); ++j2) {
405 Matrix *M = Matrix_Alloc(k1+k2+d+d, 1+k1+k2+d+1);
406 for (int k = 0; k < k1+k2; ++k) {
407 value_set_si(M->p[k][0], 1);
408 value_set_si(M->p[k][1+k], 1);
409 }
410 for (int k = 0; k < d; ++k) {
411 value_set_si(M->p[k1+k2+k][1+k1+k2+k], -1);
412 zz2value(term[i]->n.power[j][k], M->p[k1+k2+k][1+k1+k2+d]);
413 for (int l = 0; l < k1; ++l)
414 zz2value(term[i]->d.power[l][k], M->p[k1+k2+k][1+l]);
415 }
416 for (int k = 0; k < d; ++k) {
417 value_set_si(M->p[k1+k2+d+k][1+k1+k2+k], -1);
418 zz2value(gf->term[i2]->n.power[j2][k],
419 M->p[k1+k2+d+k][1+k1+k2+d]);
420 for (int l = 0; l < k2; ++l)
421 zz2value(gf->term[i2]->d.power[l][k],
422 M->p[k1+k2+d+k][1+k1+l]);
423 }
424 Polyhedron *P = Constraints2Polyhedron(M, MaxRays);
425 Matrix_Free(M);
426
427 QQ c = term[i]->n.coeff[j];
428 c *= gf->term[i2]->n.coeff[j2];
429 pp.add(c, P, MaxRays);
430
431 Polyhedron_Free(P);
432 }
433 }
434
435 gen_fun *t = pp.compute(MaxRays);
436 sum->add(one, t);
437 delete t;
438 }
439 }
440 Polyhedron_Free(U);
441 return sum;
442 }
443
444 void gen_fun::add_union(gen_fun *gf, unsigned MaxRays)
445 {
446 QQ one(1, 1), mone(-1, 1);
447
448 gen_fun *hp = Hadamard_product(gf, MaxRays);
449 add(one, gf);
450 add(mone, hp);
451 delete hp;
452 }
453
454 static void Polyhedron_Shift(Polyhedron *P, Vector *offset)
455 {
456 Value tmp;
457 value_init(tmp);
458 for (int i = 0; i < P->NbConstraints; ++i) {
459 Inner_Product(P->Constraint[i]+1, offset->p, P->Dimension, &tmp);
460 value_subtract(P->Constraint[i][1+P->Dimension],
461 P->Constraint[i][1+P->Dimension], tmp);
462 }
463 for (int i = 0; i < P->NbRays; ++i) {
464 if (value_notone_p(P->Ray[i][0]))
465 continue;
466 if (value_zero_p(P->Ray[i][1+P->Dimension]))
467 continue;
468 Vector_Combine(P->Ray[i]+1, offset->p, P->Ray[i]+1,
469 P->Ray[i][0], P->Ray[i][1+P->Dimension], P->Dimension);
470 }
471 value_clear(tmp);
472 }
473
474 void gen_fun::shift(const vec_ZZ& offset)
475 {
476 for (int i = 0; i < term.size(); ++i)
477 for (int j = 0; j < term[i]->n.power.NumRows(); ++j)
478 term[i]->n.power[j] += offset;
479
480 Vector *v = Vector_Alloc(offset.length());
481 zz2values(offset, v->p);
482 Polyhedron_Shift(context, v);
483 Vector_Free(v);
484 }
485
486 static void print_power(std::ostream& os, QQ& c, vec_ZZ& p,
487 unsigned int nparam, char **param_name)
488 {
489 bool first = true;
490
491 for (int i = 0; i < p.length(); ++i) {
492 if (p[i] == 0)
493 continue;
494 if (first) {
495 if (c.n == -1 && c.d == 1)
496 os << "-";
497 else if (c.n != 1 || c.d != 1) {
498 os << c.n;
499 if (c.d != 1)
500 os << " / " << c.d;
501 os << "*";
502 }
503 first = false;
504 } else
505 os << "*";
506 if (i < nparam)
507 os << param_name[i];
508 else
509 os << "x" << i;
510 if (p[i] == 1)
511 continue;
512 if (p[i] < 0)
513 os << "^(" << p[i] << ")";
514 else
515 os << "^" << p[i];
516 }
517 if (first) {
518 os << c.n;
519 if (c.d != 1)
520 os << " / " << c.d;
521 }
522 }
523
524 void gen_fun::print(std::ostream& os, unsigned int nparam, char **param_name) const
525 {
526 QQ mone(-1, 1);
527 for (int i = 0; i < term.size(); ++i) {
528 if (i != 0)
529 os << " + ";
530 os << "(";
531 for (int j = 0; j < term[i]->n.coeff.length(); ++j) {
532 if (j != 0 && term[i]->n.coeff[j].n > 0)
533 os << "+";
534 print_power(os, term[i]->n.coeff[j], term[i]->n.power[j],
535 nparam, param_name);
536 }
537 os << ")/(";
538 for (int j = 0; j < term[i]->d.power.NumRows(); ++j) {
539 if (j != 0)
540 os << " * ";
541 os << "(1";
542 print_power(os, mone, term[i]->d.power[j], nparam, param_name);
543 os << ")";
544 }
545 os << ")";
546 }
547 }
548
549 gen_fun::operator evalue *() const
550 {
551 evalue *EP = NULL;
552 evalue factor;
553 value_init(factor.d);
554 value_init(factor.x.n);
555 for (int i = 0; i < term.size(); ++i) {
556 unsigned nvar = term[i]->d.power.NumRows();
557 unsigned nparam = term[i]->d.power.NumCols();
558 Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + nparam + 1);
559 mat_ZZ& d = term[i]->d.power;
560 Polyhedron *U = context ? context : Universe_Polyhedron(nparam);
561
562 for (int j = 0; j < term[i]->n.coeff.length(); ++j) {
563 for (int r = 0; r < nparam; ++r) {
564 value_set_si(C->p[r][0], 0);
565 for (int c = 0; c < nvar; ++c) {
566 zz2value(d[c][r], C->p[r][1+c]);
567 }
568 Vector_Set(&C->p[r][1+nvar], 0, nparam);
569 value_set_si(C->p[r][1+nvar+r], -1);
570 zz2value(term[i]->n.power[j][r], C->p[r][1+nvar+nparam]);
571 }
572 for (int r = 0; r < nvar; ++r) {
573 value_set_si(C->p[nparam+r][0], 1);
574 Vector_Set(&C->p[nparam+r][1], 0, nvar + nparam + 1);
575 value_set_si(C->p[nparam+r][1+r], 1);
576 }
577 Polyhedron *P = Constraints2Polyhedron(C, 0);
578 evalue *E = barvinok_enumerate_ev(P, U, 0);
579 Polyhedron_Free(P);
580 if (EVALUE_IS_ZERO(*E)) {
581 free_evalue_refs(E);
582 free(E);
583 continue;
584 }
585 zz2value(term[i]->n.coeff[j].n, factor.x.n);
586 zz2value(term[i]->n.coeff[j].d, factor.d);
587 emul(&factor, E);
588 /*
589 Matrix_Print(stdout, P_VALUE_FMT, C);
590 char *test[] = { "A", "B", "C", "D", "E", "F", "G" };
591 print_evalue(stdout, E, test);
592 */
593 if (!EP)
594 EP = E;
595 else {
596 eadd(E, EP);
597 free_evalue_refs(E);
598 free(E);
599 }
600 }
601 Matrix_Free(C);
602 if (!context)
603 Polyhedron_Free(U);
604 }
605 value_clear(factor.d);
606 value_clear(factor.x.n);
607 return EP;
608 }
609
610 void gen_fun::coefficient(Value* params, Value* c) const
611 {
612 if (context && !in_domain(context, params)) {
613 value_set_si(*c, 0);
614 return;
615 }
616
617 evalue part;
618 value_init(part.d);
619 value_init(part.x.n);
620 evalue sum;
621 value_init(sum.d);
622 evalue_set_si(&sum, 0, 1);
623 Value tmp;
624 value_init(tmp);
625
626 for (int i = 0; i < term.size(); ++i) {
627 unsigned nvar = term[i]->d.power.NumRows();
628 unsigned nparam = term[i]->d.power.NumCols();
629 Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + 1);
630 mat_ZZ& d = term[i]->d.power;
631
632 for (int j = 0; j < term[i]->n.coeff.length(); ++j) {
633 for (int r = 0; r < nparam; ++r) {
634 value_set_si(C->p[r][0], 0);
635 for (int c = 0; c < nvar; ++c) {
636 zz2value(d[c][r], C->p[r][1+c]);
637 }
638 zz2value(term[i]->n.power[j][r], C->p[r][1+nvar]);
639 value_subtract(C->p[r][1+nvar], C->p[r][1+nvar], params[r]);
640 }
641 for (int r = 0; r < nvar; ++r) {
642 value_set_si(C->p[nparam+r][0], 1);
643 Vector_Set(&C->p[nparam+r][1], 0, nvar + 1);
644 value_set_si(C->p[nparam+r][1+r], 1);
645 }
646 Polyhedron *P = Constraints2Polyhedron(C, 0);
647 if (emptyQ(P)) {
648 Polyhedron_Free(P);
649 continue;
650 }
651 barvinok_count(P, &tmp, 0);
652 Polyhedron_Free(P);
653 if (value_zero_p(tmp))
654 continue;
655 zz2value(term[i]->n.coeff[j].n, part.x.n);
656 zz2value(term[i]->n.coeff[j].d, part.d);
657 value_multiply(part.x.n, part.x.n, tmp);
658 eadd(&part, &sum);
659 }
660 Matrix_Free(C);
661 }
662
663 assert(value_one_p(sum.d));
664 value_assign(*c, sum.x.n);
665
666 value_clear(tmp);
667 value_clear(part.d);
668 value_clear(part.x.n);
669 value_clear(sum.d);
670 value_clear(sum.x.n);
671 }
lattice point enumeration library