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 #define ALLOC(type) (type*)malloc(sizeof(type))
  12
  13 /* returns an evalue that corresponds to
  14  *
  15  * c/(*den) x_param
  16  */
  17 static evalue *term(int param, ZZ& c, Value *den = NULL)
  18 {
  19     evalue *EP = new evalue();
  20     value_init(EP->d);
  21     value_set_si(EP->d,0);
  22     EP->x.p = new_enode(polynomial, 2, param + 1);
  23     evalue_set_si(&EP->x.p->arr[0], 0, 1);
  24     value_init(EP->x.p->arr[1].x.n);
  25     if (den == NULL)
  26         value_set_si(EP->x.p->arr[1].d, 1);
  27     else
  28         value_assign(EP->x.p->arr[1].d, *den);
  29     zz2value(c, EP->x.p->arr[1].x.n);
  30     return EP;
  31 }
  32
  33 /* returns an evalue that corresponds to
  34  *
  35  *   sum_i p[i] * x_i
  36  */
  37 evalue *multi_monom(vec_ZZ& p)
  38 {
  39     evalue *X = new evalue();
  40     value_init(X->d);
  41     value_init(X->x.n);
  42     unsigned nparam = p.length()-1;
  43     zz2value(p[nparam], X->x.n);
  44     value_set_si(X->d, 1);
  45     for (int i = 0; i < nparam; ++i) {
  46         if (p[i] == 0)
  47             continue;
  48         evalue *T = term(i, p[i]);
  49         eadd(T, X);
  50         free_evalue_refs(T);
  51         delete T;
  52     }
  53     return X;
  54 }
  55
  56 /*
  57  * Check whether mapping polyhedron P on the affine combination
  58  * num yields a range that has a fixed quotient on integer
  59  * division by d
  60  * If zero is true, then we are only interested in the quotient
  61  * for the cases where the remainder is zero.
  62  * Returns NULL if false and a newly allocated value if true.
  63  */
  64 static Value *fixed_quotient(Polyhedron *P, vec_ZZ& num, Value d, bool zero)
  65 {
  66     Value* ret = NULL;
  67     int len = num.length();
  68     Matrix *T = Matrix_Alloc(2, len);
  69     zz2values(num, T->p[0]);
  70     value_set_si(T->p[1][len-1], 1);
  71     Polyhedron *I = Polyhedron_Image(P, T, P->NbConstraints);
  72     Matrix_Free(T);
  73
  74     int i;
  75     for (i = 0; i < I->NbRays; ++i)
  76         if (value_zero_p(I->Ray[i][2])) {
  77             Polyhedron_Free(I);
  78             return NULL;
  79         }
  80
  81     Value min, max;
  82     value_init(min);
  83     value_init(max);
  84     int bounded = line_minmax(I, &min, &max);
  85     assert(bounded);
  86
  87     if (zero)
  88         mpz_cdiv_q(min, min, d);
  89     else
  90         mpz_fdiv_q(min, min, d);
  91     mpz_fdiv_q(max, max, d);
  92
  93     if (value_eq(min, max)) {
  94         ret = ALLOC(Value);
  95         value_init(*ret);
  96         value_assign(*ret, min);
  97     }
  98     value_clear(min);
  99     value_clear(max);
 100     return ret;
 101 }
 102
 103 /*
 104  * Normalize linear expression coef modulo m
 105  * Removes common factor and reduces coefficients
 106  * Returns index of first non-zero coefficient or len
 107  */
 108 int normal_mod(Value *coef, int len, Value *m)
 109 {
 110     Value gcd;
 111     value_init(gcd);
 112
 113     Vector_Gcd(coef, len, &gcd);
 114     Gcd(gcd, *m, &gcd);
 115     Vector_AntiScale(coef, coef, gcd, len);
 116
 117     value_division(*m, *m, gcd);
 118     value_clear(gcd);
 119
 120     if (value_one_p(*m))
 121         return len;
 122
 123     int j;
 124     for (j = 0; j < len; ++j)
 125         mpz_fdiv_r(coef[j], coef[j], *m);
 126     for (j = 0; j < len; ++j)
 127         if (value_notzero_p(coef[j]))
 128             break;
 129
 130     return j;
 131 }
 132
 133 static bool mod_needed(Polyhedron *PD, vec_ZZ& num, Value d, evalue *E)
 134 {
 135     Value *q = fixed_quotient(PD, num, d, false);
 136
 137     if (!q)
 138         return true;
 139
 140     value_oppose(*q, *q);
 141     evalue EV;
 142     value_init(EV.d);
 143     value_set_si(EV.d, 1);
 144     value_init(EV.x.n);
 145     value_multiply(EV.x.n, *q, d);
 146     eadd(&EV, E);
 147     free_evalue_refs(&EV);
 148     value_clear(*q);
 149     free(q);
 150     return false;
 151 }
 152
 153 /* modifies f argument ! */
 154 static void ceil_mod(Value *coef, int len, Value d, ZZ& f, evalue *EP, Polyhedron *PD)
 155 {
 156     Value m;
 157     value_init(m);
 158     value_set_si(m, -1);
 159
 160     Vector_Scale(coef, coef, m, len);
 161
 162     value_assign(m, d);
 163     int j = normal_mod(coef, len, &m);
 164
 165     if (j == len) {
 166         value_clear(m);
 167         return;
 168     }
 169
 170     vec_ZZ num;
 171     values2zz(coef, num, len);
 172
 173     ZZ g;
 174     value2zz(m, g);
 175
 176     evalue tmp;
 177     value_init(tmp.d);
 178     evalue_set_si(&tmp, 0, 1);
 179
 180     int p = j;
 181     if (g % 2 == 0)
 182         while (j < len-1 && (num[j] == g/2 || num[j] == 0))
 183             ++j;
 184     if ((j < len-1 && num[j] > g/2) || (j == len-1 && num[j] >= (g+1)/2)) {
 185         for (int k = j; k < len-1; ++k)
 186             if (num[k] != 0)
 187                 num[k] = g - num[k];
 188         num[len-1] = g - 1 - num[len-1];
 189         value_assign(tmp.d, m);
 190         ZZ t = f*(g-1);
 191         zz2value(t, tmp.x.n);
 192         eadd(&tmp, EP);
 193         f = -f;
 194     }
 195
 196     if (p >= len-1) {
 197         ZZ t = num[len-1] * f;
 198         zz2value(t, tmp.x.n);
 199         value_assign(tmp.d, m);
 200         eadd(&tmp, EP);
 201     } else {
 202         evalue *E = multi_monom(num);
 203         evalue EV;
 204         value_init(EV.d);
 205
 206         if (PD && !mod_needed(PD, num, m, E)) {
 207             value_init(EV.x.n);
 208             zz2value(f, EV.x.n);
 209             value_assign(EV.d, m);
 210             emul(&EV, E);
 211             eadd(E, EP);
 212         } else {
 213             value_init(EV.x.n);
 214             value_set_si(EV.x.n, 1);
 215             value_assign(EV.d, m);
 216             emul(&EV, E);
 217             value_clear(EV.x.n);
 218             value_set_si(EV.d, 0);
 219             EV.x.p = new_enode(fractional, 3, -1);
 220             evalue_copy(&EV.x.p->arr[0], E);
 221             evalue_set_si(&EV.x.p->arr[1], 0, 1);
 222             value_init(EV.x.p->arr[2].x.n);
 223             zz2value(f, EV.x.p->arr[2].x.n);
 224             value_set_si(EV.x.p->arr[2].d, 1);
 225
 226             eadd(&EV, EP);
 227         }
 228
 229         free_evalue_refs(&EV);
 230         free_evalue_refs(E);
 231         delete E;
 232     }
 233
 234     free_evalue_refs(&tmp);
 235
 236 out:
 237     value_clear(m);
 238 }
 239
 240 #ifdef USE_MODULO
 241 static void ceil(Value *coef, int len, Value d, ZZ& f,
 242                  evalue *EP, Polyhedron *PD) {
 243     ceil_mod(coef, len, d, f, EP, PD);
 244 }
 245 #else
 246 static void ceil(Value *coef, int len, Value d, ZZ& f,
 247                  evalue *EP, Polyhedron *PD) {
 248     ceil_mod(coef, len, d, f, EP, PD);
 249     evalue_mod2table(EP, len-1);
 250 }
 251 #endif
 252
 253 evalue* bv_ceil3(Value *coef, int len, Value d, Polyhedron *P)
 254 {
 255     Vector *val = Vector_Alloc(len);
 256
 257     Value t;
 258     value_init(t);
 259     value_set_si(t, -1);
 260     Vector_Scale(coef, val->p, t, len);
 261     value_absolute(t, d);
 262
 263     vec_ZZ num;
 264     values2zz(val->p, num, len);
 265     evalue *EP = multi_monom(num);
 266
 267     evalue tmp;
 268     value_init(tmp.d);
 269     value_init(tmp.x.n);
 270     value_set_si(tmp.x.n, 1);
 271     value_assign(tmp.d, t);
 272
 273     emul(&tmp, EP);
 274
 275     ZZ one;
 276     one = 1;
 277     ceil_mod(val->p, len, t, one, EP, P);
 278     value_clear(t);
 279
 280     /* copy EP to malloc'ed evalue */
 281     evalue *E = ALLOC(evalue);
 282     *E = *EP;
 283     delete EP;
 284
 285     free_evalue_refs(&tmp);
 286     Vector_Free(val);
 287
 288     return E;
 289 }
 290
 291 void lattice_point(Value* values, Polyhedron *i, vec_ZZ& vertex)
 292 {
 293     unsigned dim = i->Dimension;
 294     if(!value_one_p(values[dim])) {
 295         Matrix* Rays = rays(i);
 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 = rays(i);
 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             mpz_cdiv_q(lambda->p[j], lambda->p[j], lambda->p[dim]);
 306         value_set_si(lambda->p[dim], 1);
 307         Vector *A = Vector_Alloc(dim+1);
 308         Vector_Matrix_Product(lambda->p, Rays, A->p);
 309         Vector_Free(lambda);
 310         Matrix_Free(Rays);
 311         values2zz(A->p, vertex, dim);
 312         Vector_Free(A);
 313     } else
 314         values2zz(values, vertex, dim);
 315 }
 316
 317 static void vertex_period(
 318                     Polyhedron *i, vec_ZZ& lambda, Matrix *T,
 319                     Value lcm, int p, Vector *val,
 320                     evalue *E, evalue* ev,
 321                     ZZ& offset)
 322 {
 323     unsigned nparam = T->NbRows - 1;
 324     unsigned dim = i->Dimension;
 325     Value tmp;
 326     ZZ nump;
 327
 328     if (p == nparam) {
 329         vec_ZZ vertex;
 330         ZZ num, l;
 331         Vector * values = Vector_Alloc(dim + 1);
 332         Vector_Matrix_Product(val->p, T, values->p);
 333         value_assign(values->p[dim], lcm);
 334         lattice_point(values->p, i, vertex);
 335         num = vertex * lambda;
 336         value2zz(lcm, l);
 337         num *= l;
 338         num += offset;
 339         value_init(ev->x.n);
 340         zz2value(num, ev->x.n);
 341         value_assign(ev->d, lcm);
 342         Vector_Free(values);
 343         return;
 344     }
 345
 346     value_init(tmp);
 347     vec_ZZ vertex;
 348     values2zz(T->p[p], vertex, dim);
 349     nump = vertex * lambda;
 350     if (First_Non_Zero(val->p, p) == -1) {
 351         value_assign(tmp, lcm);
 352         evalue *ET = term(p, nump, &tmp);
 353         eadd(ET, E);
 354         free_evalue_refs(ET);
 355         delete ET;
 356     }
 357
 358     value_assign(tmp, lcm);
 359     if (First_Non_Zero(T->p[p], dim) != -1)
 360         Vector_Gcd(T->p[p], dim, &tmp);
 361     Gcd(tmp, lcm, &tmp);
 362     if (value_lt(tmp, lcm)) {
 363         ZZ count;
 364
 365         value_division(tmp, lcm, tmp);
 366         value_set_si(ev->d, 0);
 367         ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
 368         value2zz(tmp, count);
 369         do {
 370             value_decrement(tmp, tmp);
 371             --count;
 372             ZZ new_offset = offset - count * nump;
 373             value_assign(val->p[p], tmp);
 374             vertex_period(i, lambda, T, lcm, p+1, val, E,
 375                           &ev->x.p->arr[VALUE_TO_INT(tmp)], new_offset);
 376         } while (value_pos_p(tmp));
 377     } else
 378         vertex_period(i, lambda, T, lcm, p+1, val, E, ev, offset);
 379     value_clear(tmp);
 380 }
 381
 382 /* Returns the power of (t+1) in the term of a rational generating function,
 383  * i.e., the scalar product of the actual lattice point and lambda.
 384  * The lattice point is the unique lattice point in the fundamental parallelepiped
 385  * of the unimodual cone i shifted to the parametric vertex W/lcm.
 386  *
 387  * The rows of W refer to the coordinates of the vertex
 388  * The first nparam columns are the coefficients of the parameters
 389  * and the final column is the constant term.
 390  * lcm is the common denominator of all coefficients.
 391  *
 392  * PD is the parameter domain, which, if != NULL, may be used to simply the
 393  * resulting expression.
 394  */
 395 #ifdef USE_MODULO
 396 evalue* lattice_point(
 397     Polyhedron *i, vec_ZZ& lambda, Matrix *W, Value lcm, Polyhedron *PD)
 398 {
 399     unsigned nparam = W->NbColumns - 1;
 400
 401     Matrix* Rays = rays2(i);
 402     Matrix *T = Transpose(Rays);
 403     Matrix *T2 = Matrix_Copy(T);
 404     Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
 405     int ok = Matrix_Inverse(T2, inv);
 406     assert(ok);
 407     Matrix_Free(Rays);
 408     Matrix_Free(T2);
 409     mat_ZZ vertex;
 410     matrix2zz(W, vertex, W->NbRows, W->NbColumns);
 411
 412     vec_ZZ num;
 413     num = lambda * vertex;
 414
 415     evalue *EP = multi_monom(num);
 416
 417     evalue tmp;
 418     value_init(tmp.d);
 419     value_init(tmp.x.n);
 420     value_set_si(tmp.x.n, 1);
 421     value_assign(tmp.d, lcm);
 422
 423     emul(&tmp, EP);
 424
 425     Matrix *L = Matrix_Alloc(inv->NbRows, W->NbColumns);
 426     Matrix_Product(inv, W, L);
 427
 428     mat_ZZ RT;
 429     matrix2zz(T, RT, T->NbRows, T->NbColumns);
 430     Matrix_Free(T);
 431
 432     vec_ZZ p = lambda * RT;
 433
 434     for (int i = 0; i < L->NbRows; ++i) {
 435         ceil_mod(L->p[i], nparam+1, lcm, p[i], EP, PD);
 436     }
 437
 438     Matrix_Free(L);
 439
 440     Matrix_Free(inv);
 441     free_evalue_refs(&tmp);
 442     return EP;
 443 }
 444 #else
 445 evalue* lattice_point(
 446     Polyhedron *i, vec_ZZ& lambda, Matrix *W, Value lcm, Polyhedron *PD)
 447 {
 448     Matrix *T = Transpose(W);
 449     unsigned nparam = T->NbRows - 1;
 450
 451     evalue *EP = new evalue();
 452     value_init(EP->d);
 453     evalue_set_si(EP, 0, 1);
 454
 455     evalue ev;
 456     Vector *val = Vector_Alloc(nparam+1);
 457     value_set_si(val->p[nparam], 1);
 458     ZZ offset(INIT_VAL, 0);
 459     value_init(ev.d);
 460     vertex_period(i, lambda, T, lcm, 0, val, EP, &ev, offset);
 461     Vector_Free(val);
 462     eadd(&ev, EP);
 463     free_evalue_refs(&ev);
 464
 465     Matrix_Free(T);
 466
 467     reduce_evalue(EP);
 468
 469     return EP;
 470 }
 471 #endif
 472
 473 /* returns the unique lattice point in the fundamental parallelepiped
 474  * of the unimodual cone C shifted to the parametric vertex V.
 475  *
 476  * The return values num and E_vertex are such that
 477  * coordinate i of this lattice point is equal to
 478  *
 479  *          num[i] + E_vertex[i]
 480  */
 481 void lattice_point(Param_Vertices *V, Polyhedron *C, vec_ZZ& num,
 482                    evalue **E_vertex)
 483 {
 484     unsigned nparam = V->Vertex->NbColumns - 2;
 485     unsigned dim = C->Dimension;
 486     vec_ZZ vertex;
 487     vertex.SetLength(nparam+1);
 488
 489     Value lcm, tmp;
 490     value_init(lcm);
 491     value_init(tmp);
 492     value_set_si(lcm, 1);
 493
 494     for (int j = 0; j < V->Vertex->NbRows; ++j) {
 495         value_lcm(lcm, V->Vertex->p[j][nparam+1], &lcm);
 496     }
 497
 498     if (value_notone_p(lcm)) {
 499         Matrix * mv = Matrix_Alloc(dim, nparam+1);
 500         for (int j = 0 ; j < dim; ++j) {
 501             value_division(tmp, lcm, V->Vertex->p[j][nparam+1]);
 502             Vector_Scale(V->Vertex->p[j], mv->p[j], tmp, nparam+1);
 503         }
 504
 505         Matrix* Rays = rays2(C);
 506         Matrix *T = Transpose(Rays);
 507         Matrix *T2 = Matrix_Copy(T);
 508         Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
 509         int ok = Matrix_Inverse(T2, inv);
 510         assert(ok);
 511         Matrix_Free(Rays);
 512         Matrix_Free(T2);
 513         Matrix *L = Matrix_Alloc(inv->NbRows, mv->NbColumns);
 514         Matrix_Product(inv, mv, L);
 515         Matrix_Free(inv);
 516
 517         evalue f;
 518         value_init(f.d);
 519         value_init(f.x.n);
 520
 521         ZZ one;
 522
 523         evalue *remainders[dim];
 524         for (int i = 0; i < dim; ++i) {
 525             remainders[i] = evalue_zero();
 526             one = 1;
 527             ceil(L->p[i], nparam+1, lcm, one, remainders[i], 0);
 528         }
 529         Matrix_Free(L);
 530
 531
 532         for (int i = 0; i < V->Vertex->NbRows; ++i) {
 533             values2zz(mv->p[i], vertex, nparam+1);
 534             E_vertex[i] = multi_monom(vertex);
 535             num[i] = 0;
 536
 537             value_set_si(f.x.n, 1);
 538             value_assign(f.d, lcm);
 539
 540             emul(&f, E_vertex[i]);
 541
 542             for (int j = 0; j < dim; ++j) {
 543                 if (value_zero_p(T->p[i][j]))
 544                     continue;
 545                 evalue cp;
 546                 value_init(cp.d);
 547                 evalue_copy(&cp, remainders[j]);
 548                 if (value_notone_p(T->p[i][j])) {
 549                     value_set_si(f.d, 1);
 550                     value_assign(f.x.n, T->p[i][j]);
 551                     emul(&f, &cp);
 552                 }
 553                 eadd(&cp, E_vertex[i]);
 554                 free_evalue_refs(&cp);
 555             }
 556         }
 557         for (int i = 0; i < dim; ++i) {
 558             free_evalue_refs(remainders[i]);
 559             free(remainders[i]);
 560         }
 561
 562         free_evalue_refs(&f);
 563
 564         Matrix_Free(T);
 565         Matrix_Free(mv);
 566         value_clear(lcm);
 567         value_clear(tmp);
 568         return;
 569     }
 570     value_clear(lcm);
 571     value_clear(tmp);
 572
 573     for (int i = 0; i < V->Vertex->NbRows; ++i) {
 574         /* fixed value */
 575         if (First_Non_Zero(V->Vertex->p[i], nparam) == -1) {
 576             E_vertex[i] = 0;
 577             value2zz(V->Vertex->p[i][nparam], num[i]);
 578         } else {
 579             values2zz(V->Vertex->p[i], vertex, nparam+1);
 580             E_vertex[i] = multi_monom(vertex);
 581             num[i] = 0;
 582         }
 583     }
 584 }
 585