/SqlSpatialIndexing/Lib/Cdd/IO.cs

using System;

using System.Collections.Generic;

using System.Text;

using System.IO;



namespace CddSharp

{

    public partial class Cdd

    {

        /* cddio.c:  Basic Input and Output Procedures for cddlib

           written by Komei Fukuda, fukuda@ifor.math.ethz.ch

           Version 0.94, Aug. 4, 2005

        */



        /* cddlib : C-library of the double description method for

           computing all vertices and extreme rays of the polyhedron 

           P= {x :  b - A x >= 0}.  

           Please read COPYING (GNU General Public Licence) and

           the manual cddlibman.tex for detail.

        */



#if nemkell

void SetInputFile(FILE **f,DataFileType inputfile,ErrorType *Error)

{

  int opened=0,stop,quit=0;

  int i,dotpos=0,trial=0;

  char ch;

  char *tempname;

  

  

  *Error=NoError;

  while (!opened && !quit) {

    fprintf(stderr,"\n>> Input file: ");

    scanf("%s",inputfile);

    ch=getchar();

    stop=false;

    for (i=0; i<filenamelen && !stop; i++){

      ch=inputfile[i];

      switch (ch) {

        case '.': 

          dotpos=i+1;

          break;

        case ';':  case ' ':  case '\0':  case '\n':  case '\t':     

          stop=true;

          tempname=(char*)calloc(filenamelen,sizeof(ch));

          strncpy(tempname,inputfile,i);

          strcpy(inputfile,tempname);

          free(tempname);

          break;

      }

    }

    if ( ( *f = fopen(inputfile,"r") )!= null) {

      fprintf(stderr,"input file %s is open\n",inputfile);

      opened=1;

      *Error=NoError;

    }

    else{

      fprintf(stderr,"The file %s not found\n",inputfile);

      trial++;

      if (trial>5) {

        *Error=IFileNotFound;

        quit=1;

      }

    }

  }

}



void SetWriteFileName(DataFileType inputfile, DataFileType outfile, char cflag, RepresentationType rep)

{

  char *extension;

  DataFileType ifilehead="";

  int i,dotpos;

  

  switch (cflag) {

    case 'o':

      switch (rep) {

        case Generator:

          extension=".ine"; break;     /* output file for ine data */

        case Inequality:

          extension=".ext"; break;     /* output file for ext data */

        default: 

          extension=".xxx";break;

      }

      break;



    case 'a':         /* decide for output adjacence */

      if (rep==Inequality)

        extension=".ead";       /* adjacency file for ext data */

      else

        extension=".iad";       /* adjacency file for ine data */

      break;

    case 'i':         /* decide for output incidence */

      if (rep==Inequality)

        extension=".ecd";       /* ext incidence file */

      else

        extension=".icd";       /* ine incidence file */

      break;

    case 'n':         /* decide for input incidence */

      if (rep==Inequality)

        extension=".icd";       /* ine incidence file */

      else

        extension=".ecd";       /* ext incidence file */

      break;

    case 'j':        /* decide for input adjacency */

      if (rep==Inequality)

        extension=".iad";       /* ine adjacency file */

      else

        extension=".ead";       /* ext adjacency file */

      break;

    case 'l':

      extension=".ddl";break;   /* log file */

    case 'd':

      extension=".dex";break;   /* decomposition output */

    case 'p':

      extension="sub.ine";break;  /* preprojection sub inequality file */

    case 'v':

      extension=".solved";break;  /* verify_input file */

    case 's':

      extension=".lps";break;   /* LP solution file */

    default:

      extension=".xxx";break;

  }

  dotpos=-1;

  for (i=0; i< strlen(inputfile); i++){

    if (inputfile[i]=='.') dotpos=i;

  }

  if (dotpos>1) strncpy(ifilehead, inputfile, dotpos);

  else strcpy(ifilehead,inputfile);

  if (strlen(inputfile)<=0) strcpy(ifilehead,"tempcdd");

  strcpy(outfile,ifilehead); 

  strcat(outfile,extension); 

  if (strcmp(inputfile, outfile)==0) {

    strcpy(outfile,inputfile); 

    strcat(outfile,extension); 

  }

/*  fprintf(stderr,"outfile name = %s\n",outfile);  */

}





NumberType GetNumberType(char *line)

{

  NumberType nt;



  if (strncmp(line, "integer", 7)==0) {

    nt = Integer;

  }

  else if (strncmp(line, "rational", 8)==0) {

    nt = Rational;

  }

  else if (strncmp(line, "real", 4)==0) {

    nt = Real;

  }

  else { 

    nt=Unknown;

  }

  return nt;

}



void ProcessCommandLine(FILE *f, Matrix M, char *line)

{

  char newline[linelenmax];

  int j;

  double value;



  init(value);

  if (strncmp(line, "hull", 4)==0) {

    M.representation = Generator;

  }

  if (strncmp(line, "debug", 5)==0) {

    debug = true;

//#ifdef GMPRATIONAL

    ddf_debug = ddf_TRUE;

//#endif

  }

  if (strncmp(line, "partial_enum", 12)==0 ||

       strncmp(line, "equality", 8)==0  ||

       strncmp(line, "linearity", 9)==0 ) {

    fgets(newline,linelenmax,f);    

    SetLinearity(M,newline);

  }

  if (strncmp(line, "maximize", 8)==0 ||

      strncmp(line, "minimize", 8)==0) {

    if (strncmp(line, "maximize", 8)==0) M.objective=LPmax;

    else M.objective=LPmin;

    for (j = 1; j <= M.colsize; j++) {

    if (M.numbtype==Real) {

//#if !defined(GMPRATIONAL)

        double rvalue;

        fscanf(f, "%lf", &rvalue);

        set_d(value, rvalue);

//#endif

      } else {

        fread_rational_value (f, value);

      }

      set(M.rowvec[j - 1],value);

      if (debug) {fprintf(stderr,"cost(%5ld) =",j); WriteNumber(stderr,value);}

    }  /*of j*/

  }

  clear(value);

}

#endif



        public bool AppendMatrix2Poly(Polyhedra poly, Matrix M)

        {

            bool success = false;

            Matrix Mpoly, Mnew = null;

            ErrorType err;



            if (poly != null && poly.m >= 0 && poly.d >= 0 &&

                poly.d == M.colsize && M.rowsize > 0)

            {

                Mpoly = CopyInput(poly);

                Mnew = AppendMatrix(Mpoly, M);



                poly = DDMatrix2Poly(Mnew, out err);



                if (err == ErrorType.NoError) success = true;

            }

            return success;

        }



        public Matrix MatrixCopy(Matrix M)

        {

            Matrix Mcopy = null;

            int m;

            int d;



            m = M.rowsize;

            d = M.colsize;

            if (m >= 0 && d >= 0)

            {

                Mcopy = CreateMatrix(m, d);

                CopyAmatrix(Mcopy.matrix, M.matrix, m, d);

                CopyArow(Mcopy.rowvec, M.rowvec, d);

                M.linset.CopyTo(Mcopy.linset);

                Mcopy.numbtype = M.numbtype;

                Mcopy.representation = M.representation;

                Mcopy.objective = M.objective;

            }

            return Mcopy;

        }



        public Matrix CopyMatrix(Matrix M)

        {

            return MatrixCopy(M);

        }



        public Matrix MatrixNormalizedCopy(Matrix M)

        {

            Matrix Mcopy = null;

            int m;

            int d;



            m = M.rowsize;

            d = M.colsize;

            if (m >= 0 && d >= 0)

            {

                Mcopy = CreateMatrix(m, d);

                CopyNormalizedAmatrix(Mcopy.matrix, M.matrix, m, d);

                CopyArow(Mcopy.rowvec, M.rowvec, d);

                M.linset.CopyTo(Mcopy.linset);

                Mcopy.numbtype = M.numbtype;

                Mcopy.representation = M.representation;

                Mcopy.objective = M.objective;

            }

            return Mcopy;

        }





        public Matrix MatrixAppend(Matrix M1, Matrix M2)

