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
  11 using std::cerr;
  12 using std::endl;
  13
  14 #define ALLOC(type) (type*)malloc(sizeof(type))
  15
  16 /* returns an evalue that corresponds to
  17  *
  18  * c/(*den) x_param
  19  */
  20 static evalue *term(int param, ZZ& c, Value *den = NULL)
  21 {
  22     evalue *EP = new evalue();
  23     value_init(EP->d);
  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);
  28     if (den == NULL)
  29         value_set_si(EP->x.p->arr[1].d, 1);
  30     else
  31         value_assign(EP->x.p->arr[1].d, *den);
  32     zz2value(c, EP->x.p->arr[1].x.n);
  33     return EP;
  34 }
  35
  36 /* returns an evalue that corresponds to
  37  *
  38  *   sum_i p[i] * x_i
  39  */
  40 evalue *multi_monom(vec_ZZ& p)
  41 {
  42     evalue *X = new evalue();
  43     value_init(X->d);
  44     value_init(X->x.n);
  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) {
  49         if (p[i] == 0)
  50             continue;
  51         evalue *T = term(i, p[i]);
  52         eadd(T, X);
  53         free_evalue_refs(T);
  54         delete T;
  55     }
  56     return X;
  57 }
  58
  59 /*
  60  * Check whether mapping polyhedron P on the affine combination
  61  * num yields a range that has a fixed quotient on integer
  62  * division by d
  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.
  66  */
  67 static Value *fixed_quotient(Polyhedron *P, vec_ZZ& num, Value d, bool zero)
  68 {
  69     Value* ret = NULL;
  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);
  75     Matrix_Free(T);
  76
  77     int i;
  78     for (i = 0; i < I->NbRays; ++i)
  79         if (value_zero_p(I->Ray[i][2])) {
  80             Polyhedron_Free(I);
  81             return NULL;
  82         }
  83
  84     Value min, max;
  85     value_init(min);
  86     value_init(max);
  87     int bounded = line_minmax(I, &min, &max);
  88     assert(bounded);
  89
  90     if (zero)
  91         mpz_cdiv_q(min, min, d);
  92     else
  93         mpz_fdiv_q(min, min, d);
  94     mpz_fdiv_q(max, max, d);
  95
  96     if (value_eq(min, max)) {
  97         ret = ALLOC(Value);
  98         value_init(*ret);
  99         value_assign(*ret, min);
 100     }
 101     value_clear(min);
 102     value_clear(max);
 103     return ret;
 104 }
 105
 106 /*
 107  * Normalize linear expression coef modulo m
 108  * Removes common factor and reduces coefficients
 109  * Returns index of first non-zero coefficient or len
 110  */
 111 int normal_mod(Value *coef, int len, Value *m)
 112 {
 113     Value gcd;
 114     value_init(gcd);
 115
 116     Vector_Gcd(coef, len, &gcd);
 117     Gcd(gcd, *m, &gcd);
 118     Vector_AntiScale(coef, coef, gcd, len);
 119
 120     value_division(*m, *m, gcd);
 121     value_clear(gcd);
 122
 123     if (value_one_p(*m))
 124         return len;
 125
 126     int j;
 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]))
 131             break;
 132
 133     return j;
 134 }
 135
 136 static bool mod_needed(Polyhedron *PD, vec_ZZ& num, Value d, evalue *E)
 137 {
 138     Value *q = fixed_quotient(PD, num, d, false);
 139
 140     if (!q)
 141         return true;
 142
 143     value_oppose(*q, *q);
 144     evalue EV;
 145     value_init(EV.d);
 146     value_set_si(EV.d, 1);
 147     value_init(EV.x.n);
 148     value_multiply(EV.x.n, *q, d);
 149     eadd(&EV, E);
 150     free_evalue_refs(&EV);
 151     value_clear(*q);
 152     free(q);
 153     return false;
 154 }
 155
 156 /* modifies f argument ! */
 157 static void ceil_mod(Value *coef, int len, Value d, ZZ& f, evalue *EP, Polyhedron *PD)
 158 {
 159     Value m;
 160     value_init(m);
 161     value_set_si(m, -1);
 162
 163     Vector_Scale(coef, coef, m, len);
 164
 165     value_assign(m, d);
 166     int j = normal_mod(coef, len, &m);
 167
 168     if (j == len) {
 169         value_clear(m);
 170         return;
 171     }
 172
 173     vec_ZZ num;
 174     values2zz(coef, num, len);
 175
 176     ZZ g;
 177     value2zz(m, g);
 178
 179     evalue tmp;
 180     value_init(tmp.d);
 181     evalue_set_si(&tmp, 0, 1);
 182
 183     int p = j;
 184     if (g % 2 == 0)
 185         while (j < len-1 && (num[j] == g/2 || num[j] == 0))
 186             ++j;
 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)
 189             if (num[k] != 0)
 190                 num[k] = g - num[k];
 191         num[len-1] = g - 1 - num[len-1];
 192         value_assign(tmp.d, m);
 193         ZZ t = f*(g-1);
 194         zz2value(t, tmp.x.n);
 195         eadd(&tmp, EP);
 196         f = -f;
 197     }
 198
 199     if (p >= len-1) {
 200         ZZ t = num[len-1] * f;
 201         zz2value(t, tmp.x.n);
 202         value_assign(tmp.d, m);
 203         eadd(&tmp, EP);
 204     } else {
 205         evalue *E = multi_monom(num);
 206         evalue EV;
 207         value_init(EV.d);
 208
 209         if (PD && !mod_needed(PD, num, m, E)) {
 210             value_init(EV.x.n);
 211             zz2value(f, EV.x.n);
 212             value_assign(EV.d, m);
 213             emul(&EV, E);
 214             eadd(E, EP);
 215         } else {
 216             value_init(EV.x.n);
 217             value_set_si(EV.x.n, 1);
 218             value_assign(EV.d, m);
 219             emul(&EV, E);
 220             value_clear(EV.x.n);
 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);
 228
 229             eadd(&EV, EP);
 230         }
 231
 232         free_evalue_refs(&EV);
 233         free_evalue_refs(E);
 234         delete E;
 235     }
 236
 237     free_evalue_refs(&tmp);
 238
 239 out:
 240     value_clear(m);
 241 }
 242
 243 static void ceil(Value *coef, int len, Value d, ZZ& f,
 244                  evalue *EP, Polyhedron *PD, barvinok_options *options)
 245 {
 246     ceil_mod(coef, len, d, f, EP, PD);
 247     if (options->lookup_table)
 248         evalue_mod2table(EP, len-1);
 249 }
 250
 251 evalue* bv_ceil3(Value *coef, int len, Value d, Polyhedron *P)
 252 {
 253     Vector *val = Vector_Alloc(len);
 254
 255     Value t;
 256     value_init(t);
 257     value_set_si(t, -1);
 258     Vector_Scale(coef, val->p, t, len);
 259     value_absolute(t, d);
 260
 261     vec_ZZ num;
 262     values2zz(val->p, num, len);
 263     evalue *EP = multi_monom(num);
 264
 265     evalue tmp;
 266     value_init(tmp.d);
 267     value_init(tmp.x.n);
 268     value_set_si(tmp.x.n, 1);
 269     value_assign(tmp.d, t);
 270
 271     emul(&tmp, EP);
 272
 273     ZZ one;
 274     one = 1;
 275     ceil_mod(val->p, len, t, one, EP, P);
 276     value_clear(t);
 277
 278     /* copy EP to malloc'ed evalue */
 279     evalue *E = ALLOC(evalue);
 280     *E = *EP;
 281     delete EP;
 282
 283     free_evalue_refs(&tmp);
 284     Vector_Free(val);
 285
 286     return E;
 287 }
 288
 289 void lattice_point(Value* values, const mat_ZZ& rays, vec_ZZ& vertex, int *closed)
 290 {
 291     unsigned dim = rays.NumRows();
 292     if (value_one_p(values[dim]) && !closed)
 293         values2zz(values, vertex, dim);
 294     else {
 295         Matrix* Rays = rays2matrix(rays);
 296         Matrix *inv = Matrix_Alloc(Rays->NbRows, Rays->NbColumns);
 297         int ok = Matrix_Inverse(Rays, inv);
 298         assert(ok);
 299         Matrix_Free(Rays);
 300         Rays = rays2matrix(rays);
 301         Vector *lambda = Vector_Alloc(dim+1);
 302         Vector_Matrix_Product(values, inv, lambda->p);
 303         Matrix_Free(inv);
 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]);
 307             else {
 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]);
 310             }
 311         value_set_si(lambda->p[dim], 1);
 312         Vector *A = Vector_Alloc(dim+1);
 313         Vector_Matrix_Product(lambda->p, Rays, A->p);
 314         Vector_Free(lambda);
 315         Matrix_Free(Rays);
 316         values2zz(A->p, vertex, dim);
 317         Vector_Free(A);
 318     }
 319 }
 320
 321 /* Compute the lattice points in the vertex cone at "values" with rays "rays".
 322  * The lattice points are returned in "vertex".
 323  *
 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
 326  * the cosets.
 327  * We compute
 328  * [k 1] [ d1*W  0 ] [  U'  0 ] = [k 1] T2
 329  *       [ -v   d1 ] [  0  d2 ]
 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.
 341  */
 342 void lattice_point(Value* values, const mat_ZZ& rays, mat_ZZ& vertex,
 343                    unsigned long det, int *closed)
 344 {
 345     unsigned dim = rays.NumRows();
 346     vertex.SetDims(det, dim);
 347     if (det == 1) {
 348         lattice_point(values, rays, vertex[0], closed);
 349         return;
 350     }
 351     Matrix* Rays = rays2matrix2(rays);
 352     Matrix *U, *W, *D;
 353     Smith(Rays, &U, &W, &D);
 354     Matrix_Free(Rays);
 355     Matrix_Free(U);
 356
 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);
 360     Matrix_Free(W);
 361     Value tmp;
 362     value_init(tmp);
 363     value_set_si(tmp, -1);
 364     Vector_Scale(values, T->p[dim], tmp, dim);
 365     value_clear(tmp);
 366     value_assign(T->p[dim][dim], values[dim]);
 367
 368     Rays = rays2matrix(rays);
 369     Matrix *inv = Matrix_Alloc(Rays->NbRows, Rays->NbColumns);
 370     int ok = Matrix_Inverse(Rays, inv);
 371     assert(ok);
 372     Matrix_Free(Rays);
 373
 374     Matrix *T2 = Matrix_Alloc(dim+1, dim+1);
 375     Matrix_Product(T, inv, T2);
 376     Matrix_Free(T);
 377
 378     Rays = rays2matrix(rays);
 379
 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]);
 389         }
 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]);
 394             else {
 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]);
 397             }
 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);
 407     }
 408     Vector_Free(k);
 409     Vector_Free(lambda);
 410     Vector_Free(lambda2);
 411     Matrix_Free(D);
 412     Matrix_Free(Rays);
 413     Matrix_Free(inv);
 414
 415     Matrix_Free(T2);
 416 }
 417
 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,
 422                     ZZ& offset)
 423 {
 424     unsigned nparam = T->NbRows - 1;
 425     unsigned dim = rays.NumRows();
 426     Value tmp;
 427     ZZ nump;
 428
 429     if (p == nparam) {
 430         vec_ZZ vertex;
 431         ZZ num, l;
 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;
 437         value2zz(lcm, l);
 438         num *= l;
 439         num += offset;
 440         value_init(ev->x.n);
 441         zz2value(num, ev->x.n);
 442         value_assign(ev->d, lcm);
 443         Vector_Free(values);
 444         return;
 445     }
 446
 447     value_init(tmp);
 448     vec_ZZ vertex;
 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);
 454         eadd(ET, E);
 455         free_evalue_refs(ET);
 456         delete ET;
 457     }
 458
 459     value_assign(tmp, lcm);
 460     if (First_Non_Zero(T->p[p], dim) != -1)
 461         Vector_Gcd(T->p[p], dim, &tmp);
 462     Gcd(tmp, lcm, &tmp);
 463     if (value_lt(tmp, lcm)) {
 464         ZZ count;
 465
 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);
 470         do {
 471             value_decrement(tmp, tmp);
 472             --count;
 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));
 478     } else
 479         vertex_period(rays, lambda, T, lcm, p+1, val, E, ev, offset);
 480     value_clear(tmp);
 481 }
 482
 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.
 487  *
 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.
 492  *
 493  * PD is the parameter domain, which, if != NULL, may be used to simply the
 494  * resulting expression.
 495  */
 496 static evalue* lattice_point_fractional(const mat_ZZ& rays, vec_ZZ& lambda,
 497                                         Matrix *W, Value lcm, Polyhedron *PD)
 498 {
 499     unsigned nparam = W->NbColumns - 1;
 500
 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);
 506     assert(ok);
 507     Matrix_Free(Rays);
 508     Matrix_Free(T2);
 509     mat_ZZ vertex;
 510     matrix2zz(W, vertex, W->NbRows, W->NbColumns);
 511
 512     vec_ZZ num;
 513     num = lambda * vertex;
 514
 515     evalue *EP = multi_monom(num);
 516
 517     evalue tmp;
 518     value_init(tmp.d);
 519     value_init(tmp.x.n);
 520     value_set_si(tmp.x.n, 1);
 521     value_assign(tmp.d, lcm);
 522
 523     emul(&tmp, EP);
 524
 525     Matrix *L = Matrix_Alloc(inv->NbRows, W->NbColumns);
 526     Matrix_Product(inv, W, L);
 527
 528     mat_ZZ RT;
 529     matrix2zz(T, RT, T->NbRows, T->NbColumns);
 530     Matrix_Free(T);
 531
 532     vec_ZZ p = lambda * RT;
 533
 534     for (int i = 0; i < L->NbRows; ++i) {
 535         ceil_mod(L->p[i], nparam+1, lcm, p[i], EP, PD);
 536     }
 537
 538     Matrix_Free(L);
 539
 540     Matrix_Free(inv);
 541     free_evalue_refs(&tmp);
 542     return EP;
 543 }
 544
 545 static evalue* lattice_point_table(const mat_ZZ& rays, vec_ZZ& lambda, Matrix *W,
 546                                    Value lcm, Polyhedron *PD)
 547 {
 548     Matrix *T = Transpose(W);
 549     unsigned nparam = T->NbRows - 1;
 550
 551     evalue *EP = new evalue();
 552     value_init(EP->d);
 553     evalue_set_si(EP, 0, 1);
 554
 555     evalue ev;
 556     Vector *val = Vector_Alloc(nparam+1);
 557     value_set_si(val->p[nparam], 1);
 558     ZZ offset(INIT_VAL, 0);
 559     value_init(ev.d);
 560     vertex_period(rays, lambda, T, lcm, 0, val, EP, &ev, offset);
 561     Vector_Free(val);
 562     eadd(&ev, EP);
 563     free_evalue_refs(&ev);
 564
 565     Matrix_Free(T);
 566
 567     reduce_evalue(EP);
 568
 569     return EP;
 570 }
 571
 572 evalue* lattice_point(const mat_ZZ& rays, vec_ZZ& lambda, Matrix *W,
 573                       Value lcm, Polyhedron *PD, barvinok_options *options)
 574 {
 575     if (options->lookup_table)
 576         return lattice_point_table(rays, lambda, W, lcm, PD);
 577     else
 578         return lattice_point_fractional(rays, lambda, W, lcm, PD);
 579 }
 580
 581 /* returns the unique lattice point in the fundamental parallelepiped
 582  * of the unimodual cone C shifted to the parametric vertex V.
 583  *
 584  * The return values num and E_vertex are such that
 585  * coordinate i of this lattice point is equal to
 586  *
 587  *          num[i] + E_vertex[i]
 588  */
 589 void lattice_point(Param_Vertices *V, const mat_ZZ& rays, vec_ZZ& num,
 590                    evalue **E_vertex, barvinok_options *options)
 591 {
 592     unsigned nparam = V->Vertex->NbColumns - 2;
 593     unsigned dim = rays.NumCols();
 594     vec_ZZ vertex;
 595     vertex.SetLength(nparam+1);
 596
 597     Value lcm, tmp;
 598     value_init(lcm);
 599     value_init(tmp);
 600     value_set_si(lcm, 1);
 601
 602     for (int j = 0; j < V->Vertex->NbRows; ++j) {
 603         value_lcm(lcm, V->Vertex->p[j][nparam+1], &lcm);
 604     }
 605
 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);
 611         }
 612
 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);
 618         assert(ok);
 619         Matrix_Free(Rays);
 620         Matrix_Free(T2);
 621         Matrix *L = Matrix_Alloc(inv->NbRows, mv->NbColumns);
 622         Matrix_Product(inv, mv, L);
 623         Matrix_Free(inv);
 624
 625         evalue f;
 626         value_init(f.d);
 627         value_init(f.x.n);
 628
 629         ZZ one;
 630
 631         evalue *remainders[dim];
 632         for (int i = 0; i < dim; ++i) {
 633             remainders[i] = evalue_zero();
 634             one = 1;
 635             ceil(L->p[i], nparam+1, lcm, one, remainders[i], 0, options);
 636         }
 637         Matrix_Free(L);
 638
 639
 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);
 643             num[i] = 0;
 644
 645             value_set_si(f.x.n, 1);
 646             value_assign(f.d, lcm);
 647
 648             emul(&f, E_vertex[i]);
 649
 650             for (int j = 0; j < dim; ++j) {
 651                 if (value_zero_p(T->p[i][j]))
 652                     continue;
 653                 evalue cp;
 654                 value_init(cp.d);
 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]);
 659                     emul(&f, &cp);
 660                 }
 661                 eadd(&cp, E_vertex[i]);
 662                 free_evalue_refs(&cp);
 663             }
 664         }
 665         for (int i = 0; i < dim; ++i) {
 666             free_evalue_refs(remainders[i]);
 667             free(remainders[i]);
 668         }
 669
 670         free_evalue_refs(&f);
 671
 672         Matrix_Free(T);
 673         Matrix_Free(mv);
 674         value_clear(lcm);
 675         value_clear(tmp);
 676         return;
 677     }
 678     value_clear(lcm);
 679     value_clear(tmp);
 680
 681     for (int i = 0; i < V->Vertex->NbRows; ++i) {
 682         /* fixed value */
 683         if (First_Non_Zero(V->Vertex->p[i], nparam) == -1) {
 684             E_vertex[i] = 0;
 685             value2zz(V->Vertex->p[i][nparam], num[i]);
 686         } else {
 687             values2zz(V->Vertex->p[i], vertex, nparam+1);
 688             E_vertex[i] = multi_monom(vertex);
 689             num[i] = 0;
 690         }
 691     }
 692 }
 693