lattice_point.cc

   1 #include <gmp.h>
   2 #include <NTL/mat_ZZ.h>
   3 #include <NTL/vec_ZZ.h>
   4 extern "C" {
   5 #include <polylib/polylibgmp.h>
   6 }
   7 #include <barvinok/barvinok.h>
   8 #include <barvinok/evalue.h>
   9 #include <barvinok/util.h>
  10 #include "config.h"
  11 #include "conversion.h"
  12 #include "lattice_point.h"
  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 #ifdef USE_MODULO
 244 static void ceil(Value *coef, int len, Value d, ZZ& f,
 245                  evalue *EP, Polyhedron *PD) {
 246     ceil_mod(coef, len, d, f, EP, PD);
 247 }
 248 #else
 249 static void ceil(Value *coef, int len, Value d, ZZ& f,
 250                  evalue *EP, Polyhedron *PD) {
 251     ceil_mod(coef, len, d, f, EP, PD);
 252     evalue_mod2table(EP, len-1);
 253 }
 254 #endif
 255
 256 evalue* bv_ceil3(Value *coef, int len, Value d, Polyhedron *P)
 257 {
 258     Vector *val = Vector_Alloc(len);
 259
 260     Value t;
 261     value_init(t);
 262     value_set_si(t, -1);
 263     Vector_Scale(coef, val->p, t, len);
 264     value_absolute(t, d);
 265
 266     vec_ZZ num;
 267     values2zz(val->p, num, len);
 268     evalue *EP = multi_monom(num);
 269
 270     evalue tmp;
 271     value_init(tmp.d);
 272     value_init(tmp.x.n);
 273     value_set_si(tmp.x.n, 1);
 274     value_assign(tmp.d, t);
 275
 276     emul(&tmp, EP);
 277
 278     ZZ one;
 279     one = 1;
 280     ceil_mod(val->p, len, t, one, EP, P);
 281     value_clear(t);
 282
 283     /* copy EP to malloc'ed evalue */
 284     evalue *E = ALLOC(evalue);
 285     *E = *EP;
 286     delete EP;
 287
 288     free_evalue_refs(&tmp);
 289     Vector_Free(val);
 290
 291     return E;
 292 }
 293
 294 void lattice_point(Value* values, Polyhedron *i, vec_ZZ& vertex)
 295 {
 296     unsigned dim = i->Dimension;
 297     if(!value_one_p(values[dim])) {
 298         Matrix* Rays = rays(i);
 299         Matrix *inv = Matrix_Alloc(Rays->NbRows, Rays->NbColumns);
 300         int ok = Matrix_Inverse(Rays, inv);
 301         assert(ok);
 302         Matrix_Free(Rays);
 303         Rays = rays(i);
 304         Vector *lambda = Vector_Alloc(dim+1);
 305         Vector_Matrix_Product(values, inv, lambda->p);
 306         Matrix_Free(inv);
 307         for (int j = 0; j < dim; ++j)
 308             mpz_cdiv_q(lambda->p[j], lambda->p[j], lambda->p[dim]);
 309         value_set_si(lambda->p[dim], 1);
 310         Vector *A = Vector_Alloc(dim+1);
 311         Vector_Matrix_Product(lambda->p, Rays, A->p);
 312         Vector_Free(lambda);
 313         Matrix_Free(Rays);
 314         values2zz(A->p, vertex, dim);
 315         Vector_Free(A);
 316     } else
 317         values2zz(values, vertex, dim);
 318 }
 319
 320 static void vertex_period(
 321                     Polyhedron *i, vec_ZZ& lambda, Matrix *T,
 322                     Value lcm, int p, Vector *val,
 323                     evalue *E, evalue* ev,
 324                     ZZ& offset)
 325 {
 326     unsigned nparam = T->NbRows - 1;
 327     unsigned dim = i->Dimension;
 328     Value tmp;
 329     ZZ nump;
 330
 331     if (p == nparam) {
 332         vec_ZZ vertex;
 333         ZZ num, l;
 334         Vector * values = Vector_Alloc(dim + 1);
 335         Vector_Matrix_Product(val->p, T, values->p);
 336         value_assign(values->p[dim], lcm);
 337         lattice_point(values->p, i, vertex);
 338         num = vertex * lambda;
 339         value2zz(lcm, l);
 340         num *= l;
 341         num += offset;
 342         value_init(ev->x.n);
 343         zz2value(num, ev->x.n);
 344         value_assign(ev->d, lcm);
 345         Vector_Free(values);
 346         return;
 347     }
 348
 349     value_init(tmp);
 350     vec_ZZ vertex;
 351     values2zz(T->p[p], vertex, dim);
 352     nump = vertex * lambda;
 353     if (First_Non_Zero(val->p, p) == -1) {
 354         value_assign(tmp, lcm);
 355         evalue *ET = term(p, nump, &tmp);
 356         eadd(ET, E);
 357         free_evalue_refs(ET);
 358         delete ET;
 359     }
 360
 361     value_assign(tmp, lcm);
 362     if (First_Non_Zero(T->p[p], dim) != -1)
 363         Vector_Gcd(T->p[p], dim, &tmp);
 364     Gcd(tmp, lcm, &tmp);
 365     if (value_lt(tmp, lcm)) {
 366         ZZ count;
 367
 368         value_division(tmp, lcm, tmp);
 369         value_set_si(ev->d, 0);
 370         ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
 371         value2zz(tmp, count);
 372         do {
 373             value_decrement(tmp, tmp);
 374             --count;
 375             ZZ new_offset = offset - count * nump;
 376             value_assign(val->p[p], tmp);
 377             vertex_period(i, lambda, T, lcm, p+1, val, E,
 378                           &ev->x.p->arr[VALUE_TO_INT(tmp)], new_offset);
 379         } while (value_pos_p(tmp));
 380     } else
 381         vertex_period(i, lambda, T, lcm, p+1, val, E, ev, offset);
 382     value_clear(tmp);
 383 }
 384
 385 /* Returns the power of (t+1) in the term of a rational generating function,
 386  * i.e., the scalar product of the actual lattice point and lambda.
 387  * The lattice point is the unique lattice point in the fundamental parallelepiped
 388  * of the unimodual cone i shifted to the parametric vertex W/lcm.
 389  *
 390  * The rows of W refer to the coordinates of the vertex
 391  * The first nparam columns are the coefficients of the parameters
 392  * and the final column is the constant term.
 393  * lcm is the common denominator of all coefficients.
 394  *
 395  * PD is the parameter domain, which, if != NULL, may be used to simply the
 396  * resulting expression.
 397  */
 398 #ifdef USE_MODULO
 399 evalue* lattice_point(
 400     Polyhedron *i, vec_ZZ& lambda, Matrix *W, Value lcm, Polyhedron *PD)
 401 {
 402     unsigned nparam = W->NbColumns - 1;
 403
 404     Matrix* Rays = rays2(i);
 405     Matrix *T = Transpose(Rays);
 406     Matrix *T2 = Matrix_Copy(T);
 407     Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
 408     int ok = Matrix_Inverse(T2, inv);
 409     assert(ok);
 410     Matrix_Free(Rays);
 411     Matrix_Free(T2);
 412     mat_ZZ vertex;
 413     matrix2zz(W, vertex, W->NbRows, W->NbColumns);
 414
 415     vec_ZZ num;
 416     num = lambda * vertex;
 417
 418     evalue *EP = multi_monom(num);
 419
 420     evalue tmp;
 421     value_init(tmp.d);
 422     value_init(tmp.x.n);
 423     value_set_si(tmp.x.n, 1);
 424     value_assign(tmp.d, lcm);
 425
 426     emul(&tmp, EP);
 427
 428     Matrix *L = Matrix_Alloc(inv->NbRows, W->NbColumns);
 429     Matrix_Product(inv, W, L);
 430
 431     mat_ZZ RT;
 432     matrix2zz(T, RT, T->NbRows, T->NbColumns);
 433     Matrix_Free(T);
 434
 435     vec_ZZ p = lambda * RT;
 436
 437     for (int i = 0; i < L->NbRows; ++i) {
 438         ceil_mod(L->p[i], nparam+1, lcm, p[i], EP, PD);
 439     }
 440
 441     Matrix_Free(L);
 442
 443     Matrix_Free(inv);
 444     free_evalue_refs(&tmp);
 445     return EP;
 446 }
 447 #else
 448 evalue* lattice_point(
 449     Polyhedron *i, vec_ZZ& lambda, Matrix *W, Value lcm, Polyhedron *PD)
 450 {
 451     Matrix *T = Transpose(W);
 452     unsigned nparam = T->NbRows - 1;
 453
 454     evalue *EP = new evalue();
 455     value_init(EP->d);
 456     evalue_set_si(EP, 0, 1);
 457
 458     evalue ev;
 459     Vector *val = Vector_Alloc(nparam+1);
 460     value_set_si(val->p[nparam], 1);
 461     ZZ offset(INIT_VAL, 0);
 462     value_init(ev.d);
 463     vertex_period(i, lambda, T, lcm, 0, val, EP, &ev, offset);
 464     Vector_Free(val);
 465     eadd(&ev, EP);
 466     free_evalue_refs(&ev);
 467
 468     Matrix_Free(T);
 469
 470     reduce_evalue(EP);
 471
 472     return EP;
 473 }
 474 #endif
 475
 476 /* returns the unique lattice point in the fundamental parallelepiped
 477  * of the unimodual cone C shifted to the parametric vertex V.
 478  *
 479  * The return values num and E_vertex are such that
 480  * coordinate i of this lattice point is equal to
 481  *
 482  *          num[i] + E_vertex[i]
 483  */
 484 void lattice_point(Param_Vertices *V, Polyhedron *C, vec_ZZ& num,
 485                    evalue **E_vertex)
 486 {
 487     unsigned nparam = V->Vertex->NbColumns - 2;
 488     unsigned dim = C->Dimension;
 489     vec_ZZ vertex;
 490     vertex.SetLength(nparam+1);
 491
 492     Value lcm, tmp;
 493     value_init(lcm);
 494     value_init(tmp);
 495     value_set_si(lcm, 1);
 496
 497     for (int j = 0; j < V->Vertex->NbRows; ++j) {
 498         value_lcm(lcm, V->Vertex->p[j][nparam+1], &lcm);
 499     }
 500
 501     if (value_notone_p(lcm)) {
 502         Matrix * mv = Matrix_Alloc(dim, nparam+1);
 503         for (int j = 0 ; j < dim; ++j) {
 504             value_division(tmp, lcm, V->Vertex->p[j][nparam+1]);
 505             Vector_Scale(V->Vertex->p[j], mv->p[j], tmp, nparam+1);
 506         }
 507
 508         Matrix* Rays = rays2(C);
 509         Matrix *T = Transpose(Rays);
 510         Matrix *T2 = Matrix_Copy(T);
 511         Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
 512         int ok = Matrix_Inverse(T2, inv);
 513         assert(ok);
 514         Matrix_Free(Rays);
 515         Matrix_Free(T2);
 516         Matrix *L = Matrix_Alloc(inv->NbRows, mv->NbColumns);
 517         Matrix_Product(inv, mv, L);
 518         Matrix_Free(inv);
 519
 520         evalue f;
 521         value_init(f.d);
 522         value_init(f.x.n);
 523
 524         ZZ one;
 525
 526         evalue *remainders[dim];
 527         for (int i = 0; i < dim; ++i) {
 528             remainders[i] = evalue_zero();
 529             one = 1;
 530             ceil(L->p[i], nparam+1, lcm, one, remainders[i], 0);
 531         }
 532         Matrix_Free(L);
 533
 534
 535         for (int i = 0; i < V->Vertex->NbRows; ++i) {
 536             values2zz(mv->p[i], vertex, nparam+1);
 537             E_vertex[i] = multi_monom(vertex);
 538             num[i] = 0;
 539
 540             value_set_si(f.x.n, 1);
 541             value_assign(f.d, lcm);
 542
 543             emul(&f, E_vertex[i]);
 544
 545             for (int j = 0; j < dim; ++j) {
 546                 if (value_zero_p(T->p[i][j]))
 547                     continue;
 548                 evalue cp;
 549                 value_init(cp.d);
 550                 evalue_copy(&cp, remainders[j]);
 551                 if (value_notone_p(T->p[i][j])) {
 552                     value_set_si(f.d, 1);
 553                     value_assign(f.x.n, T->p[i][j]);
 554                     emul(&f, &cp);
 555                 }
 556                 eadd(&cp, E_vertex[i]);
 557                 free_evalue_refs(&cp);
 558             }
 559         }
 560         for (int i = 0; i < dim; ++i) {
 561             free_evalue_refs(remainders[i]);
 562             free(remainders[i]);
 563         }
 564
 565         free_evalue_refs(&f);
 566
 567         Matrix_Free(T);
 568         Matrix_Free(mv);
 569         value_clear(lcm);
 570         value_clear(tmp);
 571         return;
 572     }
 573     value_clear(lcm);
 574     value_clear(tmp);
 575
 576     for (int i = 0; i < V->Vertex->NbRows; ++i) {
 577         /* fixed value */
 578         if (First_Non_Zero(V->Vertex->p[i], nparam) == -1) {
 579             E_vertex[i] = 0;
 580             value2zz(V->Vertex->p[i][nparam], num[i]);
 581         } else {
 582             values2zz(V->Vertex->p[i], vertex, nparam+1);
 583             E_vertex[i] = multi_monom(vertex);
 584             num[i] = 0;
 585         }
 586     }
 587 }
 588