iscc: call_bin_op: handle binary test operations
[barvinok.git] / lattice_point.cc
blob4b0528fae9c54634ab0e2838d440cda91c89e64c
1 #include <assert.h>
2 #include <NTL/mat_ZZ.h>
3 #include <NTL/vec_ZZ.h>
4 #include <barvinok/barvinok.h>
5 #include <barvinok/evalue.h>
6 #include <barvinok/util.h>
7 #include "config.h"
8 #include "conversion.h"
9 #include "lattice_point.h"
10 #include "param_util.h"
12 using std::cerr;
13 using std::endl;
15 #define ALLOC(type) (type*)malloc(sizeof(type))
17 /* returns an evalue that corresponds to
19 * c/(*den) x_param
21 static evalue *term(int param, ZZ& c, Value *den = NULL)
23 evalue *EP = new evalue();
24 value_init(EP->d);
25 value_set_si(EP->d,0);
26 EP->x.p = new_enode(polynomial, 2, param + 1);
27 evalue_set_si(&EP->x.p->arr[0], 0, 1);
28 value_init(EP->x.p->arr[1].x.n);
29 if (den == NULL)
30 value_set_si(EP->x.p->arr[1].d, 1);
31 else
32 value_assign(EP->x.p->arr[1].d, *den);
33 zz2value(c, EP->x.p->arr[1].x.n);
34 return EP;
37 /* returns an evalue that corresponds to
39 * sum_i p[i] * x_i
41 evalue *multi_monom(vec_ZZ& p)
43 evalue *X = ALLOC(evalue);
44 value_init(X->d);
45 value_init(X->x.n);
46 unsigned nparam = p.length()-1;
47 zz2value(p[nparam], X->x.n);
48 value_set_si(X->d, 1);
49 for (int i = 0; i < nparam; ++i) {
50 if (p[i] == 0)
51 continue;
52 evalue *T = term(i, p[i]);
53 eadd(T, X);
54 free_evalue_refs(T);
55 delete T;
57 return X;
61 * Check whether mapping polyhedron P on the affine combination
62 * num yields a range that has a fixed quotient on integer
63 * division by d
64 * If zero is true, then we are only interested in the quotient
65 * for the cases where the remainder is zero.
66 * Returns NULL if false and a newly allocated value if true.
68 static Value *fixed_quotient(Polyhedron *P, vec_ZZ& num, Value d, bool zero)
70 Value* ret = NULL;
71 int len = num.length();
72 Matrix *T = Matrix_Alloc(2, len);
73 zz2values(num, T->p[0]);
74 value_set_si(T->p[1][len-1], 1);
75 Polyhedron *I = Polyhedron_Image(P, T, P->NbConstraints);
76 Matrix_Free(T);
78 int i;
79 for (i = 0; i < I->NbRays; ++i)
80 if (value_zero_p(I->Ray[i][2])) {
81 Polyhedron_Free(I);
82 return NULL;
85 Value min, max;
86 value_init(min);
87 value_init(max);
88 int bounded = line_minmax(I, &min, &max);
89 assert(bounded);
91 if (zero)
92 mpz_cdiv_q(min, min, d);
93 else
94 mpz_fdiv_q(min, min, d);
95 mpz_fdiv_q(max, max, d);
97 if (value_eq(min, max)) {
98 ret = ALLOC(Value);
99 value_init(*ret);
100 value_assign(*ret, min);
102 value_clear(min);
103 value_clear(max);
104 return ret;
108 * Normalize linear expression coef modulo m
109 * Removes common factor and reduces coefficients
110 * Returns index of first non-zero coefficient or len
112 int normal_mod(Value *coef, int len, Value *m)
114 Value gcd;
115 value_init(gcd);
117 Vector_Gcd(coef, len, &gcd);
118 value_gcd(gcd, gcd, *m);
119 Vector_AntiScale(coef, coef, gcd, len);
121 value_division(*m, *m, gcd);
122 value_clear(gcd);
124 if (value_one_p(*m))
125 return len;
127 int j;
128 for (j = 0; j < len; ++j)
129 mpz_fdiv_r(coef[j], coef[j], *m);
130 for (j = 0; j < len; ++j)
131 if (value_notzero_p(coef[j]))
132 break;
134 return j;
137 static bool mod_needed(Polyhedron *PD, vec_ZZ& num, Value d, evalue *E)
139 Value *q = fixed_quotient(PD, num, d, false);
141 if (!q)
142 return true;
144 value_oppose(*q, *q);
145 evalue EV;
146 value_init(EV.d);
147 value_set_si(EV.d, 1);
148 value_init(EV.x.n);
149 value_multiply(EV.x.n, *q, d);
150 eadd(&EV, E);
151 free_evalue_refs(&EV);
152 value_clear(*q);
153 free(q);
154 return false;
157 /* Computes the fractional part of the affine expression specified
158 * by coef (of length nvar+1) and the denominator denom.
159 * If PD is not NULL, then it specifies additional constraints
160 * on the variables that may be used to simplify the resulting
161 * fractional part expression.
163 * Modifies coef argument !
165 evalue *fractional_part(Value *coef, Value denom, int nvar, Polyhedron *PD)
167 Value m;
168 value_init(m);
169 evalue *EP = evalue_zero();
170 int sign = 1;
172 value_assign(m, denom);
173 int j = normal_mod(coef, nvar+1, &m);
175 if (j == nvar+1) {
176 value_clear(m);
177 return EP;
180 vec_ZZ num;
181 values2zz(coef, num, nvar+1);
183 ZZ g;
184 value2zz(m, g);
186 evalue tmp;
187 value_init(tmp.d);
188 evalue_set_si(&tmp, 0, 1);
190 int p = j;
191 if (g % 2 == 0)
192 while (j < nvar && (num[j] == g/2 || num[j] == 0))
193 ++j;
194 if ((j < nvar && num[j] > g/2) || (j == nvar && num[j] >= (g+1)/2)) {
195 for (int k = j; k < nvar; ++k)
196 if (num[k] != 0)
197 num[k] = g - num[k];
198 num[nvar] = g - 1 - num[nvar];
199 value_assign(tmp.d, m);
200 ZZ t = sign*(g-1);
201 zz2value(t, tmp.x.n);
202 eadd(&tmp, EP);
203 sign = -sign;
206 if (p >= nvar) {
207 ZZ t = num[nvar] * sign;
208 zz2value(t, tmp.x.n);
209 value_assign(tmp.d, m);
210 eadd(&tmp, EP);
211 } else {
212 evalue *E = multi_monom(num);
213 evalue EV;
214 value_init(EV.d);
216 if (PD && !mod_needed(PD, num, m, E)) {
217 value_init(EV.x.n);
218 value_set_si(EV.x.n, sign);
219 value_assign(EV.d, m);
220 emul(&EV, E);
221 eadd(E, EP);
222 } else {
223 value_init(EV.x.n);
224 value_set_si(EV.x.n, 1);
225 value_assign(EV.d, m);
226 emul(&EV, E);
227 value_clear(EV.x.n);
228 value_set_si(EV.d, 0);
229 EV.x.p = new_enode(fractional, 3, -1);
230 evalue_copy(&EV.x.p->arr[0], E);
231 evalue_set_si(&EV.x.p->arr[1], 0, 1);
232 value_init(EV.x.p->arr[2].x.n);
233 value_set_si(EV.x.p->arr[2].x.n, sign);
234 value_set_si(EV.x.p->arr[2].d, 1);
236 eadd(&EV, EP);
239 free_evalue_refs(&EV);
240 evalue_free(E);
243 free_evalue_refs(&tmp);
245 out:
246 value_clear(m);
248 return EP;
251 /* Computes the ceil of the affine expression specified
252 * by coef (of length nvar+1) and the denominator denom.
253 * If PD is not NULL, then it specifies additional constraints
254 * on the variables that may be used to simplify the resulting
255 * ceil expression.
257 * Modifies coef argument !
259 evalue *ceiling(Value *coef, Value denom, int nvar, Polyhedron *PD)
261 evalue *EP, *f;
262 EP = affine2evalue(coef, denom, nvar);
263 Vector_Oppose(coef, coef, nvar+1);
264 f = fractional_part(coef, denom, nvar, PD);
265 eadd(f, EP);
266 evalue_free(f);
267 return EP;
270 static evalue *ceil(Value *coef, int len, Value d,
271 barvinok_options *options)
273 evalue *c;
275 Vector_Oppose(coef, coef, len);
276 c = fractional_part(coef, d, len-1, NULL);
277 if (options->lookup_table)
278 evalue_mod2table(c, len-1);
279 return c;
282 void lattice_point_fixed(Value *vertex, Value *vertex_res,
283 Matrix *Rays, Matrix *Rays_res,
284 Value *point)
286 unsigned dim = Rays->NbRows;
287 if (value_one_p(vertex[dim]))
288 Vector_Copy(vertex_res, point, Rays_res->NbColumns);
289 else {
290 Matrix *R2 = Matrix_Copy(Rays);
291 Matrix *inv = Matrix_Alloc(Rays->NbRows, Rays->NbColumns);
292 int ok = Matrix_Inverse(R2, inv);
293 assert(ok);
294 Matrix_Free(R2);
295 Vector *lambda = Vector_Alloc(dim);
296 Vector_Matrix_Product(vertex, inv, lambda->p);
297 Matrix_Free(inv);
298 for (int j = 0; j < dim; ++j)
299 mpz_cdiv_q(lambda->p[j], lambda->p[j], vertex[dim]);
300 Vector_Matrix_Product(lambda->p, Rays_res, point);
301 Vector_Free(lambda);
305 static Matrix *Matrix_AddRowColumn(Matrix *M)
307 Matrix *M2 = Matrix_Alloc(M->NbRows+1, M->NbColumns+1);
308 for (int i = 0; i < M->NbRows; ++i)
309 Vector_Copy(M->p[i], M2->p[i], M->NbColumns);
310 value_set_si(M2->p[M->NbRows][M->NbColumns], 1);
311 return M2;
314 #define FORALL_COSETS(det,D,i,k) \
315 do { \
316 Vector *k = Vector_Alloc(D->NbRows+1); \
317 value_set_si(k->p[D->NbRows], 1); \
318 for (unsigned long i = 0; i < det; ++i) { \
319 if (i) \
320 for (int j = D->NbRows-1; j >= 0; --j) { \
321 value_increment(k->p[j], k->p[j]); \
322 if (value_eq(k->p[j], D->p[j][j])) \
323 value_set_si(k->p[j], 0); \
324 else \
325 break; \
328 #define END_FORALL_COSETS \
331 Vector_Free(k); \
332 } while(0);
334 /* Compute the lattice points in the vertex cone at "values" with rays "rays".
335 * The lattice points are returned in "vertex".
337 * Rays has the generators as rows and so does W.
338 * We first compute { m-v, u_i^* } with m = k W, where k runs through
339 * the cosets.
340 * We compute
341 * [k 1] [ d1*W 0 ] [ U' 0 ] = [k 1] T2
342 * [ -v d1 ] [ 0 d2 ]
343 * where d1 and d2 are the denominators of v and U^{-1}=U'/d2.
344 * Then lambda = { k } (componentwise)
345 * We compute x - floor(x) = {x} = { a/b } as fdiv_r(a,b)/b
346 * For open rays/facets, we need values in (0,1] rather than [0,1),
347 * so we compute {{x}} = x - ceil(x-1) = a/b - ceil((a-b)/b)
348 * = (a - b cdiv_q(a-b,b) - b + b)/b
349 * = (cdiv_r(a,b)+b)/b
350 * Finally, we compute v + lambda * U
351 * The denominator of lambda can be d1*d2, that of lambda2 = lambda*U
352 * can be at most d1, since it is integer if v = 0.
353 * The denominator of v + lambda2 is 1.
355 * The _res variants of the input variables may have been multiplied with
356 * a (list of) nonorthogonal vector(s) and may therefore have fewer columns
357 * than their original counterparts.
359 void lattice_points_fixed(Value *vertex, Value *vertex_res,
360 Matrix *Rays, Matrix *Rays_res, Matrix *points,
361 unsigned long det)
363 unsigned dim = Rays->NbRows;
364 if (det == 1) {
365 lattice_point_fixed(vertex, vertex_res, Rays, Rays_res,
366 points->p[0]);
367 return;
369 Matrix *U, *W, *D;
370 Smith(Rays, &U, &W, &D);
371 Matrix_Free(U);
373 /* Sanity check */
374 unsigned long det2 = 1;
375 for (int i = 0 ; i < D->NbRows; ++i)
376 det2 *= mpz_get_ui(D->p[i][i]);
377 assert(det == det2);
379 Matrix *T = Matrix_Alloc(W->NbRows+1, W->NbColumns+1);
380 for (int i = 0; i < W->NbRows; ++i)
381 Vector_Scale(W->p[i], T->p[i], vertex[dim], W->NbColumns);
382 Matrix_Free(W);
383 Vector_Oppose(vertex, T->p[dim], dim);
384 value_assign(T->p[dim][dim], vertex[dim]);
386 Matrix *R2 = Matrix_AddRowColumn(Rays);
387 Matrix *inv = Matrix_Alloc(R2->NbRows, R2->NbColumns);
388 int ok = Matrix_Inverse(R2, inv);
389 assert(ok);
390 Matrix_Free(R2);
392 Matrix *T2 = Matrix_Alloc(dim+1, dim+1);
393 Matrix_Product(T, inv, T2);
394 Matrix_Free(T);
396 Vector *lambda = Vector_Alloc(dim+1);
397 Vector *lambda2 = Vector_Alloc(Rays_res->NbColumns);
398 FORALL_COSETS(det, D, i, k)
399 Vector_Matrix_Product(k->p, T2, lambda->p);
400 for (int j = 0; j < dim; ++j)
401 mpz_fdiv_r(lambda->p[j], lambda->p[j], lambda->p[dim]);
402 Vector_Matrix_Product(lambda->p, Rays_res, lambda2->p);
403 for (int j = 0; j < lambda2->Size; ++j)
404 assert(mpz_divisible_p(lambda2->p[j], inv->p[dim][dim]));
405 Vector_AntiScale(lambda2->p, lambda2->p, inv->p[dim][dim], lambda2->Size);
406 Vector_Add(lambda2->p, vertex_res, lambda2->p, lambda2->Size);
407 for (int j = 0; j < lambda2->Size; ++j)
408 assert(mpz_divisible_p(lambda2->p[j], vertex[dim]));
409 Vector_AntiScale(lambda2->p, points->p[i], vertex[dim], lambda2->Size);
410 END_FORALL_COSETS
411 Vector_Free(lambda);
412 Vector_Free(lambda2);
413 Matrix_Free(D);
414 Matrix_Free(inv);
416 Matrix_Free(T2);
419 /* Returns the power of (t+1) in the term of a rational generating function,
420 * i.e., the scalar product of the actual lattice point and lambda.
421 * The lattice point is the unique lattice point in the fundamental parallelepiped
422 * of the unimodual cone i shifted to the parametric vertex W/lcm.
424 * The rows of W refer to the coordinates of the vertex
425 * The first nparam columns are the coefficients of the parameters
426 * and the final column is the constant term.
427 * lcm is the common denominator of all coefficients.
429 static evalue **lattice_point_fractional(const mat_ZZ& rays, vec_ZZ& lambda,
430 Matrix *V,
431 unsigned long det)
433 unsigned nparam = V->NbColumns-2;
434 evalue **E = new evalue *[det];
436 Matrix* Rays = zz2matrix(rays);
437 Matrix *T = Transpose(Rays);
438 Matrix *T2 = Matrix_Copy(T);
439 Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
440 int ok = Matrix_Inverse(T2, inv);
441 assert(ok);
442 Matrix_Free(T2);
443 mat_ZZ vertex;
444 matrix2zz(V, vertex, V->NbRows, V->NbColumns-1);
446 vec_ZZ num;
447 num = lambda * vertex;
449 evalue *EP = multi_monom(num);
451 evalue_div(EP, V->p[0][nparam+1]);
453 Matrix *L = Matrix_Alloc(inv->NbRows, V->NbColumns);
454 Matrix_Product(inv, V, L);
456 mat_ZZ RT;
457 matrix2zz(T, RT, T->NbRows, T->NbColumns);
458 Matrix_Free(T);
460 vec_ZZ p = lambda * RT;
462 Value tmp;
463 value_init(tmp);
465 if (det == 1) {
466 for (int i = 0; i < L->NbRows; ++i) {
467 evalue *f;
468 Vector_Oppose(L->p[i], L->p[i], nparam+1);
469 f = fractional_part(L->p[i], V->p[i][nparam+1], nparam, NULL);
470 zz2value(p[i], tmp);
471 evalue_mul(f, tmp);
472 eadd(f, EP);
473 evalue_free(f);
475 E[0] = EP;
476 } else {
477 for (int i = 0; i < L->NbRows; ++i)
478 value_assign(L->p[i][nparam+1], V->p[i][nparam+1]);
480 Value denom;
481 value_init(denom);
482 mpz_set_ui(denom, det);
483 value_multiply(denom, L->p[0][nparam+1], denom);
485 Matrix *U, *W, *D;
486 Smith(Rays, &U, &W, &D);
487 Matrix_Free(U);
489 /* Sanity check */
490 unsigned long det2 = 1;
491 for (int i = 0 ; i < D->NbRows; ++i)
492 det2 *= mpz_get_ui(D->p[i][i]);
493 assert(det == det2);
495 Matrix_Transposition(inv);
496 Matrix *T2 = Matrix_Alloc(W->NbRows, inv->NbColumns);
497 Matrix_Product(W, inv, T2);
498 Matrix_Free(W);
500 unsigned dim = D->NbRows;
501 Vector *lambda = Vector_Alloc(dim);
503 Vector *row = Vector_Alloc(nparam+1);
504 FORALL_COSETS(det, D, i, k)
505 Vector_Matrix_Product(k->p, T2, lambda->p);
506 E[i] = ALLOC(evalue);
507 value_init(E[i]->d);
508 evalue_copy(E[i], EP);
509 for (int j = 0; j < L->NbRows; ++j) {
510 evalue *f;
511 Vector_Oppose(L->p[j], row->p, nparam+1);
512 value_addmul(row->p[nparam], L->p[j][nparam+1], lambda->p[j]);
513 f = fractional_part(row->p, denom, nparam, NULL);
514 zz2value(p[j], tmp);
515 evalue_mul(f, tmp);
516 eadd(f, E[i]);
517 evalue_free(f);
519 END_FORALL_COSETS
520 Vector_Free(row);
522 Vector_Free(lambda);
523 Matrix_Free(T2);
524 Matrix_Free(D);
526 value_clear(denom);
527 evalue_free(EP);
529 value_clear(tmp);
531 Matrix_Free(Rays);
532 Matrix_Free(L);
533 Matrix_Free(inv);
535 return E;
538 static evalue **lattice_point(const mat_ZZ& rays, vec_ZZ& lambda,
539 Param_Vertices *V,
540 unsigned long det,
541 barvinok_options *options)
543 evalue **lp = lattice_point_fractional(rays, lambda, V->Vertex, det);
544 if (options->lookup_table) {
545 for (int i = 0; i < det; ++i)
546 evalue_mod2table(lp[i], V->Vertex->NbColumns-2);
548 return lp;
551 Matrix *relative_coordinates(Param_Vertices *V, Matrix *basis)
553 unsigned nparam = V->Vertex->NbColumns - 2;
554 Matrix *T = Matrix_Copy(basis);
555 Matrix *inv = Matrix_Alloc(T->NbRows, T->NbColumns);
556 int ok = Matrix_Inverse(T, inv);
557 assert(ok);
558 Matrix_Free(T);
560 Param_Vertex_Common_Denominator(V);
561 /* temporarily ignore (common) denominator */
562 V->Vertex->NbColumns--;
563 Matrix *L = Matrix_Alloc(inv->NbRows, V->Vertex->NbColumns);
564 Matrix_Product(inv, V->Vertex, L);
565 V->Vertex->NbColumns++;
566 Matrix_Free(inv);
568 return L;
571 /* returns the unique lattice point in the fundamental parallelepiped
572 * of the unimodual cone C shifted to the parametric vertex V.
574 * The return values num and E_vertex are such that
575 * coordinate i of this lattice point is equal to
577 * num[i] + E_vertex[i]
579 void lattice_point(Param_Vertices *V, const mat_ZZ& rays, vec_ZZ& num,
580 evalue **E_vertex, barvinok_options *options)
582 unsigned nparam = V->Vertex->NbColumns - 2;
583 unsigned dim = rays.NumCols();
585 /* It doesn't really make sense to call lattice_point when dim == 0,
586 * but apparently it happens from indicator_constructor in lexmin.
588 if (!dim)
589 return;
591 vec_ZZ vertex;
592 vertex.SetLength(nparam+1);
594 Value tmp;
595 value_init(tmp);
597 assert(V->Vertex->NbRows > 0);
598 Param_Vertex_Common_Denominator(V);
600 if (value_notone_p(V->Vertex->p[0][nparam+1])) {
601 Matrix* Rays = zz2matrix(rays);
602 Matrix *T = Transpose(Rays);
603 Matrix_Free(Rays);
604 Matrix *L = relative_coordinates(V, T);
606 evalue f;
607 value_init(f.d);
608 value_init(f.x.n);
610 evalue **remainders = new evalue *[dim];
611 for (int i = 0; i < dim; ++i)
612 remainders[i] = ceil(L->p[i], nparam+1, V->Vertex->p[0][nparam+1],
613 options);
614 Matrix_Free(L);
617 for (int i = 0; i < V->Vertex->NbRows; ++i) {
618 values2zz(V->Vertex->p[i], vertex, nparam+1);
619 E_vertex[i] = multi_monom(vertex);
620 num[i] = 0;
622 value_set_si(f.x.n, 1);
623 value_assign(f.d, V->Vertex->p[0][nparam+1]);
625 emul(&f, E_vertex[i]);
627 for (int j = 0; j < dim; ++j) {
628 if (value_zero_p(T->p[i][j]))
629 continue;
630 evalue cp;
631 value_init(cp.d);
632 evalue_copy(&cp, remainders[j]);
633 if (value_notone_p(T->p[i][j])) {
634 value_set_si(f.d, 1);
635 value_assign(f.x.n, T->p[i][j]);
636 emul(&f, &cp);
638 eadd(&cp, E_vertex[i]);
639 free_evalue_refs(&cp);
642 for (int i = 0; i < dim; ++i)
643 evalue_free(remainders[i]);
644 delete [] remainders;
646 free_evalue_refs(&f);
648 Matrix_Free(T);
649 value_clear(tmp);
650 return;
652 value_clear(tmp);
654 for (int i = 0; i < V->Vertex->NbRows; ++i) {
655 /* fixed value */
656 if (First_Non_Zero(V->Vertex->p[i], nparam) == -1) {
657 E_vertex[i] = 0;
658 value2zz(V->Vertex->p[i][nparam], num[i]);
659 } else {
660 values2zz(V->Vertex->p[i], vertex, nparam+1);
661 E_vertex[i] = multi_monom(vertex);
662 num[i] = 0;
667 static int lattice_point_fixed(Param_Vertices* V, const mat_ZZ& rays,
668 vec_ZZ& lambda, term_info* term, unsigned long det)
670 unsigned nparam = V->Vertex->NbColumns - 2;
671 unsigned dim = rays.NumCols();
673 for (int i = 0; i < dim; ++i)
674 if (First_Non_Zero(V->Vertex->p[i], nparam) != -1)
675 return 0;
677 Vector *fixed = Vector_Alloc(dim+1);
678 for (int i = 0; i < dim; ++i)
679 value_assign(fixed->p[i], V->Vertex->p[i][nparam]);
680 value_assign(fixed->p[dim], V->Vertex->p[0][nparam+1]);
682 mat_ZZ vertex;
683 Matrix *points = Matrix_Alloc(det, dim);
684 Matrix* Rays = zz2matrix(rays);
685 lattice_points_fixed(fixed->p, fixed->p, Rays, Rays, points, det);
686 Matrix_Free(Rays);
687 matrix2zz(points, vertex, points->NbRows, points->NbColumns);
688 Matrix_Free(points);
689 term->E = NULL;
690 term->constant = vertex * lambda;
691 Vector_Free(fixed);
693 return 1;
696 /* Returns the power of (t+1) in the term of a rational generating function,
697 * i.e., the scalar product of the actual lattice point and lambda.
698 * The lattice point is the unique lattice point in the fundamental parallelepiped
699 * of the unimodual cone i shifted to the parametric vertex V.
701 * The result is returned in term.
703 void lattice_point(Param_Vertices* V, const mat_ZZ& rays, vec_ZZ& lambda,
704 term_info* term, unsigned long det,
705 barvinok_options *options)
707 unsigned nparam = V->Vertex->NbColumns - 2;
708 unsigned dim = rays.NumCols();
709 mat_ZZ vertex;
710 vertex.SetDims(V->Vertex->NbRows, nparam+1);
712 Param_Vertex_Common_Denominator(V);
714 if (lattice_point_fixed(V, rays, lambda, term, det))
715 return;
717 if (det != 1 || value_notone_p(V->Vertex->p[0][nparam+1])) {
718 term->E = lattice_point(rays, lambda, V, det, options);
719 return;
721 for (int i = 0; i < V->Vertex->NbRows; ++i) {
722 assert(value_one_p(V->Vertex->p[i][nparam+1])); // for now
723 values2zz(V->Vertex->p[i], vertex[i], nparam+1);
726 vec_ZZ num;
727 num = lambda * vertex;
729 int nn = 0;
730 for (int j = 0; j < nparam; ++j)
731 if (num[j] != 0)
732 ++nn;
733 if (nn >= 1) {
734 term->E = new evalue *[1];
735 term->E[0] = multi_monom(num);
736 } else {
737 term->E = NULL;
738 term->constant.SetLength(1);
739 term->constant[0] = num[nparam];