        {

            Matrix M = null;

            int i, m, m1, m2;

            int j, d, d1, d2;



            m1 = M1.rowsize;

            d1 = M1.colsize;

            m2 = M2.rowsize;

            d2 = M2.colsize;



            m = m1 + m2;

            d = d1;



            if (d1 >= 0 && d1 == d2 && m1 >= 0 && m2 >= 0)

            {

                M = CreateMatrix(m, d);

                CopyAmatrix(M.matrix, M1.matrix, m1, d);

                CopyArow(M.rowvec, M1.rowvec, d);

                for (i = 0; i < m1; i++)

                {

                    if (M1.linset.Get(i + 1)) M.linset.Set(i + 1, true);

                }

                for (i = 0; i < m2; i++)

                {

                    for (j = 0; j < d; j++)

                        M.matrix[m1 + i][j] = M2.matrix[i][j];

                    /* append the second matrix */

                    if (M2.linset.Get(i + 1)) M.linset.Set(m1 + i + 1, true);

                }

                M.numbtype = M1.numbtype;

            }

            return M;

        }



        public Matrix MatrixNormalizedSortedCopy(Matrix M, out int[] newpos)  /* 094 */

        {

            /* Sort the rows of double[][] lexicographically, and return a link to this sorted copy. 

            The vector newpos is allocated, where newpos[i] returns the new row index

            of the original row i (i=1,...,M.rowsize). */

            Matrix Mcopy = null, Mnorm = null;

            int m, i;

            int d;

            int[] roworder;



            /* if (newpos!=null) free(newpos); */

            m = M.rowsize;

            d = M.colsize;

            roworder = new int[m + 1];

            newpos = new int[m + 1];

            if (m >= 0 && d >= 0)

            {

                Mnorm = MatrixNormalizedCopy(M);

                Mcopy = CreateMatrix(m, d);

                for (i = 1; i <= m; i++) roworder[i] = i;



                RandomPermutation(roworder, m, 123);

                QuickSort(roworder, 1, m, Mnorm.matrix, d);



                PermuteCopyAmatrix(Mcopy.matrix, Mnorm.matrix, m, d, roworder);

                CopyArow(Mcopy.rowvec, M.rowvec, d);

                for (i = 1; i <= m; i++)

                {

                    if (M.linset.Get(roworder[i])) Mcopy.linset.Set(i, true);

                    newpos[roworder[i]] = i;

                }

                Mcopy.numbtype = M.numbtype;

                Mcopy.representation = M.representation;

                Mcopy.objective = M.objective;



            }



            return Mcopy;

        }



        public Matrix MatrixUniqueCopy(Matrix M, out int[] newpos)

        {

            /* Remove row duplicates, and return a link to this sorted copy. 

               Linearity rows have priority over the other rows.

               It is better to call this after sorting with MatrixNormalizedSortedCopy.

               The vector newpos is allocated, where *newpos[i] returns the new row index

               of the original row i (i=1,...,M.rowsize).  *newpos[i] is negative if the original

               row is dominated by -*newpos[i] and eliminated in the new copy.

            */

            Matrix Mcopy = null;

            int m, i, uniqrows;

            BitArray preferredrows;

            int d;

            int[] roworder;



            /* if (newpos!=null) free(newpos); */

            m = M.rowsize;

            d = M.colsize;

            preferredrows = M.linset;

            roworder = new int[m + 1];

            if (m >= 0 && d >= 0)

            {

                for (i = 1; i <= m; i++) roworder[i] = i;

                UniqueRows(roworder, 1, m, M.matrix, d, preferredrows, out uniqrows);



                Mcopy = CreateMatrix(uniqrows, d);

                PermutePartialCopyAmatrix(Mcopy.matrix, M.matrix, m, d, roworder, 1, m);

                CopyArow(Mcopy.rowvec, M.rowvec, d);

                for (i = 1; i <= m; i++)

                {

                    if (roworder[i] > 0 && M.linset.Get(i)) Mcopy.linset.Set(roworder[i], true);

                }

                Mcopy.numbtype = M.numbtype;

                Mcopy.representation = M.representation;

                Mcopy.objective = M.objective;

            }

            newpos = roworder;

            return Mcopy;

        }





        public Matrix MatrixNormalizedSortedUniqueCopy(Matrix M, out int[] newpos) /* 094 */

        {

            /* Sort and remove row duplicates, and return a link to this sorted copy. 

               Linearity rows have priority over the other rows.

               It is better to call this after sorting with MatrixNormalizedSortedCopy.

               The vector newpos is allocated, where *newpos[i] returns the new row index

               of the original row i (i=1,...,M.rowsize).  *newpos[i] is negative if the original

               row is dominated by -*newpos[i] and eliminated in the new copy.

            */

            Matrix M1 = null, M2 = null;

            int m, i;

            int d;

            int[] newpos1 = null, newpos1r = null, newpos2 = null;



            /* if (newpos!=null) free(newpos); */

            m = M.rowsize;

            d = M.colsize;

            newpos = new int[m + 1];

            newpos1r = new int[m + 1];

            if (m >= 0 && d >= 0)

            {

                M1 = MatrixNormalizedSortedCopy(M, out newpos1);

                for (i = 1; i <= m; i++) newpos1r[newpos1[i]] = i;  /* reverse of newpos1 */

                M2 = MatrixUniqueCopy(M1, out newpos2);

                M2.linset.SetAll(false);    //emptyset

                for (i = 1; i <= m; i++)

                {

                    if (newpos2[newpos1[i]] > 0)

                    {

                        //printf("newpos1[%ld]=%ld, newpos2[newpos1[%ld]]=%ld\n",i,newpos1[i], i,newpos2[newpos1[i]]);

                        if (M.linset.Get(i)) M2.linset.Set(newpos2[newpos1[i]], true);

                        newpos[i] = newpos2[newpos1[i]];

                    }

                    else

                    {

                        newpos[i] = -newpos1r[-newpos2[newpos1[i]]];

                    }

                }



            }



            return M2;

        }



        public Matrix MatrixSortedUniqueCopy(Matrix M, out int[] newpos)  /* 094 */

        {

            /* Same as MatrixNormalizedSortedUniqueCopy except that it returns a unnormalized origial data

               with original ordering.

            */

            Matrix M1 = null, M2 = null;

            int m, i, k, ii;

            int d;

            int[] newpos1 = null, newpos1r = null, newpos2 = null;



            /* if (newpos!=null) free(newpos); */

            m = M.rowsize;

            d = M.colsize;

            newpos = new int[m + 1];

            newpos1r = new int[m + 1];

            if (m >= 0 && d >= 0)

            {

                M1 = MatrixNormalizedSortedCopy(M, out newpos1);

                for (i = 1; i <= m; i++) newpos1r[newpos1[i]] = i;  /* reverse of newpos1 */

                M2 = MatrixUniqueCopy(M1, out newpos2);



                M2.linset.SetAll(false);

                for (i = 1; i <= m; i++)

                {

                    if (newpos2[newpos1[i]] > 0)

                    {

                        if (M.linset.Get(i)) M2.linset.Set(newpos2[newpos1[i]], true);

                        newpos[i] = newpos2[newpos1[i]];

                    }

                    else

                    {

                        newpos[i] = -newpos1r[-newpos2[newpos1[i]]];

                    }

                }



                ii = 0;

                M2.linset.SetAll(false);

                for (i = 1; i <= m; i++)

                {

                    k = newpos[i];

                    if (k > 0)

                    {

                        ii += 1;

                        newpos[i] = ii;

                        CopyArow(M2.matrix[ii - 1], M.matrix[i - 1], d);

                        if (M.linset.Get(i)) M2.linset.Set(ii, true);

                    }

                }



            }



            return M2;

        }



        public Matrix AppendMatrix(Matrix M1, Matrix M2)

        {

            return MatrixAppend(M1, M2);

        }



        public bool MatrixAppendTo(ref Matrix M1, Matrix M2)

