scarf.cc

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