?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
evalue_read.c: add evalue_read_from_str for reading from a string
[barvinok.git] / lattice_point.cc
blob1bab51a85e6b708dd7a03a8f28ed463f3fe46b6a
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"
11
12 using std::cerr;
13 using std::endl;
14
15 #define ALLOC(type) (type*)malloc(sizeof(type))
16
17 /* returns an evalue that corresponds to
18 *
19 * c/(*den) x_param
20 */
21 static evalue *term(int param, ZZ& c, Value *den = NULL)
22 {
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;
35 }
36
37 /* returns an evalue that corresponds to
38 *
39 * sum_i p[i] * x_i
40 */
41 evalue *multi_monom(vec_ZZ& p)
42 {
43 evalue *X = new 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;
56 }
57 return X;
58 }
59
60 /*
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.
67 */
68 static Value *fixed_quotient(Polyhedron *P, vec_ZZ& num, Value d, bool zero)
69 {
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);
77
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;
83 }
84
85 Value min, max;
86 value_init(min);
87 value_init(max);
88 int bounded = line_minmax(I, &min, &max);
89 assert(bounded);
90
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);
96
97 if (value_eq(min, max)) {
98 ret = ALLOC(Value);
99 value_init(*ret);
100 value_assign(*ret, min);
101 }
102 value_clear(min);
103 value_clear(max);
104 return ret;
105 }
106
107 /*
108 * Normalize linear expression coef modulo m
109 * Removes common factor and reduces coefficients
110 * Returns index of first non-zero coefficient or len
111 */
112 int normal_mod(Value *coef, int len, Value *m)
113 {
114 Value gcd;
115 value_init(gcd);
116
117 Vector_Gcd(coef, len, &gcd);
118 Gcd(gcd, *m, &gcd);
119 Vector_AntiScale(coef, coef, gcd, len);
120
121 value_division(*m, *m, gcd);
122 value_clear(gcd);
123
124 if (value_one_p(*m))
125 return len;
126
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;
133
134 return j;
135 }
136
137 static bool mod_needed(Polyhedron *PD, vec_ZZ& num, Value d, evalue *E)
138 {
139 Value *q = fixed_quotient(PD, num, d, false);
140
141 if (!q)
142 return true;
143
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;
155 }
156
157 /* modifies coef argument ! */
158 static void fractional_part(Value *coef, int len, Value d, ZZ f, evalue *EP,
159 Polyhedron *PD, bool up)
160 {
161 Value m;
162 value_init(m);
163
164 if (up) {
165 /* f {{ x }} = f - f { -x } */
166 zz2value(f, m);
167 evalue_add_constant(EP, m);
168 Vector_Oppose(coef, coef, len);
169 f = -f;
170 }
171
172 value_assign(m, d);
173 int j = normal_mod(coef, len, &m);
174
175 if (j == len) {
176 value_clear(m);
177 return;
178 }
179
180 vec_ZZ num;
181 values2zz(coef, num, len);
182
183 ZZ g;
184 value2zz(m, g);
185
186 evalue tmp;
187 value_init(tmp.d);
188 evalue_set_si(&tmp, 0, 1);
189
190 int p = j;
191 if (g % 2 == 0)
192 while (j < len-1 && (num[j] == g/2 || num[j] == 0))
193 ++j;
194 if ((j < len-1 && num[j] > g/2) || (j == len-1 && num[j] >= (g+1)/2)) {
195 for (int k = j; k < len-1; ++k)
196 if (num[k] != 0)
197 num[k] = g - num[k];
198 num[len-1] = g - 1 - num[len-1];
199 value_assign(tmp.d, m);
200 ZZ t = f*(g-1);
201 zz2value(t, tmp.x.n);
202 eadd(&tmp, EP);
203 f = -f;
204 }
205
206 if (p >= len-1) {
207 ZZ t = num[len-1] * f;
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);
215
216 if (PD && !mod_needed(PD, num, m, E)) {
217 value_init(EV.x.n);
218 zz2value(f, EV.x.n);
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 zz2value(f, EV.x.p->arr[2].x.n);
234 value_set_si(EV.x.p->arr[2].d, 1);
235
236 eadd(&EV, EP);
237 }
238
239 free_evalue_refs(&EV);
240 free_evalue_refs(E);
241 delete E;
242 }
243
244 free_evalue_refs(&tmp);
245
246 out:
247 value_clear(m);
248 }
249
250 static void ceil(Value *coef, int len, Value d, ZZ& f,
251 evalue *EP, barvinok_options *options)
252 {
253 Vector_Oppose(coef, coef, len);
254 fractional_part(coef, len, d, f, EP, NULL, false);
255 if (options->lookup_table)
256 evalue_mod2table(EP, len-1);
257 }
258
259 evalue* bv_ceil3(Value *coef, int len, Value d, Polyhedron *P)
260 {
261 Vector *val = Vector_Alloc(len);
262
263 Value t;
264 value_init(t);
265 value_set_si(t, -1);
266 Vector_Scale(coef, val->p, t, len);
267 value_absolute(t, d);
268
269 vec_ZZ num;
270 values2zz(val->p, num, len);
271 evalue *EP = multi_monom(num);
272
273 evalue tmp;
274 value_init(tmp.d);
275 value_init(tmp.x.n);
276 value_set_si(tmp.x.n, 1);
277 value_assign(tmp.d, t);
278
279 emul(&tmp, EP);
280
281 ZZ one;
282 one = 1;
283 Vector_Oppose(val->p, val->p, len);
284 fractional_part(val->p, len, t, one, EP, P, false);
285 value_clear(t);
286
287 /* copy EP to malloc'ed evalue */
288 evalue *E = ALLOC(evalue);
289 *E = *EP;
290 delete EP;
291
292 free_evalue_refs(&tmp);
293 Vector_Free(val);
294
295 return E;
296 }
297
298 void lattice_point(Value* values, const mat_ZZ& rays, vec_ZZ& vertex, int *closed)
299 {
300 unsigned dim = rays.NumRows();
301 if (value_one_p(values[dim]) && !closed)
302 values2zz(values, vertex, dim);
303 else {
304 Matrix* Rays = rays2matrix(rays);
305 Matrix *inv = Matrix_Alloc(Rays->NbRows, Rays->NbColumns);
306 int ok = Matrix_Inverse(Rays, inv);
307 assert(ok);
308 Matrix_Free(Rays);
309 Rays = rays2matrix(rays);
310 Vector *lambda = Vector_Alloc(dim+1);
311 Vector_Matrix_Product(values, inv, lambda->p);
312 Matrix_Free(inv);
313 for (int j = 0; j < dim; ++j)
314 if (!closed || closed[j])
315 mpz_cdiv_q(lambda->p[j], lambda->p[j], lambda->p[dim]);
316 else {
317 value_addto(lambda->p[j], lambda->p[j], lambda->p[dim]);
318 mpz_fdiv_q(lambda->p[j], lambda->p[j], lambda->p[dim]);
319 }
320 value_set_si(lambda->p[dim], 1);
321 Vector *A = Vector_Alloc(dim+1);
322 Vector_Matrix_Product(lambda->p, Rays, A->p);
323 Vector_Free(lambda);
324 Matrix_Free(Rays);
325 values2zz(A->p, vertex, dim);
326 Vector_Free(A);
327 }
328 }
329
330 #define FORALL_COSETS(det,D,i,k) \
331 do { \
332 Vector *k = Vector_Alloc(D->NbRows+1); \
333 value_set_si(k->p[D->NbRows], 1); \
334 for (unsigned long i = 0; i < det; ++i) { \
335 unsigned long _fc_val = i; \
336 for (int j = 0; j < D->NbRows; ++j) { \
337 value_set_si(k->p[j], _fc_val % mpz_get_ui(D->p[j][j]));\
338 _fc_val /= mpz_get_ui(D->p[j][j]); \
339 }
340 #define END_FORALL_COSETS \
341 } \
342 Vector_Free(k); \
343 } while(0);
344
345 /* Compute the lattice points in the vertex cone at "values" with rays "rays".
346 * The lattice points are returned in "vertex".
347 *
348 * Rays has the generators as rows and so does W.
349 * We first compute { m-v, u_i^* } with m = k W, where k runs through
350 * the cosets.
351 * We compute
352 * [k 1] [ d1*W 0 ] [ U' 0 ] = [k 1] T2
353 * [ -v d1 ] [ 0 d2 ]
354 * where d1 and d2 are the denominators of v and U^{-1}=U'/d2.
355 * Then lambda = { k } (componentwise)
356 * We compute x - floor(x) = {x} = { a/b } as fdiv_r(a,b)/b
357 * For open rays/facets, we need values in (0,1] rather than [0,1),
358 * so we compute {{x}} = x - ceil(x-1) = a/b - ceil((a-b)/b)
359 * = (a - b cdiv_q(a-b,b) - b + b)/b
360 * = (cdiv_r(a,b)+b)/b
361 * Finally, we compute v + lambda * U
362 * The denominator of lambda can be d1*d2, that of lambda2 = lambda*U
363 * can be at most d1, since it is integer if v = 0.
364 * The denominator of v + lambda2 is 1.
365 */
366 void lattice_point(Value* values, const mat_ZZ& rays, mat_ZZ& vertex,
367 unsigned long det, int *closed)
368 {
369 unsigned dim = rays.NumRows();
370 vertex.SetDims(det, dim);
371 if (det == 1) {
372 lattice_point(values, rays, vertex[0], closed);
373 return;
374 }
375 Matrix* Rays = zz2matrix(rays);
376 Matrix *U, *W, *D;
377 Smith(Rays, &U, &W, &D);
378 Matrix_Free(Rays);
379 Matrix_Free(U);
380
381 /* Sanity check */
382 unsigned long det2 = 1;
383 for (int i = 0 ; i < D->NbRows; ++i)
384 det2 *= mpz_get_ui(D->p[i][i]);
385 assert(det == det2);
386
387 Matrix *T = Matrix_Alloc(W->NbRows+1, W->NbColumns+1);
388 for (int i = 0; i < W->NbRows; ++i)
389 Vector_Scale(W->p[i], T->p[i], values[dim], W->NbColumns);
390 Matrix_Free(W);
391 Value tmp;
392 value_init(tmp);
393 value_set_si(tmp, -1);
394 Vector_Scale(values, T->p[dim], tmp, dim);
395 value_clear(tmp);
396 value_assign(T->p[dim][dim], values[dim]);
397
398 Rays = rays2matrix(rays);
399 Matrix *inv = Matrix_Alloc(Rays->NbRows, Rays->NbColumns);
400 int ok = Matrix_Inverse(Rays, inv);
401 assert(ok);
402 Matrix_Free(Rays);
403
404 Matrix *T2 = Matrix_Alloc(dim+1, dim+1);
405 Matrix_Product(T, inv, T2);
406 Matrix_Free(T);
407
408 Rays = rays2matrix(rays);
409
410 Vector *lambda = Vector_Alloc(dim+1);
411 Vector *lambda2 = Vector_Alloc(dim+1);
412 FORALL_COSETS(det, D, i, k)
413 Vector_Matrix_Product(k->p, T2, lambda->p);
414 for (int j = 0; j < dim; ++j)
415 if (!closed || closed[j])
416 mpz_fdiv_r(lambda->p[j], lambda->p[j], lambda->p[dim]);
417 else {
418 mpz_cdiv_r(lambda->p[j], lambda->p[j], lambda->p[dim]);
419 value_addto(lambda->p[j], lambda->p[j], lambda->p[dim]);
420 }
421 Vector_Matrix_Product(lambda->p, Rays, lambda2->p);
422 for (int j = 0; j < dim; ++j)
423 assert(mpz_divisible_p(lambda2->p[j], inv->p[dim][dim]));
424 Vector_AntiScale(lambda2->p, lambda2->p, inv->p[dim][dim], dim+1);
425 Vector_Add(lambda2->p, values, lambda2->p, dim);
426 for (int j = 0; j < dim; ++j)
427 assert(mpz_divisible_p(lambda2->p[j], values[dim]));
428 Vector_AntiScale(lambda2->p, lambda2->p, values[dim], dim+1);
429 values2zz(lambda2->p, vertex[i], dim);
430 END_FORALL_COSETS
431 Vector_Free(lambda);
432 Vector_Free(lambda2);
433 Matrix_Free(D);
434 Matrix_Free(Rays);
435 Matrix_Free(inv);
436
437 Matrix_Free(T2);
438 }
439
440 /* Returns the power of (t+1) in the term of a rational generating function,
441 * i.e., the scalar product of the actual lattice point and lambda.
442 * The lattice point is the unique lattice point in the fundamental parallelepiped
443 * of the unimodual cone i shifted to the parametric vertex W/lcm.
444 *
445 * The rows of W refer to the coordinates of the vertex
446 * The first nparam columns are the coefficients of the parameters
447 * and the final column is the constant term.
448 * lcm is the common denominator of all coefficients.
449 */
450 static evalue **lattice_point_fractional(const mat_ZZ& rays, vec_ZZ& lambda,
451 Matrix *V,
452 unsigned long det, int *closed)
453 {
454 unsigned nparam = V->NbColumns-2;
455 evalue **E = new evalue *[det];
456
457 Matrix* Rays = zz2matrix(rays);
458 Matrix *T = Transpose(Rays);
459 Matrix *T2 = Matrix_Copy(T);
460 Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
461 int ok = Matrix_Inverse(T2, inv);
462 assert(ok);
463 Matrix_Free(T2);
464 mat_ZZ vertex;
465 matrix2zz(V, vertex, V->NbRows, V->NbColumns-1);
466
467 vec_ZZ num;
468 num = lambda * vertex;
469
470 evalue *EP = multi_monom(num);
471
472 evalue_div(EP, V->p[0][nparam+1]);
473
474 Matrix *L = Matrix_Alloc(inv->NbRows, V->NbColumns);
475 Matrix_Product(inv, V, L);
476
477 mat_ZZ RT;
478 matrix2zz(T, RT, T->NbRows, T->NbColumns);
479 Matrix_Free(T);
480
481 vec_ZZ p = lambda * RT;
482
483 if (det == 1) {
484 for (int i = 0; i < L->NbRows; ++i) {
485 Vector_Oppose(L->p[i], L->p[i], nparam+1);
486 fractional_part(L->p[i], nparam+1, V->p[i][nparam+1], p[i],
487 EP, NULL, closed && !closed[i]);
488 }
489 E[0] = EP;
490 } else {
491 for (int i = 0; i < L->NbRows; ++i)
492 value_assign(L->p[i][nparam+1], V->p[i][nparam+1]);
493
494 Value denom;
495 value_init(denom);
496 mpz_set_ui(denom, det);
497 value_multiply(denom, L->p[0][nparam+1], denom);
498
499 Matrix *U, *W, *D;
500 Smith(Rays, &U, &W, &D);
501 Matrix_Free(U);
502
503 /* Sanity check */
504 unsigned long det2 = 1;
505 for (int i = 0 ; i < D->NbRows; ++i)
506 det2 *= mpz_get_ui(D->p[i][i]);
507 assert(det == det2);
508
509 Matrix_Transposition(inv);
510 Matrix *T2 = Matrix_Alloc(W->NbRows, inv->NbColumns);
511 Matrix_Product(W, inv, T2);
512 Matrix_Free(W);
513
514 unsigned dim = D->NbRows;
515 Vector *lambda = Vector_Alloc(dim);
516
517 Vector *row = Vector_Alloc(nparam+1);
518 FORALL_COSETS(det, D, i, k)
519 Vector_Matrix_Product(k->p, T2, lambda->p);
520 E[i] = new evalue();
521 value_init(E[i]->d);
522 evalue_copy(E[i], EP);
523 for (int j = 0; j < L->NbRows; ++j) {
524 Vector_Oppose(L->p[j], row->p, nparam+1);
525 value_addmul(row->p[nparam], L->p[j][nparam+1], lambda->p[j]);
526 fractional_part(row->p, nparam+1, denom, p[j],
527 E[i], NULL, closed && !closed[j]);
528 }
529 END_FORALL_COSETS
530 Vector_Free(row);
531
532 Vector_Free(lambda);
533 Matrix_Free(T2);
534 Matrix_Free(D);
535
536 value_clear(denom);
537 free_evalue_refs(EP);
538 delete EP;
539 }
540
541 Matrix_Free(Rays);
542 Matrix_Free(L);
543 Matrix_Free(inv);
544
545 return E;
546 }
547
548 static evalue **lattice_point(const mat_ZZ& rays, vec_ZZ& lambda,
549 Param_Vertices *V,
550 unsigned long det, int *closed,
551 barvinok_options *options)
552 {
553 evalue **lp = lattice_point_fractional(rays, lambda, V->Vertex, det, closed);
554 if (options->lookup_table) {
555 for (int i = 0; i < det; ++i)
556 evalue_mod2table(lp[i], V->Vertex->NbColumns-2);
557 }
558 return lp;
559 }
560
561 /* returns the unique lattice point in the fundamental parallelepiped
562 * of the unimodual cone C shifted to the parametric vertex V.
563 *
564 * The return values num and E_vertex are such that
565 * coordinate i of this lattice point is equal to
566 *
567 * num[i] + E_vertex[i]
568 */
569 void lattice_point(Param_Vertices *V, const mat_ZZ& rays, vec_ZZ& num,
570 evalue **E_vertex, barvinok_options *options)
571 {
572 unsigned nparam = V->Vertex->NbColumns - 2;
573 unsigned dim = rays.NumCols();
574 vec_ZZ vertex;
575 vertex.SetLength(nparam+1);
576
577 Value lcm, tmp;
578 value_init(lcm);
579 value_init(tmp);
580 value_set_si(lcm, 1);
581
582 for (int j = 0; j < V->Vertex->NbRows; ++j) {
583 value_lcm(lcm, V->Vertex->p[j][nparam+1], &lcm);
584 }
585
586 if (value_notone_p(lcm)) {
587 Matrix * mv = Matrix_Alloc(dim, nparam+1);
588 for (int j = 0 ; j < dim; ++j) {
589 value_division(tmp, lcm, V->Vertex->p[j][nparam+1]);
590 Vector_Scale(V->Vertex->p[j], mv->p[j], tmp, nparam+1);
591 }
592
593 Matrix* Rays = zz2matrix(rays);
594 Matrix *T = Transpose(Rays);
595 Matrix *T2 = Matrix_Copy(T);
596 Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
597 int ok = Matrix_Inverse(T2, inv);
598 assert(ok);
599 Matrix_Free(Rays);
600 Matrix_Free(T2);
601 Matrix *L = Matrix_Alloc(inv->NbRows, mv->NbColumns);
602 Matrix_Product(inv, mv, L);
603 Matrix_Free(inv);
604
605 evalue f;
606 value_init(f.d);
607 value_init(f.x.n);
608
609 ZZ one;
610
611 evalue *remainders[dim];
612 for (int i = 0; i < dim; ++i) {
613 remainders[i] = evalue_zero();
614 one = 1;
615 ceil(L->p[i], nparam+1, lcm, one, remainders[i], options);
616 }
617 Matrix_Free(L);
618
619
620 for (int i = 0; i < V->Vertex->NbRows; ++i) {
621 values2zz(mv->p[i], vertex, nparam+1);
622 E_vertex[i] = multi_monom(vertex);
623 num[i] = 0;
624
625 value_set_si(f.x.n, 1);
626 value_assign(f.d, lcm);
627
628 emul(&f, E_vertex[i]);
629
630 for (int j = 0; j < dim; ++j) {
631 if (value_zero_p(T->p[i][j]))
632 continue;
633 evalue cp;
634 value_init(cp.d);
635 evalue_copy(&cp, remainders[j]);
636 if (value_notone_p(T->p[i][j])) {
637 value_set_si(f.d, 1);
638 value_assign(f.x.n, T->p[i][j]);
639 emul(&f, &cp);
640 }
641 eadd(&cp, E_vertex[i]);
642 free_evalue_refs(&cp);
643 }
644 }
645 for (int i = 0; i < dim; ++i) {
646 free_evalue_refs(remainders[i]);
647 free(remainders[i]);
648 }
649
650 free_evalue_refs(&f);
651
652 Matrix_Free(T);
653 Matrix_Free(mv);
654 value_clear(lcm);
655 value_clear(tmp);
656 return;
657 }
658 value_clear(lcm);
659 value_clear(tmp);
660
661 for (int i = 0; i < V->Vertex->NbRows; ++i) {
662 /* fixed value */
663 if (First_Non_Zero(V->Vertex->p[i], nparam) == -1) {
664 E_vertex[i] = 0;
665 value2zz(V->Vertex->p[i][nparam], num[i]);
666 } else {
667 values2zz(V->Vertex->p[i], vertex, nparam+1);
668 E_vertex[i] = multi_monom(vertex);
669 num[i] = 0;
670 }
671 }
672 }
673
674 static int lattice_point_fixed(Param_Vertices* V, const mat_ZZ& rays,
675 vec_ZZ& lambda, term_info* term, unsigned long det, int *closed)
676 {
677 unsigned nparam = V->Vertex->NbColumns - 2;
678 unsigned dim = rays.NumCols();
679
680 for (int i = 0; i < dim; ++i)
681 if (First_Non_Zero(V->Vertex->p[i], nparam) != -1)
682 return 0;
683
684 Vector *fixed = Vector_Alloc(dim+1);
685 for (int i = 0; i < dim; ++i)
686 value_assign(fixed->p[i], V->Vertex->p[i][nparam]);
687 value_assign(fixed->p[dim], V->Vertex->p[0][nparam+1]);
688
689 mat_ZZ vertex;
690 lattice_point(fixed->p, rays, vertex, det, closed);
691 term->E = NULL;
692 term->constant = vertex * lambda;
693 Vector_Free(fixed);
694
695 return 1;
696 }
697
698 /* Returns the power of (t+1) in the term of a rational generating function,
699 * i.e., the scalar product of the actual lattice point and lambda.
700 * The lattice point is the unique lattice point in the fundamental parallelepiped
701 * of the unimodual cone i shifted to the parametric vertex V.
702 *
703 * The result is returned in term.
704 */
705 void lattice_point(Param_Vertices* V, const mat_ZZ& rays, vec_ZZ& lambda,
706 term_info* term, unsigned long det, int *closed,
707 barvinok_options *options)
708 {
709 unsigned nparam = V->Vertex->NbColumns - 2;
710 unsigned dim = rays.NumCols();
711 mat_ZZ vertex;
712 vertex.SetDims(V->Vertex->NbRows, nparam+1);
713
714 Param_Vertex_Common_Denominator(V);
715
716 if (lattice_point_fixed(V, rays, lambda, term, det, closed))
717 return;
718
719 if (det != 1 || closed || value_notone_p(V->Vertex->p[0][nparam+1])) {
720 term->E = lattice_point(rays, lambda, V, det, closed, options);
721 return;
722 }
723 for (int i = 0; i < V->Vertex->NbRows; ++i) {
724 assert(value_one_p(V->Vertex->p[i][nparam+1])); // for now
725 values2zz(V->Vertex->p[i], vertex[i], nparam+1);
726 }
727
728 vec_ZZ num;
729 num = lambda * vertex;
730
731 int nn = 0;
732 for (int j = 0; j < nparam; ++j)
733 if (num[j] != 0)
734 ++nn;
735 if (nn >= 1) {
736 term->E = new evalue *[1];
737 term->E[0] = multi_monom(num);
738 } else {
739 term->E = NULL;
740 term->constant.SetLength(1);
741 term->constant[0] = num[nparam];
742 }
743 }
lattice point enumeration library