        {

            Matrix M = null;

            int i, m, m1, m2;

            int j, d, d1, d2;

            bool success = false;



            m1 = M1.rowsize;

            d1 = M1.colsize;

            m2 = M2.rowsize;

            d2 = M2.colsize;



            m = m1 + m2;

            d = d1;



            if (d1 >= 0 && d1 == d2 && m1 >= 0 && m2 >= 0)

            {

                M = CreateMatrix(m, d);

                CopyAmatrix(M.matrix, M1.matrix, m1, d);

                CopyArow(M.rowvec, M1.rowvec, d);

                for (i = 0; i < m1; i++)

                {

                    if (M1.linset.Get(i + 1)) M.linset.Set(i + 1, true);

                }

                for (i = 0; i < m2; i++)

                {

                    for (j = 0; j < d; j++)

                        M.matrix[m1 + i][j] = M2.matrix[i][j];

                    /* append the second matrix */

                    if (M2.linset.Get(i + 1)) M.linset.Set(m1 + i + 1, true);

                }

                M.numbtype = M1.numbtype;



                M1 = M;

                success = true;

            }

            return success;

        }



        public bool MatrixRowRemove(Matrix M, int r) /* 092 */

        {

            int i, m;

            int d;

            bool success = false;



            m = M.rowsize;

            d = M.colsize;



            if (r >= 1 && r <= m)

            {

                M.rowsize = m - 1;



                M.linset.Set(r, false);

                /* slide the row headers */

                for (i = r; i < m; i++)

                {

                    M.matrix[i - 1] = M.matrix[i];

                    if (M.linset.Get(i + 1))

                    {

                        M.linset.Set(i + 1, false);

                        M.linset.Set(i, true);

                    }

                }

                success = true;

            }

            return success;

        }



        public bool MatrixRowRemove2(Matrix M, int r, out int[] newpos) /* 094 */

        {

            int i, m;

            int d;

            bool success = false;

            int[] roworder;



            newpos = null;



            m = M.rowsize;

            d = M.colsize;



            if (r >= 1 && r <= m)

            {

                roworder = new int[m + 1];

                M.rowsize = m - 1;



                M.linset.Set(r, false);

                /* slide the row headers */

                for (i = 1; i < r; i++) roworder[i] = i;

                roworder[r] = 0; /* meaning it is removed */

                for (i = r; i < m; i++)

                {

                    M.matrix[i - 1] = M.matrix[i];

                    roworder[i + 1] = i;

                    if (M.linset.Get(i + 1))

                    {

                        M.linset.Set(i + 1, false);

                        M.linset.Set(i, true);

                    }

                }

                success = true;

            }

            return success;

        }



        public Matrix MatrixSubmatrix(Matrix M, BitArray delset) /* 092 */

        {

            Matrix Msub = null;

            int i, isub = 1, m, msub;

            int d;



            m = M.rowsize;

            d = M.colsize;

            msub = m;

            if (m >= 0 && d >= 0)

            {

                for (i = 1; i <= m; i++)

                {

                    if (delset.Get(i)) msub -= 1;

                }

                Msub = CreateMatrix(msub, d);

                for (i = 1; i <= m; i++)

                {

                    if (!delset.Get(i))

                    {

                        CopyArow(Msub.matrix[isub - 1], M.matrix[i - 1], d);

                        if (M.linset.Get(i))

                        {

                            Msub.linset.Set(isub, true);

                        }

                        isub++;

                    }

                }

                CopyArow(Msub.rowvec, M.rowvec, d);

                Msub.numbtype = M.numbtype;

                Msub.representation = M.representation;

                Msub.objective = M.objective;

            }

            return Msub;

        }



        public Matrix MatrixSubmatrix2(Matrix M, BitArray delset, out int[] newpos) /* 092 */

        { /* returns a pointer to a new matrix which is a submatrix of M with rows in delset

  removed.  *newpos[i] returns the position of the original row i in the new matrix.

  It is -1 if and only if it is deleted.

  */



            Matrix Msub = null;

            int i, isub = 1, m, msub;

            int d;

            int[] roworder;



            newpos = null;



            m = M.rowsize;

            d = M.colsize;

            msub = m;

            if (m >= 0 && d >= 0)

            {

                roworder = new int[m + 1];

                for (i = 1; i <= m; i++)

                {

                    if (delset.Get(i)) msub -= 1;

                }

                Msub = CreateMatrix(msub, d);

                for (i = 1; i <= m; i++)

                {

                    if (delset.Get(i))

                    {

                        roworder[i] = 0; /* zero means the row i is removed */

                    }

                    else

                    {

                        CopyArow(Msub.matrix[isub - 1], M.matrix[i - 1], d);

                        if (M.linset.Get(i))

                        {

                            Msub.linset.Set(isub, true);

                        }

                        roworder[i] = isub;

                        isub++;

                    }

                }

                newpos = roworder;

                CopyArow(Msub.rowvec, M.rowvec, d);

                Msub.numbtype = M.numbtype;

                Msub.representation = M.representation;

                Msub.objective = M.objective;

            }

            return Msub;

        }





        public Matrix MatrixSubmatrix2L(Matrix M, BitArray delset, out int[] newpos) /* 094 */

        { /* This is same as MatrixSubmatrix2 except that the linearity rows will be shifted up

     so that they are at the top of the matrix.

  */

            Matrix Msub = null;

            int i, iL, iI, m, msub;

            int d;

            int[] roworder;



            newpos = null;



            m = M.rowsize;

            d = M.colsize;

            msub = m;

            if (m >= 0 && d >= 0)

            {

                roworder = new int[m + 1];

                for (i = 1; i <= m; i++)

                {

                    if (delset.Get(i)) msub -= 1;

                }

                Msub = CreateMatrix(msub, d);

                iL = 1; iI = M.linset.Count + 1;  /* starting positions */

                for (i = 1; i <= m; i++)

                {

                    if (delset.Get(i))

                    {

                        roworder[i] = 0; /* zero means the row i is removed */

                    }

                    else

                    {

                        if (M.linset.Get(i))

                        {

                            CopyArow(Msub.matrix[iL - 1], M.matrix[i - 1], d);

                            Msub.linset.Set(i, false);

                            Msub.linset.Set(iL, true);

                            roworder[i] = iL;

                            iL += 1;

                        }

                        else

                        {

                            CopyArow(Msub.matrix[iI - 1], M.matrix[i - 1], d);

                            roworder[i] = iI;

                            iI += 1;

                        }

                    }

                }

                newpos = roworder;

                CopyArow(Msub.rowvec, M.rowvec, d);

                Msub.numbtype = M.numbtype;

                Msub.representation = M.representation;

                Msub.objective = M.objective;

            }

            return Msub;

        }



        public int MatrixRowsRemove(ref Matrix M, BitArray delset) /* 094 */

        {

            Matrix Msub = null;

            int success;



            Msub = MatrixSubmatrix(M, delset);



            M = Msub;

            success = 1;

            return success;

        }



        public int MatrixRowsRemove2(ref Matrix M, BitArray delset, out int[] newpos) /* 094 */

        {

            Matrix Msub = null;

            int success;



            Msub = MatrixSubmatrix2(M, delset, out newpos);



            M = Msub;

            success = 1;

            return success;

        }



        public int MatrixShiftupLinearity(ref Matrix M, out int[] newpos) /* 094 */

        {

            Matrix Msub = null;

            int success;

            BitArray delset;



            delset = new BitArray(M.rowsize);  /* emptyset */

            Msub = MatrixSubmatrix2L(M, delset, out newpos);



            M = Msub;



            success = 1;

            return success;

        }



        public Polyhedra CreatePolyhedraData(int m, int d)

        {

            int i;

            Polyhedra poly = null;



            poly = new Polyhedra();

            poly.child = null; /* this links the homogenized cone data */

            poly.m = m;

            poly.d = d;

            poly.n = -1;  /* the size of output is not known */

            poly.m_alloc = m + 2; /* the allocated row size of matrix A */

            poly.d_alloc = d;   /* the allocated col size of matrix A */

            poly.numbtype = NumberType.Real;

            InitializeAmatrix(poly.m_alloc, poly.d_alloc, out poly.A);

            InitializeArow(d, out poly.c);           /* cost vector */

            poly.representation = RepresentationType.Inequality;

            poly.homogeneous = false;



            poly.EqualityIndex = new int[m + 2];

            /* size increased to m+2 in 092b because it is used by the child cone, 

               This is a bug fix suggested by Thao Dang. */

            /* ith component is 1 if it is equality, -1 if it is strict inequality, 0 otherwise. */

            for (i = 0; i <= m + 1; i++) poly.EqualityIndex[i] = 0;



            poly.IsEmpty = -1;  /* initially set to -1, neither true nor false, meaning unknown */



            poly.NondegAssumed = false;

            poly.InitBasisAtBottom = false;

            poly.RestrictedEnumeration = false;

            poly.RelaxedEnumeration = false;



            poly.AincGenerated = false;  /* Ainc is a set array to store the input incidence. */



            return poly;

        }



