lattice_point.cc

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