lattice_point.cc

   1 #include <assert.h>
   2 #include <NTL/mat_ZZ.h>
   3 #include <NTL/vec_ZZ.h>
   4 #include <barvinok/polylib.h>
   5 #include <barvinok/options.h>
   6 #include <barvinok/evalue.h>
   7 #include <barvinok/util.h>
   8 #include "config.h"
   9 #include "conversion.h"
  10 #include "lattice_point.h"
  11 #include "param_util.h"
  12
  13 using std::cerr;
  14 using std::endl;
  15
  16 #define ALLOC(type) (type*)malloc(sizeof(type))
  17
  18 /* returns an evalue that corresponds to
  19  *
  20  * c/(*den) x_param
  21  */
  22 static evalue *term(int param, ZZ& c, Value *den = NULL)
  23 {
  24     evalue *EP = new evalue();
  25     value_init(EP->d);
  26     value_set_si(EP->d,0);
  27     EP->x.p = new_enode(polynomial, 2, param + 1);
  28     evalue_set_si(&EP->x.p->arr[0], 0, 1);
  29     value_init(EP->x.p->arr[1].x.n);
  30     if (den == NULL)
  31         value_set_si(EP->x.p->arr[1].d, 1);
  32     else
  33         value_assign(EP->x.p->arr[1].d, *den);
  34     zz2value(c, EP->x.p->arr[1].x.n);
  35     return EP;
  36 }
  37
  38 /* returns an evalue that corresponds to
  39  *
  40  *   sum_i p[i] * x_i
  41  */
  42 evalue *multi_monom(vec_ZZ& p)
  43 {
  44     evalue *X = ALLOC(evalue);
  45     value_init(X->d);
  46     value_init(X->x.n);
  47     unsigned nparam = p.length()-1;
  48     zz2value(p[nparam], X->x.n);
  49     value_set_si(X->d, 1);
  50     for (int i = 0; i < nparam; ++i) {
  51         if (p[i] == 0)
  52             continue;
  53         evalue *T = term(i, p[i]);
  54         eadd(T, X);
  55         free_evalue_refs(T);
  56         delete T;
  57     }
  58     return X;
  59 }
  60
  61 /*
  62  * Check whether mapping polyhedron P on the affine combination
  63  * num yields a range that has a fixed quotient on integer
  64  * division by d
  65  * If zero is true, then we are only interested in the quotient
  66  * for the cases where the remainder is zero.
  67  * Returns NULL if false and a newly allocated value if true.
  68  */
  69 static Value *fixed_quotient(Polyhedron *P, vec_ZZ& num, Value d, bool zero)
  70 {
  71     Value* ret = NULL;
  72     int len = num.length();
  73     Matrix *T = Matrix_Alloc(2, len);
  74     zz2values(num, T->p[0]);
  75     value_set_si(T->p[1][len-1], 1);
  76     Polyhedron *I = Polyhedron_Image(P, T, P->NbConstraints);
  77     Matrix_Free(T);
  78
  79     int i;
  80     for (i = 0; i < I->NbRays; ++i)
  81         if (value_zero_p(I->Ray[i][2])) {
  82             Polyhedron_Free(I);
  83             return NULL;
  84         }
  85
  86     Value min, max;
  87     value_init(min);
  88     value_init(max);
  89     int bounded = line_minmax(I, &min, &max);
  90     assert(bounded);
  91
  92     if (zero)
  93         mpz_cdiv_q(min, min, d);
  94     else
  95         mpz_fdiv_q(min, min, d);
  96     mpz_fdiv_q(max, max, d);
  97
  98     if (value_eq(min, max)) {
  99         ret = ALLOC(Value);
 100         value_init(*ret);
 101         value_assign(*ret, min);
 102     }
 103     value_clear(min);
 104     value_clear(max);
 105     return ret;
 106 }
 107
 108 /*
 109  * Normalize linear expression coef modulo m
 110  * Removes common factor and reduces coefficients
 111  * Returns index of first non-zero coefficient or len
 112  */
 113 int normal_mod(Value *coef, int len, Value *m)
 114 {
 115     Value gcd;
 116     value_init(gcd);
 117
 118     Vector_Gcd(coef, len, &gcd);
 119     value_gcd(gcd, gcd, *m);
 120     Vector_AntiScale(coef, coef, gcd, len);
 121
 122     value_division(*m, *m, gcd);
 123     value_clear(gcd);
 124
 125     if (value_one_p(*m))
 126         return len;
 127
 128     int j;
 129     for (j = 0; j < len; ++j)
 130         mpz_fdiv_r(coef[j], coef[j], *m);
 131     for (j = 0; j < len; ++j)
 132         if (value_notzero_p(coef[j]))
 133             break;
 134
 135     return j;
 136 }
 137
 138 static bool mod_needed(Polyhedron *PD, vec_ZZ& num, Value d, evalue *E)
 139 {
 140     Value *q = fixed_quotient(PD, num, d, false);
 141
 142     if (!q)
 143         return true;
 144
 145     value_oppose(*q, *q);
 146     evalue EV;
 147     value_init(EV.d);
 148     value_set_si(EV.d, 1);
 149     value_init(EV.x.n);
 150     value_multiply(EV.x.n, *q, d);
 151     eadd(&EV, E);
 152     free_evalue_refs(&EV);
 153     value_clear(*q);
 154     free(q);
 155     return false;
 156 }
 157
 158 /* Computes the fractional part of the affine expression specified
 159  * by coef (of length nvar+1) and the denominator denom.
 160  * If PD is not NULL, then it specifies additional constraints
 161  * on the variables that may be used to simplify the resulting
 162  * fractional part expression.
 163  *
 164  * Modifies coef argument !
 165  */
 166 evalue *fractional_part(Value *coef, Value denom, int nvar, Polyhedron *PD)
 167 {
 168     Value m;
 169     value_init(m);
 170     evalue *EP = evalue_zero();
 171     int sign = 1;
 172
 173     value_assign(m, denom);
 174     int j = normal_mod(coef, nvar+1, &m);
 175
 176     if (j == nvar+1) {
 177         value_clear(m);
 178         return EP;
 179     }
 180
 181     vec_ZZ num;
 182     values2zz(coef, num, nvar+1);
 183
 184     ZZ g;
 185     value2zz(m, g);
 186
 187     evalue tmp;
 188     value_init(tmp.d);
 189     evalue_set_si(&tmp, 0, 1);
 190
 191     int p = j;
 192     if (g % 2 == 0)
 193         while (j < nvar && (num[j] == g/2 || num[j] == 0))
 194             ++j;
 195     if ((j < nvar && num[j] > g/2) || (j == nvar && num[j] >= (g+1)/2)) {
 196         for (int k = j; k < nvar; ++k)
 197             if (num[k] != 0)
 198                 num[k] = g - num[k];
 199         num[nvar] = g - 1 - num[nvar];
 200         value_assign(tmp.d, m);
 201         ZZ t = sign*(g-1);
 202         zz2value(t, tmp.x.n);
 203         eadd(&tmp, EP);
 204         sign = -sign;
 205     }
 206
 207     if (p >= nvar) {
 208         ZZ t = num[nvar] * sign;
 209         zz2value(t, tmp.x.n);
 210         value_assign(tmp.d, m);
 211         eadd(&tmp, EP);
 212     } else {
 213         evalue *E = multi_monom(num);
 214         evalue EV;
 215         value_init(EV.d);
 216
 217         if (PD && !mod_needed(PD, num, m, E)) {
 218             value_init(EV.x.n);
 219             value_set_si(EV.x.n, sign);
 220             value_assign(EV.d, m);
 221             emul(&EV, E);
 222             eadd(E, EP);
 223         } else {
 224             value_init(EV.x.n);
 225             value_set_si(EV.x.n, 1);
 226             value_assign(EV.d, m);
 227             emul(&EV, E);
 228             value_clear(EV.x.n);
 229             value_set_si(EV.d, 0);
 230             EV.x.p = new_enode(fractional, 3, -1);
 231             evalue_copy(&EV.x.p->arr[0], E);
 232             evalue_set_si(&EV.x.p->arr[1], 0, 1);
 233             value_init(EV.x.p->arr[2].x.n);
 234             value_set_si(EV.x.p->arr[2].x.n, sign);
 235             value_set_si(EV.x.p->arr[2].d, 1);
 236
 237             eadd(&EV, EP);
 238         }
 239
 240         free_evalue_refs(&EV);
 241         evalue_free(E);
 242     }
 243
 244     free_evalue_refs(&tmp);
 245
 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 unimodular 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     Matrix *T = Matrix_Copy(basis);
 554     Matrix *inv = Matrix_Alloc(T->NbRows, T->NbColumns);
 555     int ok = Matrix_Inverse(T, inv);
 556     assert(ok);
 557     Matrix_Free(T);
 558
 559     Param_Vertex_Common_Denominator(V);
 560     /* temporarily ignore (common) denominator */
 561     V->Vertex->NbColumns--;
 562     Matrix *L = Matrix_Alloc(inv->NbRows, V->Vertex->NbColumns);
 563     Matrix_Product(inv, V->Vertex, L);
 564     V->Vertex->NbColumns++;
 565     Matrix_Free(inv);
 566
 567     return L;
 568 }
 569
 570 /* returns the unique lattice point in the fundamental parallelepiped
 571  * of the unimodular cone C shifted to the parametric vertex V.
 572  *
 573  * The return values num and E_vertex are such that
 574  * coordinate i of this lattice point is equal to
 575  *
 576  *          num[i] + E_vertex[i]
 577  */
 578 void lattice_point(Param_Vertices *V, const mat_ZZ& rays, vec_ZZ& num,
 579                    evalue **E_vertex, barvinok_options *options)
 580 {
 581     unsigned nparam = V->Vertex->NbColumns - 2;
 582     unsigned dim = rays.NumCols();
 583
 584     /* It doesn't really make sense to call lattice_point when dim == 0,
 585      * but apparently it happens from indicator_constructor in lexmin.
 586      */
 587     if (!dim)
 588         return;
 589
 590     vec_ZZ vertex;
 591     vertex.SetLength(nparam+1);
 592
 593     assert(V->Vertex->NbRows > 0);
 594     Param_Vertex_Common_Denominator(V);
 595
 596     if (value_notone_p(V->Vertex->p[0][nparam+1])) {
 597         Matrix* Rays = zz2matrix(rays);
 598         Matrix *T = Transpose(Rays);
 599         Matrix_Free(Rays);
 600         Matrix *L = relative_coordinates(V, T);
 601
 602         evalue f;
 603         value_init(f.d);
 604         value_init(f.x.n);
 605
 606         evalue **remainders = new evalue *[dim];
 607         for (int i = 0; i < dim; ++i)
 608             remainders[i] = ceil(L->p[i], nparam+1, V->Vertex->p[0][nparam+1],
 609                                  options);
 610         Matrix_Free(L);
 611
 612
 613         for (int i = 0; i < V->Vertex->NbRows; ++i) {
 614             values2zz(V->Vertex->p[i], vertex, nparam+1);
 615             E_vertex[i] = multi_monom(vertex);
 616             num[i] = 0;
 617
 618             value_set_si(f.x.n, 1);
 619             value_assign(f.d, V->Vertex->p[0][nparam+1]);
 620
 621             emul(&f, E_vertex[i]);
 622
 623             for (int j = 0; j < dim; ++j) {
 624                 if (value_zero_p(T->p[i][j]))
 625                     continue;
 626                 evalue cp;
 627                 value_init(cp.d);
 628                 evalue_copy(&cp, remainders[j]);
 629                 if (value_notone_p(T->p[i][j])) {
 630                     value_set_si(f.d, 1);
 631                     value_assign(f.x.n, T->p[i][j]);
 632                     emul(&f, &cp);
 633                 }
 634                 eadd(&cp, E_vertex[i]);
 635                 free_evalue_refs(&cp);
 636             }
 637         }
 638         for (int i = 0; i < dim; ++i)
 639             evalue_free(remainders[i]);
 640         delete [] remainders;
 641
 642         free_evalue_refs(&f);
 643
 644         Matrix_Free(T);
 645         return;
 646     }
 647
 648     for (int i = 0; i < V->Vertex->NbRows; ++i) {
 649         /* fixed value */
 650         if (First_Non_Zero(V->Vertex->p[i], nparam) == -1) {
 651             E_vertex[i] = 0;
 652             value2zz(V->Vertex->p[i][nparam], num[i]);
 653         } else {
 654             values2zz(V->Vertex->p[i], vertex, nparam+1);
 655             E_vertex[i] = multi_monom(vertex);
 656             num[i] = 0;
 657         }
 658     }
 659 }
 660
 661 static int lattice_point_fixed(Param_Vertices* V, const mat_ZZ& rays,
 662     vec_ZZ& lambda, term_info* term, unsigned long det)
 663 {
 664     unsigned nparam = V->Vertex->NbColumns - 2;
 665     unsigned dim = rays.NumCols();
 666
 667     for (int i = 0; i < dim; ++i)
 668         if (First_Non_Zero(V->Vertex->p[i], nparam) != -1)
 669             return 0;
 670
 671     Vector *fixed = Vector_Alloc(dim+1);
 672     for (int i = 0; i < dim; ++i)
 673         value_assign(fixed->p[i], V->Vertex->p[i][nparam]);
 674     value_assign(fixed->p[dim], V->Vertex->p[0][nparam+1]);
 675
 676     mat_ZZ vertex;
 677     Matrix *points = Matrix_Alloc(det, dim);
 678     Matrix* Rays = zz2matrix(rays);
 679     lattice_points_fixed(fixed->p, fixed->p, Rays, Rays, points, det);
 680     Matrix_Free(Rays);
 681     matrix2zz(points, vertex, points->NbRows, points->NbColumns);
 682     Matrix_Free(points);
 683     term->E = NULL;
 684     term->constant = vertex * lambda;
 685     Vector_Free(fixed);
 686
 687     return 1;
 688 }
 689
 690 /* Returns the power of (t+1) in the term of a rational generating function,
 691  * i.e., the scalar product of the actual lattice point and lambda.
 692  * The lattice point is the unique lattice point in the fundamental parallelepiped
 693  * of the unimodular cone i shifted to the parametric vertex V.
 694  *
 695  * The result is returned in term.
 696  */
 697 void lattice_point(Param_Vertices* V, const mat_ZZ& rays, vec_ZZ& lambda,
 698     term_info* term, unsigned long det,
 699     barvinok_options *options)
 700 {
 701     unsigned nparam = V->Vertex->NbColumns - 2;
 702     mat_ZZ vertex;
 703     vertex.SetDims(V->Vertex->NbRows, nparam+1);
 704
 705     Param_Vertex_Common_Denominator(V);
 706
 707     if (lattice_point_fixed(V, rays, lambda, term, det))
 708         return;
 709
 710     if (det != 1 || value_notone_p(V->Vertex->p[0][nparam+1])) {
 711         term->E = lattice_point(rays, lambda, V, det, options);
 712         return;
 713     }
 714     for (int i = 0; i < V->Vertex->NbRows; ++i) {
 715         assert(value_one_p(V->Vertex->p[i][nparam+1]));  // for now
 716         values2zz(V->Vertex->p[i], vertex[i], nparam+1);
 717     }
 718
 719     vec_ZZ num;
 720     num = lambda * vertex;
 721
 722     int nn = 0;
 723     for (int j = 0; j < nparam; ++j)
 724         if (num[j] != 0)
 725             ++nn;
 726     if (nn >= 1) {
 727         term->E = new evalue *[1];
 728         term->E[0] = multi_monom(num);
 729     } else {
 730         term->E = NULL;
 731         term->constant.SetLength(1);
 732         term->constant[0] = num[nparam];
 733     }
 734 }