2 #include <barvinok/util.h>
4 #include "lattice_point.h"
9 /* Computes the integer points in the fundamental parallelepiped and
10 * passes them along (in num) to the counter specific (i.e., specialization
11 * specific) add_lattice_points.
13 void counter_base::handle(const mat_ZZ
& rays
, Value
*V
, const QQ
& c
,
14 unsigned long det
, barvinok_options
*options
)
16 Matrix
* Rays
= zz2matrix(rays
);
19 assert(c
.n
== 1 || c
.n
== -1);
20 int sign
= to_int(c
.n
);
22 Matrix_Vector_Product(Rays
, lambda
->p
, den
->p_Init
);
23 for (int k
= 0; k
< dim
; ++k
)
24 if (value_zero_p(den
->p_Init
[k
])) {
28 Inner_Product(lambda
->p
, V
, dim
, &tmp
);
29 lattice_points_fixed(V
, &tmp
, Rays
, den
, num
, det
);
36 add_lattice_points(sign
);
39 static void add_falling_powers(dpoly
& n
, Value degree
)
41 value_increment(n
.coeff
->p
[0], n
.coeff
->p
[0]);
42 if (n
.coeff
->Size
== 1)
45 int min
= n
.coeff
->Size
-1;
46 if (value_posz_p(degree
) && value_cmp_si(degree
, min
) < 0)
47 min
= VALUE_TO_INT(degree
);
51 value_assign(tmp
, degree
);
52 value_addto(n
.coeff
->p
[1], n
.coeff
->p
[1], tmp
);
53 for (int i
= 2; i
<= min
; ++i
) {
54 value_decrement(degree
, degree
);
55 value_multiply(tmp
, tmp
, degree
);
56 mpz_divexact_ui(tmp
, tmp
, i
);
57 value_addto(n
.coeff
->p
[i
], n
.coeff
->p
[i
], tmp
);
62 void counter::add_lattice_points(int sign
)
65 for (int k
= 0; k
< num
->NbRows
; ++k
)
66 add_falling_powers(d
, num
->p_Init
[k
]);
67 dpoly
n(dim
, den
->p_Init
[0], 1);
68 for (int k
= 1; k
< dim
; ++k
) {
69 dpoly
fact(dim
, den
->p_Init
[k
], 1);
72 d
.div(n
, count
, sign
);
78 void tcounter::setup_todd(unsigned dim
)
80 value_set_si(todd
.coeff
->p
[0], 1);
83 value_set_si(d
.coeff
->p
[dim
], 1);
84 for (int i
= dim
-1; i
>= 0; --i
)
85 mpz_mul_ui(d
.coeff
->p
[i
], d
.coeff
->p
[i
+1], i
+2);
87 todd_denom
= todd
.div(d
);
88 /* shift denominators up -> divide by (dim+1)! */
89 for (int i
= todd
.coeff
->Size
-1; i
>= 1; --i
)
90 value_assign(todd_denom
->p
[i
], todd_denom
->p
[i
-1]);
91 value_set_si(todd_denom
->p
[0], 1);
94 void tcounter::adapt_todd(dpoly
& t
, const Value c
)
96 if (t
.coeff
->Size
<= 1)
99 value_multiply(t
.coeff
->p
[1], t
.coeff
->p
[1], tmp
);
100 for (int i
= 2; i
< t
.coeff
->Size
; ++i
) {
101 value_multiply(tmp
, tmp
, c
);
102 value_multiply(t
.coeff
->p
[i
], t
.coeff
->p
[i
], tmp
);
106 void tcounter::add_powers(dpoly
& n
, const Value c
)
108 value_increment(n
.coeff
->p
[0], n
.coeff
->p
[0]);
109 if (n
.coeff
->Size
== 1)
111 value_assign(tmp
, c
);
112 value_addto(n
.coeff
->p
[1], n
.coeff
->p
[1], tmp
);
113 for (int i
= 2; i
< n
.coeff
->Size
; ++i
) {
114 value_multiply(tmp
, tmp
, c
);
115 value_addto(n
.coeff
->p
[i
], n
.coeff
->p
[i
], tmp
);
119 void tcounter::add_lattice_points(int sign
)
122 value_assign(denom
, den
->p_Init
[0]);
123 adapt_todd(t
, den
->p_Init
[0]);
124 for (int k
= 1; k
< dim
; ++k
) {
126 value_multiply(denom
, denom
, den
->p_Init
[k
]);
127 adapt_todd(fact
, den
->p_Init
[k
]);
132 for (int k
= 0; k
< num
->NbRows
; ++k
)
133 add_powers(n
, num
->p_Init
[k
]);
135 for (int i
= 0; i
< n
.coeff
->Size
; ++i
)
136 value_multiply(n
.coeff
->p
[i
], n
.coeff
->p
[i
], todd_denom
->p
[i
]);
137 value_multiply(denom
, denom
, todd_denom
->p
[todd_denom
->Size
-1]);
139 value_set_si(tmp
, 1);
140 for (int i
= 2; i
< n
.coeff
->Size
; ++i
) {
141 mpz_mul_ui(tmp
, tmp
, i
);
142 mpz_divexact(n
.coeff
->p
[i
], n
.coeff
->p
[i
], tmp
);
145 value_multiply(tmp
, t
.coeff
->p
[0], n
.coeff
->p
[n
.coeff
->Size
-1]);
146 for (int i
= 1; i
< n
.coeff
->Size
; ++i
)
147 value_addmul(tmp
, t
.coeff
->p
[i
], n
.coeff
->p
[n
.coeff
->Size
-1-i
]);
149 value_assign(mpq_numref(tcount
), tmp
);
150 value_assign(mpq_denref(tcount
), denom
);
151 mpq_canonicalize(tcount
);
153 mpq_sub(count
, count
, tcount
);
155 mpq_add(count
, count
, tcount
);