scarf.cc

   1 #include <assert.h>
   2 #include <vector>
   3 #include <barvinok/barvinok.h>
   4 #include <barvinok/util.h>
   5 #include "config.h"
   6
   7 using std::vector;
   8
   9 static Matrix *extract_matrix(Polyhedron *P, unsigned dim)
  10 {
  11     Matrix *A;
  12     int n_col;
  13
  14     n_col = 0;
  15     for (int i = 0; i < P->NbConstraints; ++i)
  16         if (value_notzero_p(P->Constraint[i][1+dim]) ||
  17             value_notzero_p(P->Constraint[i][1+dim+1]))
  18             ++n_col;
  19
  20     A = Matrix_Alloc(2, n_col+2);
  21     n_col = 0;
  22     for (int i = 0; i < P->NbConstraints; ++i) {
  23         if (value_zero_p(P->Constraint[i][1+dim]) &&
  24             value_zero_p(P->Constraint[i][1+dim+1]))
  25             continue;
  26         value_assign(A->p[0][n_col], P->Constraint[i][1+dim]);
  27         value_assign(A->p[1][n_col], P->Constraint[i][1+dim+1]);
  28         ++n_col;
  29     }
  30     value_set_si(A->p[0][n_col], 1);
  31     value_set_si(A->p[1][n_col+1], 1);
  32
  33     return A;
  34 }
  35
  36 static int lex_sign(Value *v, int len)
  37 {
  38     int first;
  39
  40     first = First_Non_Zero(v, len);
  41     return first == -1 ? 0 : value_sign(v[first]);
  42 }
  43
  44 static void set_pos(int pos[4], int actual, int wanted)
  45 {
  46     if (actual == wanted)
  47         return;
  48     int t = pos[actual];
  49     pos[actual] = pos[wanted];
  50     pos[wanted] = t;
  51 }
  52
  53 static Matrix *normalize_matrix(Matrix *A, int pos[4], int *n)
  54 {
  55     Matrix *T, *B;
  56     Value tmp, tmp2, factor;
  57     int type = -1;
  58
  59     value_init(tmp);
  60     value_init(tmp2);
  61     value_init(factor);
  62
  63     T = Matrix_Alloc(2, 2);
  64     Extended_Euclid(A->p[0][pos[0]], A->p[1][pos[0]],
  65                     &T->p[0][0], &T->p[0][1], &tmp);
  66     value_division(T->p[1][0], A->p[1][pos[0]], tmp);
  67     value_division(T->p[1][1], A->p[0][pos[0]], tmp);
  68     value_oppose(T->p[0][0], T->p[0][0]);
  69     value_oppose(T->p[0][1], T->p[0][1]);
  70     value_oppose(T->p[1][0], T->p[1][0]);
  71
  72     B = Matrix_Alloc(2, A->NbColumns);
  73     Matrix_Product(T, A, B);
  74     Matrix_Free(T);
  75
  76     /* Make zero in first position negative */
  77     if (lex_sign(B->p[1], B->NbColumns) > 0) {
  78         value_set_si(tmp, -1);
  79         Vector_Scale(B->p[1], B->p[1], tmp, B->NbColumns);
  80     }
  81
  82     /* First determine whether the matrix is of sign pattern I or II
  83      * (Theorem 1.11)
  84      */
  85     if (*n == 3) {
  86         assert(value_neg_p(B->p[1][pos[1]]));
  87         assert(value_pos_p(B->p[1][pos[2]]));
  88
  89         value_set_si(factor, 0);
  90         for (int i = 1; i <= 2; ++i) {
  91             value_pdivision(tmp, B->p[0][pos[i]], B->p[1][pos[i]]);
  92             value_increment(tmp, tmp);
  93             if (value_gt(tmp, factor))
  94                 value_assign(factor, tmp);
  95         }
  96         value_oppose(factor, factor);
  97         value_set_si(tmp, 1);
  98         Vector_Combine(B->p[0], B->p[1], B->p[0],
  99                        tmp, factor, B->NbColumns);
 100         Vector_Exchange(B->p[0], B->p[1], B->NbColumns);
 101         /* problems with three constraints are considered
 102          * to be of sign pattern II
 103          */
 104         type = 2;
 105     } else {
 106         int i;
 107         for (i = 1; i <= 3; ++i)
 108             if (value_zero_p(B->p[1][pos[i]]))
 109                 break;
 110         if (i <= 3) {
 111             /* put zero in position 3 */
 112             set_pos(pos, i, 3);
 113
 114             /* put positive one in position 1 */
 115             for (i = 1; i <= 3; ++i)
 116                 if (value_pos_p(B->p[1][pos[i]]))
 117                     break;
 118             set_pos(pos, i, 1);
 119
 120             value_set_si(factor, 0);
 121             for (int i = 1; i <= 2; ++i) {
 122                 value_pdivision(tmp, B->p[0][pos[i]], B->p[1][pos[i]]);
 123                 value_increment(tmp, tmp);
 124                 if (value_gt(tmp, factor))
 125                     value_assign(factor, tmp);
 126             }
 127             value_oppose(factor, factor);
 128             value_set_si(tmp, 1);
 129             Vector_Combine(B->p[0], B->p[1], B->p[0], tmp, factor, B->NbColumns);
 130
 131             assert(value_notzero_p(B->p[0][pos[3]]));
 132             type = value_pos_p(B->p[0][pos[3]]) ? 1 : 2;
 133         } else {
 134             int neg = 0;
 135             int sign = lex_sign(B->p[1], B->NbColumns);
 136             assert(sign < 0);
 137             for (int i = 1; i <= 3; ++i)
 138                 if (value_neg_p(B->p[1][pos[i]]))
 139                     ++neg;
 140             assert(neg == 1 || neg == 2);
 141             if (neg == 1) {
 142                 int i;
 143                 /* put negative one in position 1 */
 144                 for (i = 1; i <= 3; ++i)
 145                     if (value_neg_p(B->p[1][pos[i]]))
 146                         break;
 147                 set_pos(pos, i, 1);
 148
 149                 value_set_si(factor, 0);
 150                 for (int i = 1; i <= 3; ++i) {
 151                     value_pdivision(tmp, B->p[0][pos[i]], B->p[1][pos[i]]);
 152                     value_increment(tmp, tmp);
 153                     if (value_gt(tmp, factor))
 154                         value_assign(factor, tmp);
 155                 }
 156                 value_oppose(factor, factor);
 157                 value_set_si(tmp, 1);
 158                 Vector_Combine(B->p[0], B->p[1], B->p[0],
 159                                tmp, factor, B->NbColumns);
 160                 Vector_Exchange(B->p[0], B->p[1], B->NbColumns);
 161                 type = 1;
 162             } else {
 163                 int i;
 164                 /* put positive one in position 1 */
 165                 for (i = 1; i <= 3; ++i)
 166                     if (value_pos_p(B->p[1][pos[i]]))
 167                         break;
 168                 set_pos(pos, i, 1);
 169
 170                 value_set_si(factor, 0);
 171                 for (int i = 1; i <= 3; ++i) {
 172                     value_pdivision(tmp, B->p[0][pos[i]], B->p[1][pos[i]]);
 173                     value_increment(tmp, tmp);
 174                     if (value_gt(tmp, factor))
 175                         value_assign(factor, tmp);
 176                 }
 177                 value_oppose(factor, factor);
 178                 value_set_si(tmp, 1);
 179                 Vector_Combine(B->p[0], B->p[1], B->p[0],
 180                                tmp, factor, B->NbColumns);
 181                 type = 1;
 182             }
 183         }
 184     }
 185
 186     assert(type != -1);
 187
 188     if (type == 2) {
 189         for (;;) {
 190             value_oppose(tmp, B->p[0][pos[1]]);
 191             value_pdivision(factor, tmp, B->p[1][pos[1]]);
 192             value_oppose(tmp, B->p[1][pos[2]]);
 193             value_pdivision(tmp, tmp, B->p[0][pos[2]]);
 194             if (value_zero_p(factor) && value_zero_p(tmp))
 195                 break;
 196             assert(value_zero_p(factor) || value_zero_p(tmp));
 197             if (value_pos_p(factor)) {
 198                 value_set_si(tmp, 1);
 199                 Vector_Combine(B->p[0], B->p[1], B->p[0], tmp, factor, B->NbColumns);
 200                 if (value_zero_p(B->p[0][pos[1]])) {
 201                     /* We will deal with this later */
 202                     assert(lex_sign(B->p[0], B->NbColumns) < 0);
 203                 }
 204             } else {
 205                 value_set_si(factor, 1);
 206                 Vector_Combine(B->p[0], B->p[1], B->p[1], tmp, factor, B->NbColumns);
 207                 if (value_zero_p(B->p[1][pos[2]])) {
 208                     /* We will deal with this later */
 209                     assert(lex_sign(B->p[1], B->NbColumns) < 0);
 210                 }
 211             }
 212         }
 213     } else {
 214         int neg;
 215         int sign;
 216         bool progress = true;
 217         while (progress) {
 218             progress = false;
 219             for (int i = 0; i <= 1; ++i) {
 220                 value_set_si(factor, -1);
 221                 for (int j = 1; j <= 3; ++j) {
 222                     if (value_zero_p(B->p[1-i][pos[j]]))
 223                         continue;
 224                     value_oppose(tmp, B->p[i][pos[j]]);
 225                     value_pdivision(tmp, tmp, B->p[1-i][pos[j]]);
 226                     if (value_neg_p(factor) || value_lt(tmp, factor))
 227                         value_assign(factor, tmp);
 228                 }
 229                 if (value_pos_p(factor)) {
 230                     value_set_si(tmp, 1);
 231                     Vector_Combine(B->p[i], B->p[1-i], B->p[i], tmp, factor,
 232                                    B->NbColumns);
 233                     sign = lex_sign(B->p[i], B->NbColumns);
 234                     for (int j = 1; j <= 3; ++j) {
 235                         if (value_notzero_p(B->p[i][pos[j]]))
 236                             continue;
 237                         /* a zero is interpreted to be of sign sign */
 238                         if ((sign > 0 && value_pos_p(B->p[1-i][pos[j]])) ||
 239                             (sign < 0 && value_neg_p(B->p[1-i][pos[j]]))) {
 240                             /* the zero is of the wrong sign => back-off one */
 241                             value_set_si(tmp2, -1);
 242                             Vector_Combine(B->p[i], B->p[1-i], B->p[i], tmp, tmp2,
 243                                            B->NbColumns);
 244                             value_decrement(factor, factor);
 245                         }
 246                     }
 247                     /* We may have backed-off, so we need to check again. */
 248                     if (value_pos_p(factor))
 249                         progress = true;
 250                 }
 251             }
 252         }
 253         sign = 0;
 254         for (int i = 0; i < B->NbColumns; ++i) {
 255             value_addto(tmp, B->p[0][i], B->p[1][i]);
 256             if (value_zero_p(tmp))
 257                 continue;
 258             sign = value_neg_p(tmp) ? -1 : 1;
 259             break;
 260         }
 261         neg = 0;
 262         for (int i = 1; i <= 3; ++i) {
 263             value_addto(tmp, B->p[0][pos[i]], B->p[1][pos[i]]);
 264             if (value_neg_p(tmp) || (sign < 0 && value_zero_p(tmp)))
 265                 ++neg;
 266         }
 267         assert(neg <= 2);
 268         switch(neg) {
 269             int i;
 270         case 1:
 271             /* cases 4 and 5 in Theorem 11.1 */
 272             value_set_si(tmp, 1);
 273             Vector_Combine(B->p[0], B->p[1], B->p[1], tmp, tmp, B->NbColumns);
 274
 275             /* put positive pair in position 3 */
 276             for (i = 1; i <= 3; ++i)
 277                 if (value_pos_p(B->p[0][pos[i]]) && value_pos_p(B->p[1][pos[i]]))
 278                     break;
 279             assert(i <= 3);
 280             set_pos(pos, i, 3);
 281
 282             break;
 283         case 2:
 284             /* cases 1 and 2 in Theorem 11.1 */
 285             value_set_si(tmp, 1);
 286             Vector_Combine(B->p[0], B->p[1], B->p[0], tmp, tmp, B->NbColumns);
 287
 288             /* put positive one in position 2 */
 289             for (i = 1; i <= 3; ++i)
 290                 if (value_pos_p(B->p[0][pos[i]]))
 291                     break;
 292             assert(i <= 3);
 293             set_pos(pos, i, 2);
 294
 295             /* fourth constraint is redundant with respect to neighborhoods */
 296             *n = 3;
 297             break;
 298         case 0:
 299             /* We will deal with these later */
 300             assert(0);
 301         }
 302     }
 303
 304     value_clear(tmp);
 305     value_clear(tmp2);
 306     value_clear(factor);
 307
 308     return B;
 309 }
 310
 311 struct simplex {
 312     Value   last;       // last multiple of offset in link
 313     Vector *offset;
 314     Matrix *M;          // rows: elements different from (0,0)
 315     int     mask;
 316
 317     simplex(int d) {
 318         M = Matrix_Alloc(d, 2);
 319         offset = NULL;
 320     }
 321     simplex(int d, int mask, Value last) {
 322         M = Matrix_Alloc(d, 2);
 323         offset = Vector_Alloc(2);
 324         value_init(this->last);
 325         value_assign(this->last, last);
 326         this->mask = mask;
 327     }
 328     void transform(Matrix *T);
 329     void normalize();
 330     Polyhedron *shrunk_polyhedron(Polyhedron *P, int dim, Matrix *A,
 331                                     unsigned MaxRays);
 332     void print(FILE *out);
 333 };
 334
 335 void simplex::print(FILE *out)
 336 {
 337     if (!offset)
 338         Matrix_Print(out, P_VALUE_FMT, M);
 339     else {
 340         fprintf(out, "%d %d\n", M->NbRows, M->NbColumns);
 341         for (int j = 0; j < M->NbRows; ++j) {
 342             for (int k = 0; k < M->NbColumns; ++k)
 343                 value_print(out, P_VALUE_FMT, M->p[j][k]);
 344             if (mask & (1 << j)) {
 345                 fprintf(out, " + k * ");
 346                 for (int k = 0; k < M->NbColumns; ++k)
 347                     value_print(out, P_VALUE_FMT, offset->p[k]);
 348             }
 349             fprintf(out, "\n");
 350         }
 351         fprintf(out, "\t0 <= k <= ");
 352         value_print(out, P_VALUE_FMT, last);
 353         fprintf(out, "\n");
 354     }
 355 }
 356
 357 static bool lex_smaller(Value *v1, Value *v2, int n)
 358 {
 359     for (int i = 0; i < n; ++i)
 360         if (value_lt(v1[i], v2[i]))
 361             return true;
 362         else if (value_gt(v1[i], v2[i]))
 363             return false;
 364     return false;
 365 }
 366
 367 void simplex::transform(Matrix *T)
 368 {
 369     Matrix *M2 = M;
 370     M = Matrix_Alloc(M2->NbRows, M2->NbColumns);
 371     Matrix_Product(M2, T, M);
 372     Matrix_Free(M2);
 373
 374     if (offset) {
 375         Vector *offset2 = offset;
 376         offset = Vector_Alloc(offset2->Size);
 377         Vector_Matrix_Product(offset2->p, T, offset->p);
 378         Vector_Free(offset2);
 379     }
 380 }
 381
 382 void simplex::normalize()
 383 {
 384     int lexmin = 0;
 385     for (int i = 1; i < M->NbRows; ++i)
 386         if (lex_smaller(M->p[i], M->p[lexmin], 2))
 387             lexmin = i;
 388     if (lex_sign(M->p[lexmin], 2) < 0) {
 389         Value tmp;
 390         value_init(tmp);
 391         value_set_si(tmp, -1);
 392         Vector_Scale(M->p[lexmin], M->p[lexmin], tmp, 2);
 393         value_set_si(tmp, 1);
 394         for (int i = 0; i < M->NbRows; ++i) {
 395             if (i == lexmin)
 396                 continue;
 397             Vector_Combine(M->p[lexmin], M->p[i], M->p[i], tmp, tmp, 2);
 398         }
 399         if (offset && (mask & (1 << lexmin))) {
 400             value_set_si(tmp, -1);
 401             Vector_Scale(offset->p, offset->p, tmp, 2);
 402             mask ^= (1 << M->NbRows) - 1 - (1 << lexmin);
 403         }
 404         value_clear(tmp);
 405     }
 406 }
 407
 408 Polyhedron *simplex::shrunk_polyhedron(Polyhedron *P, int dim, Matrix *A,
 409                                         unsigned MaxRays)
 410 {
 411     Matrix *Constraints, *b;
 412     Vector *b_offset = NULL;
 413     Polyhedron *Q;
 414     Value min, min_var, tmp;
 415     value_init(tmp);
 416     value_init(min);
 417     value_init(min_var);
 418     int constant;
 419
 420     b = Matrix_Alloc(M->NbRows, A->NbColumns);
 421     Matrix_Product(M, A, b);
 422
 423     if (offset) {
 424         b_offset = Vector_Alloc(A->NbColumns);
 425         Vector_Matrix_Product(offset->p, A, b_offset->p);
 426     }
 427
 428     if (!offset)
 429         Constraints = Polyhedron2Constraints(P);
 430     else {
 431         Constraints = Matrix_Alloc(P->NbConstraints+2, P->Dimension+2+1);
 432         for (int i = 0; i < P->NbConstraints; ++i) {
 433             Vector_Copy(P->Constraint[i], Constraints->p[i], 1+dim+2);
 434             Vector_Copy(P->Constraint[i]+1+dim+2, Constraints->p[i]+1+dim+2+1,
 435                         (P->Dimension+2)-(1+dim+2));
 436         }
 437         value_set_si(Constraints->p[P->NbConstraints][0], 1);
 438         value_set_si(Constraints->p[P->NbConstraints][1+dim+2], 1);
 439         value_set_si(Constraints->p[P->NbConstraints+1][0], 1);
 440         value_set_si(Constraints->p[P->NbConstraints+1][1+dim+2], -1);
 441         value_assign(Constraints->p[P->NbConstraints+1][Constraints->NbColumns-1],
 442                      last);
 443     }
 444     constant = Constraints->NbColumns - 1;
 445
 446     for (int i = 0, j = 0; i < P->NbConstraints; ++i) {
 447         if (value_zero_p(Constraints->p[i][1+dim]) &&
 448             value_zero_p(Constraints->p[i][1+dim+1]))
 449             continue;
 450         value_set_si(min, 0);
 451         for (int k = 0; k < b->NbRows; ++k) {
 452             if (offset && (mask & (1 << k)))
 453                 continue;
 454             if (value_lt(b->p[k][j], min))
 455                 value_assign(min, b->p[k][j]);
 456         }
 457         value_set_si(min_var, 0);
 458         if (offset) {
 459             if (value_neg_p(b_offset->p[j])) {
 460                 value_oppose(min_var, b_offset->p[j]);
 461                 value_multiply(min_var, min_var, last);
 462                 value_increment(min_var, min_var);
 463             }
 464             for (int k = 0; k < b->NbRows; ++k) {
 465                 if (!(mask & (1 << k)))
 466                     continue;
 467                 if (value_lt(b->p[k][j], min_var))
 468                     value_assign(min_var, b->p[k][j]);
 469             }
 470         }
 471         if (!offset || value_pos_p(b_offset->p[j])) {
 472             if (value_le(min, min_var))
 473                 value_addto(Constraints->p[i][constant],
 474                             Constraints->p[i][constant], min);
 475             else {
 476                 value_assign(tmp, min_var);
 477                 value_addmul(tmp, last, b_offset->p[j]);
 478                 if (value_le(tmp, min)) {
 479                     value_addto(Constraints->p[i][constant],
 480                                 Constraints->p[i][constant], min_var);
 481                     value_addto(Constraints->p[i][1+dim+2],
 482                                 Constraints->p[i][1+dim+2], b_offset->p[j]);
 483                 } else {
 484                     int lastrow = Constraints->NbRows;
 485                     int cols = Constraints->NbColumns;
 486                     Matrix *C = Constraints;
 487                     Constraints = AddANullRow(Constraints);
 488                     Matrix_Free(C);
 489                     Vector_Copy(Constraints->p[i], Constraints->p[lastrow], cols);
 490                     value_addto(Constraints->p[i][constant],
 491                                 Constraints->p[i][constant], min_var);
 492                     value_addto(Constraints->p[i][1+dim+2],
 493                                 Constraints->p[i][1+dim+2], b_offset->p[j]);
 494                     value_addto(Constraints->p[lastrow][constant],
 495                                 Constraints->p[lastrow][constant], min);
 496                 }
 497             }
 498         } else {
 499             if (value_le(min_var, min)) {
 500                 value_addto(Constraints->p[i][constant],
 501                             Constraints->p[i][constant], min_var);
 502                 value_addto(Constraints->p[i][1+dim+2],
 503                             Constraints->p[i][1+dim+2], b_offset->p[j]);
 504             } else {
 505                 value_assign(tmp, min_var);
 506                 value_addmul(tmp, last, b_offset->p[j]);
 507                 if (value_le(min, tmp)) {
 508                     value_addto(Constraints->p[i][constant],
 509                                 Constraints->p[i][constant], min);
 510                 } else {
 511                     int lastrow = Constraints->NbRows;
 512                     int cols = Constraints->NbColumns;
 513                     Matrix *C = Constraints;
 514                     Constraints = AddANullRow(Constraints);
 515                     Matrix_Free(C);
 516                     Vector_Copy(Constraints->p[i], Constraints->p[lastrow], cols);
 517                     value_addto(Constraints->p[i][constant],
 518                                 Constraints->p[i][constant], min_var);
 519                     value_addto(Constraints->p[i][1+dim+2],
 520                                 Constraints->p[i][1+dim+2], b_offset->p[j]);
 521                     value_addto(Constraints->p[lastrow][constant],
 522                                 Constraints->p[lastrow][constant], min);
 523                 }
 524             }
 525         }
 526         ++j;
 527     }
 528     Q = Constraints2Polyhedron(Constraints, MaxRays);
 529
 530     if (b_offset)
 531         Vector_Free(b_offset);
 532     Matrix_Free(b);
 533     Matrix_Free(Constraints);
 534     value_clear(tmp);
 535     value_clear(min);
 536     value_clear(min_var);
 537
 538     return Q;
 539 }
 540
 541 struct scarf_complex {
 542     vector<simplex> simplices;
 543     void add(Matrix *B, int pos[4], int n);
 544     void add(Matrix *T, simplex s);
 545     void print(FILE *out);
 546     ~scarf_complex() {
 547         for (int i = 0; i < simplices.size(); ++i) {
 548             Matrix_Free(simplices[i].M);
 549             if (simplices[i].offset) {
 550                 Vector_Free(simplices[i].offset);
 551                 value_clear(simplices[i].last);
 552             }
 553         }
 554     }
 555 };
 556
 557 void scarf_complex::add(Matrix *T, simplex s)
 558 {
 559     s.transform(T);
 560     s.normalize();
 561     if (s.offset && lex_sign(s.offset->p, 2) < 0) {
 562         Value factor;
 563         Value tmp;
 564         value_init(factor);
 565         value_init(tmp);
 566         /* compute the smallest multiple (factor) of the offset that
 567          * makes on of the vertices lexico-negative.
 568          */
 569         int lexmin = -1;
 570         for (int i = 0; i < s.M->NbRows; ++i) {
 571             if (!(s.mask & (1 << i)))
 572                 continue;
 573             if (lexmin == -1 || lex_smaller(s.M->p[i], s.M->p[lexmin], 2))
 574                 lexmin = i;
 575         }
 576         if (value_zero_p(s.offset->p[0])) {
 577             if (value_pos_p(s.M->p[lexmin][0]))
 578                 value_increment(factor, s.last);
 579             else {
 580                 value_oppose(factor, s.M->p[lexmin][1]);
 581                 mpz_cdiv_q(factor, factor, s.offset->p[1]);
 582             }
 583         } else {
 584             value_oppose(factor, s.M->p[lexmin][0]);
 585             mpz_cdiv_q(factor, factor, s.offset->p[0]);
 586             if (mpz_divisible_p(s.M->p[lexmin][0], s.offset->p[0])) {
 587                 value_assign(tmp, s.M->p[lexmin][1]);
 588                 value_addmul(tmp, factor, s.offset->p[1]);
 589                 if (value_pos_p(tmp))
 590                     value_increment(factor, factor);
 591             }
 592         }
 593         if (value_le(factor, s.last)) {
 594             simplex part(s.M->NbRows, s.mask, s.last);
 595             Vector_Copy(s.offset->p, part.offset->p, 2);
 596             value_set_si(tmp, 1);
 597             for (int i = 0; i < s.M->NbRows; ++i) {
 598                 if (s.mask & (1 << i))
 599                     Vector_Combine(s.M->p[i], s.offset->p, part.M->p[i],
 600                                    tmp, factor, 2);
 601                 else
 602                     Vector_Copy(s.M->p[i], part.M->p[i], 2);
 603             }
 604             value_subtract(part.last, part.last, factor);
 605             value_decrement(s.last, factor);
 606             part.normalize();
 607             simplices.push_back(part);
 608         }
 609         value_clear(tmp);
 610         value_clear(factor);
 611     }
 612     simplices.push_back(s);
 613 }
 614
 615 void scarf_complex::add(Matrix *B, int pos[4], int n)
 616 {
 617     Matrix *T;
 618
 619     T = Matrix_Alloc(2, 2);
 620     Vector_Copy(B->p[0]+B->NbColumns-2, T->p[0], 2);
 621     Vector_Copy(B->p[1]+B->NbColumns-2, T->p[1], 2);
 622
 623     if (n == 3 || value_neg_p(B->p[0][pos[3]])) {
 624         assert(n == 3 || value_neg_p(B->p[1][pos[3]]));
 625
 626         simplex s1(1);
 627         value_set_si(s1.M->p[0][0], 0);
 628         value_set_si(s1.M->p[0][1], 1);
 629         add(T, s1);
 630
 631         simplex s2(1);
 632         value_set_si(s2.M->p[0][0], 1);
 633         value_set_si(s2.M->p[0][1], 1);
 634         add(T, s2);
 635
 636         simplex s3(1);
 637         value_set_si(s3.M->p[0][0], 1);
 638         value_set_si(s3.M->p[0][1], 0);
 639         add(T, s3);
 640
 641         simplex s4(2);
 642         value_set_si(s4.M->p[0][0], 0);
 643         value_set_si(s4.M->p[0][1], 1);
 644         value_set_si(s4.M->p[1][0], 1);
 645         value_set_si(s4.M->p[1][1], 1);
 646         add(T, s4);
 647
 648         simplex s5(2);
 649         value_set_si(s5.M->p[0][0], 1);
 650         value_set_si(s5.M->p[0][1], 0);
 651         value_set_si(s5.M->p[1][0], 1);
 652         value_set_si(s5.M->p[1][1], 1);
 653         add(T, s5);
 654     } else {
 655         Matrix *h;
 656         Vector *offset;
 657         bool initial = true;
 658         bool progress = true;
 659         Value tmp, tmp2, factor;
 660         int sign;
 661
 662         value_init(tmp);
 663         value_init(tmp2);
 664         value_init(factor);
 665
 666         assert(value_pos_p(B->p[0][pos[3]]));
 667         assert(value_pos_p(B->p[1][pos[3]]));
 668
 669         h = Matrix_Alloc(3, 2);
 670         value_set_si(h->p[0][0], 1);
 671         value_set_si(h->p[0][1], 0);
 672         value_set_si(h->p[1][0], 0);
 673         value_set_si(h->p[1][1], 1);
 674         value_set_si(h->p[2][0], 1);
 675         value_set_si(h->p[2][1], 1);
 676
 677         offset = Vector_Alloc(2);
 678
 679         while (progress) {
 680             progress = false;
 681             for (int i = 0; i <= 1; ++i) {
 682                 value_set_si(factor, -1);
 683                 for (int j = 1; j <= 2; ++j) {
 684                     if (value_zero_p(B->p[1-i][pos[j]]))
 685                         continue;
 686                     value_oppose(tmp, B->p[i][pos[j]]);
 687                     value_pdivision(tmp, tmp, B->p[1-i][pos[j]]);
 688                     if (value_neg_p(factor) || value_lt(tmp, factor))
 689                         value_assign(factor, tmp);
 690                 }
 691                 if (value_pos_p(factor)) {
 692                     value_set_si(tmp, 1);
 693                     Vector_Combine(B->p[i], B->p[1-i], B->p[i], tmp, factor,
 694                                    B->NbColumns);
 695                     sign = lex_sign(B->p[i], B->NbColumns);
 696                     for (int j = 1; j <= 2; ++j) {
 697                         if (value_notzero_p(B->p[i][pos[j]]))
 698                             continue;
 699                         /* a zero is interpreted to be of sign sign */
 700                         if ((sign > 0 && value_pos_p(B->p[1-i][pos[j]])) ||
 701                             (sign < 0 && value_neg_p(B->p[1-i][pos[j]]))) {
 702                             /* the zero is of the wrong sign => back-off one */
 703                             value_set_si(tmp2, -1);
 704                             Vector_Combine(B->p[i], B->p[1-i], B->p[i], tmp, tmp2,
 705                                            B->NbColumns);
 706                             value_decrement(factor, factor);
 707                         }
 708                     }
 709                     /* We may have backed-off, so we need to check again. */
 710                     if (value_pos_p(factor)) {
 711                         progress = true;
 712                         value_set_si(tmp, 1);
 713                         value_set_si(tmp2, -1);
 714
 715                         Vector_Combine(h->p[2], h->p[i], offset->p, tmp, tmp2, 2);
 716
 717                         if (initial) {
 718                             /* the initial simplices not in any link */
 719                             simplex l1(1);
 720                             Vector_Copy(h->p[0], l1.M->p[0], 2);
 721                             add(T, l1);
 722
 723                             simplex l2(1);
 724                             Vector_Copy(h->p[1], l2.M->p[0], 2);
 725                             add(T, l2);
 726
 727                             simplex l3(1);
 728                             Vector_Combine(h->p[0], h->p[1], l3.M->p[0],
 729                                            tmp, tmp2, 2);
 730                             add(T, l3);
 731
 732                             simplex t1(2);
 733                             Vector_Copy(h->p[0], t1.M->p[0], 2);
 734                             Vector_Copy(h->p[1], t1.M->p[1], 2);
 735                             add(T, t1);
 736
 737                             simplex t2(2);
 738                             Vector_Combine(h->p[0], h->p[1], t2.M->p[0],
 739                                            tmp, tmp2, 2);
 740                             Vector_Combine(h->p[2], h->p[1], t2.M->p[1],
 741                                            tmp, tmp2, 2);
 742                             add(T, t2);
 743                         } else {
 744                             /* update h */
 745                             Vector_Combine(h->p[i], offset->p, h->p[i],
 746                                            tmp, tmp, 2);
 747                             Vector_Combine(h->p[2], offset->p, h->p[2],
 748                                            tmp, tmp, 2);
 749                             value_decrement(factor, factor);
 750                         }
 751
 752                         simplex q(3, 0x4 | (1 << i), factor);
 753                         Vector_Copy(h->p[0], q.M->p[0], 2);
 754                         Vector_Copy(h->p[1], q.M->p[1], 2);
 755                         Vector_Copy(h->p[2], q.M->p[2], 2);
 756                         Vector_Copy(offset->p, q.offset->p, 2);
 757                         add(T, q);
 758
 759                         simplex t1(2, 0x3, factor);
 760                         Vector_Copy(h->p[i], t1.M->p[0], 2);
 761                         Vector_Copy(h->p[2], t1.M->p[1], 2);
 762                         Vector_Copy(offset->p, t1.offset->p, 2);
 763                         add(T, t1);
 764
 765                         simplex t2(2, 0x2, factor);
 766                         Vector_Copy(h->p[1-i], t2.M->p[0], 2);
 767                         Vector_Copy(h->p[2], t2.M->p[1], 2);
 768                         Vector_Copy(offset->p, t2.offset->p, 2);
 769                         add(T, t2);
 770
 771                         simplex l(1, 0x1, factor);
 772                         Vector_Copy(h->p[2], l.M->p[0], 2);
 773                         Vector_Copy(offset->p, l.offset->p, 2);
 774                         add(T, l);
 775
 776                         /* update h */
 777                         Vector_Combine(h->p[i], offset->p, h->p[i], tmp, factor, 2);
 778                         Vector_Combine(h->p[2], offset->p, h->p[2], tmp, factor, 2);
 779
 780                         initial = false;
 781                     }
 782                 }
 783             }
 784         }
 785         if (initial) {
 786             /* the initial simplices not in any link */
 787             simplex l1(1);
 788             Vector_Copy(h->p[0], l1.M->p[0], 2);
 789             add(T, l1);
 790
 791             simplex l2(1);
 792             Vector_Copy(h->p[1], l2.M->p[0], 2);
 793             add(T, l2);
 794
 795             simplex l3(1);
 796             Vector_Combine(h->p[0], h->p[1], l3.M->p[0],
 797                            tmp, tmp2, 2);
 798             add(T, l3);
 799
 800             simplex t1(2);
 801             Vector_Copy(h->p[0], t1.M->p[0], 2);
 802             Vector_Copy(h->p[1], t1.M->p[1], 2);
 803             add(T, t1);
 804
 805             simplex t2(2);
 806             Vector_Combine(h->p[0], h->p[1], t2.M->p[0],
 807                            tmp, tmp2, 2);
 808             Vector_Combine(h->p[2], h->p[1], t2.M->p[1],
 809                            tmp, tmp2, 2);
 810             add(T, t2);
 811
 812             {
 813                 /* the simplices in a link, here of length 1 */
 814                 simplex q(3);
 815                 Vector_Copy(h->p[0], q.M->p[0], 2);
 816                 Vector_Copy(h->p[1], q.M->p[1], 2);
 817                 Vector_Copy(h->p[2], q.M->p[2], 2);
 818                 add(T, q);
 819
 820                 simplex t1(2);
 821                 Vector_Copy(h->p[0], t1.M->p[0], 2);
 822                 Vector_Copy(h->p[2], t1.M->p[1], 2);
 823                 add(T, t1);
 824
 825                 simplex t2(2);
 826                 Vector_Copy(h->p[1], t2.M->p[0], 2);
 827                 Vector_Copy(h->p[2], t2.M->p[1], 2);
 828                 add(T, t2);
 829
 830                 simplex l(1);
 831                 Vector_Copy(h->p[2], l.M->p[0], 2);
 832                 add(T, l);
 833             }
 834         }
 835
 836         Vector_Free(offset);
 837         Matrix_Free(h);
 838
 839         value_clear(tmp);
 840         value_clear(tmp2);
 841         value_clear(factor);
 842     }
 843
 844     Matrix_Free(T);
 845 }
 846
 847 void scarf_complex::print(FILE *out)
 848 {
 849     for (int i = 0; i < simplices.size(); ++i)
 850         simplices[i].print(out);
 851 }
 852
 853 struct scarf_collector {
 854     virtual void add(Polyhedron *P, int sign, Polyhedron *C,
 855                      barvinok_options *options) = 0;
 856     virtual ~scarf_collector() {}
 857 };
 858
 859 static void scarf(Polyhedron *P, unsigned exist, unsigned nparam,
 860                   barvinok_options *options, scarf_collector& col)
 861 {
 862     Matrix *A, *B;
 863     int dim = P->Dimension - exist - nparam;
 864     assert(exist == 2);
 865     int pos[4];
 866     Polyhedron *U;
 867     gen_fun *gf;
 868     int n;
 869
 870     A = extract_matrix(P, dim);
 871
 872     n = A->NbColumns - 2;
 873     assert(n >= 3 && n <= 4);
 874
 875     int l = 0;
 876     for (int i = 0; i < n; ++i) {
 877         int j;
 878         for (j = 0; j < l; ++j)
 879             if (value_eq(A->p[0][pos[j]], A->p[0][i]) &&
 880                 value_eq(A->p[1][pos[j]], A->p[1][i]))
 881                 break;
 882         if (j < l)
 883             continue;
 884         pos[l++] = i;
 885     }
 886
 887     assert(l >= 3 && l <= 4);
 888     B = normalize_matrix(A, pos, &l);
 889
 890     scarf_complex scarf;
 891     scarf.add(B, pos, l);
 892
 893     U = Universe_Polyhedron(nparam);
 894     col.add(P, 0, U, options);
 895     for (int i = 0; i < scarf.simplices.size(); ++i) {
 896         Polyhedron *Q;
 897         int sign = (scarf.simplices[i].M->NbRows % 2) ? -1 : 1;
 898         Q = scarf.simplices[i].shrunk_polyhedron(P, dim, A, options->MaxRays);
 899         col.add(Q, sign, U, options);
 900         Polyhedron_Free(Q);
 901     }
 902     Polyhedron_Free(U);
 903
 904     Matrix_Free(B);
 905
 906     Matrix_Free(A);
 907 }
 908
 909 struct scarf_collector_gf : public scarf_collector {
 910     QQ c;
 911     gen_fun *gf;
 912
 913     scarf_collector_gf() {
 914         c.d = 1;
 915     }
 916     virtual void add(Polyhedron *P, int sign, Polyhedron *C,
 917                      barvinok_options *options);
 918 };
 919
 920 void scarf_collector_gf::add(Polyhedron *P, int sign, Polyhedron *C,
 921                              barvinok_options *options)
 922 {
 923     if (!sign)
 924         gf = barvinok_series_with_options(P, C, options);
 925     else {
 926         gen_fun *gf2;
 927         c.n = sign;
 928         gf2 = barvinok_series_with_options(P, C, options);
 929         gf->add(c, gf2, options);
 930         delete gf2;
 931     }
 932 }
 933
 934 gen_fun *barvinok_enumerate_scarf_series(Polyhedron *P,
 935                           unsigned exist, unsigned nparam, barvinok_options *options)
 936 {
 937     scarf_collector_gf scgf;
 938     scarf(P, exist, nparam, options, scgf);
 939     return scgf.gf;
 940 }
 941
 942 struct scarf_collector_ev : public scarf_collector {
 943     evalue mone;
 944     evalue *EP;
 945
 946     scarf_collector_ev() {
 947         value_init(mone.d);
 948         evalue_set_si(&mone, -1, 1);
 949     }
 950     ~scarf_collector_ev() {
 951         free_evalue_refs(&mone);
 952     }
 953     virtual void add(Polyhedron *P, int sign, Polyhedron *C,
 954                      barvinok_options *options);
 955 };
 956
 957 void scarf_collector_ev::add(Polyhedron *P, int sign, Polyhedron *C,
 958                              barvinok_options *options)
 959 {
 960     if (!sign)
 961         EP = barvinok_enumerate_with_options(P, C, options);
 962     else {
 963         evalue *E2;
 964         E2 = barvinok_enumerate_with_options(P, C, options);
 965         if (sign < 0)
 966             emul(&mone, E2);
 967         eadd(E2, EP);
 968         evalue_free(E2);
 969     }
 970 }
 971
 972 evalue *barvinok_enumerate_scarf(Polyhedron *P,
 973                           unsigned exist, unsigned nparam, barvinok_options *options)
 974 {
 975     scarf_collector_ev scev;
 976     scarf(P, exist, nparam, options, scev);
 977     return scev.EP;
 978 }