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