        public bool InitializeConeData(int m, int d, out Cone cone)

        {

            bool success = true;

            int j;



            cone = new Cone();



            /* INPUT: A given representation of a cone: inequality */

            cone.m = m;

            cone.d = d;

            cone.m_alloc = m + 2; /* allocated row size of matrix A */

            cone.d_alloc = d;   /* allocated col size of matrix A, B and Bsave */

            cone.numbtype = NumberType.Real;

            cone.parent = null;



            /* CONTROL: variables to control computation */

            cone.Iteration = 0;



            cone.HalfspaceOrder = RowOrderType.LexMin;



            cone.ArtificialRay = null;

            cone.FirstRay = null;

            cone.LastRay = null; /* The second description: Generator */

            cone.PosHead = null;

            cone.ZeroHead = null;

            cone.NegHead = null;

            cone.PosLast = null;

            cone.ZeroLast = null;

            cone.NegLast = null;

            cone.RecomputeRowOrder = true;

            cone.PreOrderedRun = false;

            cone.GroundSet = new BitArray(cone.m_alloc);

            cone.EqualitySet = new BitArray(cone.m_alloc);

            cone.NonequalitySet = new BitArray(cone.m_alloc);

            cone.AddedHalfspaces = new BitArray(cone.m_alloc);

            cone.WeaklyAddedHalfspaces = new BitArray(cone.m_alloc);

            cone.InitialHalfspaces = new BitArray(cone.m_alloc);

            cone.RayCount = 0;

            cone.FeasibleRayCount = 0;

            cone.WeaklyFeasibleRayCount = 0;

            cone.TotalRayCount = 0;

            cone.ZeroRayCount = 0;

            cone.EdgeCount = 0;

            cone.TotalEdgeCount = 0;

            cone.count_int = 0;

            cone.count_int_good = 0;

            cone.count_int_bad = 0;

            cone.rseed = 1;  /* random seed for random row permutation */



            InitializeBmatrix(cone.d_alloc, out cone.B);

            InitializeBmatrix(cone.d_alloc, out cone.Bsave);

            InitializeAmatrix(cone.m_alloc, cone.d_alloc, out cone.A);



            cone.Edges = new Adjacency[cone.m_alloc];

            cone.InitialRayIndex = new int[d + 1];

            cone.OrderVector = new int[cone.m_alloc + 1];



            cone.newcol = new int[cone.d + 1];

            for (j = 0; j <= cone.d; j++) cone.newcol[j] = j;  /* identity map, initially */

            cone.LinearityDim = -2; /* -2 if it is not computed */

            cone.ColReduced = false;

            cone.d_orig = d;



            /* STATES: variables to represent current state. */

            /*cone.Error;

              cone.CompStatus;

              cone.starttime;

              cone.endtime;

            */



            return success;

        }



        public Cone ConeDataLoad(Polyhedra poly)

        {

            Cone cone = null;

            int d, j;

            int m, i;



            m = poly.m;

            d = poly.d;

            if (!(poly.homogeneous) && poly.representation == RepresentationType.Inequality)

            {

                m = poly.m + 1;

            }

            poly.m1 = m;



            InitializeConeData(m, d, out cone);

            cone.representation = poly.representation;



            /* Points to the original polyhedra data, and reversely */

            cone.parent = poly;

            poly.child = cone;



            for (i = 1; i <= poly.m; i++)

                for (j = 1; j <= cone.d; j++)

                    cone.A[i - 1][j - 1] = poly.A[i - 1][j - 1];



            if (poly.representation == RepresentationType.Inequality && !(poly.homogeneous))

            {

                cone.A[m - 1][0] = 1.0;

                for (j = 2; j <= d; j++) cone.A[m - 1][j - 1] = 0.0;

            }



            return cone;

        }



#if nemkell

void SetLinearity(Matrix M, char *line)

{

  int i=0;

  int eqsize,var;

  char *next;

  const char ct[]=", ";  /* allows separators "," and " ". */



  next=strtok(line,ct);

  eqsize=atol(next); 

  while (i < eqsize && (next=strtok(null,ct))!=null) {

     var=atol(next);

     set_addelem(M.linset,var); i++;

  }

  if (i!=eqsize) {

    fprintf(stderr,"* Warning: there are inconsistencies in linearity setting.\n");

  }

  return;

}



Matrix PolyFile2Matrix (FILE *f, ErrorType *Error)

{

  Matrix M=null;

  int m_input,i;

  int d_input,j;

  RepresentationType rep=Inequality;

  double value;

  bool found=false, newformat=false, successful=false, linearity=false;

  char command[linelenmax], comsave[linelenmax], numbtype[wordlenmax];

  NumberType NT;

//#if !defined(GMPRATIONAL)

  double rvalue;

//#endif



  init(value);

  (*Error)=NoError;

  while (!found)

  {

    if (fscanf(f,"%s",command)==EOF) {

      (*Error)=ImproperInputFormat;

      goto _L99;

    }

    else {

      if (strncmp(command, "V-representation", 16)==0) {

        rep=Generator; newformat=true;

      }

      if (strncmp(command, "H-representation", 16)==0){

        rep=Inequality; newformat=true;

      }

      if (strncmp(command, "partial_enum", 12)==0 || 

          strncmp(command, "equality", 8)==0  ||

          strncmp(command, "linearity", 9)==0 ) {

        linearity=true;

        fgets(comsave,linelenmax,f);

      }

      if (strncmp(command, "begin", 5)==0) found=true;

    }

  }

  fscanf(f, "%ld %ld %s", &m_input, &d_input, numbtype);

  fprintf(stderr,"size = %ld x %ld\nNumber Type = %s\n", m_input, d_input, numbtype);

  NT=GetNumberType(numbtype);

  if (NT==Unknown) {

      (*Error)=ImproperInputFormat;

      goto _L99;

    } 

  M=CreateMatrix(m_input, d_input);

  M.representation=rep;

  M.numbtype=NT;



  for (i = 1; i <= m_input; i++) {

    for (j = 1; j <= d_input; j++) {

      if (NT==Real) {

//#if defined GMPRATIONAL

        *Error=NoRealNumberSupport;

        goto _L99;

//#else

        fscanf(f, "%lf", &rvalue);

        set_d(value, rvalue);

//#endif

      } else {

        fread_rational_value (f, value);

      }

      set(M.matrix[i-1][j - 1],value);

      if (debug) {fprintf(stderr,"a(%3ld,%5ld) = ",i,j); WriteNumber(stderr,value);}

    }  /*of j*/

  }  /*of i*/

  if (fscanf(f,"%s",command)==EOF) {

   	 (*Error)=ImproperInputFormat;

  	 goto _L99;

  }

  else if (strncmp(command, "end", 3)!=0) {

     if (debug) fprintf(stderr,"'end' missing or illegal extra data: %s\n",command);

     (*Error)=ImproperInputFormat;

     goto _L99;

  }

  

  successful=true;

  if (linearity) {

    SetLinearity(M,comsave);

  }

  while (!feof(f)) {

    fscanf(f,"%s", command);

    ProcessCommandLine(f, M, command);

    fgets(command,linelenmax,f); /* skip the CR/LF */

  } 



_L99: ;

  clear(value);

  /* if (f!=null) fclose(f); */

  return M;

}



#endif



        public Polyhedra DDMatrix2Poly(Matrix M, out ErrorType err)

        {

            int i;

            int j;

            Polyhedra poly = null;



            err = ErrorType.NoError;

            if (M.rowsize < 0 || M.colsize < 0)

            {

                err = ErrorType.NegativeMatrixSize;

                return null;

            }

            poly = CreatePolyhedraData(M.rowsize, M.colsize);

            poly.representation = M.representation;

            poly.homogeneous = true;



            for (i = 1; i <= M.rowsize; i++)

            {

                if (M.linset.Get(i))

                {

                    poly.EqualityIndex[i] = 1;

                }

                for (j = 1; j <= M.colsize; j++)

                {

                    poly.A[i - 1][j - 1] = M.matrix[i - 1][j - 1];

                    if (j == 1 && M.matrix[i - 1][j - 1] != 0.0) poly.homogeneous = false;

                }  /*of j*/

            }  /*of i*/

            DoubleDescription(poly, out err);



            return poly;

        }



