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>
8 #include "conversion.h"
9 #include "lattice_point.h"
14 #define ALLOC(type) (type*)malloc(sizeof(type))
16 /* returns an evalue that corresponds to
20 static evalue
*term(int param
, ZZ
& c
, Value
*den
= NULL
)
22 evalue
*EP
= new evalue();
24 value_set_si(EP
->d
,0);
25 EP
->x
.p
= new_enode(polynomial
, 2, param
+ 1);
26 evalue_set_si(&EP
->x
.p
->arr
[0], 0, 1);
27 value_init(EP
->x
.p
->arr
[1].x
.n
);
29 value_set_si(EP
->x
.p
->arr
[1].d
, 1);
31 value_assign(EP
->x
.p
->arr
[1].d
, *den
);
32 zz2value(c
, EP
->x
.p
->arr
[1].x
.n
);
36 /* returns an evalue that corresponds to
40 evalue
*multi_monom(vec_ZZ
& p
)
42 evalue
*X
= new evalue();
45 unsigned nparam
= p
.length()-1;
46 zz2value(p
[nparam
], X
->x
.n
);
47 value_set_si(X
->d
, 1);
48 for (int i
= 0; i
< nparam
; ++i
) {
51 evalue
*T
= term(i
, p
[i
]);
60 * Check whether mapping polyhedron P on the affine combination
61 * num yields a range that has a fixed quotient on integer
63 * If zero is true, then we are only interested in the quotient
64 * for the cases where the remainder is zero.
65 * Returns NULL if false and a newly allocated value if true.
67 static Value
*fixed_quotient(Polyhedron
*P
, vec_ZZ
& num
, Value d
, bool zero
)
70 int len
= num
.length();
71 Matrix
*T
= Matrix_Alloc(2, len
);
72 zz2values(num
, T
->p
[0]);
73 value_set_si(T
->p
[1][len
-1], 1);
74 Polyhedron
*I
= Polyhedron_Image(P
, T
, P
->NbConstraints
);
78 for (i
= 0; i
< I
->NbRays
; ++i
)
79 if (value_zero_p(I
->Ray
[i
][2])) {
87 int bounded
= line_minmax(I
, &min
, &max
);
91 mpz_cdiv_q(min
, min
, d
);
93 mpz_fdiv_q(min
, min
, d
);
94 mpz_fdiv_q(max
, max
, d
);
96 if (value_eq(min
, max
)) {
99 value_assign(*ret
, min
);
107 * Normalize linear expression coef modulo m
108 * Removes common factor and reduces coefficients
109 * Returns index of first non-zero coefficient or len
111 int normal_mod(Value
*coef
, int len
, Value
*m
)
116 Vector_Gcd(coef
, len
, &gcd
);
118 Vector_AntiScale(coef
, coef
, gcd
, len
);
120 value_division(*m
, *m
, gcd
);
127 for (j
= 0; j
< len
; ++j
)
128 mpz_fdiv_r(coef
[j
], coef
[j
], *m
);
129 for (j
= 0; j
< len
; ++j
)
130 if (value_notzero_p(coef
[j
]))
136 static bool mod_needed(Polyhedron
*PD
, vec_ZZ
& num
, Value d
, evalue
*E
)
138 Value
*q
= fixed_quotient(PD
, num
, d
, false);
143 value_oppose(*q
, *q
);
146 value_set_si(EV
.d
, 1);
148 value_multiply(EV
.x
.n
, *q
, d
);
150 free_evalue_refs(&EV
);
156 /* modifies f argument ! */
157 static void ceil_mod(Value
*coef
, int len
, Value d
, ZZ
& f
, evalue
*EP
, Polyhedron
*PD
)
163 Vector_Scale(coef
, coef
, m
, len
);
166 int j
= normal_mod(coef
, len
, &m
);
174 values2zz(coef
, num
, len
);
181 evalue_set_si(&tmp
, 0, 1);
185 while (j
< len
-1 && (num
[j
] == g
/2 || num
[j
] == 0))
187 if ((j
< len
-1 && num
[j
] > g
/2) || (j
== len
-1 && num
[j
] >= (g
+1)/2)) {
188 for (int k
= j
; k
< len
-1; ++k
)
191 num
[len
-1] = g
- 1 - num
[len
-1];
192 value_assign(tmp
.d
, m
);
194 zz2value(t
, tmp
.x
.n
);
200 ZZ t
= num
[len
-1] * f
;
201 zz2value(t
, tmp
.x
.n
);
202 value_assign(tmp
.d
, m
);
205 evalue
*E
= multi_monom(num
);
209 if (PD
&& !mod_needed(PD
, num
, m
, E
)) {
212 value_assign(EV
.d
, m
);
217 value_set_si(EV
.x
.n
, 1);
218 value_assign(EV
.d
, m
);
221 value_set_si(EV
.d
, 0);
222 EV
.x
.p
= new_enode(fractional
, 3, -1);
223 evalue_copy(&EV
.x
.p
->arr
[0], E
);
224 evalue_set_si(&EV
.x
.p
->arr
[1], 0, 1);
225 value_init(EV
.x
.p
->arr
[2].x
.n
);
226 zz2value(f
, EV
.x
.p
->arr
[2].x
.n
);
227 value_set_si(EV
.x
.p
->arr
[2].d
, 1);
232 free_evalue_refs(&EV
);
237 free_evalue_refs(&tmp
);
243 static void ceil(Value
*coef
, int len
, Value d
, ZZ
& f
,
244 evalue
*EP
, Polyhedron
*PD
, barvinok_options
*options
)
246 ceil_mod(coef
, len
, d
, f
, EP
, PD
);
247 if (options
->lookup_table
)
248 evalue_mod2table(EP
, len
-1);
251 evalue
* bv_ceil3(Value
*coef
, int len
, Value d
, Polyhedron
*P
)
253 Vector
*val
= Vector_Alloc(len
);
258 Vector_Scale(coef
, val
->p
, t
, len
);
259 value_absolute(t
, d
);
262 values2zz(val
->p
, num
, len
);
263 evalue
*EP
= multi_monom(num
);
268 value_set_si(tmp
.x
.n
, 1);
269 value_assign(tmp
.d
, t
);
275 ceil_mod(val
->p
, len
, t
, one
, EP
, P
);
278 /* copy EP to malloc'ed evalue */
279 evalue
*E
= ALLOC(evalue
);
283 free_evalue_refs(&tmp
);
289 void lattice_point(Value
* values
, const mat_ZZ
& rays
, vec_ZZ
& vertex
, int *closed
)
291 unsigned dim
= rays
.NumRows();
292 if (value_one_p(values
[dim
]) && !closed
)
293 values2zz(values
, vertex
, dim
);
295 Matrix
* Rays
= rays2matrix(rays
);
296 Matrix
*inv
= Matrix_Alloc(Rays
->NbRows
, Rays
->NbColumns
);
297 int ok
= Matrix_Inverse(Rays
, inv
);
300 Rays
= rays2matrix(rays
);
301 Vector
*lambda
= Vector_Alloc(dim
+1);
302 Vector_Matrix_Product(values
, inv
, lambda
->p
);
304 for (int j
= 0; j
< dim
; ++j
)
305 if (!closed
|| closed
[j
])
306 mpz_cdiv_q(lambda
->p
[j
], lambda
->p
[j
], lambda
->p
[dim
]);
308 value_addto(lambda
->p
[j
], lambda
->p
[j
], lambda
->p
[dim
]);
309 mpz_fdiv_q(lambda
->p
[j
], lambda
->p
[j
], lambda
->p
[dim
]);
311 value_set_si(lambda
->p
[dim
], 1);
312 Vector
*A
= Vector_Alloc(dim
+1);
313 Vector_Matrix_Product(lambda
->p
, Rays
, A
->p
);
316 values2zz(A
->p
, vertex
, dim
);
321 /* Compute the lattice points in the vertex cone at "values" with rays "rays".
322 * The lattice points are returned in "vertex".
324 * Rays has the generators as rows and so does W.
325 * We first compute { m-v, u_i^* } with m = k W, where k runs through
328 * [k 1] [ d1*W 0 ] [ U' 0 ] = [k 1] T2
330 * where d1 and d2 are the denominators of v and U^{-1}=U'/d2.
331 * Then lambda = { k } (componentwise)
332 * We compute x - floor(x) = {x} = { a/b } as fdiv_r(a,b)/b
333 * For open rays/facets, we need values in (0,1] rather than [0,1),
334 * so we compute {{x}} = x - ceil(x-1) = a/b - ceil((a-b)/b)
335 * = (a - b cdiv_q(a-b,b) - b + b)/b
336 * = (cdiv_r(a,b)+b)/b
337 * Finally, we compute v + lambda * U
338 * The denominator of lambda can be d1*d2, that of lambda2 = lambda*U
339 * can be at most d1, since it is integer if v = 0.
340 * The denominator of v + lambda2 is 1.
342 void lattice_point(Value
* values
, const mat_ZZ
& rays
, mat_ZZ
& vertex
,
343 unsigned long det
, int *closed
)
345 unsigned dim
= rays
.NumRows();
346 vertex
.SetDims(det
, dim
);
348 lattice_point(values
, rays
, vertex
[0], closed
);
351 Matrix
* Rays
= rays2matrix2(rays
);
353 Smith(Rays
, &U
, &W
, &D
);
357 Matrix
*T
= Matrix_Alloc(W
->NbRows
+1, W
->NbColumns
+1);
358 for (int i
= 0; i
< W
->NbRows
; ++i
)
359 Vector_Scale(W
->p
[i
], T
->p
[i
], values
[dim
], W
->NbColumns
);
363 value_set_si(tmp
, -1);
364 Vector_Scale(values
, T
->p
[dim
], tmp
, dim
);
366 value_assign(T
->p
[dim
][dim
], values
[dim
]);
368 Rays
= rays2matrix(rays
);
369 Matrix
*inv
= Matrix_Alloc(Rays
->NbRows
, Rays
->NbColumns
);
370 int ok
= Matrix_Inverse(Rays
, inv
);
374 Matrix
*T2
= Matrix_Alloc(dim
+1, dim
+1);
375 Matrix_Product(T
, inv
, T2
);
378 Rays
= rays2matrix(rays
);
380 Vector
*k
= Vector_Alloc(dim
+1);
381 value_set_si(k
->p
[dim
], 1);
382 Vector
*lambda
= Vector_Alloc(dim
+1);
383 Vector
*lambda2
= Vector_Alloc(dim
+1);
384 for (unsigned long i
= 0; i
< det
; ++i
) {
385 unsigned long val
= i
;
386 for (int j
= 0; j
< dim
; ++j
) {
387 value_set_si(k
->p
[j
], val
% mpz_get_ui(D
->p
[j
][j
]));
388 val
/= mpz_get_ui(D
->p
[j
][j
]);
390 Vector_Matrix_Product(k
->p
, T2
, lambda
->p
);
391 for (int j
= 0; j
< dim
; ++j
)
392 if (!closed
|| closed
[j
])
393 mpz_fdiv_r(lambda
->p
[j
], lambda
->p
[j
], lambda
->p
[dim
]);
395 mpz_cdiv_r(lambda
->p
[j
], lambda
->p
[j
], lambda
->p
[dim
]);
396 value_addto(lambda
->p
[j
], lambda
->p
[j
], lambda
->p
[dim
]);
398 Vector_Matrix_Product(lambda
->p
, Rays
, lambda2
->p
);
399 for (int j
= 0; j
< dim
; ++j
)
400 assert(mpz_divisible_p(lambda2
->p
[j
], inv
->p
[dim
][dim
]));
401 Vector_AntiScale(lambda2
->p
, lambda2
->p
, inv
->p
[dim
][dim
], dim
+1);
402 Vector_Add(lambda2
->p
, values
, lambda2
->p
, dim
);
403 for (int j
= 0; j
< dim
; ++j
)
404 assert(mpz_divisible_p(lambda2
->p
[j
], values
[dim
]));
405 Vector_AntiScale(lambda2
->p
, lambda2
->p
, values
[dim
], dim
+1);
406 values2zz(lambda2
->p
, vertex
[i
], dim
);
410 Vector_Free(lambda2
);
418 static void vertex_period(
419 const mat_ZZ
& rays
, vec_ZZ
& lambda
, Matrix
*T
,
420 Value lcm
, int p
, Vector
*val
,
421 evalue
*E
, evalue
* ev
,
424 unsigned nparam
= T
->NbRows
- 1;
425 unsigned dim
= rays
.NumRows();
432 Vector
* values
= Vector_Alloc(dim
+ 1);
433 Vector_Matrix_Product(val
->p
, T
, values
->p
);
434 value_assign(values
->p
[dim
], lcm
);
435 lattice_point(values
->p
, rays
, vertex
, NULL
);
436 num
= vertex
* lambda
;
441 zz2value(num
, ev
->x
.n
);
442 value_assign(ev
->d
, lcm
);
449 values2zz(T
->p
[p
], vertex
, dim
);
450 nump
= vertex
* lambda
;
451 if (First_Non_Zero(val
->p
, p
) == -1) {
452 value_assign(tmp
, lcm
);
453 evalue
*ET
= term(p
, nump
, &tmp
);
455 free_evalue_refs(ET
);
459 value_assign(tmp
, lcm
);
460 if (First_Non_Zero(T
->p
[p
], dim
) != -1)
461 Vector_Gcd(T
->p
[p
], dim
, &tmp
);
463 if (value_lt(tmp
, lcm
)) {
466 value_division(tmp
, lcm
, tmp
);
467 value_set_si(ev
->d
, 0);
468 ev
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
469 value2zz(tmp
, count
);
471 value_decrement(tmp
, tmp
);
473 ZZ new_offset
= offset
- count
* nump
;
474 value_assign(val
->p
[p
], tmp
);
475 vertex_period(rays
, lambda
, T
, lcm
, p
+1, val
, E
,
476 &ev
->x
.p
->arr
[VALUE_TO_INT(tmp
)], new_offset
);
477 } while (value_pos_p(tmp
));
479 vertex_period(rays
, lambda
, T
, lcm
, p
+1, val
, E
, ev
, offset
);
483 /* Returns the power of (t+1) in the term of a rational generating function,
484 * i.e., the scalar product of the actual lattice point and lambda.
485 * The lattice point is the unique lattice point in the fundamental parallelepiped
486 * of the unimodual cone i shifted to the parametric vertex W/lcm.
488 * The rows of W refer to the coordinates of the vertex
489 * The first nparam columns are the coefficients of the parameters
490 * and the final column is the constant term.
491 * lcm is the common denominator of all coefficients.
493 * PD is the parameter domain, which, if != NULL, may be used to simply the
494 * resulting expression.
496 static evalue
* lattice_point_fractional(const mat_ZZ
& rays
, vec_ZZ
& lambda
,
497 Matrix
*W
, Value lcm
, Polyhedron
*PD
)
499 unsigned nparam
= W
->NbColumns
- 1;
501 Matrix
* Rays
= rays2matrix2(rays
);
502 Matrix
*T
= Transpose(Rays
);
503 Matrix
*T2
= Matrix_Copy(T
);
504 Matrix
*inv
= Matrix_Alloc(T2
->NbRows
, T2
->NbColumns
);
505 int ok
= Matrix_Inverse(T2
, inv
);
510 matrix2zz(W
, vertex
, W
->NbRows
, W
->NbColumns
);
513 num
= lambda
* vertex
;
515 evalue
*EP
= multi_monom(num
);
520 value_set_si(tmp
.x
.n
, 1);
521 value_assign(tmp
.d
, lcm
);
525 Matrix
*L
= Matrix_Alloc(inv
->NbRows
, W
->NbColumns
);
526 Matrix_Product(inv
, W
, L
);
529 matrix2zz(T
, RT
, T
->NbRows
, T
->NbColumns
);
532 vec_ZZ p
= lambda
* RT
;
534 for (int i
= 0; i
< L
->NbRows
; ++i
) {
535 ceil_mod(L
->p
[i
], nparam
+1, lcm
, p
[i
], EP
, PD
);
541 free_evalue_refs(&tmp
);
545 static evalue
* lattice_point_table(const mat_ZZ
& rays
, vec_ZZ
& lambda
, Matrix
*W
,
546 Value lcm
, Polyhedron
*PD
)
548 Matrix
*T
= Transpose(W
);
549 unsigned nparam
= T
->NbRows
- 1;
551 evalue
*EP
= new evalue();
553 evalue_set_si(EP
, 0, 1);
556 Vector
*val
= Vector_Alloc(nparam
+1);
557 value_set_si(val
->p
[nparam
], 1);
558 ZZ
offset(INIT_VAL
, 0);
560 vertex_period(rays
, lambda
, T
, lcm
, 0, val
, EP
, &ev
, offset
);
563 free_evalue_refs(&ev
);
572 evalue
* lattice_point(const mat_ZZ
& rays
, vec_ZZ
& lambda
, Matrix
*W
,
573 Value lcm
, Polyhedron
*PD
, barvinok_options
*options
)
575 if (options
->lookup_table
)
576 return lattice_point_table(rays
, lambda
, W
, lcm
, PD
);
578 return lattice_point_fractional(rays
, lambda
, W
, lcm
, PD
);
581 /* returns the unique lattice point in the fundamental parallelepiped
582 * of the unimodual cone C shifted to the parametric vertex V.
584 * The return values num and E_vertex are such that
585 * coordinate i of this lattice point is equal to
587 * num[i] + E_vertex[i]
589 void lattice_point(Param_Vertices
*V
, const mat_ZZ
& rays
, vec_ZZ
& num
,
590 evalue
**E_vertex
, barvinok_options
*options
)
592 unsigned nparam
= V
->Vertex
->NbColumns
- 2;
593 unsigned dim
= rays
.NumCols();
595 vertex
.SetLength(nparam
+1);
600 value_set_si(lcm
, 1);
602 for (int j
= 0; j
< V
->Vertex
->NbRows
; ++j
) {
603 value_lcm(lcm
, V
->Vertex
->p
[j
][nparam
+1], &lcm
);
606 if (value_notone_p(lcm
)) {
607 Matrix
* mv
= Matrix_Alloc(dim
, nparam
+1);
608 for (int j
= 0 ; j
< dim
; ++j
) {
609 value_division(tmp
, lcm
, V
->Vertex
->p
[j
][nparam
+1]);
610 Vector_Scale(V
->Vertex
->p
[j
], mv
->p
[j
], tmp
, nparam
+1);
613 Matrix
* Rays
= rays2matrix2(rays
);
614 Matrix
*T
= Transpose(Rays
);
615 Matrix
*T2
= Matrix_Copy(T
);
616 Matrix
*inv
= Matrix_Alloc(T2
->NbRows
, T2
->NbColumns
);
617 int ok
= Matrix_Inverse(T2
, inv
);
621 Matrix
*L
= Matrix_Alloc(inv
->NbRows
, mv
->NbColumns
);
622 Matrix_Product(inv
, mv
, L
);
631 evalue
*remainders
[dim
];
632 for (int i
= 0; i
< dim
; ++i
) {
633 remainders
[i
] = evalue_zero();
635 ceil(L
->p
[i
], nparam
+1, lcm
, one
, remainders
[i
], 0, options
);
640 for (int i
= 0; i
< V
->Vertex
->NbRows
; ++i
) {
641 values2zz(mv
->p
[i
], vertex
, nparam
+1);
642 E_vertex
[i
] = multi_monom(vertex
);
645 value_set_si(f
.x
.n
, 1);
646 value_assign(f
.d
, lcm
);
648 emul(&f
, E_vertex
[i
]);
650 for (int j
= 0; j
< dim
; ++j
) {
651 if (value_zero_p(T
->p
[i
][j
]))
655 evalue_copy(&cp
, remainders
[j
]);
656 if (value_notone_p(T
->p
[i
][j
])) {
657 value_set_si(f
.d
, 1);
658 value_assign(f
.x
.n
, T
->p
[i
][j
]);
661 eadd(&cp
, E_vertex
[i
]);
662 free_evalue_refs(&cp
);
665 for (int i
= 0; i
< dim
; ++i
) {
666 free_evalue_refs(remainders
[i
]);
670 free_evalue_refs(&f
);
681 for (int i
= 0; i
< V
->Vertex
->NbRows
; ++i
) {
683 if (First_Non_Zero(V
->Vertex
->p
[i
], nparam
) == -1) {
685 value2zz(V
->Vertex
->p
[i
][nparam
], num
[i
]);
687 values2zz(V
->Vertex
->p
[i
], vertex
, nparam
+1);
688 E_vertex
[i
] = multi_monom(vertex
);