genfun.cc

   1 #include <iostream>
   2 #include <assert.h>
   3 #include <genfun.h>
   4 #include <barvinok.h>
   5 #include <barvinok2.h>
   6
   7 using std::cout;
   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 static int lex_cmp(mat_ZZ& a, mat_ZZ& b)
  20 {
  21     assert(a.NumCols() == b.NumCols());
  22     int alen = a.NumRows();
  23     int blen = b.NumRows();
  24     int len = alen < blen ? alen : blen;
  25
  26     for (int i = 0; i < len; ++i) {
  27         int s = lex_cmp(a[i], b[i]);
  28         if (s)
  29             return s;
  30     }
  31     return alen-blen;
  32 }
  33
  34 void gen_fun::add(ZZ& cn, ZZ& cd, vec_ZZ& num, mat_ZZ& den)
  35 {
  36     if (cn == 0)
  37         return;
  38
  39     short_rat * r = new short_rat;
  40     r->n.coeff.SetDims(1, 2);
  41     r->n.coeff[0][0] = cn;
  42     r->n.coeff[0][1] = cd;
  43     r->n.power.SetDims(1, num.length());
  44     r->n.power[0] = num;
  45     r->d.power = den;
  46
  47     for (int i = 0; i < r->d.power.NumRows(); ++i) {
  48         int j;
  49         for (j = 0; j < r->d.power.NumCols(); ++j)
  50             if (r->d.power[i][j] != 0)
  51                 break;
  52         if (r->d.power[i][j] < 0) {
  53             r->d.power[i] = -r->d.power[i];
  54             r->n.coeff[0][0] = -r->n.coeff[0][0];
  55             r->n.power[0] += r->d.power[i];
  56         }
  57     }
  58
  59     for (int i = 0; i < term.size(); ++i)
  60         if (lex_cmp(term[i]->d.power, r->d.power) == 0) {
  61             int len = term[i]->n.coeff.NumRows();
  62             int dim = term[i]->n.power.NumCols();
  63             term[i]->n.coeff.SetDims(len+1, 2);
  64             term[i]->n.power.SetDims(len+1, dim);
  65             term[i]->n.coeff[len] = r->n.coeff[0];
  66             term[i]->n.power[len] = r->n.power[0];
  67             delete r;
  68             return;
  69         }
  70
  71     term.push_back(r);
  72 }
  73
  74 static void print_power(vec_ZZ& c, vec_ZZ& p,
  75                         unsigned int nparam, char **param_name)
  76 {
  77     bool first = true;
  78
  79     for (int i = 0; i < p.length(); ++i) {
  80         if (p[i] == 0)
  81             continue;
  82         if (first) {
  83             if (c[0] == -1 && c[1] == 1)
  84                 cout << "-";
  85             else if (c[0] != 1 || c[1] != 1) {
  86                 cout << c[0];
  87                 if (c[1] != 1)
  88                     cout << " / " << c[1];
  89                 cout << " * ";
  90             }
  91             first = false;
  92         } else
  93             cout << " * ";
  94         if (i < nparam)
  95             cout << param_name[i];
  96         else
  97             cout << "x" << i;
  98         if (p[i] == 1)
  99             continue;
 100         if (p[i] < 0)
 101             cout << "^(" << p[i] << ")";
 102         else
 103             cout << "^" << p[i];
 104     }
 105     if (first) {
 106         cout << c[0];
 107         if (c[1] != 1)
 108             cout << " / " << c[1];
 109     }
 110 }
 111
 112 void gen_fun::print(unsigned int nparam, char **param_name)
 113 {
 114     vec_ZZ mone;
 115     mone.SetLength(2);
 116     mone[0] = -1;
 117     mone[1] = 1;
 118     for (int i = 0; i < term.size(); ++i) {
 119         if (i != 0)
 120             cout << " + ";
 121         cout << "(";
 122         for (int j = 0; j < term[i]->n.coeff.NumRows(); ++j) {
 123             if (j != 0 && term[i]->n.coeff[j][0] > 0)
 124                 cout << "+";
 125             print_power(term[i]->n.coeff[j], term[i]->n.power[j],
 126                         nparam, param_name);
 127         }
 128         cout << ")/(";
 129         for (int j = 0; j < term[i]->d.power.NumRows(); ++j) {
 130             if (j != 0)
 131                 cout << " * ";
 132             cout << "(1";
 133             print_power(mone, term[i]->d.power[j],
 134                         nparam, param_name);
 135             cout << ")";
 136         }
 137         cout << ")";
 138     }
 139 }
 140
 141 gen_fun::operator evalue *()
 142 {
 143     evalue *EP = NULL;
 144     evalue factor;
 145     value_init(factor.d);
 146     value_init(factor.x.n);
 147     for (int i = 0; i < term.size(); ++i) {
 148         unsigned nvar = term[i]->d.power.NumRows();
 149         unsigned nparam = term[i]->d.power.NumCols();
 150         Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + nparam + 1);
 151         mat_ZZ& d = term[i]->d.power;
 152         Polyhedron *U = Universe_Polyhedron(nparam);
 153
 154         for (int j = 0; j < term[i]->n.coeff.NumRows(); ++j) {
 155             for (int r = 0; r < nparam; ++r) {
 156                 value_set_si(C->p[r][0], 0);
 157                 for (int c = 0; c < nvar; ++c) {
 158                     zz2value(d[c][r], C->p[r][1+c]);
 159                 }
 160                 Vector_Set(&C->p[r][1+nvar], 0, nparam);
 161                 value_set_si(C->p[r][1+nvar+r], -1);
 162                 zz2value(term[i]->n.power[j][r], C->p[r][1+nvar+nparam]);
 163             }
 164             for (int r = 0; r < nvar; ++r) {
 165                 value_set_si(C->p[nparam+r][0], 1);
 166                 Vector_Set(&C->p[nparam+r][1], 0, nvar + nparam + 1);
 167                 value_set_si(C->p[nparam+r][1+r], 1);
 168             }
 169             Polyhedron *P = Constraints2Polyhedron(C, 0);
 170             evalue *E = barvinok_enumerate_ev(P, U, 0);
 171             Polyhedron_Free(P);
 172             zz2value(term[i]->n.coeff[j][0], factor.x.n);
 173             zz2value(term[i]->n.coeff[j][1], factor.d);
 174             emul(&factor, E);
 175             /*
 176             Matrix_Print(stdout, P_VALUE_FMT, C);
 177             char *test[] = { "A", "B", "C", "D", "E", "F", "G" };
 178             print_evalue(stdout, E, test);
 179             */
 180             if (!EP)
 181                 EP = E;
 182             else {
 183                 eadd(E, EP);
 184                 free_evalue_refs(E);
 185             }
 186         }
 187         Matrix_Free(C);
 188         Polyhedron_Free(U);
 189     }
 190     value_clear(factor.d);
 191     value_clear(factor.x.n);
 192     return EP;
 193 }