        public Polyhedra DDMatrix2Poly2(Matrix M, RowOrderType horder, out ErrorType err)

        {

            int i;

            int j;

            Polyhedra poly = null;



            err = ErrorType.NoError;

            if (M.rowsize < 0 || M.colsize < 0)

            {

                err = ErrorType.NegativeMatrixSize;

                return null;

            }

            poly = CreatePolyhedraData(M.rowsize, M.colsize);

            poly.representation = M.representation;

            poly.homogeneous = true;



            for (i = 1; i <= M.rowsize; i++)

            {

                if (M.linset.Get(i))

                {

                    poly.EqualityIndex[i] = 1;

                }

                for (j = 1; j <= M.colsize; j++)

                {

                    poly.A[i - 1][j - 1] = M.matrix[i - 1][j - 1];

                    if (j == 1 && M.matrix[i - 1][j - 1] != 0.0) poly.homogeneous = false;

                }  /*of j*/

            }  /*of i*/

            DoubleDescription2(poly, horder, out err);



            return poly;

        }



        void MatrixIntegerFilter(Matrix M)

        {   /* setting an almost integer to the integer. */

            int i;

            int j;

            double x;



            x = 0.0;

            for (i = 0; i < M.rowsize; i++)

                for (j = 0; j < M.colsize; j++)

                {

                    SnapToInteger(out x, M.matrix[i][j]);

                    M.matrix[i][j] = x;

                }



        }



        void CopyRay(double[] a, int d_origsize, Ray RR,

          RepresentationType rep, int[] reducedcol)

        {

            int j, j1;

            double b;



            b = 0.0;

            for (j = 1; j <= d_origsize; j++)

            {

                j1 = reducedcol[j];

                if (j1 > 0)

                {

                    a[j - 1] = RR.RRay[j1 - 1];

                    /* the original column j is mapped to j1, and thus

                       copy the corresponding component */

                }

                else

                {

                    a[j - 1] = 0.0;

                    /* original column is redundant and removed for computation */

                }

            }



            b = a[0];

            if (rep == RepresentationType.Generator && b != 0)

            {

                a[0] = 1.0;

                for (j = 2; j <= d_origsize; j++)

                    a[j - 1] /= b;    /* normalization for generators */

            }

        }





        void WriteRay(TextWriter f, int d_origsize, Ray RR, RepresentationType rep, int[] reducedcol)

        {

            int j;

            double[] a;





            InitializeArow(d_origsize + 1, out a);



            CopyRay(a, d_origsize, RR, rep, reducedcol);

            for (j = 0; j < d_origsize; j++) WriteNumber(f, a[j]);

            f.WriteLine();

        }



        void WriteArow(TextWriter f, double[] a, int d)

        {

            int j;



            for (j = 0; j < d; j++) WriteNumber(f, a[j]);

            f.WriteLine();

        }



        void WriteAmatrix(TextWriter f, double[][] A, int rowmax, int colmax)

        {

            int i, j;



            if (A == null)

            {

                f.WriteLine("WriteAmatrix: The requested matrix is empty");

                return;

            }

            f.WriteLine("begin");

#if GMPRATIONAL

  fprintf(f, " %ld %ld rational\n",rowmax, colmax);

#else

            f.WriteLine(" {0} {1} real", rowmax, colmax);

#endif

            for (i = 1; i <= rowmax; i++)

            {

                for (j = 1; j <= colmax; j++)

                {

                    WriteNumber(f, A[i - 1][j - 1]);

                }

                f.WriteLine();

            }

            f.WriteLine("end");



        }



        void WriteBmatrix(TextWriter f, int d_size, double[][] B)

        {

            int j1, j2;



            if (B == null)

            {

                f.WriteLine("WriteBmatrix: The requested matrix is empty");

                return;

            }

            for (j1 = 0; j1 < d_size; j1++)

            {

                for (j2 = 0; j2 < d_size; j2++)

                {

                    WriteNumber(f, B[j1][j2]);

                }  /*of j2*/

                f.WriteLine();

            }  /*of j1*/

            f.WriteLine();



        }

#if nemkell

void WriteSetFamily(TextWriter f, SetFamily F)

{

  int i;



  if (F==null){

    f.WriteLine("WriteSetFamily: The requested family is empty");

    return;

  }

  f.WriteLine("begin");

  f.WriteLine("  {0}    {1}", F.famsize, F.setsize);

  for (i=0; i<F.famsize; i++) {

    f.Write(" {0} {1} : ", i+1, F.set[i].Count);

    //*****set_fwrite(f, F.set[i]);

  }

  f.WriteLine("end");



}



void WriteSetFamilyCompressed(FILE *f, SetFamily F)

{

  int i,card;



  if (F==null){

    fprintf(f, "WriteSetFamily: The requested family is empty\n");

    goto _L99;

  }

  fprintf(f,"begin\n");

  fprintf(f,"  %ld    %ld\n", F.famsize, F.setsize);

  for (i=0; i<F.famsize; i++) {

    card=set_card(F.set[i]);

    if (F.setsize - card >= card){

      fprintf(f, " %ld %ld : ", i+1, card);

      set_fwrite(f, F.set[i]);

    } else {

      fprintf(f, " %ld %ld : ", i+1, -card);

      set_fwrite_compl(f, F.set[i]);

    }

  }

  fprintf(f,"end\n");

_L99:;

}

#endif

        void WriteMatrix(TextWriter f, Matrix M)

        {

            int i, linsize;



            if (M == null)

            {

                f.WriteLine("WriteMatrix: The requested matrix is empty");

                return;

            }

            switch (M.representation)

            {

                case RepresentationType.Inequality:

                    f.WriteLine("H-representation"); break;

                case RepresentationType.Generator:

                    f.WriteLine("V-representation"); break;

                case RepresentationType.Unspecified:

                    break;

            }

            linsize = M.linset.Count;

            if (linsize > 0)

            {

                f.Write("linearity {0} ", linsize);

                for (i = 1; i <= M.rowsize; i++)

                    if (M.linset.Get(i)) f.Write(" {0}", i);

                f.WriteLine();

            }

            WriteAmatrix(f, M.matrix, M.rowsize, M.colsize);

            if (M.objective != LPObjectiveType.LPnone)

            {

                if (M.objective == LPObjectiveType.LPmax)

                    f.WriteLine("maximize");

                else

                    f.WriteLine("minimize");

                WriteArow(f, M.rowvec, M.colsize);

            }



        }



        void WriteLP(TextWriter f, LP lp)

        {

            if (lp == null)

            {

                f.WriteLine("WriteLP: The requested lp is empty");

                return;

            }

            f.WriteLine("H-representation");



            WriteAmatrix(f, lp.A, (lp.m) - 1, lp.d);

            if (lp.objective != LPObjectiveType.LPnone)

            {

                if (lp.objective == LPObjectiveType.LPmax)

                    f.WriteLine("maximize");

                else

                    f.WriteLine("minimize");

                WriteArow(f, lp.A[lp.objrow - 1], lp.d);

            }



        }





        void SnapToInteger(out double y, double x)

        {

            /* this is broken.  It does nothing. */

            y = x;

        }





        void WriteReal(TextWriter f, double x)

        {

            /*

          int ix1,ix2,ix;

          double ax;



          ax=x;  

          ix1= (int) (Math.Abs(ax) * 10000.0 + 0.5);

          ix2= (int) (Math.Abs(ax) + 0.5);

          ix2= ix2*10000;

          if ( ix1 == ix2) {

            if (x > 0) {

              ix = (int)(ax + 0.5);

            } else {

              ix = (int)(-ax + 0.5);

              ix = -ix;

            }

            f.Write(" {0}", ix);    //%2ld

          } else

            f.Write(" {0}.9E",ax);   //% .9E

             * */



            f.Write(x);

        }



        void WriteNumber(TextWriter f, double x)

