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 = 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;
  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         evalue_free(E);
 241     }
 242
 243     free_evalue_refs(&tmp);
 244
 245 out:
 246     value_clear(m);
 247
 248     return EP;
 249 }
 250
 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.
 256  *
 257  * Modifies coef argument !
 258  */
 259 evalue *ceiling(Value *coef, Value denom, int nvar, Polyhedron *PD)
 260 {
 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;
 268 }
 269
 270 static evalue *ceil(Value *coef, int len, Value d,
 271                     barvinok_options *options)
 272 {
 273     evalue *c;
 274
 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;
 280 }
 281
 282 void lattice_point_fixed(Value *vertex, Value *vertex_res,
 283                          Matrix *Rays, Matrix *Rays_res,
 284                          Value *point)
 285 {
 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);
 302     }
 303 }
 304
 305 static Matrix *Matrix_AddRowColumn(Matrix *M)
 306 {
 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;
 312 }
 313
 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;                                          \
 326                 }                                                       \
 327             {
 328 #define END_FORALL_COSETS                                               \
 329             }                                                           \
 330         }                                                               \
 331         Vector_Free(k);                                                 \
 332     } while(0);
 333
 334 /* Compute the lattice points in the vertex cone at "values" with rays "rays".
 335  * The lattice points are returned in "vertex".
 336  *
 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.
 354  *
 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.
 358  */
 359 void lattice_points_fixed(Value *vertex, Value *vertex_res,
 360                           Matrix *Rays, Matrix *Rays_res, Matrix *points,
 361                           unsigned long det)
 362 {
 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;
 368     }
 369     Matrix *U, *W, *D;
 370     Smith(Rays, &U, &W, &D);
 371     Matrix_Free(U);
 372
 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);
 378
 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]);
 385
 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);
 391
 392     Matrix *T2 = Matrix_Alloc(dim+1, dim+1);
 393     Matrix_Product(T, inv, T2);
 394     Matrix_Free(T);
 395
 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);
 415
 416     Matrix_Free(T2);
 417 }
 418
 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.
 423  *
 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.
 428  */
 429 static evalue **lattice_point_fractional(const mat_ZZ& rays, vec_ZZ& lambda,
 430                                          Matrix *V,
 431                                          unsigned long det)
 432 {
 433     unsigned nparam = V->NbColumns-2;
 434     evalue **E = new evalue *[det];
 435
 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);
 445
 446     vec_ZZ num;
 447     num = lambda * vertex;
 448
 449     evalue *EP = multi_monom(num);
 450
 451     evalue_div(EP, V->p[0][nparam+1]);
 452
 453     Matrix *L = Matrix_Alloc(inv->NbRows, V->NbColumns);
 454     Matrix_Product(inv, V, L);
 455
 456     mat_ZZ RT;
 457     matrix2zz(T, RT, T->NbRows, T->NbColumns);
 458     Matrix_Free(T);
 459
 460     vec_ZZ p = lambda * RT;
 461
 462     Value tmp;
 463     value_init(tmp);
 464
 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);
 474         }
 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]);
 479
 480         Value denom;
 481         value_init(denom);
 482         mpz_set_ui(denom, det);
 483         value_multiply(denom, L->p[0][nparam+1], denom);
 484
 485         Matrix *U, *W, *D;
 486         Smith(Rays, &U, &W, &D);
 487         Matrix_Free(U);
 488
 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);
 494
 495         Matrix_Transposition(inv);
 496         Matrix *T2 = Matrix_Alloc(W->NbRows, inv->NbColumns);
 497         Matrix_Product(W, inv, T2);
 498         Matrix_Free(W);
 499
 500         unsigned dim = D->NbRows;
 501         Vector *lambda = Vector_Alloc(dim);
 502
 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);
 518             }
 519         END_FORALL_COSETS
 520         Vector_Free(row);
 521
 522         Vector_Free(lambda);
 523         Matrix_Free(T2);
 524         Matrix_Free(D);
 525
 526         value_clear(denom);
 527         evalue_free(EP);
 528     }
 529     value_clear(tmp);
 530
 531     Matrix_Free(Rays);
 532     Matrix_Free(L);
 533     Matrix_Free(inv);
 534
 535     return E;
 536 }
 537
 538 static evalue **lattice_point(const mat_ZZ& rays, vec_ZZ& lambda,
 539                               Param_Vertices *V,
 540                               unsigned long det,
 541                               barvinok_options *options)
 542 {
 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);
 547     }
 548     return lp;
 549 }
 550
 551 Matrix *relative_coordinates(Param_Vertices *V, Matrix *basis)
 552 {
 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);
 559
 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);
 567
 568     return L;
 569 }
 570
 571 /* returns the unique lattice point in the fundamental parallelepiped
 572  * of the unimodual cone C shifted to the parametric vertex V.
 573  *
 574  * The return values num and E_vertex are such that
 575  * coordinate i of this lattice point is equal to
 576  *
 577  *          num[i] + E_vertex[i]
 578  */
 579 void lattice_point(Param_Vertices *V, const mat_ZZ& rays, vec_ZZ& num,
 580                    evalue **E_vertex, barvinok_options *options)
 581 {
 582     unsigned nparam = V->Vertex->NbColumns - 2;
 583     unsigned dim = rays.NumCols();
 584
 585     /* It doesn't really make sense to call lattice_point when dim == 0,
 586      * but apparently it happens from indicator_constructor in lexmin.
 587      */
 588     if (!dim)
 589         return;
 590
 591     vec_ZZ vertex;
 592     vertex.SetLength(nparam+1);
 593
 594     Value tmp;
 595     value_init(tmp);
 596
 597     assert(V->Vertex->NbRows > 0);
 598     Param_Vertex_Common_Denominator(V);
 599
 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);
 605
 606         evalue f;
 607         value_init(f.d);
 608         value_init(f.x.n);
 609
 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);
 615
 616
 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;
 621
 622             value_set_si(f.x.n, 1);
 623             value_assign(f.d, V->Vertex->p[0][nparam+1]);
 624
 625             emul(&f, E_vertex[i]);
 626
 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);
 637                 }
 638                 eadd(&cp, E_vertex[i]);
 639                 free_evalue_refs(&cp);
 640             }
 641         }
 642         for (int i = 0; i < dim; ++i)
 643             evalue_free(remainders[i]);
 644         delete [] remainders;
 645
 646         free_evalue_refs(&f);
 647
 648         Matrix_Free(T);
 649         value_clear(tmp);
 650         return;
 651     }
 652     value_clear(tmp);
 653
 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;
 663         }
 664     }
 665 }
 666
 667 static int lattice_point_fixed(Param_Vertices* V, const mat_ZZ& rays,
 668     vec_ZZ& lambda, term_info* term, unsigned long det)
 669 {
 670     unsigned nparam = V->Vertex->NbColumns - 2;
 671     unsigned dim = rays.NumCols();
 672
 673     for (int i = 0; i < dim; ++i)
 674         if (First_Non_Zero(V->Vertex->p[i], nparam) != -1)
 675             return 0;
 676
 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]);
 681
 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);
 692
 693     return 1;
 694 }
 695
 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.
 700  *
 701  * The result is returned in term.
 702  */
 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)
 706 {
 707     unsigned nparam = V->Vertex->NbColumns - 2;
 708     unsigned dim = rays.NumCols();
 709     mat_ZZ vertex;
 710     vertex.SetDims(V->Vertex->NbRows, nparam+1);
 711
 712     Param_Vertex_Common_Denominator(V);
 713
 714     if (lattice_point_fixed(V, rays, lambda, term, det))
 715         return;
 716
 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;
 720     }
 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);
 724     }
 725
 726     vec_ZZ num;
 727     num = lambda * vertex;
 728
 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];
 740     }
 741 }