text: // Copyright (C) 2002, International Business Machines
  // Corporation and others.  All Rights Reserved.
  // This code is licensed under the terms of the Eclipse Public License (EPL).
  
  #include <cstdlib>
  #include <cstdio>
  #include <cmath>
  #include <cfloat>
  #include <cassert>
  #include <iostream>
  //#define CGL_DEBUG 1
  #ifdef NDEBUG
  //#undef NDEBUG
  #endif
  #include "CoinPragma.hpp"
  #include "CoinHelperFunctions.hpp"
  #include "CoinPackedVector.hpp"
  #include "CoinPackedMatrix.hpp"
  #include "CoinIndexedVector.hpp"
  #include "OsiRowCutDebugger.hpp"
  #ifndef USE_CGL_RATIONAL
  #define USE_CGL_RATIONAL 0
  #endif
  // -1 no rational, 0 as before, 1 use code from CglGmi
  #if USE_CGL_RATIONAL>0
  #if USE_CGL_RATIONAL<=10
  #error "If USE_CGL_RATIONAL>0 must be at least 10 (maybe not more than 1000)" 
  #endif
  #include "CoinRational.hpp"
  #endif
  #define CGL_HAS_CLP_GOMORY
  #ifdef CGL_HAS_CLP_GOMORY
  #include "OsiClpSolverInterface.hpp"
  #endif
  #include "CoinFactorization.hpp"
  #undef CLP_OSL
  #if COINUTILS_BIGINDEX_IS_INT
  #define CLP_OSL 1
  #if CLP_OSL!=1&&CLP_OSL!=3
  #undef CLP_OSL
  #else
  #include "CoinOslFactorization.hpp"
  #endif
  #endif
  #include "CoinWarmStartBasis.hpp"
  #include "CglGomory.hpp"
  #include "CoinFinite.hpp"
  #ifdef CGL_DEBUG_GOMORY
  int gomory_try=CGL_DEBUG_GOMORY;
  #endif
  //-------------------------------------------------------------------
  // Generate Gomory cuts
  //------------------------------------------------------------------- 
  void CglGomory::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
  			     const CglTreeInfo info)
  {
  #ifdef CGL_DEBUG_GOMORY
    gomory_try++;
  #endif
    // Get basic problem information
    int numberColumns=si.getNumCols(); 
    
    // get integer variables and basis
    char * intVar = new char[numberColumns];
    int i;
    CoinWarmStart * warmstart = si.getWarmStart();
    CoinWarmStartBasis* warm =
      dynamic_cast<CoinWarmStartBasis*>(warmstart);
    const double * colUpper = si.getColUpper();
    const double * colLower = si.getColLower();
    //#define CLP_INVESTIGATE2
  #ifndef CLP_INVESTIGATE2
    if ((info.options&16)!=0)
  #endif
       printf("%d %d %d\n",info.inTree,info.options,info.pass);
    for (i=0;i<numberColumns;i++) {
      if (si.isInteger(i)) {
        if (colUpper[i]>colLower[i]+0.5) {
  	if (fabs(colUpper[i]-1.0)<1.0e-12&&fabs(colLower[i])<1.0e-12) {
  	  intVar[i]=1; //0-1
  	} else  if (colLower[i]>=0.0) {
  	  intVar[i] = 2; // other
  	} else {
  	  // negative bounds - I am not sure works
  	  intVar[i] = 3;
  	}
        } else {
  	intVar[i] = 4;
        }
      } else {
        intVar[i]=0;
      }
    }
    const OsiSolverInterface * useSolver=&si;
  #ifdef CGL_HAS_CLP_GOMORY
    double * objective = NULL;
    OsiClpSolverInterface * clpSolver = dynamic_cast<OsiClpSolverInterface *>(originalSolver_);
    int numberOriginalRows = -1;
    if (clpSolver) {
      useSolver = originalSolver_;
      assert (gomoryType_);
      // check simplex is plausible
      if (!clpSolver->getNumRows()||numberColumns!=clpSolver->getNumCols()) {
        delete originalSolver_;
        originalSolver_=si.clone();
        clpSolver = dynamic_cast<OsiClpSolverInterface *>(originalSolver_);
        assert (clpSolver);
        useSolver = originalSolver_;
      }
      ClpSimplex * simplex = clpSolver->getModelPtr();
      numberOriginalRows = simplex->numberRows();
      int numberRows = si.getNumRows();
      assert (numberOriginalRows<=numberRows);
      // only do if different (unless type 2x)
      int gomoryType = gomoryType_%10;
      int whenToDo = gomoryType_/10;
      if (whenToDo==2 ||(numberRows>numberOriginalRows && whenToDo==1
  		       && (info.options&512)==0) ||
  	((info.options&1024)!=0 && (info.options&512)==0
  	 && numberTimesStalled_<3)) {
        // bounds
        memcpy(simplex->columnLower(),colLower,numberColumns*sizeof(double));
        memcpy(simplex->columnUpper(),colUpper,numberColumns*sizeof(double));
        double * obj = simplex->objective();
        objective = CoinCopyOfArray(obj,numberColumns);
        const double * pi = si.getRowPrice();
        const CoinPackedMatrix * rowCopy = si.getMatrixByRow();
        const int * column = rowCopy->getIndices();
        const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
        const int * rowLength = rowCopy->getVectorLengths(); 
        const double * rowElements = rowCopy->getElements();
        const double * rowLower = si.getRowLower();
        const double * rowUpper = si.getRowUpper();
        int numberCopy;
        int numberAdd;
        double * rowLower2 = NULL;
        double * rowUpper2 = NULL;
        int * column2 = NULL;
        CoinBigIndex * rowStart2 = NULL;
        double * rowElements2 = NULL;
        char * copy = new char [numberRows-numberOriginalRows];
        memset(copy,0,numberRows-numberOriginalRows);
        if (gomoryType==2) {
  	numberCopy=0;
  	numberAdd=0;
  	for (int iRow=numberOriginalRows;iRow<numberRows;iRow++) {
  	  bool simple = true;
  	  for (CoinBigIndex k=rowStart[iRow];
  	       k<rowStart[iRow]+rowLength[iRow];k++) {
  	    double value = rowElements[k];
  	    if (value!=floor(value+0.5)) {
  	      simple=false;
  	      break;
  	    }
  	  }
  	  if (simple) {
  	    numberCopy++;
  	    numberAdd+=rowLength[iRow];
  	    copy[iRow-numberOriginalRows]=1;
  	  }
  	}
  	if (numberCopy) {
  	  //printf("Using %d rows out of %d\n",numberCopy,numberRows-numberOriginalRows);
  	  rowLower2 = new double [numberCopy];
  	  rowUpper2 = new double [numberCopy];
  	  rowStart2 = new CoinBigIndex [numberCopy+1];
  	  rowStart2[0]=0;
  	  column2 = new int [numberAdd];
  	  rowElements2 = new double [numberAdd];
  	}
        }
        numberCopy=0;
        numberAdd=0;
        const double * rowSolution = si.getRowActivity();
        double offset=0.0;
        for (int iRow=numberOriginalRows;iRow<numberRows;iRow++) {
  	if (!copy[iRow-numberOriginalRows]) {
  	  double value = pi[iRow];
  	  offset += rowSolution[iRow]*value;
  	  for (CoinBigIndex k=rowStart[iRow];
  	       k<rowStart[iRow]+rowLength[iRow];k++) {
  	    int iColumn=column[k];
  	    obj[iColumn] -= value*rowElements[k];
  	  }
  	} else {
  	  rowLower2[numberCopy]=rowLower[iRow];
  	  rowUpper2[numberCopy]=rowUpper[iRow];
  	  for (CoinBigIndex k=rowStart[iRow];
  	       k<rowStart[iRow]+rowLength[iRow];k++) {
  	    column2[numberAdd]=column[k];
  	    rowElements2[numberAdd++]=rowElements[k];
  	  }
  	  numberCopy++;
  	  rowStart2[numberCopy]=numberAdd;
  	}
        }
  #if 0
        CoinThreadRandom randomNumberGenerator;
        const double * solution = si.getColSolution();
        for (int i=0;i<numberColumns;i++) {
  	if (intVar[i]==1) {
  	  double randomNumber = randomNumberGenerator.randomDouble();
  	  //randomNumber = 0.001*floor(randomNumber*1000.0);
  	  if (solution[i]>0.5)
  	    obj[i] -= randomNumber*0.001*fabs(obj[i]);
  	  else
  	    obj[i] += randomNumber*0.001*fabs(obj[i]);
  	}
        }
  #endif
        if (numberCopy) {
  	clpSolver->addRows(numberCopy,
  			   rowStart2,column2,rowElements2,
  			   rowLower2,rowUpper2);
  	delete [] rowLower2 ;
  	delete [] rowUpper2 ;
  	delete [] column2 ;
  	delete [] rowStart2 ;
  	delete [] rowElements2 ;
        }
        delete [] copy;
        memcpy(simplex->primalColumnSolution(),si.getColSolution(),
  	     numberColumns*sizeof(double));
        warm->resize(numberOriginalRows,numberColumns);
        clpSolver->setBasis(*warm);
        delete warm;
        simplex->setDualObjectiveLimit(COIN_DBL_MAX);
        simplex->setLogLevel(0);
        simplex->primal(1);
        // check basis
        int numberTotal=simplex->numberRows()+simplex->numberColumns();
        int superbasic=0;
        for (int i=0;i<numberTotal;i++) {
  	if (simplex->getStatus(i)==ClpSimplex::superBasic)
  	  superbasic++;
        }
        if (superbasic) {
  	//printf("%d superbasic!\n",superbasic);
  	simplex->dual();
  	superbasic=0;
  	for (int i=0;i<numberTotal;i++) {
  	  if (simplex->getStatus(i)==ClpSimplex::superBasic)
  	    superbasic++;
  	}
  	assert (!superbasic);
        }
        //printf("Trying - %d its status %d objs %g %g - with offset %g\n",
        //     simplex->numberIterations(),simplex->status(),
        //     simplex->objectiveValue(),si.getObjValue(),simplex->objectiveValue()+offset);
        //simplex->setLogLevel(0);
        warm=simplex->getBasis();
        warmstart=warm;
        if (simplex->status()) {
  	//printf("BAD status %d\n",simplex->status());
  	//clpSolver->writeMps("clp");
  	//si.writeMps("si");
  	delete [] objective;
  	objective=NULL;
  	useSolver=&si;
        }
      } else {
        // don't do
        delete warmstart;
        warmstart=NULL;
        if ((info.options&1024)==0)
  	numberTimesStalled_=0;
      }
    }
  #endif
  #ifdef CGL_DEBUG
    const OsiRowCutDebugger * debugger = si.getRowCutDebugger();
    if (debugger&&!debugger->onOptimalPath(si))
      debugger = NULL;
  #else
    const OsiRowCutDebugger * debugger = NULL;
  #endif
    int numberRowCutsBefore = cs.sizeRowCuts();
  
    if (warmstart)
      generateCuts(debugger, cs, *useSolver->getMatrixByCol(), 
  		 *useSolver->getMatrixByRow(),
  		 useSolver->getColSolution(),
  		 useSolver->getColLower(), useSolver->getColUpper(), 
  		 useSolver->getRowLower(), useSolver->getRowUpper(),
  		 intVar,warm,info);
  #ifdef CGL_HAS_CLP_GOMORY
    if (objective) {
      ClpSimplex * simplex = clpSolver->getModelPtr();
      memcpy(simplex->objective(),objective,numberColumns*sizeof(double));
      delete [] objective;
      // take out locally useless cuts
      const double * solution = si.getColSolution();
      double primalTolerance = 1.0e-7;
      int numberRowCutsAfter = cs.sizeRowCuts();
      for (int k = numberRowCutsAfter - 1; k >= numberRowCutsBefore; k--) {
        const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
        double sum = 0.0;
        int n = thisCut->row().getNumElements();
        const int * column = thisCut->row().getIndices();
        const double * element = thisCut->row().getElements();
        assert (n);
        for (int i = 0; i < n; i++) {
  	double value = element[i];
  	sum += value * solution[column[i]];
        }
        if (sum > thisCut->ub() + primalTolerance) {
  	sum = sum - thisCut->ub();
        } else if (sum < thisCut->lb() - primalTolerance) {
  	sum = thisCut->lb() - sum;
        } else {
  	sum = 0.0;
        }
        if (!sum) {
  	// take out
  	cs.eraseRowCut(k);
        }
      }
  #ifdef CLP_INVESTIGATE2
      printf("OR %p pass %d inTree %c - %d cuts (but %d deleted)\n",
         originalSolver_,info.pass,info.inTree?'Y':'N',
         numberRowCutsAfter-numberRowCutsBefore,
         numberRowCutsAfter-cs.sizeRowCuts());
  #endif
    }
  #endif
  
    delete warmstart;
    delete [] intVar;
    if ((!info.inTree&&((info.options&4)==4||((info.options&8)&&!info.pass)))
        ||(info.options&16)!=0) {
      int limit = maximumLengthOfCutInTree();
      int numberRowCutsAfter = cs.sizeRowCuts();
      for (int i=numberRowCutsBefore;i<numberRowCutsAfter;i++) {
        int length = cs.rowCutPtr(i)->row().getNumElements();
        if (length<=limit)
  	cs.rowCutPtr(i)->setGloballyValid();
      }
    }
    if ((gomoryType_%10)==2) {
      // back to original
      assert(clpSolver);
      int numberRows = clpSolver->getNumRows();
      if (numberRows>numberOriginalRows) {
        int numberDelete = numberRows-numberOriginalRows;
        int * delRow = new int [numberDelete];
        for (int i=0;i<numberDelete;i++)
  	delRow[i]=i+numberOriginalRows;
        clpSolver->deleteRows(numberDelete,delRow);
        delete [] delRow;
      }
    }
  }
  
  // Returns value - floor but allowing for small errors
  inline double above_integer(double value) {
    double value2=floor(value);
    double value3=floor(value+0.5);
    if (fabs(value3-value)<1.0e-9*(fabs(value3)+1.0))
      return 0.0;
    return value-value2;
  }
  #if USE_CGL_RATIONAL <= 0
  //-------------------------------------------------------------------
  // Returns the greatest common denominator of two 
  // positive integers, a and b, found using Euclid's algorithm 
  //-------------------------------------------------------------------
  static long long int gcd(long long int a, long long int b) 
  {
    long long int remainder = -1;
  #if CGL_DEBUG>1
    printf("gcd of %ld and %ld\n",a,b);
    int nLoop=0;
  #endif
    // make sure a<=b (will always remain so)
    if(a > b) {
      // Swap a and b
      long long int temp = a;
      a = b;
      b = temp;
    }
    // if zero then gcd is nonzero (zero may occur in rhs of packed)
    if (!a) {
      if (b) {
        return b;
      } else {
        printf("**** gcd given two zeros!!\n");
        abort();
      }
    }
    while (remainder) {
  
  #if CGL_DEBUG>1
      nLoop++;
      if (nLoop>50) {
        abort();
        return -1;
      }
  #endif
      remainder = b % a;
      b = a;
      a = remainder;
    }
  #if CGL_DEBUG>1
    printf("=> %d\n",b);
  #endif
    return b;
  }
  #endif
  #define GOMORY_LONG
  #ifdef GOMORY_LONG
  #define SMALL_VALUE1 1.0e-14
  #define GOMORY_INT long long int
  #else
  #define SMALL_VALUE1 1.0e-10
  #define GOMORY_INT int
  #endif
  #if USE_CGL_RATIONAL>0
  static long computeGcd(long a, long b) {
    // This is the standard Euclidean algorithm for gcd
    long remainder = 1;
    // Make sure a<=b (will always remain so)
    if (a > b) {
      // Swap a and b
      long temp = a;
      a = b;
      b = temp;
    }
    // If zero then gcd is nonzero
    if (!a) {
      if (b) {
        return b;
      } 
      else {
        printf("### WARNING: CglGMI::computeGcd() given two zeroes!\n");
        exit(1);
      }
    }
    while (remainder) {
      remainder = b % a;
      b = a;
      a = remainder;
    }
    return b;
  } /* computeGcd */
  static bool scaleCutIntegral(double* cutElem, int* cutIndex, int cutNz,
  			     double& cutRhs, double maxdelta) {
    long gcd, lcm;
    double maxscale = 1000; 
    long maxdnom = USE_CGL_RATIONAL;
    //long numerator = 0, denominator = 0;
    // Initialize gcd and lcm
    CoinRational r = CoinRational(cutRhs, maxdelta, maxdnom);
    if (r.getNumerator() != 0){
       gcd = labs(r.getNumerator());
       lcm = r.getDenominator();
    }
    else{
  #if defined GMI_TRACE_CLEAN
        printf("Cannot compute rational number, scaling procedure aborted\n");
  #endif
      return false;
    }
    for (int i = 0; i < cutNz; ++i) {
      CoinRational r = CoinRational(cutElem[i], maxdelta, maxdnom);
      if (r.getNumerator() != 0){
         gcd = computeGcd(gcd, r.getNumerator());
         lcm *= r.getDenominator()/(computeGcd(lcm,r.getDenominator()));
      }
      else{
  #if defined GMI_TRACE_CLEAN
        printf("Cannot compute rational number, scaling procedure aborted\n");
  #endif
        return false;
      } 
    }
    double scale = ((double)lcm)/((double)gcd);
    if (fabs(scale) > maxscale) {
  #if defined GMI_TRACE_CLEAN
        printf("Scaling factor too large, scaling procedure aborted\n");
  #endif
        return false;
    }
    scale = fabs(scale);
    // Looks like we have a good scaling factor; scale and return;
    for (int i = 0; i < cutNz; ++i) {
      double value = cutElem[i]*scale;
      cutElem[i] = floor(value+0.5);
      assert (fabs(cutElem[i]-value)<1.0e-9);
    }
    {
      double value = cutRhs*scale;
      cutRhs = floor(value+0.5);
      assert (fabs(cutRhs-value)<1.0e-9);
    }
    return true;
  } /* scaleCutIntegral */
  #elif USE_CGL_RATIONAL == 0
  //-------------------------------------------------------------------
  // Returns the nearest rational with denominator < maxDenominator
  //-------------------------------------------------------------------
  typedef struct {
    int numerator;
    int denominator;
  } Rational;
  typedef struct {
    GOMORY_INT numerator;
    GOMORY_INT denominator;
  } RationalLong;
  inline Rational nearestRational(double value, int maxDenominator) 
  {
    RationalLong tryThis;
    RationalLong tryA;
    RationalLong tryB;
    double integerPart;
  
  #if CGL_DEBUG>1
    printf("Rational of %g is ",value);
  #endif
    int nLoop=0;
  
    tryA.numerator=0;
    tryA.denominator=1;
    tryB.numerator=1;
    tryB.denominator=0;
    Rational returnRational;
  
    if (fabs(value)<SMALL_VALUE1) {
      returnRational.numerator=static_cast<int>(tryA.numerator);
      returnRational.denominator=static_cast<int>(tryA.denominator);
      return returnRational;
    }
    integerPart = floor(value);
    value -= integerPart;
    tryThis.numerator = tryB.numerator* static_cast<GOMORY_INT> (integerPart) + tryA.numerator;
    tryThis.denominator = tryB.denominator* static_cast<GOMORY_INT> (integerPart) + tryA.denominator;
    tryA = tryB;
    tryB = tryThis;
  
    while (value>SMALL_VALUE1 && tryB.denominator <=maxDenominator) {
      nLoop++;
      if (nLoop>50) {
        Rational bad;
        bad.numerator=-1;
        bad.denominator=-1;
  #if CGL_DEBUG>1
        printf(" *** bad rational\n");
  #endif
        return bad;
      }
      value = 1.0/value;
      integerPart = floor(value+SMALL_VALUE1);
      value -= integerPart;
      tryThis.numerator = tryB.numerator* static_cast<GOMORY_INT> (integerPart) + tryA.numerator;
      tryThis.denominator = tryB.denominator* static_cast<GOMORY_INT>(integerPart) + tryA.denominator;
      tryA = tryB;
      tryB = tryThis;
    }
    if (tryB.denominator <= maxDenominator) {
  #if CGL_DEBUG>1
      printf("%lld/%lld\n",tryB.numerator,tryB.denominator);
  #endif
      if (tryB.denominator<COIN_INT_MAX) {
        returnRational.numerator=static_cast<int>(tryB.numerator);
        returnRational.denominator=static_cast<int>(tryB.denominator);
      } else {
        returnRational.numerator=-1;
        returnRational.denominator=-1;
  #if CGL_DEBUG>1
        printf(" *** bad rational\n");
  #endif
      }
    } else {
  #if CGL_DEBUG>1
      printf("%lld/%lld\n",tryA.numerator,tryA.denominator);
  #endif
      if (tryA.denominator<COIN_INT_MAX) {
        returnRational.numerator=static_cast<int>(tryA.numerator);
        returnRational.denominator=static_cast<int>(tryA.denominator);
      } else {
        returnRational.numerator=-1;
        returnRational.denominator=-1;
  #if CGL_DEBUG>1
        printf(" *** bad rational\n");
  #endif
      }
    }
    return returnRational;
  }
  #endif
  // Does actual work - returns number of cuts
  int
  CglGomory::generateCuts( 
  #ifdef CGL_DEBUG
  			const OsiRowCutDebugger * debugger, 
  #else
  			const OsiRowCutDebugger * , 
  #endif
                           OsiCuts & cs,
                           const CoinPackedMatrix & columnCopy,
                           const CoinPackedMatrix & rowCopy,
                           const double * colsol,
                           const double * colLower, const double * colUpper,
                           const double * rowLower, const double * rowUpper,
  			 const char * intVar,
                           const CoinWarmStartBasis* warm,
                           const CglTreeInfo info)
  {
    int infoOptions=info.options;
    bool globalCuts = (infoOptions&16)!=0;
    double testFixed = (!globalCuts) ? 1.0e-8 : -1.0;
    // get what to look at
    double away = info.inTree ? away_ : CoinMin(away_,awayAtRoot_);
    int numberRows=columnCopy.getNumRows();
    int numberColumns=columnCopy.getNumCols(); 
    CoinBigIndex numberElements=columnCopy.getNumElements();
    // Allow bigger length on initial matrix (if special setting)
    //if (limit==512&&!info.inTree&&!info.pass)
    //limit=1024;
    // Start of code to create a factorization from warm start (A) ====
    // check factorization is okay
    CoinFactorization factorization;
  #ifdef CLP_OSL
    CoinOslFactorization * factorization2=NULL;
    if (alternateFactorization_) {
      factorization2 = new CoinOslFactorization();
    }
  #endif
    // We can either set increasing rows so ...IsBasic gives pivot row
    // or we can just increment iBasic one by one
    // for now let ...iBasic give pivot row
    int status=-100;
    // probably could use pivotVariables from OsiSimplexModel
    int * rowIsBasic = new int[numberRows];
    int * columnIsBasic = new int[numberColumns];
    int i;
    int numberBasic=0;
    for (i=0;i<numberRows;i++) {
      if (warm->getArtifStatus(i) == CoinWarmStartBasis::basic) {
        rowIsBasic[i]=1;
        numberBasic++;
      } else {
        rowIsBasic[i]=-1;
      }
    }
    for (i=0;i<numberColumns;i++) {
      if (warm->getStructStatus(i) == CoinWarmStartBasis::basic) {
        columnIsBasic[i]=1;
        numberBasic++;
      } else {
        columnIsBasic[i]=-1;
      }
    }
    //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
    while (status<-98) {
  #ifdef CLP_OSL
      if (!alternateFactorization_) {
  #endif
        status=factorization.factorize(columnCopy,
  				     rowIsBasic, columnIsBasic);
        if (status==-99) factorization.areaFactor(factorization.areaFactor() * 2.0);
  #ifdef CLP_OSL
      } else {
        double areaFactor=1.0;
        status=factorization2->factorize(columnCopy,
  				      rowIsBasic, columnIsBasic,areaFactor);
        if (status==-99) areaFactor *= 2.0;
      }
  #endif
    } 
    if (status) {
  #ifdef COIN_DEVELOP
      std::cout<<"Bad factorization of basis - status "<<status<<std::endl;
  #endif
      delete [] rowIsBasic;
      delete [] columnIsBasic;
      return -1;
    }
    // End of creation of factorization (A) ====
    
  #ifdef CLP_OSL
    double relaxation = !alternateFactorization_ ? factorization.conditionNumber() :
      factorization2->conditionNumber();
  #else
    double relaxation = factorization.conditionNumber();
  #endif
    // if very small be a bit more careful
    if (relaxation<1.0e-10)
      relaxation=1.0/sqrt(relaxation);
  #ifdef COIN_DEVELOP_z
    if (relaxation>1.0e49)
      printf("condition %g\n",relaxation);
  #endif
    //printf("condition %g %g\n",relaxation,conditionNumberMultiplier_);
    relaxation *= conditionNumberMultiplier_;
    double bounds[2]={-COIN_DBL_MAX,0.0};
    int iColumn,iRow;
  
    const int * column = rowCopy.getIndices();
    const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
    const int * rowLength = rowCopy.getVectorLengths(); 
    const double * rowElements = rowCopy.getElements();
    const int * row = columnCopy.getIndices();
    const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
    const int * columnLength = columnCopy.getVectorLengths(); 
    const double * columnElements = columnCopy.getElements();
  
    // we need to do book-keeping for variables at ub
    double tolerance = 1.0e-7;
    bool * swap= new bool [numberColumns];
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      if (columnIsBasic[iColumn]<0&&
  	colUpper[iColumn]-colsol[iColumn]<=tolerance) {
        swap[iColumn]=true;
      } else {
        swap[iColumn]=false;
      }
    }
  
    // get row activities (could use solver but lets do here )
    double * rowActivity = new double [numberRows];
    CoinFillN(rowActivity,numberRows,0.0);
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      double value = colsol[iColumn];
      CoinBigIndex k;
      for (k=columnStart[iColumn];k<columnStart[iColumn]+columnLength[iColumn];k++) {
        iRow = row[k];
        rowActivity[iRow] += columnElements[k]*value;
      }
    }
    /* we need to mark rows:
       0) equality
       1) slack at lb (activity at ub)
       2) slack at ub (activity at lb)
       and 4 bit is set if slack must be integer
  
    */
    int * rowType = new int[numberRows];
    for (iRow=0;iRow<numberRows;iRow++) {
      if (rowIsBasic[iRow]<0&&rowUpper[iRow]>rowLower[iRow]+1.0e-7) {
        int type=0;
        double rhs=0.0;
        if (rowActivity[iRow]>=rowUpper[iRow]-1.0e-7) {
  	type=1;
  	rhs=rowUpper[iRow];
        } else if (rowActivity[iRow]<=rowLower[iRow]+1.0e-7) {
  	type=2;
  	rhs=rowLower[iRow];
        } else {
  	// probably large rhs
  	if (rowActivity[iRow]-rowLower[iRow]<
  	    rowUpper[iRow]-rowActivity[iRow])
  	  rowType[iRow]=2;
  	else
  	  rowType[iRow]=1;
  #ifdef CGL_DEBUG
  	assert (CoinMin(rowUpper[iRow]-rowActivity[iRow],
  		    rowActivity[iRow]-rowUpper[iRow])<1.0e-5);
  	//abort();
          continue;
  #else
  	continue;
  #endif
        }
        if (above_integer(rhs)<1.0e-10) {
  	// could be integer slack
  	bool allInteger=true;
  	CoinBigIndex k;
  	for (k=rowStart[iRow];
  	     k<rowStart[iRow]+rowLength[iRow];k++) {
  	  int iColumn=column[k];
  	  if (!intVar[iColumn]||above_integer(rowElements[k])>1.0e-10) {
  	    // not integer slacks
  	    allInteger=false;
  	    break;
  	  }
  	}
  	if (allInteger) {
  	  type |= 4;
  	}
        }
        rowType[iRow]=type;
      } else {
        // row is equality or basic
        rowType[iRow]=0;
      }
    }
  
    // Start of code to create work arrays for factorization (B) ====
    // two vectors for updating (one is work)
    CoinIndexedVector work;
    CoinIndexedVector array;
    // make sure large enough
    work.reserve(numberRows);
    array.reserve(numberRows);
    int * arrayRows = array.getIndices();
    double * arrayElements = array.denseVector();
    // End of code to create work arrays (B) ====
  
    int numberAdded=0;
    // we also need somewhere to accumulate cut
    CoinIndexedVector cutVector;
    cutVector.reserve(numberColumns+1);
    int * cutIndex = cutVector.getIndices();
    double * cutElement = cutVector.denseVector(); 
    // and for packed form (as not necessarily in order)
    // also space for sort
    bool doSorted = (infoOptions&256)!=0;
    int lengthArray = static_cast<int>(numberColumns+1+((numberColumns+1)*sizeof(int))/sizeof(double));
    if (doSorted)
      lengthArray+=numberColumns;
    double * packed = new double[lengthArray]; 
    double * sort = packed+numberColumns+1;
    int * which = reinterpret_cast<int *>(doSorted ? (sort+numberColumns): (sort));
    double tolerance1=1.0e-6;
    double tolerance2=0.9;
  #if USE_CGL_RATIONAL <= 0
    double tolerance3=1.0e-10;
  #endif
    double tolerance6=1.0e-6;
    double tolerance9=1.0e-4;
  #define MORE_GOMORY_CUTS 1
  #ifdef CLP_INVESTIGATE2
    int saveLimit = info.inTree ? 50 : 1000;
  #else
  #if MORE_GOMORY_CUTS==2||MORE_GOMORY_CUTS==3
    int saveLimit;
  #endif  
  #endif  
    // get limit on length of cut
    int limit = 0;
    if (!limit_)
      dynamicLimitInTree_ = CoinMax(50,numberColumns/40);
    if (!info.inTree) {
      limit = limitAtRoot_;
      if (!info.pass) {
        tolerance1=1.0;
        tolerance2=1.0e-2;
  #if USE_CGL_RATIONAL <= 0
        //tolerance3=1.0e-6;
  #endif
        tolerance6=1.0e-7;
        tolerance9=1.0e-5;
        if (!limit)
  	limit=numberColumns;
      } else {
        if((infoOptions&32)==0/*&&numberTimesStalled_<3*/) {
  	if (!limit) {
  	  if(numberElements>8*numberColumns)
  	    limit=numberColumns;
  	  else
  	    limit = CoinMax(1000,numberRows/4);
  	}
        } else {
  	//limit=numberColumns;
  	//numberTimesStalled_++;
        } 
      }
    } else {
      limit = limit_;
      if (!limit) {
        if (!info.pass) 
  	limit = dynamicLimitInTree_;
        else
  	limit=50;
      }
    }
    // If big - allow for rows
    if (limit>=numberColumns)
      limit += numberRows;
  #ifdef CLP_INVESTIGATE2
    if (limit>saveLimit&&!info.inTree&&(infoOptions&512)==0) 
      printf("Gomory limit changed from %d to %d, inTree %c, pass %d, r %d,c %d,e %d\n",
  	   saveLimit,limit,info.inTree ? 'Y' : 'N',info.pass,
  	   numberRows,numberColumns,numberElements);
  #endif
    int nCandidates=0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      // This returns pivot row for columns or -1 if not basic (C) ====
      int iBasic=columnIsBasic[iColumn];
      if (iBasic>=0&&intVar[iColumn]) {
        double reducedValue=above_integer(colsol[iColumn]);
        //printf("col %d bas %d val %.18g\n",iColumn,iBasic,colsol[iColumn]);
        if(intVar[iColumn]&&reducedValue<1.0-away&&reducedValue>away) {
  	if (doSorted)
  	  sort[nCandidates]=fabs(0.5-reducedValue);
  	which[nCandidates++]=iColumn;
        }
      }
    }
    CoinBigIndex nTotalEls=COIN_INT_MAX;
    if (doSorted) {
      CoinSort_2(sort,sort+nCandidates,which);
      CoinBigIndex nElsNow = columnCopy.getNumElements();
      CoinBigIndex nAdd;
      CoinBigIndex nAdd2;
      CoinBigIndex nReasonable;
      int depth=info.level;
      if (depth<2) {
        nAdd=10000;
        if (info.pass>0)
  	nAdd = CoinMin(nAdd,nElsNow+2*numberRows);
        nAdd2 = 5*numberColumns;
        nReasonable = CoinMax(nAdd2,nElsNow/8+nAdd);
        if (!depth&&!info.pass) {
  	// allow more
  	nAdd += nElsNow/2;
  	nAdd2 += nElsNow/2;
  	nReasonable += nElsNow/2;
  	limit=numberRows+numberColumns;
        }
      } else {
        nAdd = 200;
        nAdd2 = 2*numberColumns;
        nReasonable = CoinMax(nAdd2,nElsNow/8+nAdd);
      }
      nTotalEls=nReasonable;
    }
  #ifdef MORE_GOMORY_CUTS
    CoinBigIndex saveTotalEls=nTotalEls;
  #endif
  #if MORE_GOMORY_CUTS==2||MORE_GOMORY_CUTS==3
    saveLimit=limit;
    if (doSorted)
      limit=numberRows+numberColumns;
    OsiCuts longCuts;
  #endif
  #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
    OsiCuts secondaryCuts;
  #endif
    for (int kColumn=0;kColumn<nCandidates;kColumn++) {
      if (nTotalEls<=0)
        break;  // Got enough
      iColumn=which[kColumn];
      double reducedValue=above_integer(colsol[iColumn]);;
      // This returns pivot row for columns or -1 if not basic (C) ====
      int iBasic=columnIsBasic[iColumn];
      double ratio=reducedValue/(1.0-reducedValue);
      if (iBasic>=0) {
    // Debug code below computes tableau column of basic ====
        int j;
  #ifdef CGL_DEBUG
        {
  	// put column into array
  	array.setVector(columnLength[iColumn],row+columnStart[iColumn],
  			columnElements+columnStart[iColumn]);
  	// get column in tableau
  #ifdef CLP_OSL
  	if (!alternateFactorization_)
  #endif
  	  factorization.updateColumn ( &work, &array );
  #ifdef CLP_OSL
  	else
  	  factorization2->updateColumn ( &work, &array );
  #endif
  	int nn=0;
  	int numberInArray=array.getNumElements();
  	for (j=0;j<numberInArray;j++) {
  	  int indexValue=arrayRows[j];
  	  double value=arrayElements[indexValue];
  	  if (fabs(value)>1.0e-5) {
  	    assert (fabs(value-1.0)<1.0e-7);
  	    assert (indexValue==iBasic);
  	    nn++;
  	  }
  	}
  	assert (nn==1);
  	array.clear();
  	//work.checkClear();
        }
  #endif
        array.clear();
        assert(intVar[iColumn]&&reducedValue<1.0-away&&reducedValue>away);
        {
  #ifdef CGL_DEBUG
  	//cutVector.checkClear();
  #endif
  	// get row of tableau
  	double one =1.0;
  	array.setVector(1,&iBasic,&one);
  	int numberNonInteger=0;
  	//Code below computes tableau row ====
  	// get pi
  #ifdef CLP_OSL
  	if (!alternateFactorization_)
  #endif
  	  factorization.updateColumnTranspose ( &work, &array );
  #ifdef CLP_OSL
  	else
  	  factorization2->updateColumnTranspose ( &work, &array );
  #endif
  	int numberInArray=array.getNumElements();
  #ifdef CGL_DEBUG
  	// check pivot on iColumn
  	{
  	  double value=0.0;
  	  int k;
  	  // add in row of tableau
  	  for (k=columnStart[iColumn];
  	       k<columnStart[iColumn]+columnLength[iColumn];k++) {
  	    iRow = row[k];
  	    value += columnElements[k]*arrayElements[iRow];
  	  }
  	  // should be 1
  	  assert (fabs(value-1.0) < 1.0e-7);
  	}
  #endif
  	double largestFactor=0.0;
  	for (j=0;j<numberInArray;j++) {
  	  int indexValue=arrayRows[j];
  	  double value=arrayElements[indexValue];
  	  largestFactor = CoinMax(largestFactor,fabs(value));
  	}
  	//reducedValue=colsol[iColumn];
  	// coding from pg 130 of Wolsey 
  	// adjustment to rhs
  	double rhs=0.0;
  	int number=0;
  	// number of terms (so includes modifications and cancellations)
  	int numberCoefficients=0;
  #ifdef CGL_DEBUG_GOMORY
  	    if (!gomory_try)
  	      printf("start for basic column %d\n",iColumn);
  #endif
  	// columns
  	for (j=0;j<numberColumns;j++) {
  	  if (columnIsBasic[j]<0&&colUpper[j]>colLower[j]+testFixed) {
  	    double value=0.0;
  	    CoinBigIndex k;
  	    // add in row of tableau
  	    for (k=columnStart[j];k<columnStart[j]+columnLength[j];k++) {
  	      iRow = row[k];
  	      double value2 = columnElements[k]*arrayElements[iRow];
  	      largestFactor = CoinMax(largestFactor,fabs(value2));
  	      value += value2;
  	    }
  #ifdef CGL_DEBUG_GOMORY
  	    if (!gomory_try&&value)
  	      printf("col %d alpha %g colsol %g swap %c bounds %g %g\n",
  		     j,value,colsol[j],swap[j] ? 'Y' : 'N',
  		     colLower[j],colUpper[j]);
  #endif
  	    // value is entry in tableau row end (C) ====
  	    if (fabs(value)<1.0e-16) {
  	      // small value
  	      continue;
  	    } else {
  	      // left in to stop over compilation?
  	      //if (iColumn==-52) printf("for basic %d, column %d has alpha %g, colsol %g\n",
  	      //		      iColumn,j,value,colsol[j]);
  #if CGL_DEBUG>1
  	      if (iColumn==52) printf("for basic %d, column %d has alpha %g, colsol %g\n",
  				      iColumn,j,value,colsol[j]);
  #endif
  	      // deal with bounds
  	      if (swap[j]) {
  		//reducedValue -= value*colUpper[j];
  		// negate
  		value = - value;
  	      } else {
  		//reducedValue -= value*colLower[j];
  	      }
  #if CGL_DEBUG>1
  	      if (iColumn==52) printf("%d value %g reduced %g int %d rhs %g swap %d\n",
  				      j,value,reducedValue,intVar[j],rhs,swap[j]);
  #endif
  	      double coefficient;
  	      if (intVar[j]) {
  		// integer
  		coefficient = above_integer(value);
  		if (coefficient > reducedValue) {
  		  coefficient = ratio * (1.0-coefficient);
  		} 
  	      } else {
  		// continuous
  		numberNonInteger++;
  		if (value > 0.0) {
  		  coefficient = value;
  		} else {
  		  //??? sign wrong in book
  		  coefficient = -ratio*value;
  		}
  	      }
  	      if (swap[j]) {
  		// negate
  		coefficient = - coefficient;
  		rhs += colUpper[j]*coefficient;
  	      } else {
  		rhs += colLower[j]*coefficient;
  	      }
  	      if (fabs(coefficient)>= COIN_INDEXED_TINY_ELEMENT) {
  		cutElement[j] = coefficient;
  		numberCoefficients++;
  		cutIndex[number++]=j;
  		// If too many - break from loop
  		if (number>limit) 
  		  break;
  	      }
  	    }
  	  } else {
  	    // basic
  	    continue;
  	  }
  	}
  	cutVector.setNumElements(number);
  	// If too many - just clear vector and skip
  	if (number>limit) {
  	  cutVector.clear();
  	  continue;
  	}
  	//check will be cut
  	//reducedValue=above_integer(reducedValue);
  	rhs += reducedValue;
  	double violation = reducedValue;
  #ifdef CGL_DEBUG
  	std::cout<<"cut has violation of "<<violation
  		 <<" value "<<colsol[iColumn]<<std::endl;
  #endif
  	// now do slacks part
  	for (j=0;j<numberInArray;j++) {
  	  iRow=arrayRows[j];
  	  double value = arrayElements[iRow];
  	  int type=rowType[iRow];
  	  if (type&&fabs(value)>=1.0e-16) {
  	    if ((type&1)==0) {
  	      // negate to get correct coefficient
  	      value = - value;
  	    }
  	    double coefficient;
  	    if ((type&4)!=0) {
  	      // integer
  	      coefficient = above_integer(value);
  	      if (coefficient > reducedValue) {
  		coefficient = ratio * (1.0-coefficient);
  	      } 
  	    } else {
  	      numberNonInteger++;
  	      // continuous
  	      if (value > 0.0) {
  		coefficient = value;
  	      } else {
  		coefficient = -ratio*value;
  	      }
  	    }
  	    if ((type&1)!=0) {
  	      // slack at ub - treat as +1.0
  	      rhs -= coefficient*rowUpper[iRow];
  	    } else {
  	      // negate yet again ?
  	      coefficient = - coefficient;
  	      rhs -= coefficient*rowLower[iRow];
  	    }
  	    CoinBigIndex k;
  	    for (k=rowStart[iRow];
  		 k<rowStart[iRow]+rowLength[iRow];k++) {
  	      int jColumn=column[k];
  	      double value=rowElements[k];
  	      double oldValue=cutElement[jColumn];
  	      cutVector.quickAdd(jColumn,-coefficient*value);
  	      numberCoefficients++;
  	      if (!intVar[jColumn]&&!oldValue) 
  		numberNonInteger++;
  	    }
  	  }
  	}
  	//check again and pack down
  	// also change signs
  	// also zero cutElement
  	double sum=0.0;
  	rhs = - rhs;
  	int n = cutVector.getNumElements();
  	//assert (n); can be 0!
  	// If too many - just clear vector and skip
  	if (n>limit) {
  	  cutVector.clear();
  	  continue;
  	}
  #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
  	double violation2=violation;
  #endif
  	number=0;
  	numberNonInteger=0;
  	double sumCoefficients=0.0;
  	for (j=0;j<n;j++) {
  	  int jColumn =cutIndex[j];
  	  double value=-cutElement[jColumn];
  	  sumCoefficients += fabs(value);
  	  if (fabs(colsol[jColumn])>10.0)
  	    sumCoefficients += 2.0*fabs(value);
  	  cutElement[jColumn]=0.0;
  	  if (fabs(value)>1.0e-8) {
  	    sum+=value*colsol[jColumn];
  	    packed[number]=value;
  	    cutIndex[number++]=jColumn;
  	    if (!intVar[jColumn])
  	      numberNonInteger++;
            } else {
  #define LARGE_BOUND 1.0e20
              // small - adjust rhs if rhs reasonable
              if (value>0.0&&colLower[jColumn]>-LARGE_BOUND) {
                rhs -= value*colLower[jColumn];
  #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
  	      // weaken violation
  	      violation2 -= fabs(value*(colsol[jColumn]-colLower[jColumn]));
  #endif
              } else if (value<0.0&&colUpper[jColumn]<LARGE_BOUND) {
                rhs -= value*colUpper[jColumn];
  #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
  	      // weaken violation
  	      violation2 -= fabs(value*(colsol[jColumn]-colUpper[jColumn]));
  #endif
              } else if (fabs(value)>1.0e-13) {
                // take anyway
                sum+=value*colsol[jColumn];
                packed[number]=value;
                cutIndex[number++]=jColumn;
  	      if (!intVar[jColumn])
  		numberNonInteger++;
              } 
            }
  	}
  	// Final test on number
  	//if (number>limit)
  	//continue;
  	// say zeroed out
  	cutVector.setNumElements(0);
  	bool accurate2=false;
  	double difference=fabs((sum-rhs)-violation);
  	double useTolerance;
  	if (tolerance1>0.99) {
  	  // use absolute
  	  useTolerance = tolerance;
  	} else {
  	  double rhs2=CoinMax(fabs(rhs),10.0);
  	  useTolerance=rhs2*0.1*tolerance1;
  	}
  	bool accurate = (difference<useTolerance);
  #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
  	double difference2=fabs((sum-rhs)-violation2);
  #if MORE_GOMORY_CUTS==1
  	if (difference2<useTolerance&&doSorted) 
  #else
  	if (difference2<useTolerance&&doSorted&&number<saveLimit) 
  #endif
  	  accurate2=true;
  #endif
  	if (sum >rhs+tolerance2*away&&
  	    (accurate||accurate2)) {
  	  //#ifdef CGL_DEBUG
  #ifdef CGL_DEBUG
  #if CGL_DEBUG<=1
  	  if (number<=10) {
  #endif
  	    std::cout<<"col "<<iColumn;
  	    for (j=0;j<number;j++) {
  	      std::cout<<" ["<<cutIndex[j]<<","<<packed[j]<<"]";
  	    }
  	    std::cout<<" <= "<<rhs<<std::endl;
  #if CGL_DEBUG<=1
  	  }
  #endif
  #endif
  	  bool cleanedCut=numberNonInteger>0;
  	  bool dontRelax = false;
  #define TRY7 1
  #if TRY7==2
  	  double rhsBeforeRelax=rhs;
  #endif
  	  if ((!numberNonInteger||number<6)&&number&&USE_CGL_RATIONAL>=0) {
  	    // pretend not integer
  	    numberNonInteger=0;
  #if USE_CGL_RATIONAL==0
  #ifdef CGL_DEBUG
  	    assert (sizeof(Rational)==sizeof(double));
  #endif
  	    Rational * cleaned = reinterpret_cast<Rational *> (cutElement);
  	    int * xInt = reinterpret_cast<int *> (cutElement);
  	    // cut should have an integer slack so try and simplify
  	    // add in rhs and put in cutElements (remember to zero later)
  	    cutIndex[number]=numberColumns+1;
  	    packed[number]=rhs;
  	    int numberNonSmall=0;
  	    long long int lcm = 1;
  	    
  	    for (j=0;j<number+1;j++) {
  	      double value=above_integer(fabs(packed[j]));
  	      if (fabs(value)<tolerance3) {
  		// very close to integer
  		continue;
  	      } else {
  		numberNonSmall++;
  	      }
  	      
  	      cleaned[j]=nearestRational(value,100000);
  	      if (cleaned[j].denominator<0) {
  		// bad rational
  		lcm=-1;
  		break;
  	      }
  	      long long int thisGcd = gcd(lcm,cleaned[j].denominator);
  	      // may need to check for overflow
  	      lcm /= thisGcd;
  	      if (static_cast<double>(lcm)*cleaned[j].denominator<
  		  1.0e18) {
  		lcm *= cleaned[j].denominator;
  	      } else {
  		lcm=-1;
  		break;
  	      }
  	    }
  	    if (lcm>0&&numberNonSmall) {
  	      double multiplier = static_cast<double>(lcm);
  	      cleanedCut=true;
  	      int nOverflow = 0; 
  	      for (j=0; j<number+1;j++) {
  		double value = fabs(packed[j]);
  		double dxInt = value*multiplier;
  		xInt[j]= static_cast<int> (dxInt+0.5); 
  #if CGL_DEBUG>1
  		printf("%g => %g   \n",value,dxInt);
  #endif
  		if (dxInt>1.0e9||fabs(dxInt-xInt[j])> 1.0e-8) {
  		  nOverflow++;
  		  break;
  		}
  	      }
  	      
  	      if (nOverflow){
  #ifdef CGL_DEBUG
  		printf("Gomory Scaling: Warning: Overflow detected \n");
  #endif
  		numberNonInteger=-1;
  	      } else {
  		
  		// find greatest common divisor of the elements
  		j=0;
  		while (!xInt[j])
  		  j++; // skip zeros
  		long long int thisGcd = gcd(xInt[j],xInt[j+1]);
  		j++;
  		for (;j<number+1;j++) {
  		  if (xInt[j])
  		    thisGcd = gcd(thisGcd,xInt[j]);
  		}
  #if 0
  		// Check nothing too illegal - FIX this
  		for (j=0;j<number+1;j++) {
  		  double old = lcm*packed[j];
  		  int newOne;
  		  if (old>0.0)
  		    newOne=xInt[j]/thisGcd;
  		  else
  		    newOne=-xInt[j]/thisGcd;
  		  if (fabs(((double) newOne)-old)>
  		      1.0e-10*(fabs(newOne)+1.0)) {
  		    // say no good - first see if happens
  		    printf("Fix this test 456 - just skip\n");
  		    abort();
  		  }
  		} 
  #endif		  
  #if CGL_DEBUG>1
  		printf("The gcd of xInt is %i\n",thisGcd);    
  #endif
  		
  		dontRelax=true;
  		// construct new cut by dividing through by gcd and 
  		double minMultiplier=1.0e100;
  		double maxMultiplier=0.0;
  		for (j=0; j<number+1;j++) {
  		  double old=packed[j];
  		  if (old>0.0) {
  		    packed[j]=static_cast<double>(xInt[j]/thisGcd);
  		  } else {
  		    packed[j]=-static_cast<double>(xInt[j]/thisGcd);
  		  }
  #if CGL_DEBUG>1
  		  printf("%g => %g   \n",old,packed[j]);
  #endif
  		  if (packed[j]) {
  		    if (fabs(packed[j])>maxMultiplier*fabs(old))
  		      maxMultiplier = packed[j]/old;
  		    if (fabs(packed[j])<minMultiplier*fabs(old))
  		      minMultiplier = packed[j]/old;
  		  }
  		}
  		rhs = packed[number];
  #ifdef CGL_DEBUG
  		printf("min, max multipliers - %g, %g\n",
  		       minMultiplier,maxMultiplier);
  #endif
  		assert(maxMultiplier/minMultiplier>0.9999&&maxMultiplier/minMultiplier<1.0001);
  	      }
  	    } else if (!numberNonSmall) {
  	      dontRelax=true;
  	    }
  	    // erase cutElement
  	    CoinFillN(cutElement,number+1,0.0);
  #elif USE_CGL_RATIONAL>0
  	    // Use coding from CglGmi
  	    double maxdelta = 1.0e-12; //param.getEPS();
  	    if (0) {
  	      double total=0.0;
  	      for (j=0;j<number;j++) {
  		int jColumn =cutIndex[j];
  		double value=packed[j];
  		std::cout<<"("<<jColumn<<","<<value<<","<<colsol[jColumn]
  			 <<") ";
  		total += value*colsol[jColumn];
  	      }
  	      std::cout<<std::endl;
  	      printf("before %g <= %g <= %g\n",bounds[0],total,rhs);
  	    }
  	    dontRelax=scaleCutIntegral(packed,cutIndex,number,
  				       rhs,maxdelta);
  	    if (dontRelax&&false) {
  	      double total=0.0;
  	      for (j=0;j<number;j++) {
  		int jColumn =cutIndex[j];
  		double value=packed[j];
  		std::cout<<"("<<jColumn<<","<<value<<","<<colsol[jColumn]
  			 <<") ";
  		total += value*colsol[jColumn];
  	      }
  	      std::cout<<std::endl;
  	      printf("after %g <= %g <= %g\n",bounds[0],total,rhs);
  	    }
  #endif
  	  }
  	  if (!dontRelax) {
  	    // relax rhs a tiny bit
  	    //#define CGL_GOMORY_OLD_RELAX
  #ifndef CGL_GOMORY_OLD_RELAX
  #if 0
  	    double rhs2=rhs;
  	    rhs2 += 1.0e-8;
  	    // relax if lots of elements for mixed gomory
  	    if (number>=20) {
  	      rhs2  += 1.0e-7*(static_cast<double> (number/20));
  	    }
  #endif
  	    rhs += 1.0e-7;
  	    if (numberCoefficients>=10||true) {
  	      rhs  += 1.0e-7*sumCoefficients+1.0e-8*numberCoefficients;
  	    }
  #if 0
  	    if (numberCoefficients>number*3)
  	    printf("old rhs %.18g new %.18g - n,nNon,nC,sumC %d,%d,%d %g\n",
  		   rhs2,rhs,number,numberNonInteger,numberCoefficients,
  		   sumCoefficients);
  #endif
  #else
  	    rhs += 1.0e-8;
  	    // relax if lots of elements for mixed gomory
  	    if (number>=20) {
  	      rhs  += 1.0e-7*(static_cast<double> (number/20));
  	    }
  #endif
  	  }
  	  // Take off tiny elements
  	  // for first pass reject
  #ifndef CGL_GOMORY_TINY_ELEMENT
  #define CGL_GOMORY_TINY_ELEMENT 1.0e-12
  #endif
  	  {
  	    int i,number2=number;
  	    number=0;
  	    double largest=0.0;
  	    double smallest=1.0e30;
  	    for (i=0;i<number2;i++) {
  	      double value=fabs(packed[i]);
  	      if (value<CGL_GOMORY_TINY_ELEMENT) {
  		int iColumn = cutIndex[i];
  		if (colUpper[iColumn]-colLower[iColumn]<LARGE_BOUND) {
  		  if (intVar[iColumn])
  		    numberNonInteger--;
  		  // weaken cut
  		  if (packed[i]>0.0) 
  		    rhs -= value*colLower[iColumn];
  		  else
  		    rhs += value*colUpper[iColumn];
  		} else {
  		  // throw away
  		  number=limit+1;
  		  numberNonInteger=1;
  		  break;
  		}
  	      } else {
  		int iColumn = cutIndex[i];
  		if (colUpper[iColumn]!=colLower[iColumn]||globalCuts) {
  		  value=fabs(value);
  		  largest=CoinMax(largest,value);
  		  smallest=CoinMin(smallest,value);
  		  cutIndex[number]=cutIndex[i];
  		  packed[number++]=packed[i];
  		} else {
  		  // fixed so subtract out
  		  rhs -= packed[i]*colLower[iColumn];
  		  if (intVar[iColumn])
  		    numberNonInteger--;
  		}
  	      }
  	    }
  	    if (largest>1.0e10*smallest||(number>20&&smallest<number*1.0e-6)||
  		numberNonInteger<-10||!number) {
  	      //assert (number); // debug this - looks as if it could be
  	      number=limit+1; //reject
  	      numberNonInteger=1;
  	    } else if (largest>1.0e9*smallest) {
  #ifdef CLP_INVESTIGATE2
  	      printf("WOuld reject %g %g ratio %g\n",smallest,largest,
  		     smallest/largest);
  #endif
  #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
  	      accurate=false;
  #endif
  	    } else {
  #define PRINT_NUMBER 0
  #if PRINT_NUMBER
  	      if (number==PRINT_NUMBER) {
  		printf("==========\n<= %.18g ",rhs);
  		for (int i=0;i<PRINT_NUMBER;i++)
  		  printf("%.18g ",packed[i]);
  		printf("\n");
  	      }
  #endif
  	      if (number>limit)  
  		continue;
  #if TRY7==1
  	      // Just scale - unless small integers
  	      bool dontScale=true;
  	      if (fabs(rhs-floor(rhs+0.5))>1.0e-14)
  		dontScale=false;
  	      for (int i=0;i<number;i++) {
  		if (fabs(packed[i]-floor(packed[i]+0.5))>1.0e-14||
  		    fabs(packed[i])>1.0e4)
  		  dontScale=false;
  	      }
  	      if (!dontScale) {
  		double multiplier = 1.0/sqrt(largest*smallest);
  		for (int i=0;i<number;i++)
  		  packed[i] *= multiplier;
  		rhs *= multiplier;
  	      }
  #if PRINT_NUMBER
  	      if (number==PRINT_NUMBER) {
  		printf("multiplier %g %g %g\n",
  		       multiplier,smallest,largest);
  	      }
  #endif
  #elif TRY7==2
  	      // Look at ratio
  	      double scaleFactor=fabs(rhsBeforeRelax);
  	      for (int i=0;i<number;i++) {
  		double value=packed[i];
  		double ratio = fabs(rhs/value);
  		double nearest=floor(ratio+0.5);
  		if (fabs(ratio-nearest)<1.0e-6 && nearest >= 1.0) {
  		  int iColumn=cutIndex[i];
  		  if (intVar[iColumn]) {
  		    if (colLower[iColumn]>=0.0) {
  		      if (value>0.0) {
  			// better if smaller
  			if (ratio>nearest) {
  			  packed[i]=scaleFactor/nearest;
  			}			  
  		      } else {
  			// better if larger
  			if (ratio<nearest) {
  			  packed[i]=-scaleFactor/nearest;
  			}			  
  		      }
  		      //if (value!=packed[i])
  		      //printf("column %d rhs %.18g element %g ratio %.18g - new element %.18g\n",
  		      //       iColumn,rhs,value,ratio,packed[i]);
  		      //assert (fabs(value-packed[i])<1.0e-4);
  		    } else {
  		      //printf("column %d rhs %g element %g ratio %.18g - bounds %g,%g\n",
  			   //  iColumn,rhs,value,ratio,colLower[iColumn],colUpper[iColumn]);
  		    }
  		  }
  		}
  	      }
  #endif
  #if PRINT_NUMBER
  	      if (number==PRINT_NUMBER) {
  		printf("after %.18g ",rhs);
  		for (int i=0;i<PRINT_NUMBER;i++)
  		  printf("%.18g ",packed[i]);
  		printf("\n");
  	      }
  #endif
  	    }
  	  }
  	  if (number<limit||!numberNonInteger) {
  	    bounds[1]=rhs;
  	    if (number>50&&numberNonInteger)
  	      bounds[1] = rhs+tolerance6+1.0e-8*fabs(rhs); // weaken
  #if GOMORY_RELAX_NUMBER
  	    else if (number>GOMORY_RELAX_NUMBER&&numberNonInteger>1)
  	      bounds[1] = rhs+tolerance6+1.0e-8*fabs(rhs); // weaken
  #endif
  	    // if close to integer - round up
  	    double nearest=floor(bounds[1]+0.5);
  	    if (bounds[1]<nearest&&bounds[1]>nearest-1.0e-4)
  	      bounds[1]=nearest;
  	    double test = CoinMin(largestFactor*largestFactorMultiplier_,
  				  relaxation);
  	    if (number>5&&numberNonInteger&&test>1.0e-20) {
  #ifdef CLP_INVESTIGATE2
  	      printf("relaxing rhs by %g - largestFactor was %g, rel %g\n",
  	         CoinMin(test*fabs(rhs),tolerance9),largestFactor,relaxation);
  #endif
  	      //bounds[1] = CoinMax(bounds[1],
  	      //		  rhs+CoinMin(test*fabs(rhs),tolerance9)); // weaken
  	      bounds[1] = bounds[1]+CoinMin(test*fabs(rhs),tolerance9); // weaken
  	    }
  #ifdef MORE_GOMORY_CUTS
  	    if (accurate) {
  #else
  	    {
  #endif
  	      // tidy
  	      if (!cleanedCut) {
  		double range=0.0;
  		for (int k=0;k<number;k++) {
  		  int iColumn=cutIndex[k];
  		  double thisRange=CoinMin(colUpper[iColumn]-colLower[iColumn],1000.0);
  		  range += thisRange;
  		}
  		// see if close to integer
  		bool close=fabs(bounds[1]-floor(bounds[1]+0.5))*range<1.0e-6;
  		if (close) {
  		  for (int k=0;k<number;k++) {
  		    if(fabs(packed[k]-floor(packed[k]+0.5))*range>1.0e-6) {
  		      close=false;
  		      break;
  		    }
  		  }
  		  if (close) {
  #ifdef PRINT_MORE
  		    bool printIt=false;
  		    for (int k=0;k<number;k++) {
  		      if(fabs(packed[k]-floor(packed[k]+0.5))>1.0e-12) {
  		        printIt=true;
  		        break;
  		      }
  		    }
  		    if (printIt) {
  		      printf("yy %.18g >= ",bounds[1]);
    		      for (int k=0;k<number;k++) {
  		        printf("(%d,%.18g) ",cutIndex[k],packed[k]);
  		      }
  		      printf("\n");
  		    }
  #endif
  		    bounds[1]=floor(bounds[1]+0.5);
  		    for (int k=0;k<number;k++) {
  		      packed[k]=floor(packed[k]+0.5);
  		    }
  		  }
  		}
  	      }
  	      OsiRowCut rc;
  	      rc.setRow(number,cutIndex,packed,false);
  	      rc.setLb(bounds[0]);
  	      rc.setUb(bounds[1]);
  #if 0
  	      double total=0.0;
  	      for (j=0;j<number;j++) {
  		int jColumn =cutIndex[j];
  		double value=packed[j];
  		std::cout<<"("<<jColumn<<","<<value<<","<<colsol[jColumn]
  			 <<") ";
  		total += value*colsol[jColumn];
  	      }
  	      std::cout<<std::endl;
  	      printf("above CUT %g <= %g <= %g\n",bounds[0],total,bounds[1]);
  #endif	      
  #ifdef CGL_DEBUG
  	      if (debugger) {
  		assert(!debugger->invalidCut(rc));
  		if(debugger->invalidCut(rc))
  		  abort();
  	      }
  #endif
  #if MORE_GOMORY_CUTS<2
  	      nTotalEls -= number;
  	      cs.insertIfNotDuplicate(rc);
  #else
  	      if(number<saveLimit) {
  		nTotalEls -= number;
  		cs.insertIfNotDuplicate(rc);
  	      } else {
  		longCuts.insertIfNotDuplicate(rc);
  	      }
  #endif
  	      //printf("nTot %d kCol %d iCol %d ibasic %d\n",
  	      //     nTotalEls,kColumn,iColumn,iBasic);
  	      numberAdded++;
  #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
  	    } else if (accurate2) {
  	      OsiRowCut rc;
  	      rc.setRow(number,cutIndex,packed,false);
  	      rc.setLb(bounds[0]);
  	      rc.setUb(bounds[1]);   
  	      secondaryCuts.insertIfNotDuplicate(rc);
  #endif
  	    }
  	  } else {
  #ifdef CGL_DEBUG
  	    std::cout<<"cut has "<<number<<" entries - skipped"
  		     <<std::endl;
  	    if (!number)
  	      std::cout<<"******* Empty cut - infeasible"<<std::endl;
  #endif
  	  }
  	} else {
  	  // why dropped?
  #ifdef CGL_DEBUG
  	  std::cout<<"first violation "<<violation<<" now "
  		   <<sum-rhs<<" why?, rhs= "<<rhs<<std::endl;
  	  
  	  for (j=0;j<number;j++) {
  	    int jColumn =cutIndex[j];
  	    double value=packed[j];
  	    std::cout<<"("<<jColumn<<","<<value<<","<<colsol[jColumn]
  		     <<") ";
  	  }
  	  std::cout<<std::endl;
  	  //abort();
  #endif
  	}
        }
      } else {
        // not basic
  #if CGL_DEBUG>1
        {
  	// put column into array
  	array.setVector(columnLength[iColumn],row+columnStart[iColumn],
  			columnElements+columnStart[iColumn]);
  	// get column in tableau
  #ifdef CLP_OSL
  	if (!alternateFactorization_)
  #endif
  	  factorization.updateColumn ( &work, &array );
  #ifdef CLP_OSL
  	else
  	  factorization2->updateColumn ( &work, &array );
  #endif
  	int numberInArray=array.getNumElements();
  	printf("non-basic %d\n",iColumn);
  	for (int j=0;j<numberInArray;j++) {
  	  int indexValue=arrayRows[j];
  	  double value=arrayElements[indexValue];
  	  if (fabs(value)>1.0e-6) {
  	    printf("%d %g\n",indexValue,value);
  	  }
  	}
        }
  #endif
      }
    }
  #ifdef CLP_OSL
    delete factorization2;
  #endif
  
    delete [] rowActivity;
    delete [] swap;
    delete [] rowType;
    delete [] packed;
    delete [] rowIsBasic;
    delete [] columnIsBasic;
  #ifdef MORE_GOMORY_CUTS
  #if MORE_GOMORY_CUTS==1
    int numberInaccurate = secondaryCuts.sizeRowCuts();
  #ifdef CLP_INVESTIGATE2
    int numberOrdinary = numberAdded-numberInaccurate;
    if (!info.inTree&&(infoOptions&512)==0) 
      printf("Gomory has %d ordinary and %d less accurate cuts(%d els)\n",
  	   numberOrdinary,numberInaccurate,saveTotalEls-nTotalEls);
  #endif
  #elif MORE_GOMORY_CUTS==2
    int numberLong = longCuts.sizeRowCuts();
  #ifdef CLP_INVESTIGATE2
    int numberOrdinary = numberAdded-numberLong;
    if (!info.inTree&&(infoOptions&512)==0) 
      printf("Gomory has %d ordinary and %d long cuts(%d els)\n",
  	   numberOrdinary,numberLong,saveTotalEls-nTotalEls);
  #endif
  #elif MORE_GOMORY_CUTS==3
    int numberLong = longCuts.sizeRowCuts();
    int numberInaccurate = secondaryCuts.sizeRowCuts();
  #ifdef CLP_INVESTIGATE2
    int numberOrdinary = numberAdded-numberLong-numberInaccurate;
    if (!info.inTree&&(infoOptions&512)==0) 
      printf("Gomory has %d ordinary, %d long and %d less accurate cuts(%d els)\n",
  	   numberOrdinary,numberLong,numberInaccurate,saveTotalEls-nTotalEls);
  #endif
  #endif
    if (doSorted&&limit<numberColumns) {
      // Just half
      nTotalEls -= saveTotalEls/2;
  #if MORE_GOMORY_CUTS==2||MORE_GOMORY_CUTS==3
      while (nTotalEls>0) {
        for (int i=0;i<numberLong;i++) {
  	nTotalEls -= longCuts.rowCutPtr(i)->row().getNumElements();
  	cs.insertIfNotDuplicate(longCuts.rowCut(i));
  	numberAdded ++;
  	if (nTotalEls<=0)
  	  break;
        }
        break;
      }
  #endif
  #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
      while (nTotalEls>0) {
        for (int i=0;i<numberInaccurate;i++) {
  	nTotalEls -= secondaryCuts.rowCutPtr(i)->row().getNumElements();
  	cs.insertIfNotDuplicate(secondaryCuts.rowCut(i));
  	numberAdded ++;
  	if (nTotalEls<=0)
  	  break;
        }
        break;
      }
  #endif
    }
  #else
  #ifdef CLP_INVESTIGATE2
    if (!info.inTree&&(infoOptions&512)==0) 
      printf("Gomory added %d cuts(%d els)\n",numberAdded,saveTotalEls-nTotalEls);
  #endif
  #endif
    return numberAdded;
  }
  // Limit stuff
  void CglGomory::setLimit(int limit)
  {
    if (limit>=0)
      limit_=limit; 
  }
  int CglGomory::getLimit() const
  {
    return limit_;
  }
  // Limit stuff at root
  void CglGomory::setLimitAtRoot(int limit)
  {
    if (limit>=0)
      limitAtRoot_=limit;
  }
  int CglGomory::getLimitAtRoot() const
  {
    return limitAtRoot_;
  }
  // Return maximum length of cut in tree
  int 
  CglGomory::maximumLengthOfCutInTree() const
  {
    if (limit_)
      return limit_;
    else
      return dynamicLimitInTree_;
  }
  
  // Away stuff
  void CglGomory::setAway(double value)
  {
    if (value>0.0&&value<=0.5)
      away_=value;
  }
  double CglGomory::getAway() const
  {
    return away_;
  }
  
  // Away stuff at root
  void CglGomory::setAwayAtRoot(double value)
  {
    if (value>0.0&&value<=0.5)
      awayAtRoot_=value;
  }
  double CglGomory::getAwayAtRoot() const
  {
    return awayAtRoot_;
  }
  
  // ConditionNumberMultiplier stuff
  void CglGomory::setConditionNumberMultiplier(double value)
  {
    if (value>=0.0)
      conditionNumberMultiplier_=value;
  }
  double CglGomory::getConditionNumberMultiplier() const
  {
    return conditionNumberMultiplier_;
  }
  
  // LargestFactorMultiplier stuff
  void CglGomory::setLargestFactorMultiplier(double value)
  {
    if (value>=0.0)
      largestFactorMultiplier_=value;
  }
  double CglGomory::getLargestFactorMultiplier() const
  {
    return largestFactorMultiplier_;
  }
  
  //-------------------------------------------------------------------
  // Default Constructor 
  //-------------------------------------------------------------------
  CglGomory::CglGomory ()
  :
  CglCutGenerator(),
  away_(0.05),
  awayAtRoot_(0.05),
  conditionNumberMultiplier_(1.0e-18),
  largestFactorMultiplier_(1.0e-13),
  originalSolver_(NULL),
  limit_(50),
  limitAtRoot_(0),
  dynamicLimitInTree_(-1),
  numberTimesStalled_(0),
  alternateFactorization_(0),
  gomoryType_(0)
  {
  
  }
  
  //-------------------------------------------------------------------
  // Copy constructor 
  //-------------------------------------------------------------------
  CglGomory::CglGomory (const CglGomory & source) :
    CglCutGenerator(source),
    away_(source.away_),
    awayAtRoot_(source.awayAtRoot_),
    conditionNumberMultiplier_(source.conditionNumberMultiplier_),
    largestFactorMultiplier_(source.largestFactorMultiplier_),
    originalSolver_(NULL),
    limit_(source.limit_),
    limitAtRoot_(source.limitAtRoot_),
    dynamicLimitInTree_(source.dynamicLimitInTree_),
    numberTimesStalled_(source.numberTimesStalled_),
    alternateFactorization_(source.alternateFactorization_),
    gomoryType_(source.gomoryType_)
  { 
    if (source.originalSolver_)
      originalSolver_ = source.originalSolver_->clone();
  }
  
  //-------------------------------------------------------------------
  // Clone
  //-------------------------------------------------------------------
  CglCutGenerator *
  CglGomory::clone() const
  {
    return new CglGomory(*this);
  }
  
  //-------------------------------------------------------------------
  // Destructor 
  //-------------------------------------------------------------------
  CglGomory::~CglGomory ()
  {
    delete originalSolver_;
  }
  
  //----------------------------------------------------------------
  // Assignment operator 
  //-------------------------------------------------------------------
  CglGomory &
  CglGomory::operator=(const CglGomory& rhs)
  {
    if (this != &rhs) {
      CglCutGenerator::operator=(rhs);
      away_=rhs.away_;
      awayAtRoot_=rhs.awayAtRoot_;
      conditionNumberMultiplier_ = rhs.conditionNumberMultiplier_;
      largestFactorMultiplier_ = rhs.largestFactorMultiplier_;
      limit_=rhs.limit_;
      limitAtRoot_=rhs.limitAtRoot_;
      dynamicLimitInTree_ = rhs.dynamicLimitInTree_;
      numberTimesStalled_ = rhs.numberTimesStalled_;
      alternateFactorization_=rhs.alternateFactorization_; 
      gomoryType_ = rhs.gomoryType_;
      delete originalSolver_;
      if (rhs.originalSolver_)
        originalSolver_ = rhs.originalSolver_->clone();
      else
        originalSolver_=NULL;
    }
    return *this;
  }
  // This can be used to refresh any information
  void 
  CglGomory::refreshSolver(OsiSolverInterface * solver)
  {
    int numberColumns=solver->getNumCols(); 
    const double * colUpper = solver->getColUpper();
    const double * colLower = solver->getColLower();
    canDoGlobalCuts_ = true;
    if (originalSolver_) {
      delete originalSolver_;
      originalSolver_ = solver->clone();
    }
    for (int i=0;i<numberColumns;i++) {
      if (solver->isInteger(i)) {
        if (colUpper[i]>colLower[i]+1.0) {
  	canDoGlobalCuts_ = false;
  	break;
        }
      }
    }
  }
  // Pass in a copy of original solver (clone it)
  void 
  CglGomory::passInOriginalSolver(OsiSolverInterface * solver)
  {
    delete originalSolver_;
    if (solver) {
      if (!gomoryType_)
        gomoryType_=1;
      originalSolver_ = solver->clone();
    } else {
      gomoryType_=0;
      originalSolver_=NULL;
    }
  }
  // Does actual work - returns number of cuts
  int
  CglGomory::generateCuts( const OsiRowCutDebugger * debugger, 
                           OsiCuts & cs,
                           const CoinPackedMatrix & columnCopy,
                           const double * colsol,
                           const double * colLower, const double * colUpper,
                           const double * rowLower, const double * rowUpper,
  			 const char * intVar,
                           const CoinWarmStartBasis* warm,
                           const CglTreeInfo info)
  {
    CoinPackedMatrix rowCopy;
    rowCopy.reverseOrderedCopyOf(columnCopy);
    return generateCuts( debugger, cs, columnCopy, rowCopy,
  		       colsol, colLower, colUpper,
  		       rowLower, rowUpper, intVar, warm, info);
  }
  // Create C++ lines to get to current state
  std::string
  CglGomory::generateCpp( FILE * fp) 
  {
    CglGomory other;
    fprintf(fp,"0#include \"CglGomory.hpp\"\n");
    fprintf(fp,"3  CglGomory gomory;\n");
    if (limit_!=other.limit_)
      fprintf(fp,"3  gomory.setLimit(%d);\n",limit_);
    else
      fprintf(fp,"4  gomory.setLimit(%d);\n",limit_);
    if (limitAtRoot_!=other.limitAtRoot_)
      fprintf(fp,"3  gomory.setLimitAtRoot(%d);\n",limitAtRoot_);
    else
      fprintf(fp,"4  gomory.setLimitAtRoot(%d);\n",limitAtRoot_);
    if (away_!=other.away_)
      fprintf(fp,"3  gomory.setAway(%g);\n",away_);
    else
      fprintf(fp,"4  gomory.setAway(%g);\n",away_);
    if (awayAtRoot_!=other.awayAtRoot_)
      fprintf(fp,"3  gomory.setAwayAtRoot(%g);\n",awayAtRoot_);
    else
      fprintf(fp,"4  gomory.setAwayAtRoot(%g);\n",awayAtRoot_);
    if (getAggressiveness()!=other.getAggressiveness())
      fprintf(fp,"3  gomory.setAggressiveness(%d);\n",getAggressiveness());
    else
      fprintf(fp,"4  gomory.setAggressiveness(%d);\n",getAggressiveness());
    return "gomory";
  }
  