        {

#if GMPRATIONAL

  mpz_t zn,zd;



  mpz_init(zn); mpz_init(zd);

  mpq_canonicalize(x);

  mpq_get_num(zn,x);

  mpq_get_den(zd,x);

  fprintf(f," ");

  if (mpz_sgn(zn)==0){

    fprintf(f,"0");

  } else if (mpz_cmp_ui(zd,1U)==0){

    mpz_out_str(f,10,zn);

  } else {

    mpz_out_str(f,10,zn);fprintf(f,"/");mpz_out_str(f,10,zd);

  }

  mpz_clear(zn); mpz_clear(zd);

#else

            WriteReal(f, x);

#endif

        }



#if nemkell

void WriteIncidence(FILE *f, Polyhedra poly)

{

  SetFamily I;



  switch (poly.representation) {

    case Inequality:

       fprintf(f, "ecd_file: Incidence of generators and inequalities\n");

      break;

    case Generator:

       fprintf(f, "icd_file: Incidence of inequalities and generators\n");

      break;



    default:

      break;

  }

  I=CopyIncidence(poly);

  WriteSetFamilyCompressed(f,I);

  FreeSetFamily(I);

}



void WriteAdjacency(FILE *f, Polyhedra poly)

{

  SetFamily A;



  switch (poly.representation) {

    case Inequality:

       fprintf(f, "ead_file: Adjacency of generators\n");

      break;

    case Generator:

       fprintf(f, "iad_file: Adjacency of inequalities\n");

      break;



    default:

      break;

  }

  A=CopyAdjacency(poly);

  WriteSetFamilyCompressed(f,A);

  FreeSetFamily(A);

}

#endif



        void ComputeAinc(Polyhedra poly)

        {

            /* This generates the input incidence array poly.Ainc, and

               two sets: poly.Ared, poly.Adom. 

            */

            int k;

            int i, m1;

            int j;

            bool redundant;

            Matrix M = null;

            double sum, temp;



            sum = 0.0; temp = 0.0;

            if (poly.AincGenerated == true) return;



            M = CopyOutput(poly);

            poly.n = M.rowsize;

            m1 = poly.m1;

            /* this number is same as poly.m, except when

               poly is given by nonhomogeneous inequalty:

               !(poly.homogeneous) && poly.representation==Inequality,

               it is poly.m+1.   See ConeDataLoad.

            */

            poly.Ainc = new BitArray[m1];

            for (i = 1; i <= m1; i++) poly.Ainc[i - 1] = new BitArray(poly.n);

            poly.Ared = new BitArray(m1);

            poly.Adom = new BitArray(m1);



            for (k = 1; k <= poly.n; k++)

            {

                for (i = 1; i <= poly.m; i++)

                {

                    sum = 0.0;

                    for (j = 1; j <= poly.d; j++)

                    {

                        temp = poly.A[i - 1][j - 1] * M.matrix[k - 1][j - 1];

                        sum += temp;

                    }

                    if (sum == 0.0)

                    {

                        poly.Ainc[i - 1].Set(k, true);

                    }

                }

                if (!(poly.homogeneous) && poly.representation == RepresentationType.Inequality)

                {

                    if (M.matrix[k - 1][0] == 0.0)

                    {

                        poly.Ainc[m1 - 1].Set(k, true);  /* added infinity inequality (1,0,0,...,0) */

                    }

                }

            }



            for (i = 1; i <= m1; i++)

            {

                if (poly.Ainc[i - 1].Count == M.rowsize)

                {

                    poly.Adom.Set(i, true);

                }

            }

            for (i = m1; i >= 1; i--)

            {

                if (poly.Ainc[i - 1].Count == 0)

                {

                    redundant = true;

                    poly.Ared.Set(i, true);

                }

                else

                {

                    redundant = false;

                    for (k = 1; k <= m1; k++)

                    {

                        if (k != i && !poly.Ared.Get(k) && !poly.Adom.Get(k) &&

                            poly.Ainc[i - 1].IsSubsetOf(poly.Ainc[k - 1]))

                        {    //****

                            if (!redundant)

                            {

                                redundant = true;

                            }

                            poly.Ared.Set(i, true);

                        }

                    }

                }

            }



            poly.AincGenerated = true;





        }



        bool InputAdjacentQ(Polyhedra poly,

          int i1, int i2)

        /* Before calling this function, RedundantSet must be 

           a set of row indices whose removal results in a minimal

           nonredundant system to represent the input polyhedron,

           DominantSet must be the set of row indices which are

           active at every extreme points/rays.

        */

