src/baseparsuif/suifmath/named_sc_fm.cc

   1 /* file "named_sc_fm.cc" */
   2
   3 /*  Copyright (c) 1994 Stanford University
   4
   5     All rights reserved.
   6
   7     This software is provided under the terms described in
   8     the "suif_copyright.h" include file. */
   9
  10 #include <suif_copyright.h>
  11
  12 #define _MODULE_ "libsuifmath.a"
  13 #pragma implementation "named_sc_fm.h"
  14
  15 #include <stdio.h>
  16 #include <stdarg.h>
  17 #include <setjmp.h>
  18 #include <suif1.h>
  19 #if 0
  20 #include <sys/times.h>
  21 #endif
  22 #include <builder.h>
  23 #include "suifmath.h"
  24
  25
  26 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  27  * Print debug information                                             *
  28  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
  29 //#define DEBUG_PRINT
  30 //#define DEBUG_PRINT2
  31 //#define DEBUG_PRINT3
  32
  33 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  34  * turn statistics gathering and printing on                           *
  35  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
  36 //#define STAT_ON
  37
  38 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  39  * When watch is turned-on, for large problems that run forever,       *
  40  * incremental status messages will be displayed.                      *
  41  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
  42 //#define WATCH_ON
  43
  44 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  45  * Use hand-tuned code in computationally intensive areas.             *
  46  * These were developed after profiling and obsurving the              *
  47  * bottlenecks.                                                        *
  48  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
  49 #define USE_SPEED_HACK
  50
  51
  52 #define DIM_M_ST        64
  53 #define DIM_M_OFF       16
  54
  55
  56
  57
  58 int named_sc_fm::cm = 0;
  59 int named_sc_fm::cn = 0;
  60 int named_sc_fm::cp = 0;
  61 sc_fm_constraint ** named_sc_fm::L = NULL;
  62 unsigned * named_sc_fm::magic_list = NULL;
  63
  64 jmp_buf * _named_sc_fm_longjmp_env = NULL;
  65
  66
  67 char * get_name(name_table_entry & nte);
  68
  69 #ifdef STAT_ON
  70 unsigned stat_fc_cons                = 0;
  71 unsigned stat_fc_recalc              = 0;
  72 unsigned stat_fc_div_by_gcd          = 0;
  73 unsigned stat_fc_inverse             = 0;
  74 unsigned stat_fc_is_only_pos_const   = 0;
  75 unsigned stat_fc_is_coef_eq_1        = 0;
  76 unsigned stat_fr_cons                = 0;
  77 unsigned stat_fr_init                = 0;
  78 unsigned stat_fr_set                 = 0;
  79 unsigned stat_fm_cons                = 0;
  80 unsigned stat_fm_extend_block        = 0;
  81 unsigned stat_fm_step0               = 0;
  82 unsigned stat_fm_step1               = 0;
  83 unsigned stat_fm_step2               = 0;
  84 unsigned stat_fm_step3               = 0;
  85 unsigned stat_fm_step4               = 0;
  86 unsigned stat_fm_lin_step            = 0;
  87 unsigned stat_fm_bounds0             = 0;
  88 unsigned stat_fm_bounds1             = 0;
  89 unsigned stat_fm_bounds2             = 0;
  90 unsigned stat_fm_bounds3             = 0;
  91 unsigned stat_fm_is_already0         = 0;
  92 unsigned long stat_fm_is_already0x   = 0;
  93 unsigned stat_fm_is_already1         = 0;
  94 unsigned stat_fm_is_already2         = 0;
  95 unsigned stat_fm_sort                = 0;
  96 unsigned stat_fm_swap                = 0;
  97 unsigned stat_fm_remove              = 0;
  98 unsigned stat_fm_get                 = 0;
  99 unsigned stat_fm_check_valid         = 0;
 100
 101 static void print_stat();
 102 #endif
 103
 104 #ifdef STAT_ON
 105 #define CALL_STAT(NM)   NM++;
 106 #else
 107 #define CALL_STAT(NM)
 108 #endif
 109
 110
 111
 112
 113
 114 #if 0
 115 static long get_time()
 116 {
 117     char buff[512];
 118     long t = times((tms *)buff);
 119     return t;
 120 }
 121 #endif
 122
 123
 124
 125
 126 /***************************************************************************
 127  * Find a good value to set the num ineq dimension.                        *
 128  *                                                                         *
 129  *   x - wanted, pr - old value,  st - minimum size, off - incriment size  *
 130  ***************************************************************************/
 131 static int next_up(int x, int pr, int st, int off)
 132 {
 133     // Dont make the block very small compared to the last block.
 134     if(pr - off > x) return pr - off;
 135     int v = st;
 136     while(v < x) v += off;
 137     return v;
 138 }
 139
 140
 141 /***************************************************************************
 142  *   given the number of rows columns and planes, guesstimate a value for  *
 143  *   maximum number of inequalities that will result in doing fm-solve.    *
 144  ***************************************************************************/
 145 static int nrows_in_blk(int p, int r, int c)
 146 {
 147     return MAX(p*r*c, r*4);
 148 }
 149
 150
 151 /***************************************************************************
 152  *  return the sign of all the non-zero elements of a constraint.  If all  *
 153  *  elements were zero or elements have different signs return 0.          *
 154  ***************************************************************************/
 155 static int sign_constraint(constraint & c) {
 156     int sgn = 0;
 157     for(int i=0; i<c.n(); i++)
 158         if(c[i]) {
 159             if(sgn == 0)
 160                 sgn = (c[i]>0)?1:-1;
 161             else if(sgn*c[i]<0)
 162                 return 0;
 163         }
 164     return sgn;
 165 }
 166
 167
 168 /***************************************************************************
 169  *  Is a + b == 0 ?                                                        *
 170  ***************************************************************************/
 171 static boolean is_add_eq_zero(sc_fm_constraint &a, sc_fm_constraint & b)
 172 {
 173     assert(a.n() == b.n());
 174     assert(a.p() == b.p());
 175     for(int ip = 0; ip<a.p(); ip++)
 176         for(int in = 0; in<a.n(); in++)
 177             if(a[ip][in] + b[ip][in]) return FALSE;
 178     return TRUE;
 179 }
 180
 181
 182 #ifdef DEBUG_PRINT
 183 static boolean chk_is_in(int ix,
 184                          named_symcoeff_ineq & all,
 185                          named_symcoeff_ineq & curr)
 186 {
 187     assert(all.n() == curr.n());
 188     assert(all.p() == curr.p());
 189     assert((ix>=0)&&(ix<all.m()));
 190     for(int im=0; im<curr.m(); im++) {
 191         int found = TRUE;
 192         for(int ip=0; (ip<curr.p())&&(found); ip++)
 193             for(int in=0; in<curr.n(); in++)
 194                 if(all[ip][ix][in] != curr[ip][im][in]) found = FALSE;
 195         if(found) return TRUE;
 196     }
 197     return FALSE;
 198 }
 199 #endif
 200
 201
 202 /* ##################################################
 203    #####   sc_fm_constraint                    #####
 204    ################################################## */
 205
 206
 207 /***************************************************************************
 208  *  create a fm_constraint. Number of planes and columns are given         *
 209  ***************************************************************************/
 210 sc_fm_constraint::sc_fm_constraint(int in, int ip)
 211 {
 212     CALL_STAT(stat_fc_cons);
 213     nn = in;
 214     pp = ip;
 215     assert((in>=0)&&(ip>=0));
 216     L = (in*ip)?(new int[in*ip]):NULL;
 217     sgn = (in)?(new int[in]):NULL;
 218     prnk = -1;
 219     crnk = -1;
 220     magic_num = 0;
 221 }
 222
 223
 224 /***************************************************************************
 225  *   Delete                                                                *
 226  ***************************************************************************/
 227 sc_fm_constraint::~sc_fm_constraint()
 228 {
 229     if(L)   delete[] L;
 230     if(sgn) delete[] sgn;
 231 }
 232
 233
 234 /***************************************************************************
 235  *  return sign of fm_constraint.                                          *
 236  ***************************************************************************/
 237 int sc_fm_constraint::sign(int in)
 238 {
 239     // jump out of solver if sign is bad
 240     if ((sgn[in] == BAD_SIGN) && _named_sc_fm_longjmp_env)
 241         longjmp(*_named_sc_fm_longjmp_env, 1);
 242
 243     assert(sgn[in]!=BAD_SIGN);
 244     return sgn[in];
 245 }
 246
 247
 248 /***************************************************************************
 249  *  Need to call when constraint may have changed.                         *
 250  *  Recalculate the cached results.                                        *
 251  ***************************************************************************/
 252 void sc_fm_constraint::recalc()
 253 {
 254     CALL_STAT(stat_fc_recalc);
 255     crnk = -1;
 256     prnk = -1;
 257     magic_num = 0;
 258
 259     for(int i=0; i<n(); i++)
 260         recalc(i);
 261 }
 262
 263
 264 /***************************************************************************
 265  *  called for each column, Will find the                                  *
 266  *    1) sign of the column.                                               *
 267  *    2) contribute to calculating the maximum rank of the columns         *
 268  *    3) contribute to calculating the maximum rank of the planes          *
 269  *    4) contribute to calculating the unique magic number                 *
 270  ***************************************************************************/
 271 void sc_fm_constraint::recalc(int in)
 272 {
 273     int * s = &sgn[in];
 274     *s = 0;
 275     boolean bad = FALSE;
 276     for(int ip=0; ip<p(); ip++) {
 277         int c = (*this)[ip][in];
 278
 279         // this simple function will 'somewhat' garuntee that two
 280         // different systems will have two different magic-numbers
 281         // most of the time.
 282         // note that the constant term is disregarded when calculating
 283         // the magic_num.
 284         // magic_num ^= (in>0)*c*(41*ip + 31*in + 11*(ip+7)*(in+41) +
 285         //                        (0x1101 << (MAX(0,18-ip))) +
 286         //                        (0x1001 << (MAX(0,20-in))) +
 287         //                        c*(ip+5)*(in+31) +
 288         //                        7*c*c);
 289
 290         if(c > 0) {
 291             magic_num ^= (in>0)*c*(41*ip + 31*in + 11*(ip+7)*(in+41) +
 292                                    (0x1101 << (MAX(0,18-ip))) +
 293                                    (0x1001 << (MAX(0,20-in))) +
 294                                    c*(ip+5)*(in+31) +
 295                                    7*c*c);
 296             if(*s == -1) bad = TRUE;
 297             *s = 1;
 298             prnk = MAX(prnk, ip);
 299             crnk = in;
 300         } else if(c < 0) {
 301             magic_num ^= (in>0)*c*(41*ip + 31*in + 11*(ip+7)*(in+41) +
 302                                    (0x1101 << (MAX(0,18-ip))) +
 303                                    (0x1001 << (MAX(0,20-in))) +
 304                                    c*(ip+5)*(in+31) +
 305                                    7*c*c);
 306             if(*s == 1) bad = TRUE;
 307             *s = -1;
 308             prnk = MAX(prnk, ip);
 309             crnk = in;
 310         }
 311     }
 312     if(bad) {
 313         *s = BAD_SIGN;
 314     }
 315 }
 316
 317
 318 /***************************************************************************
 319  * Divide the constraint by its common-sub-expressions.                    *
 320  *                                                                         *
 321  *  1) first find the integer gcd value. if it is > 1 divide the entire    *
 322  *     constraint by the gcd.                                              *
 323  *  2) Check if a common-sub-expression (>0) exist for all the non-zero    *
 324  *     terms. If so divide by that common-sub-expression.                  *
 325  ***************************************************************************/
 326 void sc_fm_constraint::div_by_gcd()
 327 {
 328     CALL_STAT(stat_fc_div_by_gcd);
 329
 330     boolean redo = FALSE;
 331     // check for a simple integer multiple
 332     int g = -1;
 333     for(int ip=0; (ip<p())&&(g!=1); ip++) {
 334         for(int in=0; (in<n())&&(g!=1); in++)
 335             if((*this)[ip][in]) {
 336                 if(g == -1)
 337                     g = ABS((*this)[ip][in]);
 338                 else
 339                     g = gcd(g, ABS((*this)[ip][in]));
 340             }
 341     }
 342
 343     if(g==-1) return;
 344     if(g==0) return;
 345
 346     if(g != 1) {
 347         for(int ip=0; ip<p(); ip++) {
 348             for(int in=0; in<n(); in++)
 349                 (*this)[ip][in] /= g;
 350         }
 351         redo = TRUE;
 352     }
 353
 354
 355     // check for a multiple of a linear expression
 356     constraint nz(p());
 357     constraint mnz(p());
 358     constraint tmp(p());
 359     boolean foundnz = FALSE;
 360     boolean matched = TRUE;
 361     for(int in=0; (in<n())&&(matched); in++)
 362         if(!is_zero(in)) {
 363             if(!foundnz) {
 364                 get(nz, in);
 365                 if(nz.rank() == 0) return;      // not a lin expr.
 366                 int rg = nz.row_gcd();
 367                 nz /= rg;
 368
 369                 // the sign of all the non-zero elements has to be the same
 370                 int sgn = sign_constraint(nz);
 371                 if(sgn == 0) return;
 372                 nz *= sgn;                      // nz is always positive
 373
 374                 mnz = nz;
 375                 mnz *= -1;                      // nmz is -nz (always neg.)
 376                 foundnz = TRUE;
 377             } else {
 378                 get(tmp, in);
 379                 int rg = tmp.row_gcd();
 380                 tmp /= rg;
 381                 if((tmp != nz)&&(tmp != mnz))
 382                     matched = FALSE;
 383             }
 384         }
 385
 386     if(foundnz && matched) {
 387         for(int in=0; in<n(); in++) {
 388             if(!is_zero(in)) {
 389                 get(tmp, in);
 390                 int rg = tmp.row_gcd();
 391                 tmp /= rg;
 392                 if(tmp == nz) {
 393                     for(int ip=1; ip<p(); ip++)
 394                         (*this)[ip][in] = 0;
 395                     (*this)[0][in] = rg;
 396                 } else if(tmp == mnz) {
 397                     for(int ip=1; ip<p(); ip++)
 398                         (*this)[ip][in] = 0;
 399                     (*this)[0][in] = -rg;
 400                 }  else
 401                     assert(0);
 402             }
 403         }
 404         redo = TRUE;
 405     }
 406
 407     if(redo) recalc();
 408 }
 409
 410
 411 /***************************************************************************
 412  *  Find the inverse of the inequality; negate everything and substract 1  *
 413  ***************************************************************************/
 414 void sc_fm_constraint::inverse()
 415 {
 416     CALL_STAT(stat_fc_inverse);
 417
 418     for(int ip=0; ip<p(); ip++)
 419         for(int in=0; in<n(); in++)
 420             (*this)[ip][in] *= -1;
 421     (*this)[0][0] += -1;
 422     recalc();
 423 }
 424
 425
 426 /***************************************************************************
 427  * Check if only non-zero values are the constants and those are also      *
 428  * positive constants.  If so, this inequality does not carry any          *
 429  * information, thus can be eliminated.                                    *
 430  ***************************************************************************/
 431 boolean sc_fm_constraint::is_only_pos_const()
 432 {
 433     CALL_STAT(stat_fc_is_only_pos_const);
 434
 435     if(rank() > 0) return FALSE;
 436
 437     int ip;
 438     for(ip=0; ip<p(); ip++)
 439         for(int in=1; in<n(); in++)
 440             if((*this)[ip][in]) return FALSE;
 441
 442     int cnt=0;
 443     for(ip=1; ip<p(); ip++)
 444         if((*this)[ip][0]<0)
 445             return FALSE;
 446         else if((*this)[ip][0] > 0)
 447             cnt += (*this)[ip][0];
 448
 449     // since all plane vars V >= 1 we can assume aV1+bV2 >= a+b
 450     // thus c + aV1+bV2 >=0 -> c+a+b >= 0
 451     if((*this)[0][0] < -cnt) return FALSE;
 452
 453     return TRUE;
 454 }
 455
 456 /***************************************************************************
 457  *
 458  ***************************************************************************/
 459 boolean sc_fm_constraint::is_coef_eq_1(int j)
 460 {
 461     CALL_STAT(stat_fc_is_coef_eq_1);
 462
 463     if(ABS((*this)[0][j]) > 1) return FALSE;
 464     for(int ip=1; ip<p(); ip++)
 465         if((*this)[ip][j]) return FALSE;
 466     return TRUE;
 467 }
 468
 469 /***************************************************************************
 470  *  Check if a given column of the ineqality is zero.  i.e. a variable is  *
 471  *  not present in the inequality.                                         *
 472  ***************************************************************************/
 473 boolean sc_fm_constraint::is_zero(int in)
 474 {
 475     for(int ip=0; ip<p(); ip++)
 476         if((*this)[ip][in]) return FALSE;
 477     return TRUE;
 478 }
 479
 480
 481 /***************************************************************************
 482  *  Get the coefficient of a given variable (or the constant value)        *
 483  ***************************************************************************/
 484 void sc_fm_constraint::get(constraint &c, int in)
 485 {
 486     assert(c.n() == p());
 487     for(int ip=0; ip<p(); ip++)
 488         c[ip] = (*this)[ip][in];
 489 }
 490
 491
 492
 493
 494 /***************************************************************************
 495  *  Error handling                                                         *
 496  ***************************************************************************/
 497 void sc_fm_constraint::error()
 498 {
 499     if(_named_sc_fm_longjmp_env)
 500         longjmp(*_named_sc_fm_longjmp_env, 1);
 501     assert(0);
 502 }
 503
 504
 505 /***************************************************************************
 506  *  print the information in a fm_constraint
 507  ***************************************************************************/
 508 void sc_fm_constraint::print(FILE *fp)
 509 {
 510     fprintf(fp,"rank=%d, plane_rank=%d,  magic=%d\n",
 511            rank(),
 512            plane_rank(),
 513            magic_num);
 514     for(int in=0; in<n(); in++) {
 515         if(sgn[in] == BAD_SIGN)
 516             fprintf(fp,"xx: ");
 517         else
 518             fprintf(fp,"%2d: ", sgn[in]);
 519         for(int ip=0; ip<p(); ip++)
 520             fprintf(fp,"%2d ", (*this)[ip][in]);
 521         fprintf(fp,"\n");
 522     }
 523 }
 524
 525 /* ##################################################
 526    #####   sc_fm_results                       #####
 527    ################################################## */
 528
 529 sc_fm_results::sc_fm_results()
 530 {
 531     CALL_STAT(stat_fr_cons);
 532
 533     mm = 0;
 534     sz = 0;
 535     L = NULL;
 536 }
 537
 538
 539 sc_fm_results::~sc_fm_results()
 540 {
 541     if(L) delete[] L;
 542 }
 543
 544
 545 void sc_fm_results::init(int im, int s)
 546 {
 547     CALL_STAT(stat_fr_init);
 548
 549     if((im==m())&&(s==sz)) return;
 550
 551     if(L) delete[] L;
 552
 553     mm = im;
 554     sz = s;
 555
 556     L = (m()*sz)?(new int[m()*sz]):NULL;
 557 }
 558
 559
 560 void sc_fm_results::set(int im, sc_fm_constraint & c)
 561 {
 562     CALL_STAT(stat_fr_set);
 563
 564     assert(sz == c.p()*c.n());
 565     int * dat = (*this)[im];
 566     for(int ip=0; ip<c.p(); ip++)
 567         for(int in=0; in<c.n(); in++)
 568             *dat++ = c[ip][in];
 569 }
 570
 571
 572
 573 /* ##################################################
 574    #####   named_sc_fm                         #####
 575    ################################################## */
 576
 577 /***************************************************************************
 578  *  Delete the FM-solver                                                   *
 579  ***************************************************************************/
 580 named_sc_fm::~named_sc_fm()
 581 {
 582 #ifdef STAT_ON
 583     long endtime = get_time();
 584     endtime -= stat_sttime;
 585     double t = ((double)endtime)/(double)60.0;
 586     fprintf(stderr, "@@@ %5.2f @@@(%d,%d,%d)-%d/%d->(%d(%d),%d,%d)@@@\n",
 587             t,
 588             stat_stm, stat_stn, stat_stp,
 589             stat_nexp, stat_ninit,
 590             stat_ninitex, stat_maxm, stat_stn, stat_stp);
 591 #endif
 592
 593     if(results) delete[] results;
 594 }
 595
 596
 597 /***************************************************************************
 598  *   Resize the static data block                                          *
 599  ***************************************************************************/
 600 void named_sc_fm::init_block(int newm)
 601 {
 602     CALL_STAT(stat_fm_cons);
 603
 604 #ifdef STAT_ON
 605     stat_ninit++;
 606     stat_ninitex = MAX(newm, stat_ninitex);
 607 #endif
 608
 609     results = (n())?(new sc_fm_results[n()]):NULL;
 610
 611     if((newm <= cm) &&
 612        (n() == cn) &&
 613        (p() == cp)) {
 614
 615         for(int im=0; im<cm; im++) {
 616             for(int ip=0; ip<cp; ip++)
 617                 for(int in=0; in<cn; in++)
 618                     (*L[im])[ip][in] = 0;
 619             magic_list[im] = 0;
 620         }
 621         return;
 622     }
 623
 624     if(cm > 0) {
 625         for(int i=0; i<cm; i++) delete L[i];
 626         if(L) delete L;
 627         if(magic_list) delete[] magic_list;
 628     }
 629
 630     cm = next_up(newm, cm, DIM_M_ST, DIM_M_OFF);
 631     cn = n();
 632     cp = p();
 633
 634 #ifdef STAT_ON
 635     stat_maxm = MAX(stat_maxm, cm);
 636     stat_nexp++;
 637 #endif
 638
 639     L = (cm)?(new p_sc_fm_constraint[cm]):NULL;
 640     magic_list = (cm)?(new unsigned[cm]):NULL;
 641     for(int i=0; i<cm; i++) {
 642         L[i] = new sc_fm_constraint(cn, cp);
 643         magic_list[i] = 0;
 644         assert(L[i]);
 645     }
 646     xm = 0;
 647 }
 648
 649 #ifdef WATCH_ON
 650 #define WATCH_CNT       1000
 651 #define WATCH_PAIR_CALL 5000
 652 #define WATCH_SING_CALL  500
 653
 654 int watch_next = WATCH_CNT;
 655 int watch_step1 = WATCH_PAIR_CALL;
 656 int watch_step2 = WATCH_SING_CALL;
 657 int watch_ns1 = 0;
 658 int watch_ns2 = 0;
 659 #endif
 660
 661 /***************************************************************************
 662  *  Increase the size of the static data block. The block need to          *
 663  *  represent newm-inequalilites.                                          *
 664  ***************************************************************************/
 665 void named_sc_fm::extend_block(int newm)
 666 {
 667     CALL_STAT(stat_fm_extend_block);
 668
 669 #ifdef STAT_ON
 670     stat_ninit++;
 671     stat_ninitex = MAX(newm, stat_ninitex);
 672 #endif
 673
 674     if(newm<=cm)return;
 675
 676     int oldcm  = cm;
 677     cm = next_up(newm, cm, DIM_M_ST, DIM_M_OFF);
 678
 679 #ifdef STAT_ON
 680     stat_maxm = MAX(stat_maxm, cm);
 681     stat_nexp++;
 682 #endif
 683 #ifdef WATCH_ON
 684     if(cm > watch_next) {
 685         fprintf(stderr, "-watch--cm=%d--\n", cm);
 686         watch_next += WATCH_CNT;
 687     }
 688 #endif
 689
 690
 691     sc_fm_constraint ** newL = (cm)?(new p_sc_fm_constraint[cm]):NULL;
 692     unsigned * newml = (cm)?(new unsigned[cm]):NULL;
 693
 694     int i;
 695     for(i=0; i<oldcm; i++) {
 696         newL[i] = L[i];
 697         newml[i] = magic_list[i];
 698     }
 699     for(; i<cm; i++) {
 700         newL[i] = new sc_fm_constraint(cn, cp);
 701         assert(newL[i]);
 702         newml[i] = 0;
 703     }
 704
 705     if(L) delete L;
 706     L = newL;
 707
 708     if(magic_list) delete[] magic_list;
 709     magic_list = newml;
 710
 711     for(int x=0; x<cm; x++) assert(L[x]);
 712 }
 713
 714
 715
 716 /***************************************************************************
 717  *   Initialize the FM-problem.                                            *
 718  *                                                                         *
 719  * Initialize the size of the problem and the static data block using the  *
 720  * input data-set.  Calculate the cached results.                          *
 721  ***************************************************************************/
 722 void named_sc_fm::init(named_symcoeff_ineq & x, boolean k)
 723 {
 724
 725     keep = k;
 726     cutoffpt = 0;
 727     nm_p.init(x.planes());
 728     nm_c.init(x.cols());
 729
 730 #ifdef STAT_ON
 731     stat_maxm = cm;
 732     stat_stm = x.m();
 733     stat_stn = n();
 734     stat_stp = p();
 735     stat_nexp = 0;
 736     stat_ninit = stat_ninitex = 0;
 737
 738     stat_sttime = get_time();
 739 #endif
 740
 741     init_block(nrows_in_blk(p(), x.m(), n()));
 742
 743     xm = x.m();
 744     for(int im=0; im<m(); im++) {
 745         for(int ip=0; ip<p(); ip++)
 746             for(int in=0; in<n(); in++)
 747                 (*this)[im][ip][in] = x[ip][im][in];
 748         (*this)[im].recalc();
 749         magic_list[im] = (*this)[im].magic_num;
 750     }
 751 }
 752
 753 /***************************************************************************
 754  *   write back to a named_symcoeff_ineq with internal data                 *
 755  ***************************************************************************/
 756 named_symcoeff_ineq * named_sc_fm::internal_get()
 757 {
 758     CALL_STAT(stat_fm_get);
 759
 760     named_symcoeff_ineq * ret = new named_symcoeff_ineq();
 761     ret->init(p(), m(), n());
 762     ret->planes().init(nm_p);
 763     ret->cols().init(nm_c);
 764
 765     for(int im=0; im<m(); im++)
 766         for(int ip=0; ip<p(); ip++)
 767             for(int in=0; in<n(); in++)
 768                 (*ret)[ip][im][in] = (*this)[im][ip][in];
 769
 770     return ret;
 771 }
 772
 773
 774 /***************************************************************************
 775  *   update the named_symcoeff_ineq using the current result                *
 776  ***************************************************************************/
 777 named_symcoeff_ineq * named_sc_fm::get()
 778 {
 779     assert(p() == nm_p.n());
 780     assert(n() == nm_c.n());
 781     named_symcoeff_ineq * ret = new named_symcoeff_ineq();
 782     get(ret);
 783     return ret;
 784 }
 785
 786
 787 /***************************************************************************
 788  *   write back to a named_symcoeff_ineq of the current result              *
 789  ***************************************************************************/
 790 void named_sc_fm::get(named_symcoeff_ineq * ret)
 791 {
 792     CALL_STAT(stat_fm_get);
 793
 794     assert(ret);
 795     assert(nm_p.n() == p());
 796     assert(nm_c.n() == n());
 797
 798     int totm = 0;
 799     for(int i=0; i<n(); i++)
 800         totm += results[i].m();
 801
 802     ret->init(p(), totm, n());
 803     ret->planes().init(nm_p);
 804     ret->cols().init(nm_c);
 805
 806     int currm = 0;
 807     for(int inx=0; inx<n(); inx++) {
 808         sc_fm_results & res = results[inx];
 809         for(int im=0; im<res.m(); im++) {
 810             int * dat = res[im];
 811             for(int ip=0; ip<p(); ip++)
 812                 for(int in=0; in<n(); in++)
 813                     (*ret)[ip][currm][in] = *dat++;
 814             currm++;
 815         }
 816     }
 817     assert(currm == totm);
 818 }
 819
 820
 821
 822
 823
 824 /***************************************************************************
 825  * Perform Fourier-Motzkin elimination for the entire problem.             *
 826  * Return TRUE  if a valid solution exists,                                *
 827  *        FALSE if the problem is inconsistant.                            *
 828  ***************************************************************************/
 829 boolean named_sc_fm::fm_step()
 830 {
 831     CALL_STAT(stat_fm_step0);
 832
 833 #ifdef WATCH_ON
 834     watch_ns2 = 0;
 835     watch_step2 = WATCH_SING_CALL;
 836 #endif
 837     return fm_step(0, n());
 838 }
 839
 840
 841 /***************************************************************************
 842  *  Perform FM-elimination from columns j backto i (both inclusive)        *
 843  *  Special step for the constant term.                                    *
 844  ***************************************************************************/
 845 boolean named_sc_fm::fm_step(int i, int j)
 846 {
 847     CALL_STAT(stat_fm_step1);
 848
 849     if(MAX(1, i)<j)
 850         for(int ii=0; ii<m(); ii++)
 851             if(check_valid(ii) == FALSE) return FALSE;
 852
 853     for(int x=j-1; x>=MAX(1, i); x--)
 854         if(fm_step(x) == FALSE) return FALSE;
 855     if(i==0) return fm_lin_step();
 856     return TRUE;
 857 }
 858
 859
 860 /***************************************************************************
 861  *  Perform FM-elimination of the j-th variable.                           *
 862  ***************************************************************************/
 863 boolean named_sc_fm::fm_step(int j)
 864 {
 865     CALL_STAT(stat_fm_step2);
 866
 867     assert((j>0)&&(j<n()));
 868 #ifdef WATCH_ON
 869     if(++watch_ns2 >= watch_step2) {
 870         fprintf(stderr, "-watch--#call to single step=%d--\n", watch_ns2);
 871         watch_step2 += WATCH_SING_CALL;
 872     }
 873     watch_ns1 = 0;
 874     watch_step1 = WATCH_PAIR_CALL;
 875 #endif
 876
 877     int low, up, none;
 878     low = up = none = 0;
 879
 880
 881     // find number of inequalities that have a positive coefficient for j,
 882     // negagive coefficient for j and j is not a participant
 883     int i;
 884     for(i=0; i<m(); i++) {
 885         int sg = (*this)[i].sign(j);
 886         if(sg > 0) low++;
 887         else if(sg < 0) up++;
 888         else none++;
 889     }
 890
 891     // Nothing to eliminate by FM
 892     if((low == 0)||(up == 0)) {
 893         if(!keep) {
 894             if(m()) {
 895                 int * del = new int[m()];
 896                 for(int im=0; im<m(); im++)
 897                     del[im] = ((*this)[im].sign(j) != 0);
 898                 remove(del);
 899                 delete[] del;
 900             }
 901         }
 902         return TRUE;
 903     }
 904
 905     // Resize the block to handle all the new inequalities.
 906     extend_block(m() + low*up);
 907     int oldm = m();
 908
 909
 910 #ifdef DEBUG_PRINT3
 911     fprintf(stderr,"....................\n [ ");
 912     for(int ix=0; ix<m(); ix++)
 913         fprintf(stderr,"%d ", (*this)[ix].sign(j));
 914     fprintf(stderr," ]\n");
 915 #endif
 916
 917     // Only new inequalities are the ones after the cutoffpt
 918     // Thus, we do not have to eliminate two inequalities that are both
 919     // old since that should have been already done.
 920     fm_step(j,         0, cutoffpt,  cutoffpt,     oldm);
 921     fm_step(j,  cutoffpt,     oldm,         0, cutoffpt);
 922     fm_step(j,  cutoffpt,     oldm,  cutoffpt,     oldm);
 923
 924
 925 #ifdef DEBUG_PRINT
 926     fprintf(stderr,"fm_step---%d(%s)--- [oldM=%d, #chk=%d -> newM=%d]\n",
 927            j, get_name(nm_c[j]),
 928            oldm, oldm-cutoffpt, m());
 929 #ifdef DEBUG_PRINT2
 930     if(oldm != m())
 931         print(stderr);
 932     else
 933         fprintf(stderr,"No change\n");
 934 #endif
 935 #endif
 936
 937
 938     // Check if the resulting system of inequalities is consistant
 939     for(i=oldm; i<m(); i++)
 940         if(check_valid(i) == FALSE) return FALSE;
 941
 942     if(!keep) {
 943         // remove the inequalities that participated in the FM-step
 944         if(m()) {
 945             int * del = new int[m()];
 946             for(int im=0; im<m(); im++)
 947                 del[im] = ((*this)[im].sign(j) != 0);
 948             remove(del);
 949             delete[] del;
 950         }
 951     }
 952
 953
 954     return TRUE;
 955 }
 956
 957
 958 #ifdef DEBUG_PRINT3
 959 static int tmp_var;
 960 #endif
 961
 962 /***************************************************************************
 963  *  For the given range of inequalities for coef of j>0 and coef of j<0,   *
 964  *  find the pairs, get the resulting inequality after performing the FM   *
 965  *  elimination step.  If that inquality is not already in and has         *
 966  *  information, add that into the system.                                 *
 967  ***************************************************************************/
 968 void named_sc_fm::fm_step(int j,
 969                            int la, int lb,
 970                            int ua, int ub)
 971 {
 972     CALL_STAT(stat_fm_step3);
 973
 974     if((!(la<lb))||(!(ua<ub))) return;
 975
 976     for(int l=la; l<lb; l++)
 977         if((*this)[l].sign(j) > 0)
 978             for(int u=ua; u<ub; u++)
 979                 if((*this)[u].sign(j) < 0) {
 980                     int tm = xm++;
 981 #ifdef DEBUG_PRINT3
 982                     fprintf(stderr," (%d %d)%d", l, u, tm); fflush(stderr);
 983 #endif
 984
 985                     fm_step(j, (*this)[tm], (*this)[l], (*this)[u]);
 986                     (*this)[tm].div_by_gcd();
 987                     magic_list[tm] = (*this)[tm].magic_num;
 988
 989 #ifdef DEBUG_PRINT3
 990                     fprintf(stderr,"%s", (*this)[tm].is_only_pos_const()?"C":"");
 991                     if(is_already_covered(tm)) fprintf(stderr,"I(%d)", tmp_var);
 992 #endif
 993                     if((*this)[tm].is_only_pos_const())
 994                         xm--;
 995                     else if(is_already_covered(tm))
 996                         xm--;
 997                 }
 998 #ifdef DEBUG_PRINT3
 999     fprintf(stderr,"\n");
1000 #endif
1001 }
1002
1003
1004 /***************************************************************************
1005  * Perform FM-elimination of the two inequalities A and B w.r.t. the       *
1006  * variable j and deposit the result in res.                               *
1007  * Many different cases to consider; speed-hacks provided to speed-up the  *
1008  * calculations.                                                           *
1009  ***************************************************************************/
1010 void named_sc_fm::fm_step(int j,
1011                            sc_fm_constraint & res,
1012                            sc_fm_constraint & A,
1013                            sc_fm_constraint & B)
1014 {
1015     CALL_STAT(stat_fm_step4);
1016
1017     assert((A.p() == B.p())&&(A.n() == B.n()));
1018     int as = A.sign(j);
1019     int bs = B.sign(j);
1020     assert(as * bs == -1);
1021
1022
1023 #ifdef WATCH_ON
1024     if(++watch_ns1 >= watch_step1) {
1025         fprintf(stderr, "-watch--#call to pair-wise step=%d--\n", watch_ns1);
1026         watch_step1 += WATCH_PAIR_CALL;
1027     }
1028 #endif
1029
1030
1031     // check if A[j] and B[j] if of the form ca*A[j] == cb*B[j]
1032     constraint cxa(p());
1033     constraint cxb(p());
1034     A.get(cxa, j);
1035     B.get(cxb, j);
1036     int ca = cxa[cxa.rank()];
1037     int cb = cxb[cxb.rank()];
1038     ca = ABS(ca);
1039     cb = ABS(cb);
1040
1041     cxa *= cb;
1042     cxb *= ca;
1043     constraint cc = cxa + cxb;
1044     if((cc.rank()==0)&&(cc[0] == 0)) {
1045         // res = cb*A + ca*B
1046
1047 #ifndef USE_SPEED_HACK
1048         for(int ip=0; ip<p(); ip++)
1049             for(int in=0; in<n(); in++)
1050                 res[ip][in] = cb*A[ip][in] + ca*B[ip][in];
1051 #else
1052         int * pres = &res[0][0];
1053         int * bound = &res[p()-1][n()-1];
1054         int * pa   = &A[0][0];
1055         int * pb   = &B[0][0];
1056         while(pres <= bound)
1057             *pres++ = cb * *pa++ + ca * *pb++;
1058 #endif
1059
1060
1061         assert((res.p() == A.p())&&(res.n() == A.n()));
1062         res.recalc();
1063         return;
1064     }
1065
1066     if(cxa.rank() == 0) {
1067 #ifndef USE_SPEED_HACK
1068         for(int ip=0; ip<p(); ip++)
1069             for(int in=0; in<n(); in++)
1070                 res[ip][in] = B[ip][in]*ABS(A[0][j]);
1071 #else
1072         int * pres = &res[0][0];
1073         int * bound = &res[p()-1][n()-1];
1074         int * pb   = &B[0][0];
1075         int mul = ABS(A[0][j]);
1076         while(pres <= bound)
1077             *pres++ =  *pb++ * mul;
1078 #endif
1079     } else {
1080         if(B.plane_rank() != 0) {
1081             if(_named_sc_fm_longjmp_env)
1082                 longjmp(*_named_sc_fm_longjmp_env, 1);
1083             fprintf(stderr,"Cannot reduce, too complicated (1)\n");
1084             fprintf(stderr,"j=%d\n", j);
1085             fprintf(stderr,"A\n");   A.print(stderr);
1086             fprintf(stderr,"B\n");   B.print(stderr);
1087             fprintf(stderr,"as=%d bs=%d\nca=%d cb=%d\n", as, bs, ca, cb);
1088             fprintf(stderr,"cxa = "); cxa.print(stderr);
1089             fprintf(stderr,"cxb = "); cxb.print(stderr);
1090             fprintf(stderr,"cc  = "); cc.print(stderr);
1091             error();
1092         }
1093
1094 #ifndef USE_SPEED_HACK
1095         for(int ip=0; ip<p(); ip++)
1096             for(int in=0; in<n(); in++)
1097                 res[ip][in] = B[0][in]*ABS(A[ip][j]);
1098 #else
1099         for(int ip=0; ip<p(); ip++) {
1100             int * pres = &res[ip][0];
1101             int * bound = &res[ip][n()-1];
1102             int * pb   = &B[0][0];
1103             int mul = ABS(A[ip][j]);
1104             while(pres <= bound)
1105                 *pres++ = *pb++ * mul;
1106         }
1107 #endif
1108     }
1109
1110     if(cxb.rank() == 0) {
1111 #ifndef USE_SPEED_HACK
1112         for(int ip=0; ip<p(); ip++)
1113             for(int in=0; in<n(); in++)
1114                 res[ip][in] += A[ip][in]*ABS(B[0][j]);
1115 #else
1116         int * pres = &res[0][0];
1117         int * bound = &res[p()-1][n()-1];
1118         int * pa   = &A[0][0];
1119         int mul = ABS(B[0][j]);
1120         while(pres <= bound)
1121             *pres++ +=  *pa++ * mul;
1122 #endif
1123     } else {
1124         if(A.plane_rank() != 0) {
1125             if(_named_sc_fm_longjmp_env)
1126                 longjmp(*_named_sc_fm_longjmp_env, 1);
1127             fprintf(stderr,"Cannot reduce, too complicated (2)\n");
1128             fprintf(stderr,"j=%d\n", j);
1129             fprintf(stderr,"A\n");   A.print(stderr);
1130             fprintf(stderr,"B\n");   B.print(stderr);
1131             fprintf(stderr,"as=%d bs=%d\nca=%d cb=%d\n", as, bs, ca, cb);
1132             fprintf(stderr,"cxa = "); cxa.print(stderr);
1133             fprintf(stderr,"cxb = "); cxb.print(stderr);
1134             fprintf(stderr,"cc  = "); cc.print(stderr);
1135             error();
1136         }
1137
1138 #ifndef USE_SPEED_HACK
1139         for(int ip=0; ip<p(); ip++)
1140             for(int in=0; in<n(); in++)
1141                 res[ip][in] += A[0][in]*ABS(B[ip][j]);
1142 #else
1143         for(int ip=0; ip<p(); ip++) {
1144             int * pres = &res[ip][0];
1145             int * bound = &res[ip][n()-1];
1146             int * pa   = &A[0][0];
1147             int mul = ABS(B[ip][j]);
1148             while(pres <= bound)
1149                 *pres++ += *pa++ * mul;
1150         }
1151 #endif
1152
1153     }
1154     assert((res.p() == A.p())&&(res.n() == A.n()));
1155     res.recalc();
1156 }
1157
1158
1159
1160 /***************************************************************************
1161  * Is the constant term consistant.                                        *
1162  *                                                                         *
1163  * This can be treated as a linear inequality of the symbolic constants    *
1164  * (the planes).  Since each symbolic constant is > 0, need to check all   *
1165  * the inequalities with only a constant term to see if the system is      *
1166  * consistant.                                                             *
1167  ***************************************************************************/
1168 boolean named_sc_fm::fm_lin_step()
1169 {
1170     CALL_STAT(stat_fm_lin_step);
1171
1172     if(m()==0) return TRUE;
1173
1174     int cnt = 0;
1175     int * clist = new int[m()];
1176     int i;
1177     for(i=0; i<m(); i++) {
1178         if((*this)[i].rank() == 0) {
1179             cnt++;
1180             clist[i] = 1;
1181         } else
1182             clist[i] = 0;
1183     }
1184
1185     if(cnt==0) {
1186         delete[] clist;
1187         return TRUE;
1188     }
1189
1190     lin_ineq ll(cnt, p());
1191
1192     int c = 0;
1193     for(i=0; i<m(); i++)
1194         if(clist[i]) {
1195             for(int pp=0; pp<p(); pp++)
1196                 ll[c][pp] = (*this)[i][pp][0];
1197             c++;
1198         }
1199     assert(c == cnt);
1200     delete[] clist;
1201
1202
1203     nim_op O;
1204     // add ineqs saying all plane vars >= 1
1205     lin_ineq x(p()-1, 1);
1206     for(i=0; i<p()-1; i++) x[i][0] = -1;
1207     lin_ineq lla = Compose(2, 1,
1208                            O.NIM(ll),
1209                            O.NIM(Compose(1, 2,
1210                                          O.NIM(x),
1211                                          O.NIM(Ident(p()-1)))));
1212
1213 #ifdef DEBUG_PRINT
1214     named_lin_ineq leq(nm_p, lla);
1215     fprintf(stderr,"fm_step---lin-- ineqs=%d -\n", lla.m());
1216     leq.print_exp(pet_system_nl,stderr);
1217 #endif
1218
1219     if(~lla == TRUE) return FALSE;
1220     return TRUE;
1221 }
1222
1223
1224
1225
1226
1227
1228 void named_sc_fm::fm_project()
1229 {
1230     CALL_STAT(stat_fm_bounds0);
1231
1232     fm_project(1, n());
1233 }
1234
1235
1236 void named_sc_fm::fm_project(int i, int j)
1237 {
1238     CALL_STAT(stat_fm_bounds1);
1239
1240 #ifdef DEBUG_PRINT
1241     fprintf(stderr,"@@@@@@@@@@@@@F@M@@P@R@O@J@E@C@T@@@@@@@@@@@@@@@@\n");
1242     sort();
1243     print(stderr);
1244 #endif
1245     assert((0<i)&&(i<=n()));
1246     assert((0<j)&&(j<=n())&&(i<=j));
1247 #ifdef WATCH_ON
1248     watch_ns2 = 0;
1249     watch_step2 = WATCH_SING_CALL;
1250 #endif
1251
1252     keep = TRUE;
1253
1254     fm_step(i, j);
1255
1256 #ifdef DEBUG_PRINT
1257     fprintf(stderr,"@@@@@@@@@@@@@@@@@@after fm_project(%d, %d)@@@@@@@@@@@@@@@@@@@@@\n", i, j);
1258     sort();
1259     print(stderr);
1260 #endif
1261
1262     for(int x=j-1; x>=i; x--) {
1263 #ifdef DEBUG_PRINT
1264         fprintf(stderr,"post step================\n");
1265 #endif
1266         fm_project(x);
1267
1268     }
1269 }
1270
1271
1272 void named_sc_fm::fm_project(int j)
1273 {
1274     fm_bounds_do(j, TRUE);
1275     fm_bounds_rec(j);
1276 }
1277
1278
1279
1280 /***************************************************************************
1281  * Find efficient loop bounds for the variables given by the system of     *
1282  * inequalities.                                                           *
1283  * Assumption: the system of inequalities is valid                         *
1284  ***************************************************************************/
1285 void named_sc_fm::fm_bounds()
1286 {
1287     CALL_STAT(stat_fm_bounds0);
1288
1289     fm_bounds(1, n());
1290 }
1291
1292
1293 /***************************************************************************
1294  * Find efficient loop bounds for the i-th to j-th variables in the system *
1295  * of inequalities.                                                        *
1296  * Assumption: the system of inequalities is valid                         *
1297  ***************************************************************************/
1298 void named_sc_fm::fm_bounds(int i, int j)
1299 {
1300     CALL_STAT(stat_fm_bounds1);
1301
1302 #ifdef DEBUG_PRINT
1303     fprintf(stderr,"@@@@@@@@@@@@@F@M@@B@O@U@N@D@S@@@@@@@@@@@@@@@@@@\n");
1304     sort();
1305     print(stderr);
1306 #endif
1307     assert((0<i)&&(i<=n()));
1308     assert((0<j)&&(j<=n())&&(i<=j));
1309 #ifdef WATCH_ON
1310     watch_ns2 = 0;
1311     watch_step2 = WATCH_SING_CALL;
1312 #endif
1313
1314     keep = TRUE;
1315
1316
1317     fm_step(i, j);
1318
1319 #ifdef DEBUG_PRINT
1320     fprintf(stderr,"@@@@@@@@@@@@@@@@@@after fm_step(%d, %d)@@@@@@@@@@@@@@@@@@@@@\n", i, j);
1321     sort();
1322     print(stderr);
1323 #endif
1324
1325     for(int x=j-1; x>=i; x--) {
1326 #ifdef DEBUG_PRINT
1327         fprintf(stderr,"post step================\n");
1328 #endif
1329         fm_bounds(x);
1330
1331     }
1332 }
1333
1334
1335 void named_sc_fm::fm_bounds(int j)
1336 {
1337     fm_bounds_do(j, FALSE);
1338     fm_bounds_rec(j);
1339 }
1340
1341
1342
1343 void named_sc_fm::fm_bounds_rec(int x)
1344 {
1345     int cnt = 0;
1346     for(int i=0; i< m(); i++)
1347         if((*this)[i].rank() == x)
1348             cnt++;
1349
1350     if(cnt>0) {
1351         results[x].init(cnt, n()*p());
1352
1353         int curr = 0;
1354         for(int i=m()-1; i>= 0; i--)
1355             if((*this)[i].rank() == x) {
1356                 results[x].set(curr++, (*this)[i]);
1357                 remove(i);
1358             }
1359         assert(curr == cnt);
1360     }
1361 }
1362
1363
1364 /***************************************************************************
1365  *  Find efficient loop bounds for the j-th variable.                      *
1366  ***************************************************************************/
1367 void named_sc_fm::fm_bounds_do(int j, boolean no_del)
1368 {
1369     CALL_STAT(stat_fm_bounds2);
1370
1371     assert((j>0)&&(j<n()));
1372
1373     if(m() == 0) return;
1374
1375     sort();
1376
1377 #ifdef DEBUG_PRINT
1378     fprintf(stderr,"fm_bounds %d(%s) ===============================\n", j, get_name(nm_c[j]));
1379 #endif
1380
1381 #ifdef DEBUG_PRINT
1382     named_symcoeff_ineq * tmp = internal_get();
1383 #endif
1384
1385     int low, up, none;
1386     low = up = none = 0;
1387
1388     int * s_l = new int[m()];
1389     int * s_u = new int[m()];
1390     int * s_eq= new int[m()];
1391     int * dumb = new int[m()];
1392
1393     int i;
1394     for(i=0; i<m(); i++) {
1395         s_eq[i]= 0;
1396         s_l[i] = 0;
1397         s_u[i] = 0;
1398         int s = 0;
1399         int r = (*this)[i].rank();
1400         dumb[i] = r;
1401         if(r==j) {
1402             s = (*this)[i].sign(j);
1403             s_l[i] = (s > 0);
1404             s_u[i] = (s < 0);
1405         }
1406         if(s > 0) low++;
1407         else if(s < 0) up++;
1408         else none++;
1409     }
1410
1411     int numeq = 0;
1412     int simpeq = -1;
1413     for(int l=0; l<m(); l++)
1414         if(s_l[l])
1415             for(int u=0; (u<m())&&(s_eq[l]==0); u++)
1416                 if(s_u[u] && (s_eq[u]==0)) {
1417                     if(is_add_eq_zero((*this)[l], (*this)[u])) {
1418                         numeq++;
1419                         s_eq[l] = numeq;
1420                         s_eq[u] = numeq;
1421                         if(simpeq == -1)
1422                             if((*this)[l].is_coef_eq_1(j))
1423                                 simpeq = numeq;
1424                     }
1425                 }
1426
1427     if(simpeq == -1) simpeq = 1;
1428     int eq[2];
1429     eq[0] = eq[1] = -1;
1430     int ex=0;
1431     if(numeq > 0) {
1432         // Found equality relationship
1433         // Give it presidence over other inequalities
1434         for(i=m()-1; i>=0; i--)
1435             if(s_eq[i] == simpeq) {
1436                 assert(ex<=1);
1437                 eq[ex++] = i;
1438             }
1439         assert(ex == 2);
1440     }
1441
1442     if(!no_del)
1443         fm_bounds(j, low, up, eq);
1444
1445
1446
1447 #ifdef DEBUG_PRINT
1448     named_symcoeff_ineq * tmpnew = internal_get();
1449     fprintf(stderr,"   %d--rnk--e-l-u-del--- numeq=%d simpeq=\n",
1450            j, numeq, simpeq);
1451     for(i=0; i<tmp->m(); i++) {
1452         fprintf(stderr,"%2d    %d    %d %s %s  %s   ",
1453                i,
1454                dumb[i],
1455                s_eq[i],
1456                (s_l[i])?"*":" ",
1457                (s_u[i])?"*":" ",
1458                chk_is_in(i, *tmp, *tmpnew)?" ":"X");
1459         tmp->print_exp(i,stderr); fprintf(stderr,"  >= 0\n");
1460     }
1461     delete tmp;
1462     delete tmpnew;
1463 #endif
1464
1465
1466     delete[] s_l;
1467     delete[] s_u;
1468     delete[] s_eq;
1469     delete[] dumb;
1470 }
1471
1472
1473 /***************************************************************************
1474  *  For each ineqality marked in ilist, mark it deleteable if it is not    *
1475  *  necessary.                                                             *
1476  *                                                                         *
1477  *  Invert the ineqality, and check if the system is consistant.           *
1478  *  BUGBUG: if ineq. is marked delete it should be deleted before checking *
1479  *          the next. Otherwise it might interfere.                        *
1480  ***************************************************************************/
1481 void named_sc_fm::fm_bounds(int j, int low, int up, int * eq)
1482 {
1483     CALL_STAT(stat_fm_bounds3);
1484
1485     for(int i=m()-1; i>=0; i--)
1486         if((i != eq[0])&&(i != eq[1]))  {
1487             int r = (*this)[i].rank();
1488             int s = (*this)[i].sign(j);
1489             if((r==j)&&(s)) {
1490                 (*this)[i].inverse();
1491                 magic_list[i] = (*this)[i].magic_num;
1492
1493                 swap(i, m()-1);     // so this can be compared with everyone by
1494                 // upping the cutoffpt by one
1495                 int currm = m();
1496                 cutoffpt = m()-1;
1497 #ifdef DEBUG_PRINT
1498                 fprintf(stderr,"{*** Inv %d\n", i);
1499 #endif
1500                 boolean ok = fm_step(0, n());
1501 #ifdef DEBUG_PRINT
1502                 fprintf(stderr," *** Inv %d res=%d }\n", i, ok);
1503 #endif
1504                 xm = currm;
1505                 cutoffpt = 0;
1506
1507                 swap(i, m()-1);
1508
1509                 (*this)[i].inverse();
1510                 magic_list[i] = (*this)[i].magic_num;
1511
1512                 if(!ok) {
1513                     if(s > 0) {
1514                         if(--low > 0)
1515                             remove(i);
1516                     } else {
1517                         if(--up > 0)
1518                             remove(i);
1519                     }
1520                 }
1521             }
1522         }
1523 }
1524
1525
1526
1527
1528 /***************************************************************************
1529  * Swap two ineqalities in the list.                                       *
1530  ***************************************************************************/
1531 void named_sc_fm::swap(int i, int j)
1532 {
1533     CALL_STAT(stat_fm_swap);
1534
1535     assert((i>=0)&&(i<m()));
1536     assert((j>=0)&&(j<m()));
1537
1538     if(i == j) return;
1539
1540     sc_fm_constraint * tmpl = L[i];
1541     L[i] = L[j];
1542     L[j] = tmpl;
1543
1544     unsigned tmpm = magic_list[i];
1545     magic_list[i] = magic_list[j];
1546     magic_list[j] = tmpm;
1547 }
1548
1549
1550 /***************************************************************************
1551  *   Sort the inequalities according to the column rank.                   *
1552  ***************************************************************************/
1553 void named_sc_fm::sort()
1554 {
1555     CALL_STAT(stat_fm_sort);
1556
1557     for(int i=0; i<m(); i++)
1558         for(int j=i+1; j<m(); j++)
1559             if((*this)[i].rank() > (*this)[j].rank()) {
1560                 swap(i, j);
1561             }
1562 }
1563
1564
1565 /***************************************************************************
1566  *  Check if the inequality is valid.                                      *
1567  *                                                                         *
1568  *  inequality is not valid if it has no variables and all the constant    *
1569  *  are less than zero.                                                    *
1570  ***************************************************************************/
1571 boolean named_sc_fm::check_valid(int r)
1572 {
1573     CALL_STAT(stat_fm_check_valid);
1574
1575     boolean found = FALSE;
1576     for(int d=0; d<p(); d++)
1577         if((*this)[r][d][0] < 0) found = TRUE;
1578         else if((*this)[r][d][0] > 0) return TRUE;
1579
1580     if(found) {
1581         for(int d=0; d<p(); d++)
1582             for(int c=1; c<n(); c++)
1583                 if((*this)[r][d][c]) return TRUE;
1584         return FALSE;
1585     }
1586     return TRUE;
1587 }
1588
1589
1590
1591 /***************************************************************************
1592  *  Check if the chk_this inequality is already covered by another         *
1593  *  ineqality in the system.                                               *
1594  *                                                                         *
1595  *                                                                         *
1596  *  Since this function is the most expensive in the profile, many things  *
1597  *  were done to speed this up.  A magic number is calculated for each     *
1598  *  ineqality s.t. if two ineqalities are the same => they have the same   *
1599  *  magic number.  That is used as a first cut to eliminate expensive      *
1600  *  tests.  Also all the magic numbers are kept in a seperate list so that *
1601  *  the N^2 test can be done cheaper.  The expensive part is called only   *
1602  *  when the matic numbers match.                                          *
1603  *  Still this N^2 is the most expensive.                                  *
1604  ***************************************************************************/
1605 boolean named_sc_fm::is_already_covered(int chk_this)
1606 {
1607     CALL_STAT(stat_fm_is_already0);
1608
1609
1610     sc_fm_constraint & newc = (*this)[chk_this];
1611
1612     assert(newc.n() == n());
1613     assert(newc.p() == p());
1614
1615 #ifndef USE_SPEED_HACK
1616     unsigned newc_magic = newc.magic_num;
1617     unsigned * pt_magic = &magic_list[0];
1618
1619     for(int im=0; im<m(); im++, pt_magic++)
1620         if(im != chk_this) {
1621             CALL_STAT(stat_fm_is_already0x);
1622             sc_fm_constraint & currc = (*this)[im];
1623             assert(currc.magic_num == *pt_magic);
1624             if(currc.magic_num == newc_magic) {
1625                 pt_magic++;
1626                 int ix = (pt_magic - 1 - &magic_list[0]);
1627                 pt_magic--;
1628                 assert(ix == im);
1629 #ifdef DEBUG_PRINT3
1630             tmp_var = im;
1631 #endif
1632                 if(is_already_covered(newc, currc)) return TRUE;
1633             }
1634         }
1635
1636 #else   // the SPEED_HACK
1637
1638     unsigned newc_magic = newc.magic_num;
1639     unsigned * pt_magic = &magic_list[0];
1640     unsigned * pt_end  =  &magic_list[chk_this];
1641     while(pt_magic < pt_end) {
1642         CALL_STAT(stat_fm_is_already0x);
1643         if(*pt_magic++ == newc_magic) {
1644             int im = (pt_magic - 1 - &magic_list[0]);
1645             sc_fm_constraint & currc = (*this)[im];
1646 #ifdef DEBUG_PRINT3
1647             tmp_var = im;
1648 #endif
1649             if(is_already_covered(newc, currc)) return TRUE;
1650         }
1651     }
1652     pt_magic++;
1653     pt_end  =  &magic_list[m()-1];
1654     while(pt_magic <= pt_end) {
1655         CALL_STAT(stat_fm_is_already0x);
1656         if(*pt_magic++ == newc_magic) {
1657             int im = (pt_magic - 1 - &magic_list[0]);
1658             sc_fm_constraint & currc = (*this)[im];
1659 #ifdef DEBUG_PRINT3
1660             tmp_var = im;
1661 #endif
1662             if(is_already_covered(newc, currc)) return TRUE;
1663         }
1664     }
1665 #endif  // USE_SPEED_HACK
1666
1667     return FALSE;
1668 }
1669
1670
1671 /***************************************************************************
1672  *  Check if the inequalities newc is coverd by oldc                       *
1673  *                                                                         *
1674  *  If the variable terms are different, not covered.                      *
1675  *  First compare if all the signs match. This is a cheaper test that will *
1676  *  find some non-matching pairs.If the signs match then start the final   *
1677  *  (expensive) check; that is to check all the elements.  If all the      *
1678  *  elements match, variable terms are identical.                          *
1679  *  If the coefficients of the variables are identical, check the constant *
1680  *  terms.                                                                 *
1681  ***************************************************************************/
1682 boolean named_sc_fm::is_already_covered(sc_fm_constraint &newc,
1683                                          sc_fm_constraint &oldc)
1684 {
1685     CALL_STAT(stat_fm_is_already1);
1686
1687 #ifndef USE_SPEED_HACK
1688     boolean found = TRUE;
1689     for(int in=1; (in<n())&&(found); in++)
1690         if(oldc.sign_no_chk(in) != newc.sign_no_chk(in))
1691             found = FALSE;
1692 #else
1693     int * tsg = &oldc.sgn[1];
1694     int * csg = &newc.sgn[1];
1695     int * bound1 = &oldc.sgn[n()-1];
1696     while(tsg<=bound1)
1697         if(*tsg++ != *csg++)  goto false1;
1698     goto label1;
1699  false1:
1700     return FALSE;
1701  label1:
1702 #endif
1703
1704
1705 #ifndef USE_SPEED_HACK
1706     for(int ip=0; (ip<p())&&(found); ip++)  {
1707         for(in=1; in<n(); in++)
1708             if(oldc[ip][in] != newc[ip][in])
1709                 found = FALSE;
1710     }
1711     if(found) return is_already_covered_chk_const(newc, oldc);
1712
1713     return FALSE;
1714 #else
1715     int * pt = &oldc[0][1];
1716     int * pc = &newc[0][1];
1717     int * step =  &oldc[1][0];
1718     int * bound2 = &oldc[p()-1][n()-1] + 1;
1719     while(1)
1720         if(pt == step) {
1721             if(pt == bound2)
1722                 return is_already_covered_chk_const(newc, oldc);
1723             pc++;
1724             pt++;
1725             step += n();
1726         }
1727         else if(*pc++ != *pt++) return FALSE;
1728 #endif
1729 }
1730
1731
1732 /***************************************************************************
1733  *                                                                         *
1734  *  A + oc >=0 and A + nc >= 0                                             *
1735  *  for all plains            oc <= nc    -> old covers new                *
1736  *  there exist a plane s.t.  oc > nc     -> old does not cover new        *
1737  ***************************************************************************/
1738 boolean named_sc_fm::is_already_covered_chk_const(sc_fm_constraint &newc,
1739                                                    sc_fm_constraint &oldc)
1740 {
1741     CALL_STAT(stat_fm_is_already2);
1742
1743     for(int ip=0; ip<p(); ip++)  {
1744         int nv = newc[ip][0];
1745         int ov = oldc[ip][0];
1746         if(ov > nv) return FALSE;
1747     }
1748     return TRUE;
1749 }
1750
1751
1752
1753 /***************************************************************************
1754  * remove the i-th inequality iff del[i]!=0                                *
1755  ***************************************************************************/
1756 void named_sc_fm::remove(int * del)
1757 {
1758     for(int i=m()-1; i>=0; i--)
1759         if(del[i]) remove(i);
1760 }
1761
1762 /***************************************************************************
1763  * remove the i-th inequality                                              *
1764  ***************************************************************************/
1765 void named_sc_fm::remove(int i)
1766 {
1767     CALL_STAT(stat_fm_remove);
1768
1769     assert(i<m());
1770     if(i != m()-1)
1771         swap(i, m()-1);
1772     xm--;
1773 }
1774
1775
1776 /***************************************************************************
1777  *  error handling                                                         *
1778  ***************************************************************************/
1779 void named_sc_fm::error()
1780 {
1781     if(_named_sc_fm_longjmp_env)
1782         longjmp(*_named_sc_fm_longjmp_env, 1);
1783     print(stderr);
1784     assert(0);
1785 }
1786
1787
1788 /***************************************************************************
1789  *  Print the current set of inequalities.                                 *
1790  ***************************************************************************/
1791 void named_sc_fm::print(FILE *fp)
1792 {
1793     fprintf(fp,"Internal:\n");
1794     named_symcoeff_ineq * tmpi = internal_get();
1795     tmpi->print_exp(pet_system_nl,fp);
1796     fprintf(fp,"Results:\n");
1797     named_symcoeff_ineq * tmpr = get();
1798     tmpr->print_exp(pet_system_nl,fp);
1799     delete tmpi;
1800     delete tmpr;
1801 }
1802
1803
1804 void named_sc_fm::delete_all()
1805 {
1806 #ifdef STAT_ON
1807     print_stat();
1808 #endif
1809
1810     if(cm > 0) {
1811         for(int i=0; i<cm; i++) delete L[i];
1812         if(L) delete L;
1813         if(magic_list) delete[] magic_list;
1814     }
1815     L = NULL;
1816     magic_list = NULL;
1817     cm = 0; cn = 0; cp = 0;
1818 }
1819
1820 #ifdef STAT_ON
1821 static void print_stat()
1822 {
1823     fflush(stderr);
1824     fprintf(stderr, "\nStatistics of Fourier Motzkin\n");
1825     fprintf(stderr, "=============================\n");
1826     fprintf(stderr, "\n**sc_fm_constraint**\n");
1827     fprintf(stderr, "       constructors =%5d,",  stat_fc_cons);
1828     fprintf(stderr, "             recalc =%6d,",  stat_fc_recalc);
1829     fprintf(stderr, "         div_by_gcd =%6d\n", stat_fc_div_by_gcd);
1830     fprintf(stderr, "            inverse =%5d,",  stat_fc_inverse);
1831     fprintf(stderr, "  is_only_pos_const =%6d,",  stat_fc_is_only_pos_const);
1832     fprintf(stderr, "       is_coef_eq_1 =%6d\n", stat_fc_is_coef_eq_1);
1833     fprintf(stderr, "\n**sc_fm_results**\n");
1834     fprintf(stderr, "       constructors =%5d,",  stat_fr_cons);
1835     fprintf(stderr, "               init =%5d,",  stat_fr_init);
1836     fprintf(stderr, "                set =%5d\n", stat_fr_set);
1837     fprintf(stderr, "\n**named_sc_fm**\n");
1838     fprintf(stderr, "       constructors =%5d,",  stat_fm_cons);
1839     fprintf(stderr, "       extend block =%5d,",  stat_fm_extend_block);
1840     fprintf(stderr, "                get =%5d\n", stat_fm_get);
1841     fprintf(stderr, "               sort =%5d,",  stat_fm_sort);
1842     fprintf(stderr, "               swap =%5d,",  stat_fm_swap);
1843     fprintf(stderr, "             remove =%5d\n", stat_fm_remove);
1844
1845     fprintf(stderr, "\n"
1846             "      STEP           all =%6d,    range =%6d,   single =%6d\n"
1847             "                 section =%6d,     pair =%6d,      lin =%6d\n",
1848            stat_fm_step0,
1849            stat_fm_step1,
1850            stat_fm_step2,
1851            stat_fm_step3,
1852            stat_fm_step4,
1853            stat_fm_lin_step);
1854     fprintf(stderr, "\n"
1855             "    BOUNDS           all =%6d,    range =%6d,   single =%6d       chk =%6d\n",
1856            stat_fm_bounds0,
1857            stat_fm_bounds1,
1858            stat_fm_bounds2,
1859            stat_fm_bounds3);
1860     fprintf(stderr, "\n"
1861             "ALREADY IN         magic =%6d,  (total)=%10ld,  vars =%7d,  const =%7d\n",
1862            stat_fm_is_already0,
1863            stat_fm_is_already0x,
1864            stat_fm_is_already1,
1865            stat_fm_is_already2);
1866     fprintf(stderr,
1867             "                     passed magic_num test      = %6.2f%%\n"
1868             "                     magic_num test was correct = %6.2f%%\n",
1869            (double)stat_fm_is_already1/(double)stat_fm_is_already0x*100.0,
1870            (double)stat_fm_is_already2/(double)stat_fm_is_already1*100.0);
1871
1872     fprintf(stderr, "\n");
1873
1874     stat_fc_cons                = 0;
1875     stat_fc_recalc              = 0;
1876     stat_fc_div_by_gcd          = 0;
1877     stat_fc_inverse             = 0;
1878     stat_fc_is_only_pos_const   = 0;
1879     stat_fc_is_coef_eq_1        = 0;
1880     stat_fr_cons                = 0;
1881     stat_fr_init                = 0;
1882     stat_fr_set                 = 0;
1883     stat_fm_cons                = 0;
1884     stat_fm_extend_block        = 0;
1885     stat_fm_step0               = 0;
1886     stat_fm_step1               = 0;
1887     stat_fm_step2               = 0;
1888     stat_fm_step3               = 0;
1889     stat_fm_step4               = 0;
1890     stat_fm_lin_step            = 0;
1891     stat_fm_bounds0             = 0;
1892     stat_fm_bounds1             = 0;
1893     stat_fm_bounds2             = 0;
1894     stat_fm_bounds3             = 0;
1895     stat_fm_is_already0         = 0;
1896     stat_fm_is_already0x        = 0;
1897     stat_fm_is_already1         = 0;
1898     stat_fm_is_already2         = 0;
1899     stat_fm_sort                = 0;
1900     stat_fm_swap                = 0;
1901     stat_fm_remove              = 0;
1902     stat_fm_get                 = 0;
1903     stat_fm_check_valid         = 0;
1904
1905 }
1906
1907 #endif