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 static Matrix *InsertColumn(Matrix *M, int pos)
 409 {
 410     Matrix *R;
 411     int i;
 412
 413     R = Matrix_Alloc(M->NbRows, M->NbColumns+1);
 414     for (i = 0; i < M->NbRows; ++i) {
 415         Vector_Copy(M->p[i], R->p[i], pos);
 416         Vector_Copy(M->p[i]+pos, R->p[i]+pos+1, M->NbColumns-pos);
 417     }
 418     return R;
 419 }
 420
 421 Polyhedron *simplex::shrunk_polyhedron(Polyhedron *P, int dim, Matrix *A,
 422                                         unsigned MaxRays)
 423 {
 424     Matrix *Constraints, *b;
 425     Vector *b_offset = NULL;
 426     Polyhedron *Q;
 427     Value min, min_var, tmp;
 428     value_init(tmp);
 429     value_init(min);
 430     value_init(min_var);
 431     int constant;
 432
 433     b = Matrix_Alloc(M->NbRows, A->NbColumns);
 434     Matrix_Product(M, A, b);
 435
 436     if (offset) {
 437         b_offset = Vector_Alloc(A->NbColumns);
 438         Vector_Matrix_Product(offset->p, A, b_offset->p);
 439     }
 440
 441     if (!offset)
 442         Constraints = Polyhedron2Constraints(P);
 443     else {
 444         Constraints = Matrix_Alloc(P->NbConstraints+2, P->Dimension+2+1);
 445         for (int i = 0; i < P->NbConstraints; ++i) {
 446             Vector_Copy(P->Constraint[i], Constraints->p[i], 1+dim+2);
 447             Vector_Copy(P->Constraint[i]+1+dim+2, Constraints->p[i]+1+dim+2+1,
 448                         (P->Dimension+2)-(1+dim+2));
 449         }
 450         value_set_si(Constraints->p[P->NbConstraints][0], 1);
 451         value_set_si(Constraints->p[P->NbConstraints][1+dim+2], 1);
 452         value_set_si(Constraints->p[P->NbConstraints+1][0], 1);
 453         value_set_si(Constraints->p[P->NbConstraints+1][1+dim+2], -1);
 454         value_assign(Constraints->p[P->NbConstraints+1][Constraints->NbColumns-1],
 455                      last);
 456     }
 457     constant = Constraints->NbColumns - 1;
 458
 459     for (int i = 0, j = 0; i < P->NbConstraints; ++i) {
 460         if (value_zero_p(Constraints->p[i][1+dim]) &&
 461             value_zero_p(Constraints->p[i][1+dim+1]))
 462             continue;
 463         value_set_si(min, 0);
 464         for (int k = 0; k < b->NbRows; ++k) {
 465             if (offset && (mask & (1 << k)))
 466                 continue;
 467             if (value_lt(b->p[k][j], min))
 468                 value_assign(min, b->p[k][j]);
 469         }
 470         value_set_si(min_var, 0);
 471         if (offset) {
 472             if (value_neg_p(b_offset->p[j])) {
 473                 value_oppose(min_var, b_offset->p[j]);
 474                 value_multiply(min_var, min_var, last);
 475                 value_increment(min_var, min_var);
 476             }
 477             for (int k = 0; k < b->NbRows; ++k) {
 478                 if (!(mask & (1 << k)))
 479                     continue;
 480                 if (value_lt(b->p[k][j], min_var))
 481                     value_assign(min_var, b->p[k][j]);
 482             }
 483         }
 484         if (!offset || value_pos_p(b_offset->p[j])) {
 485             if (value_le(min, min_var))
 486                 value_addto(Constraints->p[i][constant],
 487                             Constraints->p[i][constant], min);
 488             else {
 489                 value_assign(tmp, min_var);
 490                 value_addmul(tmp, last, b_offset->p[j]);
 491                 if (value_le(tmp, min)) {
 492                     value_addto(Constraints->p[i][constant],
 493                                 Constraints->p[i][constant], min_var);
 494                     value_addto(Constraints->p[i][1+dim+2],
 495                                 Constraints->p[i][1+dim+2], b_offset->p[j]);
 496                 } else {
 497                     int lastrow = Constraints->NbRows;
 498                     int cols = Constraints->NbColumns;
 499                     Matrix *C = Constraints;
 500                     Constraints = AddANullRow(Constraints);
 501                     Matrix_Free(C);
 502                     Vector_Copy(Constraints->p[i], Constraints->p[lastrow], cols);
 503                     value_addto(Constraints->p[i][constant],
 504                                 Constraints->p[i][constant], min_var);
 505                     value_addto(Constraints->p[i][1+dim+2],
 506                                 Constraints->p[i][1+dim+2], b_offset->p[j]);
 507                     value_addto(Constraints->p[lastrow][constant],
 508                                 Constraints->p[lastrow][constant], min);
 509                 }
 510             }
 511         } else {
 512             if (value_le(min_var, min)) {
 513                 value_addto(Constraints->p[i][constant],
 514                             Constraints->p[i][constant], min_var);
 515                 value_addto(Constraints->p[i][1+dim+2],
 516                             Constraints->p[i][1+dim+2], b_offset->p[j]);
 517             } else {
 518                 value_assign(tmp, min_var);
 519                 value_addmul(tmp, last, b_offset->p[j]);
 520                 if (value_le(min, tmp)) {
 521                     value_addto(Constraints->p[i][constant],
 522                                 Constraints->p[i][constant], min);
 523                 } else {
 524                     int lastrow = Constraints->NbRows;
 525                     int cols = Constraints->NbColumns;
 526                     Matrix *C = Constraints;
 527                     Constraints = AddANullRow(Constraints);
 528                     Matrix_Free(C);
 529                     Vector_Copy(Constraints->p[i], Constraints->p[lastrow], cols);
 530                     value_addto(Constraints->p[i][constant],
 531                                 Constraints->p[i][constant], min_var);
 532                     value_addto(Constraints->p[i][1+dim+2],
 533                                 Constraints->p[i][1+dim+2], b_offset->p[j]);
 534                     value_addto(Constraints->p[lastrow][constant],
 535                                 Constraints->p[lastrow][constant], min);
 536                 }
 537             }
 538         }
 539         ++j;
 540     }
 541     Q = Constraints2Polyhedron(Constraints, MaxRays);
 542
 543     if (b_offset)
 544         Vector_Free(b_offset);
 545     Matrix_Free(b);
 546     Matrix_Free(Constraints);
 547     value_clear(tmp);
 548     value_clear(min);
 549     value_clear(min_var);
 550
 551     return Q;
 552 }
 553
 554 struct scarf_complex {
 555     vector<simplex> simplices;
 556     void add(Matrix *B, int pos[4], int n);
 557     void add(Matrix *T, simplex s);
 558     void print(FILE *out);
 559     ~scarf_complex() {
 560         for (int i = 0; i < simplices.size(); ++i) {
 561             Matrix_Free(simplices[i].M);
 562             if (simplices[i].offset) {
 563                 Vector_Free(simplices[i].offset);
 564                 value_clear(simplices[i].last);
 565             }
 566         }
 567     }
 568 };
 569
 570 void scarf_complex::add(Matrix *T, simplex s)
 571 {
 572     s.transform(T);
 573     s.normalize();
 574     if (s.offset && lex_sign(s.offset->p, 2) < 0) {
 575         Value factor;
 576         Value tmp;
 577         value_init(factor);
 578         value_init(tmp);
 579         /* compute the smallest multiple (factor) of the offset that
 580          * makes on of the vertices lexico-negative.
 581          */
 582         int lexmin = -1;
 583         for (int i = 0; i < s.M->NbRows; ++i) {
 584             if (!(s.mask & (1 << i)))
 585                 continue;
 586             if (lexmin == -1 || lex_smaller(s.M->p[i], s.M->p[lexmin], 2))
 587                 lexmin = i;
 588         }
 589         if (value_zero_p(s.offset->p[0])) {
 590             if (value_pos_p(s.M->p[lexmin][0]))
 591                 value_increment(factor, s.last);
 592             else {
 593                 value_oppose(factor, s.M->p[lexmin][1]);
 594                 mpz_cdiv_q(factor, factor, s.offset->p[1]);
 595             }
 596         } else {
 597             value_oppose(factor, s.M->p[lexmin][0]);
 598             mpz_cdiv_q(factor, factor, s.offset->p[0]);
 599             if (mpz_divisible_p(s.M->p[lexmin][0], s.offset->p[0])) {
 600                 value_assign(tmp, s.M->p[lexmin][1]);
 601                 value_addmul(tmp, factor, s.offset->p[1]);
 602                 if (value_pos_p(tmp))
 603                     value_increment(factor, factor);
 604             }
 605         }
 606         if (value_le(factor, s.last)) {
 607             simplex part(s.M->NbRows, s.mask, s.last);
 608             Vector_Copy(s.offset->p, part.offset->p, 2);
 609             value_set_si(tmp, 1);
 610             for (int i = 0; i < s.M->NbRows; ++i) {
 611                 if (s.mask & (1 << i))
 612                     Vector_Combine(s.M->p[i], s.offset->p, part.M->p[i],
 613                                    tmp, factor, 2);
 614                 else
 615                     Vector_Copy(s.M->p[i], part.M->p[i], 2);
 616             }
 617             value_subtract(part.last, part.last, factor);
 618             value_decrement(s.last, factor);
 619             part.normalize();
 620             simplices.push_back(part);
 621         }
 622         value_clear(tmp);
 623         value_clear(factor);
 624     }
 625     simplices.push_back(s);
 626 }
 627
 628 void scarf_complex::add(Matrix *B, int pos[4], int n)
 629 {
 630     Matrix *T;
 631
 632     T = Matrix_Alloc(2, 2);
 633     Vector_Copy(B->p[0]+B->NbColumns-2, T->p[0], 2);
 634     Vector_Copy(B->p[1]+B->NbColumns-2, T->p[1], 2);
 635
 636     if (n == 3 || value_neg_p(B->p[0][pos[3]])) {
 637         assert(n == 3 || value_neg_p(B->p[1][pos[3]]));
 638
 639         simplex s1(1);
 640         value_set_si(s1.M->p[0][0], 0);
 641         value_set_si(s1.M->p[0][1], 1);
 642         add(T, s1);
 643
 644         simplex s2(1);
 645         value_set_si(s2.M->p[0][0], 1);
 646         value_set_si(s2.M->p[0][1], 1);
 647         add(T, s2);
 648
 649         simplex s3(1);
 650         value_set_si(s3.M->p[0][0], 1);
 651         value_set_si(s3.M->p[0][1], 0);
 652         add(T, s3);
 653
 654         simplex s4(2);
 655         value_set_si(s4.M->p[0][0], 0);
 656         value_set_si(s4.M->p[0][1], 1);
 657         value_set_si(s4.M->p[1][0], 1);
 658         value_set_si(s4.M->p[1][1], 1);
 659         add(T, s4);
 660
 661         simplex s5(2);
 662         value_set_si(s5.M->p[0][0], 1);
 663         value_set_si(s5.M->p[0][1], 0);
 664         value_set_si(s5.M->p[1][0], 1);
 665         value_set_si(s5.M->p[1][1], 1);
 666         add(T, s5);
 667     } else {
 668         Matrix *h;
 669         Vector *offset;
 670         bool initial = true;
 671         bool progress = true;
 672         Value tmp, tmp2, factor;
 673         int sign;
 674
 675         value_init(tmp);
 676         value_init(tmp2);
 677         value_init(factor);
 678
 679         assert(value_pos_p(B->p[0][pos[3]]));
 680         assert(value_pos_p(B->p[1][pos[3]]));
 681
 682         h = Matrix_Alloc(3, 2);
 683         value_set_si(h->p[0][0], 1);
 684         value_set_si(h->p[0][1], 0);
 685         value_set_si(h->p[1][0], 0);
 686         value_set_si(h->p[1][1], 1);
 687         value_set_si(h->p[2][0], 1);
 688         value_set_si(h->p[2][1], 1);
 689
 690         offset = Vector_Alloc(2);
 691
 692         while (progress) {
 693             progress = false;
 694             for (int i = 0; i <= 1; ++i) {
 695                 value_set_si(factor, -1);
 696                 for (int j = 1; j <= 2; ++j) {
 697                     if (value_zero_p(B->p[1-i][pos[j]]))
 698                         continue;
 699                     value_oppose(tmp, B->p[i][pos[j]]);
 700                     value_pdivision(tmp, tmp, B->p[1-i][pos[j]]);
 701                     if (value_neg_p(factor) || value_lt(tmp, factor))
 702                         value_assign(factor, tmp);
 703                 }
 704                 if (value_pos_p(factor)) {
 705                     value_set_si(tmp, 1);
 706                     Vector_Combine(B->p[i], B->p[1-i], B->p[i], tmp, factor,
 707                                    B->NbColumns);
 708                     sign = lex_sign(B->p[i], B->NbColumns);
 709                     for (int j = 1; j <= 2; ++j) {
 710                         if (value_notzero_p(B->p[i][pos[j]]))
 711                             continue;
 712                         /* a zero is interpreted to be of sign sign */
 713                         if ((sign > 0 && value_pos_p(B->p[1-i][pos[j]])) ||
 714                             (sign < 0 && value_neg_p(B->p[1-i][pos[j]]))) {
 715                             /* the zero is of the wrong sign => back-off one */
 716                             value_set_si(tmp2, -1);
 717                             Vector_Combine(B->p[i], B->p[1-i], B->p[i], tmp, tmp2,
 718                                            B->NbColumns);
 719                             value_decrement(factor, factor);
 720                         }
 721                     }
 722                     /* We may have backed-off, so we need to check again. */
 723                     if (value_pos_p(factor)) {
 724                         progress = true;
 725                         value_set_si(tmp, 1);
 726                         value_set_si(tmp2, -1);
 727
 728                         Vector_Combine(h->p[2], h->p[i], offset->p, tmp, tmp2, 2);
 729
 730                         if (initial) {
 731                             /* the initial simplices not in any link */
 732                             simplex l1(1);
 733                             Vector_Copy(h->p[0], l1.M->p[0], 2);
 734                             add(T, l1);
 735
 736                             simplex l2(1);
 737                             Vector_Copy(h->p[1], l2.M->p[0], 2);
 738                             add(T, l2);
 739
 740                             simplex l3(1);
 741                             Vector_Combine(h->p[0], h->p[1], l3.M->p[0],
 742                                            tmp, tmp2, 2);
 743                             add(T, l3);
 744
 745                             simplex t1(2);
 746                             Vector_Copy(h->p[0], t1.M->p[0], 2);
 747                             Vector_Copy(h->p[1], t1.M->p[1], 2);
 748                             add(T, t1);
 749
 750                             simplex t2(2);
 751                             Vector_Combine(h->p[0], h->p[1], t2.M->p[0],
 752                                            tmp, tmp2, 2);
 753                             Vector_Combine(h->p[2], h->p[1], t2.M->p[1],
 754                                            tmp, tmp2, 2);
 755                             add(T, t2);
 756                         } else {
 757                             /* update h */
 758                             Vector_Combine(h->p[i], offset->p, h->p[i],
 759                                            tmp, tmp, 2);
 760                             Vector_Combine(h->p[2], offset->p, h->p[2],
 761                                            tmp, tmp, 2);
 762                             value_decrement(factor, factor);
 763                         }
 764
 765                         simplex q(3, 0x4 | (1 << i), factor);
 766                         Vector_Copy(h->p[0], q.M->p[0], 2);
 767                         Vector_Copy(h->p[1], q.M->p[1], 2);
 768                         Vector_Copy(h->p[2], q.M->p[2], 2);
 769                         Vector_Copy(offset->p, q.offset->p, 2);
 770                         add(T, q);
 771
 772                         simplex t1(2, 0x3, factor);
 773                         Vector_Copy(h->p[i], t1.M->p[0], 2);
 774                         Vector_Copy(h->p[2], t1.M->p[1], 2);
 775                         Vector_Copy(offset->p, t1.offset->p, 2);
 776                         add(T, t1);
 777
 778                         simplex t2(2, 0x2, factor);
 779                         Vector_Copy(h->p[1-i], t2.M->p[0], 2);
 780                         Vector_Copy(h->p[2], t2.M->p[1], 2);
 781                         Vector_Copy(offset->p, t2.offset->p, 2);
 782                         add(T, t2);
 783
 784                         simplex l(1, 0x1, factor);
 785                         Vector_Copy(h->p[2], l.M->p[0], 2);
 786                         Vector_Copy(offset->p, l.offset->p, 2);
 787                         add(T, l);
 788
 789                         /* update h */
 790                         Vector_Combine(h->p[i], offset->p, h->p[i], tmp, factor, 2);
 791                         Vector_Combine(h->p[2], offset->p, h->p[2], tmp, factor, 2);
 792
 793                         initial = false;
 794                     }
 795                 }
 796             }
 797         }
 798         if (initial) {
 799             /* the initial simplices not in any link */
 800             simplex l1(1);
 801             Vector_Copy(h->p[0], l1.M->p[0], 2);
 802             add(T, l1);
 803
 804             simplex l2(1);
 805             Vector_Copy(h->p[1], l2.M->p[0], 2);
 806             add(T, l2);
 807
 808             simplex l3(1);
 809             Vector_Combine(h->p[0], h->p[1], l3.M->p[0],
 810                            tmp, tmp2, 2);
 811             add(T, l3);
 812
 813             simplex t1(2);
 814             Vector_Copy(h->p[0], t1.M->p[0], 2);
 815             Vector_Copy(h->p[1], t1.M->p[1], 2);
 816             add(T, t1);
 817
 818             simplex t2(2);
 819             Vector_Combine(h->p[0], h->p[1], t2.M->p[0],
 820                            tmp, tmp2, 2);
 821             Vector_Combine(h->p[2], h->p[1], t2.M->p[1],
 822                            tmp, tmp2, 2);
 823             add(T, t2);
 824
 825             {
 826                 /* the simplices in a link, here of length 1 */
 827                 simplex q(3);
 828                 Vector_Copy(h->p[0], q.M->p[0], 2);
 829                 Vector_Copy(h->p[1], q.M->p[1], 2);
 830                 Vector_Copy(h->p[2], q.M->p[2], 2);
 831                 add(T, q);
 832
 833                 simplex t1(2);
 834                 Vector_Copy(h->p[0], t1.M->p[0], 2);
 835                 Vector_Copy(h->p[2], t1.M->p[1], 2);
 836                 add(T, t1);
 837
 838                 simplex t2(2);
 839                 Vector_Copy(h->p[1], t2.M->p[0], 2);
 840                 Vector_Copy(h->p[2], t2.M->p[1], 2);
 841                 add(T, t2);
 842
 843                 simplex l(1);
 844                 Vector_Copy(h->p[2], l.M->p[0], 2);
 845                 add(T, l);
 846             }
 847         }
 848
 849         Vector_Free(offset);
 850         Matrix_Free(h);
 851
 852         value_clear(tmp);
 853         value_clear(tmp2);
 854         value_clear(factor);
 855     }
 856
 857     Matrix_Free(T);
 858 }
 859
 860 void scarf_complex::print(FILE *out)
 861 {
 862     for (int i = 0; i < simplices.size(); ++i)
 863         simplices[i].print(out);
 864 }
 865
 866 struct scarf_collector {
 867     virtual void add(Polyhedron *P, int sign, Polyhedron *C,
 868                      barvinok_options *options) = 0;
 869     virtual ~scarf_collector() {}
 870 };
 871
 872 static void scarf(Polyhedron *P, unsigned exist, unsigned nparam,
 873                   barvinok_options *options, scarf_collector& col)
 874 {
 875     Matrix *A, *B;
 876     int dim = P->Dimension - exist - nparam;
 877     assert(exist == 2);
 878     int pos[4];
 879     Polyhedron *U;
 880     gen_fun *gf;
 881     int n;
 882
 883     A = extract_matrix(P, dim);
 884
 885     n = A->NbColumns - 2;
 886     assert(n >= 3 && n <= 4);
 887
 888     int l = 0;
 889     for (int i = 0; i < n; ++i) {
 890         int j;
 891         for (j = 0; j < l; ++j)
 892             if (value_eq(A->p[0][pos[j]], A->p[0][i]) &&
 893                 value_eq(A->p[1][pos[j]], A->p[1][i]))
 894                 break;
 895         if (j < l)
 896             continue;
 897         pos[l++] = i;
 898     }
 899
 900     assert(l >= 3 && l <= 4);
 901     B = normalize_matrix(A, pos, &l);
 902
 903     scarf_complex scarf;
 904     scarf.add(B, pos, l);
 905
 906     U = Universe_Polyhedron(nparam);
 907     col.add(P, 0, U, options);
 908     for (int i = 0; i < scarf.simplices.size(); ++i) {
 909         Polyhedron *Q;
 910         int sign = (scarf.simplices[i].M->NbRows % 2) ? -1 : 1;
 911         Q = scarf.simplices[i].shrunk_polyhedron(P, dim, A, options->MaxRays);
 912         col.add(Q, sign, U, options);
 913         Polyhedron_Free(Q);
 914     }
 915     Polyhedron_Free(U);
 916
 917     Matrix_Free(B);
 918
 919     Matrix_Free(A);
 920 }
 921
 922 struct scarf_collector_gf : public scarf_collector {
 923     QQ c;
 924     gen_fun *gf;
 925
 926     scarf_collector_gf() {
 927         c.d = 1;
 928     }
 929     virtual void add(Polyhedron *P, int sign, Polyhedron *C,
 930                      barvinok_options *options);
 931 };
 932
 933 void scarf_collector_gf::add(Polyhedron *P, int sign, Polyhedron *C,
 934                              barvinok_options *options)
 935 {
 936     if (!sign)
 937         gf = barvinok_series_with_options(P, C, options);
 938     else {
 939         gen_fun *gf2;
 940         c.n = sign;
 941         gf2 = barvinok_series_with_options(P, C, options);
 942         gf->add(c, gf2, options);
 943         delete gf2;
 944     }
 945 }
 946
 947 gen_fun *barvinok_enumerate_scarf_series(Polyhedron *P,
 948                           unsigned exist, unsigned nparam, barvinok_options *options)
 949 {
 950     scarf_collector_gf scgf;
 951     scarf(P, exist, nparam, options, scgf);
 952     return scgf.gf;
 953 }
 954
 955 struct scarf_collector_ev : public scarf_collector {
 956     evalue mone;
 957     evalue *EP;
 958
 959     scarf_collector_ev() {
 960         value_init(mone.d);
 961         evalue_set_si(&mone, -1, 1);
 962     }
 963     ~scarf_collector_ev() {
 964         free_evalue_refs(&mone);
 965     }
 966     virtual void add(Polyhedron *P, int sign, Polyhedron *C,
 967                      barvinok_options *options);
 968 };
 969
 970 void scarf_collector_ev::add(Polyhedron *P, int sign, Polyhedron *C,
 971                              barvinok_options *options)
 972 {
 973     if (!sign)
 974         EP = barvinok_enumerate_with_options(P, C, options);
 975     else {
 976         evalue *E2;
 977         E2 = barvinok_enumerate_with_options(P, C, options);
 978         if (sign < 0)
 979             emul(&mone, E2);
 980         eadd(E2, EP);
 981         evalue_free(E2);
 982     }
 983 }
 984
 985 evalue *barvinok_enumerate_scarf(Polyhedron *P,
 986                           unsigned exist, unsigned nparam, barvinok_options *options)
 987 {
 988     scarf_collector_ev scev;
 989     scarf(P, exist, nparam, options, scev);
 990     return scev.EP;
 991 }