        {

            bool adj = true;

            int i;

            BitArray common;

            int lastn = 0;



            if (poly.AincGenerated == false) ComputeAinc(poly);

            if (lastn != poly.n)

            {

                //if (lastn >0) set_free(common);

                common = new BitArray(poly.n);

                lastn = poly.n;

            }

            if (poly.Ared.Get(i1) || poly.Ared.Get(i2))

            {

                adj = false;

                return adj;

            }

            if (poly.Adom.Get(i1) || poly.Adom.Get(i2))

            {

                // dominant inequality is considered adjacencent to all others.

                adj = true;

                return adj;

            }

            common = poly.Ainc[i1 - 1].And(poly.Ainc[i2 - 1]);       //**** intersection

            i = 0;

            while (i < poly.m1 && adj == true)

            {

                i++;

                if (i != i1 && i != i2 && !poly.Ared.Get(i) &&

                    !poly.Adom.Get(i) && common.IsSubsetOf(poly.Ainc[i - 1]))

                {

                    adj = false;

                }

            }



            return adj;

        }



#if nemkell

void WriteInputIncidence(FILE *f, Polyhedra poly)

{

  SetFamily I;



  if (poly.AincGenerated==false) ComputeAinc(poly);

  switch (poly.representation) {

  case Inequality:

    fprintf(f,"icd_file: Incidence of inequalities and generators\n");

    break;



  case Generator:

   fprintf(f,"ecd_file: Incidence of generators and inequalities\n");

    break;



  default:

    break;

  }



  I=CopyInputIncidence(poly);

  WriteSetFamilyCompressed(f,I);

  FreeSetFamily(I);

}



void WriteInputAdjacency(FILE *f, Polyhedra poly)

{

  SetFamily A;



  if (poly.AincGenerated==false){

    ComputeAinc(poly);

  }

  switch (poly.representation) {

  case Inequality:

    fprintf(f, "iad_file: Adjacency of inequalities\n");

    break;



  case Generator:

    fprintf(f, "ead_file: Adjacency of generators\n");

    break;



  default:

    break;

  }

  A=CopyInputAdjacency(poly);

  WriteSetFamilyCompressed(f,A);

  FreeSetFamily(A);

}





void WriteProgramDescription(FILE *f)

{

  fprintf(f, "* cddlib: a double description library:%s\n", DDVERSION);

  fprintf(f, "* compiled for %s arithmetic.\n", ARITHMETIC);

  fprintf(f,"* %s\n",COPYRIGHT);

}



void WriteTimes(FILE *f, time_t starttime, time_t endtime)

{

  int ptime,ptime_sec,ptime_minu, ptime_hour;

  

/* 

   ptime=difftime(endtime,starttime); 

   This function is ANSI standard, but not available sometime 

*/

  ptime=(endtime - starttime);     

 /* This is to replace the line above, but it may not give 

    correct time in seconds */ 

  ptime_hour=ptime/3600;

  ptime_minu=(ptime-ptime_hour*3600)/60;

  ptime_sec=ptime%60;

  fprintf(f,"* Computation started at %s",asctime(localtime(&starttime)));

  fprintf(f,"*             ended   at %s",asctime(localtime(&endtime)));

  fprintf(f,"* Total processor time = %ld seconds\n",ptime);

  fprintf(f,"*                      = %ld h %ld m %ld s\n",ptime_hour, ptime_minu, ptime_sec);

}



void WriteDDTimes(FILE *f, Polyhedra poly)

{

  WriteTimes(f,poly.child.starttime,poly.child.endtime);

}



void WriteLPTimes(FILE *f, LPPtr lp)

{

  WriteTimes(f,lp.starttime,lp.endtime);

}



void WriteLPStats(FILE *f)

{

  time_t currenttime;

  

  time(&currenttime);

  

  fprintf(f, "\n*--- Statistics of pivots ---\n");

//#if defined GMPRATIONAL

  fprintf(f, "* f0 = %ld (float basis finding pivots)\n",ddf_statBApivots);

  fprintf(f, "* fc = %ld (float CC pivots)\n",ddf_statCCpivots);

  fprintf(f, "* f1 = %ld (float dual simplex phase I pivots)\n",ddf_statDS1pivots);

  fprintf(f, "* f2 = %ld (float dual simplex phase II pivots)\n",ddf_statDS2pivots);

  fprintf(f, "* f3 = %ld (float anticycling CC pivots)\n",ddf_statACpivots);

  fprintf(f, "* e0 = %ld (exact basis finding pivots)\n",statBApivots);

  fprintf(f, "* ec = %ld (exact CC pivots)\n",statCCpivots);

  fprintf(f, "* e1 = %ld (exact dual simplex phase I pivots)\n",statDS1pivots);

  fprintf(f, "* e2 = %ld (exact dual simplex phase II pivots)\n",statDS2pivots);

  fprintf(f, "* e3 = %ld (exact anticycling CC pivots)\n",statACpivots);

  fprintf(f, "* e4 = %ld (exact basis verification pivots)\n",statBSpivots);

//#else

  fprintf(f, "f0 = %ld (float basis finding pivots)\n",statBApivots);

  fprintf(f, "fc = %ld (float CC pivots)\n",statCCpivots);

  fprintf(f, "f1 = %ld (float dual simplex phase I pivots)\n",statDS1pivots);

  fprintf(f, "f2 = %ld (float dual simplex phase II pivots)\n",statDS2pivots);

  fprintf(f, "f3 = %ld (float anticycling CC pivots)\n",statACpivots);

//#endif

 WriteLPMode(f);

  WriteTimes(f,statStartTime, currenttime);

}



void WriteLPMode(FILE *f)

{

  fprintf(f, "\n* LP solver: ");

  switch (choiceLPSolverDefault) {

    case DualSimplex:

      fprintf(f, "DualSimplex\n");

      break;

    case CrissCross:

      fprintf(f, "Criss-Cross\n");

      break;

    default: break;

  }

  

  fprintf(f, "* Redundancy cheking solver: ");

  switch (choiceRedcheckAlgorithm) {

    case DualSimplex:

      fprintf(f, "DualSimplex\n");

      break;

    case CrissCross:

      fprintf(f, "Criss-Cross\n");

      break;

    default: break;

  }

  

  fprintf(f, "* Lexicographic pivot: ");

  if (choiceLexicoPivotQ)  fprintf(f, " on\n"); 

  else fprintf(f, " off\n"); 



}





void WriteRunningMode(FILE *f, Polyhedra poly)

{

  if (poly.child!=null){

    fprintf(f,"* roworder: ");

    switch (poly.child.HalfspaceOrder) {



    case MinIndex:

      fprintf(f, "minindex\n");

      break;



    case MaxIndex:

      fprintf(f, "maxindex\n");

      break;



    case MinCutoff:

      fprintf(f, "mincutoff\n");

      break;



    case MaxCutoff:

      fprintf(f, "maxcutoff\n");

    break;



    case MixCutoff:

      fprintf(f, "mixcutoff\n");

      break;



    case LexMin:

      fprintf(f, "lexmin\n");

      break;



    case LexMax:

      fprintf(f, "lexmax\n");

      break;



    case RandomRow:

      fprintf(f, "random  %d\n",poly.child.rseed);

      break;



    default: break;

    }

  }

}





void WriteCompletionStatus(FILE *f, Cone cone)

{

  if (cone.Iteration<cone.m && cone.CompStatus==AllFound) {

    fprintf(f,"*Computation completed at Iteration %4ld.\n", cone.Iteration);

  } 

  if (cone.CompStatus == RegionEmpty) {

    fprintf(f,"*Computation completed at Iteration %4ld because the region found empty.\n",cone.Iteration);

  }   

}



void WritePolyFile(FILE *f, Polyhedra poly)

{

  WriteAmatrix(f,poly.A,poly.m,poly.d);

}





void WriteErrorMessages(FILE *f, ErrorType Error)

{

  switch (Error) {



  case DimensionTooLarge:

    fprintf(f, "*Input Error: Input matrix is too large:\n");

    fprintf(f, "*Please increase MMAX and/or NMAX in the source code and recompile.\n");

    break;



  case IFileNotFound:

    fprintf(f, "*Input Error: Specified input file does not exist.\n");

    break;



  case OFileNotOpen:

    fprintf(f, "*Output Error: Specified output file cannot be opened.\n");

    break;



  case NegativeMatrixSize:

    fprintf(f, "*Input Error: Input matrix has a negative size:\n");

    fprintf(f, "*Please check rowsize or colsize.\n");

    break;



  case ImproperInputFormat:

    fprintf(f,"*Input Error: Input format is not correct.\n");

    fprintf(f,"*Format:\n");

    fprintf(f," begin\n");

    fprintf(f,"   m   n  NumberType(real, rational or integer)\n");

    fprintf(f,"   b  -A\n");

    fprintf(f," end\n");

    break;



  case EmptyVrepresentation:

    fprintf(f, "*Input Error: V-representation is empty:\n");

    fprintf(f, "*cddlib does not accept this trivial case for which output can be any inconsistent system.\n");

    break;



  case EmptyHrepresentation:

    fprintf(f, "*Input Error: H-representation is empty.\n");

    break;



  case EmptyRepresentation:

    fprintf(f, "*Input Error: Representation is empty.\n");

    break;



  case NoLPObjective:

    fprintf(f, "*LP Error: No LP objective (max or min) is set.\n");

    break;



  case NoRealNumberSupport:

    fprintf(f, "*LP Error: The binary (with GMP Rational) does not support Real number input.\n");

    fprintf(f, "         : Use a binary compiled without -DGMPRATIONAL option.\n");

    break;



 case NotAvailForH:

    fprintf(f, "*Error: A function is called with H-rep which does not support an H-representation.\n");

    break;



 case NotAvailForV:

    fprintf(f, "*Error: A function is called with V-rep which does not support an V-representation.\n");

    break;



 case CannotHandleLinearity:

    fprintf(f, "*Error: The function called cannot handle linearity.\n");

    break;



 case RowIndexOutOfRange:

    fprintf(f, "*Error: Specified row index is out of range\n");

    break;



 case ColIndexOutOfRange:

    fprintf(f, "*Error: Specified column index is out of range\n");

    break;



 case LPCycling:

    fprintf(f, "*Error: Possibly an LP cycling occurs.  Use the Criss-Cross method.\n");

    break;

    

 case NumericallyInconsistent:

    fprintf(f, "*Error: Numerical inconsistency is found.  Use the GMP exact arithmetic.\n");

    break;

    

  case NoError:

    fprintf(f,"*No Error found.\n");

    break;

  }

}

#endif



        SetFamily CopyIncidence(Polyhedra poly)

        {

            SetFamily F = null;

            int k;

            int i;



            if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return F;

            if (poly.AincGenerated == false) ComputeAinc(poly);

            F = CreateSetFamily(poly.n, poly.m1);

            for (i = 1; i <= poly.m1; i++)

                for (k = 1; k <= poly.n; k++)

                    if (poly.Ainc[i - 1].Get(k)) F.set[k - 1].Set(i, true);



            return F;

        }



        SetFamily CopyInputIncidence(Polyhedra poly)

        {

            int i;

            SetFamily F = null;



            if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return F;

            if (poly.AincGenerated == false) ComputeAinc(poly);

            F = CreateSetFamily(poly.m1, poly.n);

            for (i = 0; i < poly.m1; i++)

            {

                poly.Ainc[i].CopyTo(F.set[i]);

            }



            return F;

        }



        SetFamily CopyAdjacency(Polyhedra poly)

