bfcounter.cc

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