scarf.cc

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