lattice_point.cc

   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     value_gcd(gcd, gcd, *m);
 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 /* 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.
 162  *
 163  * Modifies coef argument !
 164  */
 165 evalue *fractional_part(Value *coef, Value denom, int nvar, Polyhedron *PD)
 166 {
 167     Value m;
 168     value_init(m);
 169     evalue *EP = evalue_zero();
 170     int sign = 1;
 171
 172     value_assign(m, denom);
 173     int j = normal_mod(coef, nvar+1, &m);
 174
 175     if (j == nvar+1) {
 176         value_clear(m);
 177         return EP;
 178     }
 179
 180     vec_ZZ num;
 181     values2zz(coef, num, nvar+1);
 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 < 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;
 204     }
 205
 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);
 215
 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);
 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     return EP;
 250 }
 251
 252 static evalue *ceil(Value *coef, int len, Value d,
 253                     barvinok_options *options)
 254 {
 255     evalue *c;
 256
 257     Vector_Oppose(coef, coef, len);
 258     c = fractional_part(coef, d, len-1, NULL);
 259     if (options->lookup_table)
 260         evalue_mod2table(c, len-1);
 261     return c;
 262 }
 263
 264 evalue* bv_ceil3(Value *coef, int len, Value d, Polyhedron *P)
 265 {
 266     Vector *val = Vector_Alloc(len);
 267
 268     Value t;
 269     value_init(t);
 270     value_set_si(t, -1);
 271     Vector_Scale(coef, val->p, t, len);
 272     value_absolute(t, d);
 273
 274     vec_ZZ num;
 275     values2zz(val->p, num, len);
 276     evalue *EP = multi_monom(num);
 277
 278     evalue tmp;
 279     value_init(tmp.d);
 280     value_init(tmp.x.n);
 281     value_set_si(tmp.x.n, 1);
 282     value_assign(tmp.d, t);
 283
 284     emul(&tmp, EP);
 285
 286     Vector_Oppose(val->p, val->p, len);
 287     evalue *f = fractional_part(val->p, t, len-1, P);
 288     value_clear(t);
 289
 290     eadd(f, EP);
 291     evalue_free(f);
 292
 293     /* copy EP to malloc'ed evalue */
 294     evalue *E = ALLOC(evalue);
 295     *E = *EP;
 296     delete EP;
 297
 298     free_evalue_refs(&tmp);
 299     Vector_Free(val);
 300
 301     return E;
 302 }
 303
 304 void lattice_point_fixed(Value *vertex, Value *vertex_res,
 305                          Matrix *Rays, Matrix *Rays_res,
 306                          Value *point)
 307 {
 308     unsigned dim = Rays->NbRows;
 309     if (value_one_p(vertex[dim]))
 310         Vector_Copy(vertex_res, point, Rays_res->NbColumns);
 311     else {
 312         Matrix *R2 = Matrix_Copy(Rays);
 313         Matrix *inv = Matrix_Alloc(Rays->NbRows, Rays->NbColumns);
 314         int ok = Matrix_Inverse(R2, inv);
 315         assert(ok);
 316         Matrix_Free(R2);
 317         Vector *lambda = Vector_Alloc(dim);
 318         Vector_Matrix_Product(vertex, inv, lambda->p);
 319         Matrix_Free(inv);
 320         for (int j = 0; j < dim; ++j)
 321             mpz_cdiv_q(lambda->p[j], lambda->p[j], vertex[dim]);
 322         Vector_Matrix_Product(lambda->p, Rays_res, point);
 323         Vector_Free(lambda);
 324     }
 325 }
 326
 327 static Matrix *Matrix_AddRowColumn(Matrix *M)
 328 {
 329     Matrix *M2 = Matrix_Alloc(M->NbRows+1, M->NbColumns+1);
 330     for (int i = 0; i < M->NbRows; ++i)
 331         Vector_Copy(M->p[i], M2->p[i], M->NbColumns);
 332     value_set_si(M2->p[M->NbRows][M->NbColumns], 1);
 333     return M2;
 334 }
 335
 336 #define FORALL_COSETS(det,D,i,k)                                        \
 337     do {                                                                \
 338         Vector *k = Vector_Alloc(D->NbRows+1);                          \
 339         value_set_si(k->p[D->NbRows], 1);                               \
 340         for (unsigned long i = 0; i < det; ++i) {                       \
 341             if (i)                                                      \
 342                 for (int j = D->NbRows-1; j >= 0; --j) {                \
 343                     value_increment(k->p[j], k->p[j]);                  \
 344                     if (value_eq(k->p[j], D->p[j][j]))                  \
 345                         value_set_si(k->p[j], 0);                       \
 346                     else                                                \
 347                         break;                                          \
 348                 }                                                       \
 349             {
 350 #define END_FORALL_COSETS                                               \
 351             }                                                           \
 352         }                                                               \
 353         Vector_Free(k);                                                 \
 354     } while(0);
 355
 356 /* Compute the lattice points in the vertex cone at "values" with rays "rays".
 357  * The lattice points are returned in "vertex".
 358  *
 359  * Rays has the generators as rows and so does W.
 360  * We first compute { m-v, u_i^* } with m = k W, where k runs through
 361  * the cosets.
 362  * We compute
 363  * [k 1] [ d1*W  0 ] [  U'  0 ] = [k 1] T2
 364  *       [ -v   d1 ] [  0  d2 ]
 365  * where d1 and d2 are the denominators of v and U^{-1}=U'/d2.
 366  * Then  lambda = { k } (componentwise)
 367  * We compute x - floor(x) = {x} = { a/b } as fdiv_r(a,b)/b
 368  * For open rays/facets, we need values in (0,1] rather than [0,1),
 369  * so we compute {{x}} = x - ceil(x-1) = a/b - ceil((a-b)/b)
 370  *                     = (a - b cdiv_q(a-b,b) - b + b)/b
 371  *                     = (cdiv_r(a,b)+b)/b
 372  * Finally, we compute v + lambda * U
 373  * The denominator of lambda can be d1*d2, that of lambda2 = lambda*U
 374  * can be at most d1, since it is integer if v = 0.
 375  * The denominator of v + lambda2 is 1.
 376  *
 377  * The _res variants of the input variables may have been multiplied with
 378  * a (list of) nonorthogonal vector(s) and may therefore have fewer columns
 379  * than their original counterparts.
 380  */
 381 void lattice_points_fixed(Value *vertex, Value *vertex_res,
 382                           Matrix *Rays, Matrix *Rays_res, Matrix *points,
 383                           unsigned long det)
 384 {
 385     unsigned dim = Rays->NbRows;
 386     if (det == 1) {
 387         lattice_point_fixed(vertex, vertex_res, Rays, Rays_res,
 388                             points->p[0]);
 389         return;
 390     }
 391     Matrix *U, *W, *D;
 392     Smith(Rays, &U, &W, &D);
 393     Matrix_Free(U);
 394
 395     /* Sanity check */
 396     unsigned long det2 = 1;
 397     for (int i = 0 ; i < D->NbRows; ++i)
 398         det2 *= mpz_get_ui(D->p[i][i]);
 399     assert(det == det2);
 400
 401     Matrix *T = Matrix_Alloc(W->NbRows+1, W->NbColumns+1);
 402     for (int i = 0; i < W->NbRows; ++i)
 403         Vector_Scale(W->p[i], T->p[i], vertex[dim], W->NbColumns);
 404     Matrix_Free(W);
 405     Vector_Oppose(vertex, T->p[dim], dim);
 406     value_assign(T->p[dim][dim], vertex[dim]);
 407
 408     Matrix *R2 = Matrix_AddRowColumn(Rays);
 409     Matrix *inv = Matrix_Alloc(R2->NbRows, R2->NbColumns);
 410     int ok = Matrix_Inverse(R2, inv);
 411     assert(ok);
 412     Matrix_Free(R2);
 413
 414     Matrix *T2 = Matrix_Alloc(dim+1, dim+1);
 415     Matrix_Product(T, inv, T2);
 416     Matrix_Free(T);
 417
 418     Vector *lambda = Vector_Alloc(dim+1);
 419     Vector *lambda2 = Vector_Alloc(Rays_res->NbColumns);
 420     FORALL_COSETS(det, D, i, k)
 421         Vector_Matrix_Product(k->p, T2, lambda->p);
 422         for (int j = 0; j < dim; ++j)
 423             mpz_fdiv_r(lambda->p[j], lambda->p[j], lambda->p[dim]);
 424         Vector_Matrix_Product(lambda->p, Rays_res, lambda2->p);
 425         for (int j = 0; j < lambda2->Size; ++j)
 426             assert(mpz_divisible_p(lambda2->p[j], inv->p[dim][dim]));
 427         Vector_AntiScale(lambda2->p, lambda2->p, inv->p[dim][dim], lambda2->Size);
 428         Vector_Add(lambda2->p, vertex_res, lambda2->p, lambda2->Size);
 429         for (int j = 0; j < lambda2->Size; ++j)
 430             assert(mpz_divisible_p(lambda2->p[j], vertex[dim]));
 431         Vector_AntiScale(lambda2->p, points->p[i], vertex[dim], lambda2->Size);
 432     END_FORALL_COSETS
 433     Vector_Free(lambda);
 434     Vector_Free(lambda2);
 435     Matrix_Free(D);
 436     Matrix_Free(inv);
 437
 438     Matrix_Free(T2);
 439 }
 440
 441 /* Returns the power of (t+1) in the term of a rational generating function,
 442  * i.e., the scalar product of the actual lattice point and lambda.
 443  * The lattice point is the unique lattice point in the fundamental parallelepiped
 444  * of the unimodual cone i shifted to the parametric vertex W/lcm.
 445  *
 446  * The rows of W refer to the coordinates of the vertex
 447  * The first nparam columns are the coefficients of the parameters
 448  * and the final column is the constant term.
 449  * lcm is the common denominator of all coefficients.
 450  */
 451 static evalue **lattice_point_fractional(const mat_ZZ& rays, vec_ZZ& lambda,
 452                                          Matrix *V,
 453                                          unsigned long det)
 454 {
 455     unsigned nparam = V->NbColumns-2;
 456     evalue **E = new evalue *[det];
 457
 458     Matrix* Rays = zz2matrix(rays);
 459     Matrix *T = Transpose(Rays);
 460     Matrix *T2 = Matrix_Copy(T);
 461     Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
 462     int ok = Matrix_Inverse(T2, inv);
 463     assert(ok);
 464     Matrix_Free(T2);
 465     mat_ZZ vertex;
 466     matrix2zz(V, vertex, V->NbRows, V->NbColumns-1);
 467
 468     vec_ZZ num;
 469     num = lambda * vertex;
 470
 471     evalue *EP = multi_monom(num);
 472
 473     evalue_div(EP, V->p[0][nparam+1]);
 474
 475     Matrix *L = Matrix_Alloc(inv->NbRows, V->NbColumns);
 476     Matrix_Product(inv, V, L);
 477
 478     mat_ZZ RT;
 479     matrix2zz(T, RT, T->NbRows, T->NbColumns);
 480     Matrix_Free(T);
 481
 482     vec_ZZ p = lambda * RT;
 483
 484     Value tmp;
 485     value_init(tmp);
 486
 487     if (det == 1) {
 488         for (int i = 0; i < L->NbRows; ++i) {
 489             evalue *f;
 490             Vector_Oppose(L->p[i], L->p[i], nparam+1);
 491             f = fractional_part(L->p[i], V->p[i][nparam+1], nparam, NULL);
 492             zz2value(p[i], tmp);
 493             evalue_mul(f, tmp);
 494             eadd(f, EP);
 495             evalue_free(f);
 496         }
 497         E[0] = EP;
 498     } else {
 499         for (int i = 0; i < L->NbRows; ++i)
 500             value_assign(L->p[i][nparam+1], V->p[i][nparam+1]);
 501
 502         Value denom;
 503         value_init(denom);
 504         mpz_set_ui(denom, det);
 505         value_multiply(denom, L->p[0][nparam+1], denom);
 506
 507         Matrix *U, *W, *D;
 508         Smith(Rays, &U, &W, &D);
 509         Matrix_Free(U);
 510
 511         /* Sanity check */
 512         unsigned long det2 = 1;
 513         for (int i = 0 ; i < D->NbRows; ++i)
 514             det2 *= mpz_get_ui(D->p[i][i]);
 515         assert(det == det2);
 516
 517         Matrix_Transposition(inv);
 518         Matrix *T2 = Matrix_Alloc(W->NbRows, inv->NbColumns);
 519         Matrix_Product(W, inv, T2);
 520         Matrix_Free(W);
 521
 522         unsigned dim = D->NbRows;
 523         Vector *lambda = Vector_Alloc(dim);
 524
 525         Vector *row = Vector_Alloc(nparam+1);
 526         FORALL_COSETS(det, D, i, k)
 527             Vector_Matrix_Product(k->p, T2, lambda->p);
 528             E[i] = new evalue();
 529             value_init(E[i]->d);
 530             evalue_copy(E[i], EP);
 531             for (int j = 0; j < L->NbRows; ++j) {
 532                 evalue *f;
 533                 Vector_Oppose(L->p[j], row->p, nparam+1);
 534                 value_addmul(row->p[nparam], L->p[j][nparam+1], lambda->p[j]);
 535                 f = fractional_part(row->p, denom, nparam, NULL);
 536                 zz2value(p[j], tmp);
 537                 evalue_mul(f, tmp);
 538                 eadd(f, E[i]);
 539                 evalue_free(f);
 540             }
 541         END_FORALL_COSETS
 542         Vector_Free(row);
 543
 544         Vector_Free(lambda);
 545         Matrix_Free(T2);
 546         Matrix_Free(D);
 547
 548         value_clear(denom);
 549         free_evalue_refs(EP);
 550         delete EP;
 551     }
 552     value_clear(tmp);
 553
 554     Matrix_Free(Rays);
 555     Matrix_Free(L);
 556     Matrix_Free(inv);
 557
 558     return E;
 559 }
 560
 561 static evalue **lattice_point(const mat_ZZ& rays, vec_ZZ& lambda,
 562                               Param_Vertices *V,
 563                               unsigned long det,
 564                               barvinok_options *options)
 565 {
 566     evalue **lp = lattice_point_fractional(rays, lambda, V->Vertex, det);
 567     if (options->lookup_table) {
 568         for (int i = 0; i < det; ++i)
 569             evalue_mod2table(lp[i], V->Vertex->NbColumns-2);
 570     }
 571     return lp;
 572 }
 573
 574 /* returns the unique lattice point in the fundamental parallelepiped
 575  * of the unimodual cone C shifted to the parametric vertex V.
 576  *
 577  * The return values num and E_vertex are such that
 578  * coordinate i of this lattice point is equal to
 579  *
 580  *          num[i] + E_vertex[i]
 581  */
 582 void lattice_point(Param_Vertices *V, const mat_ZZ& rays, vec_ZZ& num,
 583                    evalue **E_vertex, barvinok_options *options)
 584 {
 585     unsigned nparam = V->Vertex->NbColumns - 2;
 586     unsigned dim = rays.NumCols();
 587
 588     /* It doesn't really make sense to call lattice_point when dim == 0,
 589      * but apparently it happens from indicator_constructor in lexmin.
 590      */
 591     if (!dim)
 592         return;
 593
 594     vec_ZZ vertex;
 595     vertex.SetLength(nparam+1);
 596
 597     Value tmp;
 598     value_init(tmp);
 599
 600     assert(V->Vertex->NbRows > 0);
 601     Param_Vertex_Common_Denominator(V);
 602
 603     if (value_notone_p(V->Vertex->p[0][nparam+1])) {
 604         Matrix* Rays = zz2matrix(rays);
 605         Matrix *T = Transpose(Rays);
 606         Matrix *T2 = Matrix_Copy(T);
 607         Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
 608         int ok = Matrix_Inverse(T2, inv);
 609         assert(ok);
 610         Matrix_Free(Rays);
 611         Matrix_Free(T2);
 612         /* temporarily ignore (common) denominator */
 613         V->Vertex->NbColumns--;
 614         Matrix *L = Matrix_Alloc(inv->NbRows, V->Vertex->NbColumns);
 615         Matrix_Product(inv, V->Vertex, L);
 616         V->Vertex->NbColumns++;
 617         Matrix_Free(inv);
 618
 619         evalue f;
 620         value_init(f.d);
 621         value_init(f.x.n);
 622
 623         evalue *remainders[dim];
 624         for (int i = 0; i < dim; ++i)
 625             remainders[i] = ceil(L->p[i], nparam+1, V->Vertex->p[0][nparam+1],
 626                                  options);
 627         Matrix_Free(L);
 628
 629
 630         for (int i = 0; i < V->Vertex->NbRows; ++i) {
 631             values2zz(V->Vertex->p[i], vertex, nparam+1);
 632             E_vertex[i] = multi_monom(vertex);
 633             num[i] = 0;
 634
 635             value_set_si(f.x.n, 1);
 636             value_assign(f.d, V->Vertex->p[0][nparam+1]);
 637
 638             emul(&f, E_vertex[i]);
 639
 640             for (int j = 0; j < dim; ++j) {
 641                 if (value_zero_p(T->p[i][j]))
 642                     continue;
 643                 evalue cp;
 644                 value_init(cp.d);
 645                 evalue_copy(&cp, remainders[j]);
 646                 if (value_notone_p(T->p[i][j])) {
 647                     value_set_si(f.d, 1);
 648                     value_assign(f.x.n, T->p[i][j]);
 649                     emul(&f, &cp);
 650                 }
 651                 eadd(&cp, E_vertex[i]);
 652                 free_evalue_refs(&cp);
 653             }
 654         }
 655         for (int i = 0; i < dim; ++i)
 656             evalue_free(remainders[i]);
 657
 658         free_evalue_refs(&f);
 659
 660         Matrix_Free(T);
 661         value_clear(tmp);
 662         return;
 663     }
 664     value_clear(tmp);
 665
 666     for (int i = 0; i < V->Vertex->NbRows; ++i) {
 667         /* fixed value */
 668         if (First_Non_Zero(V->Vertex->p[i], nparam) == -1) {
 669             E_vertex[i] = 0;
 670             value2zz(V->Vertex->p[i][nparam], num[i]);
 671         } else {
 672             values2zz(V->Vertex->p[i], vertex, nparam+1);
 673             E_vertex[i] = multi_monom(vertex);
 674             num[i] = 0;
 675         }
 676     }
 677 }
 678
 679 static int lattice_point_fixed(Param_Vertices* V, const mat_ZZ& rays,
 680     vec_ZZ& lambda, term_info* term, unsigned long det)
 681 {
 682     unsigned nparam = V->Vertex->NbColumns - 2;
 683     unsigned dim = rays.NumCols();
 684
 685     for (int i = 0; i < dim; ++i)
 686         if (First_Non_Zero(V->Vertex->p[i], nparam) != -1)
 687             return 0;
 688
 689     Vector *fixed = Vector_Alloc(dim+1);
 690     for (int i = 0; i < dim; ++i)
 691         value_assign(fixed->p[i], V->Vertex->p[i][nparam]);
 692     value_assign(fixed->p[dim], V->Vertex->p[0][nparam+1]);
 693
 694     mat_ZZ vertex;
 695     Matrix *points = Matrix_Alloc(det, dim);
 696     Matrix* Rays = zz2matrix(rays);
 697     lattice_points_fixed(fixed->p, fixed->p, Rays, Rays, points, det);
 698     Matrix_Free(Rays);
 699     matrix2zz(points, vertex, points->NbRows, points->NbColumns);
 700     Matrix_Free(points);
 701     term->E = NULL;
 702     term->constant = vertex * lambda;
 703     Vector_Free(fixed);
 704
 705     return 1;
 706 }
 707
 708 /* Returns the power of (t+1) in the term of a rational generating function,
 709  * i.e., the scalar product of the actual lattice point and lambda.
 710  * The lattice point is the unique lattice point in the fundamental parallelepiped
 711  * of the unimodual cone i shifted to the parametric vertex V.
 712  *
 713  * The result is returned in term.
 714  */
 715 void lattice_point(Param_Vertices* V, const mat_ZZ& rays, vec_ZZ& lambda,
 716     term_info* term, unsigned long det,
 717     barvinok_options *options)
 718 {
 719     unsigned nparam = V->Vertex->NbColumns - 2;
 720     unsigned dim = rays.NumCols();
 721     mat_ZZ vertex;
 722     vertex.SetDims(V->Vertex->NbRows, nparam+1);
 723
 724     Param_Vertex_Common_Denominator(V);
 725
 726     if (lattice_point_fixed(V, rays, lambda, term, det))
 727         return;
 728
 729     if (det != 1 || value_notone_p(V->Vertex->p[0][nparam+1])) {
 730         term->E = lattice_point(rays, lambda, V, det, options);
 731         return;
 732     }
 733     for (int i = 0; i < V->Vertex->NbRows; ++i) {
 734         assert(value_one_p(V->Vertex->p[i][nparam+1]));  // for now
 735         values2zz(V->Vertex->p[i], vertex[i], nparam+1);
 736     }
 737
 738     vec_ZZ num;
 739     num = lambda * vertex;
 740
 741     int nn = 0;
 742     for (int j = 0; j < nparam; ++j)
 743         if (num[j] != 0)
 744             ++nn;
 745     if (nn >= 1) {
 746         term->E = new evalue *[1];
 747         term->E[0] = multi_monom(num);
 748     } else {
 749         term->E = NULL;
 750         term->constant.SetLength(1);
 751         term->constant[0] = num[nparam];
 752     }
 753 }