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 static void ceil(Value *coef, int len, Value d, ZZ& f,
 241                  evalue *EP, Polyhedron *PD, barvinok_options *options)
 242 {
 243     ceil_mod(coef, len, d, f, EP, PD);
 244     if (options->lookup_table)
 245         evalue_mod2table(EP, len-1);
 246 }
 247
 248 evalue* bv_ceil3(Value *coef, int len, Value d, Polyhedron *P)
 249 {
 250     Vector *val = Vector_Alloc(len);
 251
 252     Value t;
 253     value_init(t);
 254     value_set_si(t, -1);
 255     Vector_Scale(coef, val->p, t, len);
 256     value_absolute(t, d);
 257
 258     vec_ZZ num;
 259     values2zz(val->p, num, len);
 260     evalue *EP = multi_monom(num);
 261
 262     evalue tmp;
 263     value_init(tmp.d);
 264     value_init(tmp.x.n);
 265     value_set_si(tmp.x.n, 1);
 266     value_assign(tmp.d, t);
 267
 268     emul(&tmp, EP);
 269
 270     ZZ one;
 271     one = 1;
 272     ceil_mod(val->p, len, t, one, EP, P);
 273     value_clear(t);
 274
 275     /* copy EP to malloc'ed evalue */
 276     evalue *E = ALLOC(evalue);
 277     *E = *EP;
 278     delete EP;
 279
 280     free_evalue_refs(&tmp);
 281     Vector_Free(val);
 282
 283     return E;
 284 }
 285
 286 void lattice_point(Value* values, const mat_ZZ& rays, vec_ZZ& vertex, int *closed)
 287 {
 288     unsigned dim = rays.NumRows();
 289     if (value_one_p(values[dim]) && !closed)
 290         values2zz(values, vertex, dim);
 291     else {
 292         Matrix* Rays = rays2matrix(rays);
 293         Matrix *inv = Matrix_Alloc(Rays->NbRows, Rays->NbColumns);
 294         int ok = Matrix_Inverse(Rays, inv);
 295         assert(ok);
 296         Matrix_Free(Rays);
 297         Rays = rays2matrix(rays);
 298         Vector *lambda = Vector_Alloc(dim+1);
 299         Vector_Matrix_Product(values, inv, lambda->p);
 300         Matrix_Free(inv);
 301         for (int j = 0; j < dim; ++j)
 302             if (!closed || closed[j])
 303                 mpz_cdiv_q(lambda->p[j], lambda->p[j], lambda->p[dim]);
 304             else {
 305                 value_addto(lambda->p[j], lambda->p[j], lambda->p[dim]);
 306                 mpz_fdiv_q(lambda->p[j], lambda->p[j], lambda->p[dim]);
 307             }
 308         value_set_si(lambda->p[dim], 1);
 309         Vector *A = Vector_Alloc(dim+1);
 310         Vector_Matrix_Product(lambda->p, Rays, A->p);
 311         Vector_Free(lambda);
 312         Matrix_Free(Rays);
 313         values2zz(A->p, vertex, dim);
 314         Vector_Free(A);
 315     }
 316 }
 317
 318 static void vertex_period(
 319                     const mat_ZZ& rays, vec_ZZ& lambda, Matrix *T,
 320                     Value lcm, int p, Vector *val,
 321                     evalue *E, evalue* ev,
 322                     ZZ& offset)
 323 {
 324     unsigned nparam = T->NbRows - 1;
 325     unsigned dim = rays.NumRows();
 326     Value tmp;
 327     ZZ nump;
 328
 329     if (p == nparam) {
 330         vec_ZZ vertex;
 331         ZZ num, l;
 332         Vector * values = Vector_Alloc(dim + 1);
 333         Vector_Matrix_Product(val->p, T, values->p);
 334         value_assign(values->p[dim], lcm);
 335         lattice_point(values->p, rays, vertex, NULL);
 336         num = vertex * lambda;
 337         value2zz(lcm, l);
 338         num *= l;
 339         num += offset;
 340         value_init(ev->x.n);
 341         zz2value(num, ev->x.n);
 342         value_assign(ev->d, lcm);
 343         Vector_Free(values);
 344         return;
 345     }
 346
 347     value_init(tmp);
 348     vec_ZZ vertex;
 349     values2zz(T->p[p], vertex, dim);
 350     nump = vertex * lambda;
 351     if (First_Non_Zero(val->p, p) == -1) {
 352         value_assign(tmp, lcm);
 353         evalue *ET = term(p, nump, &tmp);
 354         eadd(ET, E);
 355         free_evalue_refs(ET);
 356         delete ET;
 357     }
 358
 359     value_assign(tmp, lcm);
 360     if (First_Non_Zero(T->p[p], dim) != -1)
 361         Vector_Gcd(T->p[p], dim, &tmp);
 362     Gcd(tmp, lcm, &tmp);
 363     if (value_lt(tmp, lcm)) {
 364         ZZ count;
 365
 366         value_division(tmp, lcm, tmp);
 367         value_set_si(ev->d, 0);
 368         ev->x.p = new_enode(periodic, VALUE_TO_INT(tmp), p+1);
 369         value2zz(tmp, count);
 370         do {
 371             value_decrement(tmp, tmp);
 372             --count;
 373             ZZ new_offset = offset - count * nump;
 374             value_assign(val->p[p], tmp);
 375             vertex_period(rays, lambda, T, lcm, p+1, val, E,
 376                           &ev->x.p->arr[VALUE_TO_INT(tmp)], new_offset);
 377         } while (value_pos_p(tmp));
 378     } else
 379         vertex_period(rays, lambda, T, lcm, p+1, val, E, ev, offset);
 380     value_clear(tmp);
 381 }
 382
 383 /* Returns the power of (t+1) in the term of a rational generating function,
 384  * i.e., the scalar product of the actual lattice point and lambda.
 385  * The lattice point is the unique lattice point in the fundamental parallelepiped
 386  * of the unimodual cone i shifted to the parametric vertex W/lcm.
 387  *
 388  * The rows of W refer to the coordinates of the vertex
 389  * The first nparam columns are the coefficients of the parameters
 390  * and the final column is the constant term.
 391  * lcm is the common denominator of all coefficients.
 392  *
 393  * PD is the parameter domain, which, if != NULL, may be used to simply the
 394  * resulting expression.
 395  */
 396 static evalue* lattice_point_fractional(const mat_ZZ& rays, vec_ZZ& lambda,
 397                                         Matrix *W, Value lcm, Polyhedron *PD)
 398 {
 399     unsigned nparam = W->NbColumns - 1;
 400
 401     Matrix* Rays = rays2matrix2(rays);
 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
 445 static evalue* lattice_point_table(const mat_ZZ& rays, vec_ZZ& lambda, Matrix *W,
 446                                    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(rays, 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
 472 evalue* lattice_point(const mat_ZZ& rays, vec_ZZ& lambda, Matrix *W,
 473                       Value lcm, Polyhedron *PD, barvinok_options *options)
 474 {
 475     if (options->lookup_table)
 476         return lattice_point_table(rays, lambda, W, lcm, PD);
 477     else
 478         return lattice_point_fractional(rays, lambda, W, lcm, PD);
 479 }
 480
 481 /* returns the unique lattice point in the fundamental parallelepiped
 482  * of the unimodual cone C shifted to the parametric vertex V.
 483  *
 484  * The return values num and E_vertex are such that
 485  * coordinate i of this lattice point is equal to
 486  *
 487  *          num[i] + E_vertex[i]
 488  */
 489 void lattice_point(Param_Vertices *V, const mat_ZZ& rays, vec_ZZ& num,
 490                    evalue **E_vertex, barvinok_options *options)
 491 {
 492     unsigned nparam = V->Vertex->NbColumns - 2;
 493     unsigned dim = rays.NumCols();
 494     vec_ZZ vertex;
 495     vertex.SetLength(nparam+1);
 496
 497     Value lcm, tmp;
 498     value_init(lcm);
 499     value_init(tmp);
 500     value_set_si(lcm, 1);
 501
 502     for (int j = 0; j < V->Vertex->NbRows; ++j) {
 503         value_lcm(lcm, V->Vertex->p[j][nparam+1], &lcm);
 504     }
 505
 506     if (value_notone_p(lcm)) {
 507         Matrix * mv = Matrix_Alloc(dim, nparam+1);
 508         for (int j = 0 ; j < dim; ++j) {
 509             value_division(tmp, lcm, V->Vertex->p[j][nparam+1]);
 510             Vector_Scale(V->Vertex->p[j], mv->p[j], tmp, nparam+1);
 511         }
 512
 513         Matrix* Rays = rays2matrix2(rays);
 514         Matrix *T = Transpose(Rays);
 515         Matrix *T2 = Matrix_Copy(T);
 516         Matrix *inv = Matrix_Alloc(T2->NbRows, T2->NbColumns);
 517         int ok = Matrix_Inverse(T2, inv);
 518         assert(ok);
 519         Matrix_Free(Rays);
 520         Matrix_Free(T2);
 521         Matrix *L = Matrix_Alloc(inv->NbRows, mv->NbColumns);
 522         Matrix_Product(inv, mv, L);
 523         Matrix_Free(inv);
 524
 525         evalue f;
 526         value_init(f.d);
 527         value_init(f.x.n);
 528
 529         ZZ one;
 530
 531         evalue *remainders[dim];
 532         for (int i = 0; i < dim; ++i) {
 533             remainders[i] = evalue_zero();
 534             one = 1;
 535             ceil(L->p[i], nparam+1, lcm, one, remainders[i], 0, options);
 536         }
 537         Matrix_Free(L);
 538
 539
 540         for (int i = 0; i < V->Vertex->NbRows; ++i) {
 541             values2zz(mv->p[i], vertex, nparam+1);
 542             E_vertex[i] = multi_monom(vertex);
 543             num[i] = 0;
 544
 545             value_set_si(f.x.n, 1);
 546             value_assign(f.d, lcm);
 547
 548             emul(&f, E_vertex[i]);
 549
 550             for (int j = 0; j < dim; ++j) {
 551                 if (value_zero_p(T->p[i][j]))
 552                     continue;
 553                 evalue cp;
 554                 value_init(cp.d);
 555                 evalue_copy(&cp, remainders[j]);
 556                 if (value_notone_p(T->p[i][j])) {
 557                     value_set_si(f.d, 1);
 558                     value_assign(f.x.n, T->p[i][j]);
 559                     emul(&f, &cp);
 560                 }
 561                 eadd(&cp, E_vertex[i]);
 562                 free_evalue_refs(&cp);
 563             }
 564         }
 565         for (int i = 0; i < dim; ++i) {
 566             free_evalue_refs(remainders[i]);
 567             free(remainders[i]);
 568         }
 569
 570         free_evalue_refs(&f);
 571
 572         Matrix_Free(T);
 573         Matrix_Free(mv);
 574         value_clear(lcm);
 575         value_clear(tmp);
 576         return;
 577     }
 578     value_clear(lcm);
 579     value_clear(tmp);
 580
 581     for (int i = 0; i < V->Vertex->NbRows; ++i) {
 582         /* fixed value */
 583         if (First_Non_Zero(V->Vertex->p[i], nparam) == -1) {
 584             E_vertex[i] = 0;
 585             value2zz(V->Vertex->p[i][nparam], num[i]);
 586         } else {
 587             values2zz(V->Vertex->p[i], vertex, nparam+1);
 588             E_vertex[i] = multi_monom(vertex);
 589             num[i] = 0;
 590         }
 591     }
 592 }
 593