bfcounter.cc

   1 #include <vector>
   2 #include "bfcounter.h"
   3 #include "lattice_point.h"
   4
   5 using std::vector;
   6 using std::cerr;
   7 using std::endl;
   8
   9 static int lex_cmp(vec_ZZ& a, vec_ZZ& b)
  10 {
  11     assert(a.length() == b.length());
  12
  13     for (int j = 0; j < a.length(); ++j)
  14         if (a[j] != b[j])
  15             return a[j] < b[j] ? -1 : 1;
  16     return 0;
  17 }
  18
  19 void bf_base::add_term(bfc_term_base *t, vec_ZZ& num_orig, vec_ZZ& extra_num)
  20 {
  21     vec_ZZ num;
  22     int d = num_orig.length();
  23     num.SetLength(d-1);
  24     for (int l = 0; l < d-1; ++l)
  25         num[l] = num_orig[l+1] + extra_num[l];
  26
  27     add_term(t, num);
  28 }
  29
  30 void bf_base::add_term(bfc_term_base *t, vec_ZZ& num)
  31 {
  32     int len = t->terms.NumRows();
  33     int i, r;
  34     for (i = 0; i < len; ++i) {
  35         r = lex_cmp(t->terms[i], num);
  36         if (r >= 0)
  37             break;
  38     }
  39     if (i == len || r > 0) {
  40         t->terms.SetDims(len+1, num.length());
  41         insert_term(t, i);
  42         t->terms[i] = num;
  43     } else {
  44         // i < len && r == 0
  45         update_term(t, i);
  46     }
  47 }
  48
  49 bfc_term_base* bf_base::find_bfc_term(bfc_vec& v, int *powers, int len)
  50 {
  51     bfc_vec::iterator i;
  52     for (i = v.begin(); i != v.end(); ++i) {
  53         int j;
  54         for (j = 0; j < len; ++j)
  55             if ((*i)->powers[j] != powers[j])
  56                 break;
  57         if (j == len)
  58             return (*i);
  59         if ((*i)->powers[j] > powers[j])
  60             break;
  61     }
  62
  63     bfc_term_base* t = new_bf_term(len);
  64     v.insert(i, t);
  65     memcpy(t->powers, powers, len * sizeof(int));
  66
  67     return t;
  68 }
  69
  70 void bf_base::reduce(mat_ZZ& factors, bfc_vec& v, barvinok_options *options)
  71 {
  72     assert(v.size() > 0);
  73     unsigned d = factors.NumCols();
  74
  75     if (d == lower)
  76         return base(factors, v);
  77
  78     bf_reducer bfr(factors, v, this);
  79
  80     bfr.reduce(options);
  81
  82     if (bfr.vn.size() > 0)
  83         reduce(bfr.nfactors, bfr.vn, options);
  84 }
  85
  86 int bf_base::setup_factors(const mat_ZZ& rays, mat_ZZ& factors,
  87                                     bfc_term_base* t, int s)
  88 {
  89     factors.SetDims(dim, dim);
  90
  91     int r;
  92
  93     for (r = 0; r < dim; ++r)
  94         t->powers[r] = 1;
  95
  96     for (r = 0; r < dim; ++r) {
  97         factors[r] = rays[r];
  98         int k;
  99         for (k = 0; k < dim; ++k)
 100             if (factors[r][k] != 0)
 101                 break;
 102         if (factors[r][k] < 0) {
 103             factors[r] = -factors[r];
 104             for (int i = 0; i < t->terms.NumRows(); ++i)
 105                 t->terms[i] += factors[r];
 106             s = -s;
 107         }
 108     }
 109
 110     return s;
 111 }
 112
 113 void bf_base::handle(const mat_ZZ& rays, Value *vertex, const QQ& c,
 114                      unsigned long det, barvinok_options *options)
 115 {
 116     bfc_term* t = new bfc_term(dim);
 117     vector< bfc_term_base * > v;
 118     v.push_back(t);
 119
 120     Matrix *points = Matrix_Alloc(det, dim);
 121     Matrix* Rays = zz2matrix(rays);
 122     lattice_points_fixed(vertex, vertex, Rays, Rays, points, det);
 123     Matrix_Free(Rays);
 124     matrix2zz(points, t->terms, points->NbRows, points->NbColumns);
 125     Matrix_Free(points);
 126
 127     // the elements of factors are always lexpositive
 128     mat_ZZ   factors;
 129     int s = setup_factors(rays, factors, t, 1);
 130
 131     t->c.SetLength(t->terms.NumRows());
 132
 133     for (int i = 0; i < t->c.length(); ++i) {
 134         t->c[i].n = s * c.n;
 135         t->c[i].d = c.d;
 136     }
 137
 138     reduce(factors, v, options);
 139 }
 140
 141 bfc_term_base* bfcounter_base::new_bf_term(int len)
 142 {
 143     bfc_term* t = new bfc_term(len);
 144     t->c.SetLength(0);
 145     return t;
 146 }
 147
 148 void bfcounter_base::set_factor(bfc_term_base *t, int k, int change)
 149 {
 150     bfc_term* bfct = static_cast<bfc_term *>(t);
 151     c = bfct->c[k];
 152     if (change)
 153         c.n = -c.n;
 154 }
 155
 156 void bfcounter_base::set_factor(bfc_term_base *t, int k, mpq_t &f, int change)
 157 {
 158     bfc_term* bfct = static_cast<bfc_term *>(t);
 159     value2zz(mpq_numref(f), c.n);
 160     value2zz(mpq_denref(f), c.d);
 161     c *= bfct->c[k];
 162     if (change)
 163         c.n = -c.n;
 164 }
 165
 166 void bfcounter_base::set_factor(bfc_term_base *t, int k, const QQ& c_factor,
 167                                 int change)
 168 {
 169     bfc_term* bfct = static_cast<bfc_term *>(t);
 170     c = bfct->c[k];
 171     c *= c_factor;
 172     if (change)
 173         c.n = -c.n;
 174 }
 175
 176 void bfcounter_base::insert_term(bfc_term_base *t, int i)
 177 {
 178     bfc_term* bfct = static_cast<bfc_term *>(t);
 179     int len = t->terms.NumRows()-1;     // already increased by one
 180
 181     bfct->c.SetLength(len+1);
 182     for (int j = len; j > i; --j) {
 183         bfct->c[j] = bfct->c[j-1];
 184         t->terms[j] = t->terms[j-1];
 185     }
 186     bfct->c[i] = c;
 187 }
 188
 189 void bfcounter_base::update_term(bfc_term_base *t, int i)
 190 {
 191     bfc_term* bfct = static_cast<bfc_term *>(t);
 192
 193     bfct->c[i] += c;
 194 }
 195
 196 void bf_reducer::compute_extra_num(int i)
 197 {
 198     clear(extra_num);
 199     changes = 0;
 200     no_param = 0;           // r from text
 201     only_param = 0;         // k-r-s from text
 202     total_power = 0;        // k from text
 203
 204     for (int j = 0; j < nf; ++j) {
 205         if (v[i]->powers[j] == 0)
 206             continue;
 207
 208         total_power += v[i]->powers[j];
 209         if (factors[j][0] == 0) {
 210             only_param += v[i]->powers[j];
 211             continue;
 212         }
 213
 214         if (old2new[j] == -1)
 215             no_param += v[i]->powers[j];
 216         else
 217             extra_num += -sign[j] * v[i]->powers[j] * nfactors[old2new[j]];
 218         changes += v[i]->powers[j];
 219     }
 220 }
 221
 222 void bf_reducer::update_powers(const std::vector<int>& powers)
 223 {
 224     for (int l = 0; l < nnf; ++l)
 225         npowers[l] = bpowers[l];
 226
 227     l_extra_num = extra_num;
 228     l_changes = changes;
 229
 230     for (int l = 0; l < powers.size(); ++l) {
 231         int n = powers[l];
 232         if (n == 0)
 233             continue;
 234         assert(old2new[l] != -1);
 235
 236         npowers[old2new[l]] += n;
 237         // interpretation of sign has been inverted
 238         // since we inverted the power for specialization
 239         if (sign[l] == 1) {
 240             l_extra_num += n * nfactors[old2new[l]];
 241             l_changes += n;
 242         }
 243     }
 244 }
 245
 246
 247 void bf_reducer::compute_reduced_factors()
 248 {
 249     unsigned nf = factors.NumRows();
 250     unsigned d = factors.NumCols();
 251     nnf = 0;
 252     nfactors.SetDims(nnf, d-1);
 253
 254     for (int i = 0; i < nf; ++i) {
 255         int j;
 256         int s = 1;
 257         for (j = 0; j < nnf; ++j) {
 258             int k;
 259             for (k = 1; k < d; ++k)
 260                 if (factors[i][k] != 0 || nfactors[j][k-1] != 0)
 261                     break;
 262             if (k < d && factors[i][k] == -nfactors[j][k-1])
 263                 s = -1;
 264             for (; k < d; ++k)
 265                 if (factors[i][k] != s * nfactors[j][k-1])
 266                     break;
 267             if (k == d)
 268                 break;
 269         }
 270         old2new[i] = j;
 271         if (j == nnf) {
 272             int k;
 273             for (k = 1; k < d; ++k)
 274                 if (factors[i][k] != 0)
 275                     break;
 276             if (k < d) {
 277                 if (factors[i][k] < 0)
 278                     s = -1;
 279                 nfactors.SetDims(++nnf, d-1);
 280                 for (int k = 1; k < d; ++k)
 281                     nfactors[j][k-1] = s * factors[i][k];
 282             } else
 283                 old2new[i] = -1;
 284         }
 285         sign[i] = s;
 286     }
 287     npowers = new int[nnf];
 288     bpowers = new int[nnf];
 289 }
 290
 291 void bf_reducer::reduce(barvinok_options *options)
 292 {
 293     compute_reduced_factors();
 294
 295     Value tmp;
 296     value_init(tmp);
 297     for (int i = 0; i < v.size(); ++i) {
 298         compute_extra_num(i);
 299
 300         if (no_param == 0) {
 301             vec_ZZ extra_num;
 302             extra_num.SetLength(d-1);
 303             int changes = 0;
 304             int *npowers = new int[nnf];
 305             for (int k = 0; k < nnf; ++k)
 306                 npowers[k] = 0;
 307             for (int k = 0; k < nf; ++k) {
 308                 assert(old2new[k] != -1);
 309                 npowers[old2new[k]] += v[i]->powers[k];
 310                 if (sign[k] == -1) {
 311                     extra_num += v[i]->powers[k] * nfactors[old2new[k]];
 312                     changes += v[i]->powers[k];
 313                 }
 314             }
 315
 316             bfc_term_base * t = bf->find_bfc_term(vn, npowers, nnf);
 317             for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
 318                 bf->set_factor(v[i], k, changes % 2);
 319                 bf->add_term(t, v[i]->terms[k], extra_num);
 320             }
 321             delete [] npowers;
 322         } else {
 323             // powers of "constant" part
 324             for (int k = 0; k < nnf; ++k)
 325                 bpowers[k] = 0;
 326             for (int k = 0; k < nf; ++k) {
 327                 if (factors[k][0] != 0)
 328                     continue;
 329                 assert(old2new[k] != -1);
 330                 bpowers[old2new[k]] += v[i]->powers[k];
 331                 if (sign[k] == -1) {
 332                     extra_num += v[i]->powers[k] * nfactors[old2new[k]];
 333                     changes += v[i]->powers[k];
 334                 }
 335             }
 336
 337             int j;
 338             for (j = 0; j < nf; ++j)
 339                 if (old2new[j] == -1 && v[i]->powers[j] > 0)
 340                     break;
 341
 342             zz2value(factors[j][0], tmp);
 343             dpoly D(no_param, tmp, 1);
 344             for (int k = 1; k < v[i]->powers[j]; ++k) {
 345                 dpoly fact(no_param, tmp, 1);
 346                 D *= fact;
 347             }
 348             for ( ; ++j < nf; )
 349                 if (old2new[j] == -1)  {
 350                     zz2value(factors[j][0], tmp);
 351                     for (int k = 0; k < v[i]->powers[j]; ++k) {
 352                         dpoly fact(no_param, tmp, 1);
 353                         D *= fact;
 354                     }
 355                 }
 356
 357             if (no_param + only_param == total_power &&
 358                     bf->constant_vertex(d)) {
 359                 bfc_term_base * t = NULL;
 360                 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
 361                     zz2value(v[i]->terms[k][0], tmp);
 362                     dpoly n(no_param, tmp);
 363                     mpq_set_si(bf->tcount, 0, 1);
 364                     n.div(D, bf->tcount, 1);
 365
 366                     if (value_zero_p(mpq_numref(bf->tcount)))
 367                         continue;
 368
 369                     if (!t)
 370                          t = bf->find_bfc_term(vn, bpowers, nnf);
 371                     bf->set_factor(v[i], k, bf->tcount, changes % 2);
 372                     bf->add_term(t, v[i]->terms[k], extra_num);
 373                 }
 374             } else {
 375                 for (int j = 0; j < v[i]->terms.NumRows(); ++j) {
 376                     zz2value(v[i]->terms[j][0], tmp);
 377                     dpoly n(no_param, tmp);
 378
 379                     dpoly_r * r = 0;
 380                     if (no_param + only_param == total_power)
 381                         r = new dpoly_r(n, nf);
 382                     else
 383                         for (int k = 0; k < nf; ++k) {
 384                             if (v[i]->powers[k] == 0)
 385                                 continue;
 386                             if (factors[k][0] == 0 || old2new[k] == -1)
 387                                 continue;
 388
 389                             zz2value(factors[k][0], tmp);
 390                             dpoly pd(no_param-1, tmp, 1);
 391
 392                             for (int l = 0; l < v[i]->powers[k]; ++l) {
 393                                 int q;
 394                                 for (q = 0; q < k; ++q)
 395                                     if (old2new[q] == old2new[k] &&
 396                                         sign[q] == sign[k])
 397                                             break;
 398
 399                                 if (r == 0)
 400                                     r = new dpoly_r(n, pd, q, nf);
 401                                 else {
 402                                     dpoly_r *nr = new dpoly_r(r, pd, q, nf);
 403                                     delete r;
 404                                     r = nr;
 405                                 }
 406                             }
 407                         }
 408
 409                     dpoly_r *rc = r->div(D);
 410                     delete r;
 411                     QQ factor;
 412                     factor.d = rc->denom;
 413
 414                     if (bf->constant_vertex(d)) {
 415                         dpoly_r_term_list& final = rc->c[rc->len-1];
 416
 417                         dpoly_r_term_list::iterator k;
 418                         for (k = final.begin(); k != final.end(); ++k) {
 419                             if ((*k)->coeff == 0)
 420                                 continue;
 421
 422                             update_powers((*k)->powers);
 423
 424                             bfc_term_base * t = bf->find_bfc_term(vn, npowers, nnf);
 425                             factor.n = (*k)->coeff;
 426                             bf->set_factor(v[i], j, factor, l_changes % 2);
 427                             bf->add_term(t, v[i]->terms[j], l_extra_num);
 428                         }
 429                     } else
 430                         bf->cum(this, v[i], j, rc, options);
 431
 432                     delete rc;
 433                 }
 434             }
 435         }
 436         delete v[i];
 437     }
 438     value_clear(tmp);
 439 }