        {

            Ray RayPtr1, RayPtr2;

            SetFamily F = null;

            int pos1, pos2;

            int lstart, k, n;

            BitArray linset, allset;

            bool adj;



            if (poly.n == 0 && poly.homogeneous && poly.representation == RepresentationType.Inequality)

            {

                n = 1; /* the origin (the unique vertex) should be output. */

            }

            else n = poly.n;

            linset = new BitArray(n);

            allset = new BitArray(n);

            if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return F;

            F = CreateSetFamily(n, n);

            if (n <= 0) return F;

            poly.child.LastRay.Next = null;

            for (RayPtr1 = poly.child.FirstRay, pos1 = 1; RayPtr1 != null;

                          RayPtr1 = RayPtr1.Next, pos1++)

            {

                for (RayPtr2 = poly.child.FirstRay, pos2 = 1; RayPtr2 != null;

                                RayPtr2 = RayPtr2.Next, pos2++)

                {

                    if (RayPtr1 != RayPtr2)

                    {

                        CheckAdjacency(poly.child, RayPtr1, RayPtr2, out adj);

                        if (adj)

                        {

                            F.set[pos1 - 1].Set(pos2, true);

                        }

                    }

                }

            }

            lstart = poly.n - poly.ldim + 1;

            allset = allset.Not();  /* allset is set to the ground set. */

            for (k = lstart; k <= poly.n; k++)

            {

                linset.Set(k, true);     /* linearity set */

                allset.CopyTo(F.set[k - 1]);  /* linearity generator is adjacent to all */

            }

            for (k = 1; k < lstart; k++)

            {

                F.set[k - 1] = F.set[k - 1].Or(linset); // union

                /* every generator is adjacent to all linearity generators */

            }



            return F;

        }



        SetFamily CopyInputAdjacency(Polyhedra poly)

        {

            int i, j;

            SetFamily F = null;



            if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return F;

            if (poly.AincGenerated == false) ComputeAinc(poly);

            F = CreateSetFamily(poly.m1, poly.m1);

            for (i = 1; i <= poly.m1; i++)

            {

                for (j = 1; j <= poly.m1; j++)

                {

                    if (i != j && InputAdjacentQ(poly, i, j))

                    {

                        F.set[i - 1].Set(j, true);

                    }

                }

            }



            return F;

        }



        Matrix CopyOutput(Polyhedra poly)

        {

            Ray Ray;

            Matrix M = null;

            int i = 0, total;

            int j, j1;

            double b;

            RepresentationType outputrep = RepresentationType.Inequality;

            bool outputorigin = false;



            b = 0.0;

            total = poly.child.LinearityDim + poly.child.FeasibleRayCount;



            if (poly.child.d <= 0 || poly.child.newcol[1] == 0) total = total - 1;

            if (poly.representation == RepresentationType.Inequality) outputrep = RepresentationType.Generator;

            if (total == 0 && poly.homogeneous && poly.representation == RepresentationType.Inequality)

            {

                total = 1;

                outputorigin = true;

                /* the origin (the unique vertex) should be output. */

            }

            if (poly.child == null || poly.child.CompStatus != CompStatusType.AllFound) return M;



            M = CreateMatrix(total, poly.d);

            Ray = poly.child.FirstRay;

            while (Ray != null)

            {

                if (Ray.feasible)

                {

                    CopyRay(M.matrix[i], poly.d, Ray, outputrep, poly.child.newcol);

                    i++;  /* 086 */

                }

                Ray = Ray.Next;

            }

            for (j = 2; j <= poly.d; j++)

            {

                if (poly.child.newcol[j] == 0)

                {

                    /* original column j is dependent on others and removed for the cone */

                    b = poly.child.Bsave[0][j - 1];

                    if (outputrep == RepresentationType.Generator && b > 0.0)

                    {

                        M.matrix[i][0] = 1.0;  /* Normalize */

                        for (j1 = 1; j1 < poly.d; j1++)

                            M.matrix[i][j1] = poly.child.Bsave[j1][j - 1] / b;

                    }

                    else

                    {

                        for (j1 = 0; j1 < poly.d; j1++)

                            M.matrix[i][j1] = poly.child.Bsave[j1][j - 1];

                    }

                    M.linset.Set(i + 1, true);

                    i++;

                }

            }

            if (outputorigin)

            {

                /* output the origin for homogeneous H-polyhedron with no rays. */

                M.matrix[0][0] = 1.0;

                for (j = 1; j < poly.d; j++)

                {

                    M.matrix[0][j] = 0.0;

                }

            }

            MatrixIntegerFilter(M);

            if (poly.representation == RepresentationType.Inequality)

                M.representation = RepresentationType.Generator;

            else

                M.representation = RepresentationType.Inequality;





            return M;

        }



        Matrix CopyInput(Polyhedra poly)

        {

            Matrix M = null;

            int i;



            M = CreateMatrix(poly.m, poly.d);

            CopyAmatrix(M.matrix, poly.A, poly.m, poly.d);



            for (i = 1; i <= poly.m; i++)

                if (poly.EqualityIndex[i] == 1) M.linset.Set(i, true);

            MatrixIntegerFilter(M);

            if (poly.representation == RepresentationType.Generator)

                M.representation = RepresentationType.Generator;

            else

                M.representation = RepresentationType.Inequality;

            return M;

        }



        public Matrix CopyGenerators(Polyhedra poly)

        {

            Matrix M = null;



            if (poly.representation == RepresentationType.Generator)

            {

                M = CopyInput(poly);

            }

            else

            {

                M = CopyOutput(poly);

            }

            return M;

        }



        Matrix CopyInequalities(Polyhedra poly)

        {

            Matrix M = null;



            if (poly.representation == RepresentationType.Inequality)

            {

                M = CopyInput(poly);

            }

            else

            {

                M = CopyOutput(poly);

            }

            return M;

        }



#if nemkell

/****************************************************************************************/

/*  rational number (a/b) read is taken from Vinci by Benno Bueeler and Andreas Enge    */

/****************************************************************************************/

void sread_rational_value (char *s, double value)

   /* reads a rational value from the specified string "s" and assigns it to "value"    */

   

{

   char     *numerator_s=null, *denominator_s=null, *position;

   int      sign = 1;

   double   numerator, denominator;

//#if defined GMPRATIONAL

   mpz_t znum, zden;

//#else

   double rvalue;

//#endif

  

   /* determine the sign of the number */

   numerator_s = s;

   if (s [0] == '-')

   {  sign = -1;

      numerator_s++;

   }

   else if (s [0] == '+')

      numerator_s++;

      

   /* look for a sign '/' and eventually split the number in numerator and denominator */

   position = strchr (numerator_s, '/');

   if (position != null)

   {  *position = '\0'; /* terminates the numerator */

      denominator_s = position + 1;

   };



   /* determine the floating point values of numerator and denominator */

   numerator=atol (numerator_s);

  

   if (position != null)

   { 

     denominator=atol (denominator_s);  

   }

   else denominator = 1;



/* 

   fprintf(stderr,"\nrational_read: numerator %f\n",numerator);

   fprintf(stderr,"rational_read: denominator %f\n",denominator);

   fprintf(stderr,"rational_read: sign %d\n",sign); 

*/



//#if defined GMPRATIONAL

   mpz_init_set_str(znum,numerator_s,10);

   if (sign<0) mpz_neg(znum,znum);

   mpz_init(zden); mpz_set_ui(zden,1);

   if (denominator_s!=null) mpz_init_set_str(zden,denominator_s,10);

   mpq_set_num(value,znum); mpq_set_den(value,zden);

   mpq_canonicalize(value);

   mpz_clear(znum); mpz_clear(zden);

   /* num=(int)sign * numerator; */

   /* den=(unsigned int) denominator; */

   /* mpq_set_si(value, num, den); */

//#elif defined GMPFLOAT

   rvalue=sign * numerator/ (signed int) denominator;

   mpf_set_d(value, rvalue);

//#else

   rvalue=sign * numerator/ (signed int) denominator;

   dset_d(value, rvalue);

//#endif

   if (debug) {

     fprintf(stderr,"rational_read: "); 

     WriteNumber(stderr,value); fprintf(stderr,"\n");

   }

}

   



void fread_rational_value (FILE *f, double value)

   /* reads a rational value from the specified file "f" and assigns it to "value"      */

   

{

   char     number_s [255];

   double rational_value;

   

   init(rational_value);

   fscanf(f, "%s ", number_s);

   sread_rational_value (number_s, rational_value);

   set(value,rational_value);

   clear(rational_value);

}

   

/****************************************************************************************/





/* end of cddio.c */

#endif



    }

}