if (bflag==NULL || mlast!=lp->m){     if (mlast!=lp->m && mlast>0) {       free(bflag);   /* called previously with different lp->m */       free(OrderVector);     }     bflag=(long *) calloc(lp->m+1,sizeof(long*));     OrderVector=(long *)calloc(lp->m+1,sizeof(long*));      /* initialize only for the first time or when a larger space is needed */          mlast=lp->m;  }  /* Initializing control variables. */  dd_ComputeRowOrderVector2(lp->m,lp->d,lp->A,OrderVector,dd_MinIndex,rseed);  lp->re=0; lp->se=0; pivots1=0;  dd_ResetTableau(lp->m,lp->d,lp->B,lp->nbindex,bflag,lp->objrow,lp->rhscol);  dd_FindLPBasis(lp->m,lp->d,lp->A,lp->B,OrderVector,lp->equalityset,      lp->nbindex,bflag,lp->objrow,lp->rhscol,&s,&found,&(lp->LPS),&pivots0);  lp->pivots[0]=pivots0;  if (!found){     lp->se=s;     goto _L99;     /* No LP basis is found, and thus Inconsistent.       Output the evidence column. */  }  stop=dd_FALSE;  do {   /* Criss-Cross Method */#if !defined GMPRATIONAL    if (pivots1>maxpivots) {      *err=dd_LPCycling;      fprintf(stderr,"max number %ld of pivots performed by the criss-cross method. Most likely due to the floating-point arithmetics error.\n", maxpivots);      goto _L99;  /* failure due to max no. of pivots performed */    }#endif    dd_SelectCrissCrossPivot(lp->m,lp->d,lp->A,lp->B,bflag,       lp->objrow,lp->rhscol,&r,&s,&chosen,&(lp->LPS));    if (chosen) {      dd_GaussianColumnPivot2(lp->m,lp->d,lp->A,lp->B,lp->nbindex,bflag,r,s);      pivots1++;    } else {      switch (lp->LPS){        case dd_Inconsistent: lp->re=r;        case dd_DualInconsistent: lp->se=s;        default: break;      }      stop=dd_TRUE;    }  } while(!stop);  _L99:  lp->pivots[1]=pivots1;  dd_statCCpivots+=pivots1;  dd_SetSolutions(lp->m,lp->d,lp->A,lp->B,   lp->objrow,lp->rhscol,lp->LPS,&(lp->optvalue),lp->sol,lp->dsol,lp->posset_extra,lp->nbindex,lp->re,lp->se,bflag);  free(nbtemp);}void dd_SetSolutions(dd_rowrange m_size,dd_colrange d_size,   dd_Amatrix A,dd_Bmatrix T,   dd_rowrange objrow,dd_colrange rhscol,dd_LPStatusType LPS,   mytype *optvalue,dd_Arow sol,dd_Arow dsol,dd_rowset posset, dd_colindex nbindex,   dd_rowrange re,dd_colrange se,dd_rowindex bflag)/* Assign the solution vectors to sol,dsol,*optvalue after solvingthe LP.*/{  dd_rowrange i;  dd_colrange j;  mytype x,sw;  int localdebug=dd_FALSE;    dd_init(x); dd_init(sw);  switch (LPS){  case dd_Optimal:    for (j=1;j<=d_size; j++) {      dd_set(sol[j-1],T[j-1][rhscol-1]);      dd_TableauEntry(&x,m_size,d_size,A,T,objrow,j);      dd_neg(dsol[j-1],x);      dd_TableauEntry(optvalue,m_size,d_size,A,T,objrow,rhscol);      if (localdebug) {fprintf(stderr,"dsol[%ld]= ",nbindex[j]); dd_WriteNumber(stderr, dsol[j-1]); }    }    for (i=1; i<=m_size; i++) {      if (bflag[i]==-1) {  /* i is a basic variable */        dd_TableauEntry(&x,m_size,d_size,A,T,i,rhscol);        if (dd_Positive(x)) set_addelem(posset, i);      }    }    break;  case dd_Inconsistent:    if (localdebug) fprintf(stderr,"DualSimplexSolve: LP is inconsistent.\n");    for (j=1;j<=d_size; j++) {      dd_set(sol[j-1],T[j-1][rhscol-1]);      dd_TableauEntry(&x,m_size,d_size,A,T,re,j);      dd_neg(dsol[j-1],x);      if (localdebug) {fprintf(stderr,"dsol[%ld]= ",nbindex[j]); dd_WriteNumber(stderr,dsol[j-1]);}    }    break;  case dd_DualInconsistent:    for (j=1;j<=d_size; j++) {      dd_set(sol[j-1],T[j-1][se-1]);      dd_TableauEntry(&x,m_size,d_size,A,T,objrow,j);      dd_neg(dsol[j-1],x);      if (localdebug) {fprintf(stderr,"dsol[%ld]= \n",nbindex[j]);dd_WriteNumber(stderr,dsol[j-1]);}    }    if (localdebug) printf( "DualSimplexSolve: LP is dual inconsistent.\n");    break;  case dd_StrucDualInconsistent:    dd_TableauEntry(&x,m_size,d_size,A,T,objrow,se);    if (dd_Positive(x)) dd_set(sw,dd_one);    else dd_neg(sw,dd_one);    for (j=1;j<=d_size; j++) {      dd_mul(sol[j-1],sw,T[j-1][se-1]);      dd_TableauEntry(&x,m_size,d_size,A,T,objrow,j);      dd_neg(dsol[j-1],x);      if (localdebug) {fprintf(stderr,"dsol[%ld]= ",nbindex[j]);dd_WriteNumber(stderr,dsol[j-1]);}    }    if (localdebug) fprintf(stderr,"DualSimplexSolve: LP is dual inconsistent.\n");    break;  default:break;  }  dd_clear(x); dd_clear(sw);}void dd_RandomPermutation2(dd_rowindex OV,long t,unsigned int seed){  long k,j,ovj;  double u,xk,r,rand_max=(double) RAND_MAX;  int localdebug=dd_FALSE;  srand(seed);  for (j=t; j>1 ; j--) {    r=rand();    u=r/rand_max;    xk=(double)(j*u +1);    k=(long)xk;    if (localdebug) fprintf(stderr,"u=%g, k=%ld, r=%g, randmax= %g\n",u,k,r,rand_max);    ovj=OV[j];    OV[j]=OV[k];    OV[k]=ovj;    if (localdebug) fprintf(stderr,"row %ld is exchanged with %ld\n",j,k);   }}void dd_ComputeRowOrderVector2(dd_rowrange m_size,dd_colrange d_size,dd_Amatrix A,    dd_rowindex OV,dd_RowOrderType ho,unsigned int rseed){  long i,itemp;    OV[0]=0;   switch (ho){  case dd_MaxIndex:    for(i=1; i<=m_size; i++) OV[i]=m_size-i+1;    break;  case dd_LexMin:    for(i=1; i<=m_size; i++) OV[i]=i;    dd_QuickSort(OV,1,m_size,A,d_size);   break;  case dd_LexMax:    for(i=1; i<=m_size; i++) OV[i]=i;    dd_QuickSort(OV,1,m_size,A,d_size);    for(i=1; i<=m_size/2;i++){   /* just reverse the order */      itemp=OV[i];      OV[i]=OV[m_size-i+1];      OV[m_size-i+1]=itemp;    }    break;  case dd_RandomRow:    for(i=1; i<=m_size; i++) OV[i]=i;    if (rseed<=0) rseed=1;    dd_RandomPermutation2(OV,m_size,rseed);    break;  case dd_MinIndex:     for(i=1; i<=m_size; i++) OV[i]=i;    break;  default:     for(i=1; i<=m_size; i++) OV[i]=i;    break; }}void dd_SelectPreorderedNext2(dd_rowrange m_size,dd_colrange d_size,    rowset excluded,dd_rowindex OV,dd_rowrange *hnext){  dd_rowrange i,k;    *hnext=0;  for (i=1; i<=m_size && *hnext==0; i++){    k=OV[i];    if (!set_member(k,excluded)) *hnext=k ;  }}#ifdef GMPRATIONALddf_LPObjectiveType Obj2Obj(dd_LPObjectiveType obj){   ddf_LPObjectiveType objf=ddf_LPnone;   switch (obj) {   case dd_LPnone: objf=ddf_LPnone; break;   case dd_LPmax: objf=ddf_LPmax; break;   case dd_LPmin: objf=ddf_LPmin; break;   }   return objf;}ddf_LPPtr dd_LPgmp2LPf(dd_LPPtr lp){  dd_rowrange i;  dd_colrange j;  ddf_LPType *lpf;  double val;  dd_boolean localdebug=dd_FALSE;  if (localdebug) fprintf(stderr,"Converting a GMP-LP to a float-LP.\n");    lpf=ddf_CreateLPData(Obj2Obj(lp->objective), ddf_Real, lp->m, lp->d);  lpf->Homogeneous = lp->Homogeneous;  lpf->eqnumber=lp->eqnumber;  /* this records the number of equations */  for (i = 1; i <= lp->m; i++) {    if (set_member(i, lp->equalityset)) set_addelem(lpf->equalityset,i);              /* it is equality. Its reversed row will not be in this set */      for (j = 1; j <= lp->d; j++) {        val=mpq_get_d(lp->A[i-1][j-1]);        ddf_set_d(lpf->A[i-1][j-1],val);      }  /*of j*/  }  /*of i*/  return lpf;}#endifdd_boolean dd_LPSolve(dd_LPPtr lp,dd_LPSolverType solver,dd_ErrorType *err)/* The current version of dd_LPSolve that solves an LP with floating-arithmetics firstand then with the specified arithimetics if it is GMP.When LP is inconsistent then *re returns the evidence row.When LP is dual-inconsistent then *se returns the evidence column.*/{  int i;  dd_boolean found=dd_FALSE;#ifdef GMPRATIONAL  ddf_LPPtr lpf;  ddf_ErrorType errf;  dd_boolean LPScorrect=dd_FALSE;  dd_boolean localdebug=dd_FALSE;#endif  *err=dd_NoError;  lp->solver=solver;     time(&lp->starttime);#ifndef GMPRATIONAL  switch (lp->solver) {    case dd_CrissCross:      dd_CrissCrossSolve(lp,err);      break;    case dd_DualSimplex:      dd_DualSimplexSolve(lp,err);      break;  }#else  lpf=dd_LPgmp2LPf(lp);  switch (lp->solver) {    case dd_CrissCross:      ddf_CrissCrossSolve(lpf,&errf);    /* First, run with double float. */	  if (errf==ddf_NoError){   /* 094a:  fix for a bug reported by Dima Pasechnik */        dd_BasisStatus(lpf,lp, &LPScorrect);    /* Check the basis. */	  } else {LPScorrect=dd_FALSE;}      if (!LPScorrect) {         if (localdebug) printf("BasisStatus: the current basis is NOT verified with GMP. Rerun with GMP.\n");         dd_CrissCrossSolve(lp,err);  /* Rerun with GMP if fails. */      } else {         if (localdebug) printf("BasisStatus: the current basis is verified with GMP. The LP Solved.\n");      }      break;    case dd_DualSimplex:      ddf_DualSimplexSolve(lpf,&errf);    /* First, run with double float. */	  if (errf==ddf_NoError){   /* 094a:  fix for a bug reported by Dima Pasechnik */        dd_BasisStatus(lpf,lp, &LPScorrect);    /* Check the basis. */	  } else {LPScorrect=dd_FALSE;}      if (!LPScorrect){         if (localdebug) printf("BasisStatus: the current basis is NOT verified with GMP. Rerun with GMP.\n");         dd_DualSimplexSolve(lp,err);  /* Rerun with GMP if fails. */         if (localdebug){            printf("*total number pivots = %ld (ph0 = %ld, ph1 = %ld, ph2 = %ld, ph3 = %ld, ph4 = %ld)\n",               lp->total_pivots,lp->pivots[0],lp->pivots[1],lp->pivots[2],lp->pivots[3],lp->pivots[4]);            ddf_WriteLPResult(stdout, lpf, errf);            dd_WriteLP(stdout, lp);         }      } else {         if (localdebug) printf("BasisStatus: the current basis is verified with GMP. The LP Solved.\n");      }      break;  }  ddf_FreeLPData(lpf);#endif  time(&lp->endtime);  lp->total_pivots=0;  for (i=0; i<=4; i++) lp->total_pivots+=lp->pivots[i];  if (*err==dd_NoError) found=dd_TRUE;  return found;}dd_boolean dd_LPSolve0(dd_LPPtr lp,dd_LPSolverType solver,dd_ErrorType *err)/* The original version of dd_LPSolve that solves an LP with specified arithimetics.When LP is inconsistent then *re returns the evidence row.When LP is dual-inconsistent then *se returns the evidence column.*/{  int i;  dd_boolean found=dd_FALSE;  *err=dd_NoError;  lp->solver=solver;  time(&lp->starttime);  switch (lp->solver) {    case dd_CrissCross:      dd_CrissCrossSolve(lp,err);      break;    case dd_DualSimplex:      dd_DualSimplexSolve(lp,err);      break;  }  time(&lp->endtime);  lp->total_pivots=0;  for (i=0; i<=4; i++) lp->total_pivots+=lp->pivots[i];  if (*err==dd_NoError) found=dd_TRUE;  return found;}dd_LPPtr dd_MakeLPforInteriorFinding(dd_LPPtr lp)/* Delete the objective row,   add an extra column with -1's to the matrix A,   add an extra row with (bceil, 0,...,0,-1),   add an objective row with (0,...,0,1), and    rows & columns, and change m_size and d_size accordingly, to output new_A.  This sets up the LP:  maximize      x_{d+1}  s.t.    A x + x_{d+1}  <=  b                x_{d+1}  <=  bm * bmax,  where bm is set to 2 by default, and bmax=max{1, b[1],...,b[m_size]}.  Note that the equalitions (linearity) in the input lp will be ignored.*/{  dd_rowrange m;  dd_colrange d;  dd_NumberType numbtype;  dd_LPObjectiveType obj;  dd_LPType *lpnew;  dd_rowrange i;   dd_colrange j;  mytype bm,bmax,bceil;  int localdebug=dd_FALSE;  dd_init(bm); dd_init(bmax); dd_init(bceil);  dd_add(bm,dd_one,dd_one); dd_set(bmax,dd_one);  numbtype=lp->numbtype;  m=lp->m+1;  d=lp->d+1;  obj=dd_LPmax;  lpnew=dd_CreateLPData(obj, numbtype, m, d);  for (i=1; i<=lp->m; i++) {    if (dd_Larger(lp->A[i-1][lp->rhscol-1],bmax))       dd_set(bmax,lp->A[i-1][lp->rhscol-1]);  }  dd_mul(bceil,bm,bmax);  if (