--- a
+++ b/feasible_joint_stiffness/lrslib-062/lrslib.c
@@ -0,0 +1,5701 @@
+/* lrslib.c     library code for lrs                     */
+
+/* modified by Gary Roumanis for multithread plrs compatability */
+/* truncate needs mod to supress last pivot */
+/* need to add a test for non-degenerate pivot step in reverse I guess */
+/* Copyright: David Avis 2005,2011 avis@cs.mcgill.ca         */
+
+/* This program is free software; you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation; either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program; if not, write to the Free Software
+   Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA  02110-1335, USA.
+ */
+
+
+
+#include <stdio.h>
+#include <string.h>
+#include "lrslib.h"
+
+/* Globals; these need to be here, rather than lrslib.h, so they are
+   not multiply defined. */
+
+FILE *lrs_cfp;			/* output file for checkpoint information       */
+FILE *lrs_ifp;			/* input file pointer       */
+FILE *lrs_ofp;			/* output file pointer      */
+
+
+static unsigned long dict_count, dict_limit, cache_tries, cache_misses;
+
+/* Variables and functions global to this file only */
+static long lrs_checkpoint_seconds = 0;
+
+static long lrs_global_count = 0;	/* Track how many lrs_dat records are 
+					   allocated */
+
+static lrs_dat_p *lrs_global_list[MAX_LRS_GLOBALS + 1];
+
+static lrs_dic *new_lrs_dic (long m, long d, long m_A);
+
+
+static void cache_dict (lrs_dic ** D_p, lrs_dat * global, long i, long j);
+static long check_cache (lrs_dic ** D_p, lrs_dat * global, long *i_p, long *j_p);
+static void save_basis (lrs_dic * D, lrs_dat * Q);
+
+static void lrs_dump_state ();
+
+static void pushQ (lrs_dat * global, long m, long d, long m_A);
+
+#ifdef TIMES
+static void ptimes ();
+static double get_time();
+#endif
+
+
+/*******************************/
+/* signals handling            */
+/*******************************/
+#ifdef SIGNALS
+static void checkpoint ();
+static void die_gracefully ();
+static void setup_signals ();
+static void timecheck ();
+#endif
+
+/*******************************/
+/* functions  for external use */
+/*******************************/
+
+/*******************************************************/
+/* lrs_main is driver for lrs.c does H/V enumeration   */
+/* showing function calls intended for public use      */
+/*******************************************************/
+long
+lrs_main (int argc, char *argv[])
+
+{
+	
+	lrs_dic *P;			/* structure for holding current dictionary and indices */
+	lrs_dat *Q;			/* structure for holding static problem data            */
+
+	lrs_mp_vector output;		/* holds one line of output; ray,vertex,facet,linearity */
+	lrs_mp_matrix Lin;		/* holds input linearities if any are found             */
+	long col;			/* output column index for dictionary                   */
+	long startcol = 0;
+	long prune = FALSE;		/* if TRUE, getnextbasis will prune tree and backtrack  */
+
+/* global variables lrs_ifp and lrs_ofp are file pointers for input and output   */
+/* they default to stdin and stdout, but may be overidden by command line parms. */
+
+
+/***************************************************
+ Step 0: 
+  Do some global initialization that should only be done once,
+  no matter how many lrs_dat records are allocated. db
+
+***************************************************/
+
+#ifdef PLRS
+	if (!lrs_mp_init (ZERO, stdin, stdout))  /* initialize arithmetic */
+		exit(1);
+#else
+  if ( !lrs_init ("\n*lrs:"))
+    return 1;
+  printf("\n%s",AUTHOR);
+#endif
+
+/*********************************************************************************/
+/* Step 1: Allocate lrs_dat, lrs_dic and set up the problem                      */
+/*********************************************************************************/
+#ifdef PLRS
+Q = lrs_alloc_dat ("");	/* allocate and init structure for static problem data */
+
+	std::ifstream input_file;
+	input_file.open(argv[0]); /* Open input file */
+	plrs_read_dat(Q, input_file);	/* read first part of problem data to get dimensions and problem type: H- or V- input representation     */
+	
+	P = lrs_alloc_dic (Q);	/* allocate and initialize lrs_dic                     */
+	if (P == NULL)
+		return 1;
+	plrs_read_dic (P, Q, input_file); /* read remainder of input to setup P and Q            */
+#else
+  Q = lrs_alloc_dat ("LRS globals");	/* allocate and init structure for static problem data */
+
+  if (Q == NULL)
+    return 1;
+
+  if (!lrs_read_dat (Q, argc, argv))	/* read first part of problem data to get dimensions   */
+    return 1;                   	/* and problem type: H- or V- input representation     */
+
+  P = lrs_alloc_dic (Q);	/* allocate and initialize lrs_dic                     */
+  if (P == NULL)
+    return 1;
+
+  if (!lrs_read_dic (P, Q))	/* read remainder of input to setup P and Q            */
+    return 1;
+#endif
+
+  output = lrs_alloc_mp_vector (Q->n);	/* output holds one line of output from dictionary     */
+
+
+/*********************************************************************************/
+/* Step 2: Find a starting cobasis from default of specified order               */
+/*         P is created to hold  active dictionary data and may be cached        */
+/*         Lin is created if necessary to hold linearity space                   */
+/*         Print linearity space if any, and retrieve output from first dict.    */
+/*********************************************************************************/
+
+  if (!lrs_getfirstbasis (&P, Q, &Lin, FALSE))
+    return 1;
+
+  /* Pivot to a starting dictionary                      */
+  /* There may have been column redundancy               */
+  /* If so the linearity space is obtained and redundant */
+  /* columns are removed. User can access linearity space */
+  /* from lrs_mp_matrix Lin dimensions nredundcol x d+1  */
+
+
+
+  if (Q->homogeneous && Q->hull)
+    startcol++;			/* col zero not treated as redundant   */
+
+	if(!Q->restart)
+		for (col = startcol; col < Q->nredundcol; col++)	/* print linearity space               */
+			lrs_printoutput (Q, Lin[col]);	/* Array Lin[][] holds the coeffs.     */
+
+
+
+/*********************************************************************************/
+/* Step 3: Terminate if lponly option set, otherwise initiate a reverse          */
+/*         search from the starting dictionary. Get output for each new dict.    */
+/*********************************************************************************/
+
+
+
+  /* We initiate reverse search from this dictionary       */
+  /* getting new dictionaries until the search is complete */
+  /* User can access each output line from output which is */
+  /* vertex/ray/facet from the lrs_mp_vector output         */
+  /* prune is TRUE if tree should be pruned at current node */
+  prune=lrs_checkbound(P,Q);
+  do
+    {
+
+//2015.6.5   after maxcobases reached, generate subtrees that have not been enumerated
+
+     if ((Q->maxcobases > 0) &&  (Q->count[2] >=Q->maxcobases))
+        {
+          if(!lrs_leaf(P,Q))      
+                                       /* do not return cobases of a leaf */
+             lrs_printcobasis(P,Q,ZERO);
+            
+          prune=TRUE;
+
+        }     // if Q-> maxcobases...
+
+         for (col = 0; col <= P->d; col++)          /* print output vertex/ray if any */
+	   if (lrs_getsolution (P, Q, output, col))
+	       lrs_printoutput (Q, output);
+
+  }while (!Q->lponly && lrs_getnextbasis (&P, Q, prune));  // do ...
+
+  if (Q->lponly)
+    lrs_lpoutput(P,Q,output);
+  else
+    lrs_printtotals (P, Q);	/* print final totals, including estimates       */
+
+  lrs_clear_mp_vector(output, Q->n);
+
+/* 2015.9.16  fix memory leaks on Gcd Lcm Lin */
+  if(Q->nredundcol > 0)
+     lrs_clear_mp_matrix(Lin,Q->nredundcol,Q->n);
+  if(Q->runs > 0)
+    { 
+      free(Q->isave);
+      free(Q->jsave);
+    }
+  long savem=P->m;              /* need this to clear Q*/
+  lrs_free_dic (P,Q);           /* deallocate lrs_dic */
+  Q->m=savem;
+
+  lrs_free_dat (Q);             /* deallocate lrs_dat */
+
+#ifndef PLRS
+  	lrs_close ("lrs:");
+#endif
+
+  return 0;
+}
+/*********************************************/
+/* end of model test program for lrs library */
+/*********************************************/
+
+
+
+
+/***********************************/
+/* 		PLRS		   */
+/***********************************/
+
+#ifdef PLRS
+
+void plrs_read_dat (lrs_dat * Q, std::ifstream &input_file)
+{
+	
+	string line;
+	bool begin = false;
+
+	if(input_file.is_open()){
+		while(input_file.good()){
+			getline(input_file, line);
+			
+			if(line.find("*") == 0){
+				//Ignore lines starting with *
+			}else if (line.find("H-representation") != string::npos){
+				Q->hull = FALSE;
+			}else if(line.find("hull")!= string::npos || line.find("V-representation")!= string::npos){
+				Q->hull = TRUE;
+		   		Q->polytope = TRUE;			
+			}else if(line.find("digits")!= string::npos){
+				long dec_digits;
+				istringstream ss(line);
+				if(!(ss>>dec_digits) && !lrs_set_digits(dec_digits)){
+					printf("\nError reading digits data!\n");				
+					exit(1);				
+				}
+			}else if(line.find("nonnegative")!= string::npos){
+				 Q->nonnegative = TRUE;
+			}else if(line.find("linearity") != string::npos){
+				//Remove the following characters
+				char chars[] = "linearity";
+				for(unsigned int i = 0; i < sizeof(chars); ++i){
+					line.erase(remove(line.begin(), line.end(), chars[i]), line.end());
+				}
+				
+				plrs_readlinearity (Q, line);
+			}else if(line.find("begin")!= string::npos){
+				begin = true;
+				break;
+
+
+			}else{
+				//Q->name = line.c_str();
+			}
+		}
+
+                if(Q->hull) 
+                      Q->getvolume=TRUE;
+
+		if(!begin){
+			printf("\nNo begin line!\n");   
+			fprintf(lrs_ofp,"\nNo begin line!\n");   
+			exit(1);
+		}
+		
+		
+		getline(input_file, line);
+		istringstream ss(line);
+		string type;
+
+
+		if(!(ss >> Q->m >> Q->n >> type)){
+		printf("\nNo data in file!\n");
+			exit(1);
+		}
+
+		if(!type.find("integer") && !type.find("rational")){
+			printf("\nData type must be integer or rational!\n");
+			exit(1);
+		}
+		
+	}else{
+		printf("\nError reading input file!\n");
+		exit(1);
+	}
+
+
+	if (Q->m == 0)
+	{
+		printf("\nNo input given!\n");
+		exit(1);
+	}
+	/* inputd may be reduced in preprocessing of linearities and redund cols */
+}
+
+/* read constraint matrix and set up problem and dictionary  */
+void plrs_read_dic (lrs_dic * P, lrs_dat * Q, std::ifstream &input_file)
+{
+
+	lrs_mp Temp, mpone;
+	lrs_mp_vector oD;		/* Denom for objective function */
+
+	long i, j;
+	string line;
+
+
+	/* assign local variables to structures */
+
+	lrs_mp_matrix A;
+	lrs_mp_vector Gcd, Lcm;
+	long hull = Q->hull;
+	long m, d;
+
+	lrs_alloc_mp(Temp); lrs_alloc_mp(mpone);
+	A = P->A;
+	m = Q->m;
+	d = Q->inputd;
+
+	Gcd = Q->Gcd;
+	Lcm = Q->Lcm;
+
+	oD = lrs_alloc_mp_vector (d);
+
+
+	itomp (ONE, mpone);
+	itomp (ONE, A[0][0]);
+	itomp (ONE, Lcm[0]);
+	itomp (ONE, Gcd[0]);
+
+
+		
+	for (i = 1; i <= m; i++)	/* read in input matrix row by row                 */
+	{
+
+		itomp (ONE, Lcm[i]);	/* Lcm of denominators */
+		itomp (ZERO, Gcd[i]);	/* Gcd of numerators */
+
+		if(!input_file.good()){
+			printf("\nInput data incorrectly formatted\n");
+			exit(1);
+		}
+
+/* allow embedded CRs in multiline input for matrix rows */
+/* there must be an easier way ....  but this seems to work */
+                j=hull;
+                while (j <= d)     /* hull data copied to cols 1..d */
+                {
+                 if(!input_file.good()){
+                        printf("\nInput incorrectly formatted\n");
+                        exit(1);
+                }
+
+                getline(input_file, line);
+                istringstream ss(line);
+                const char* ptr1;
+                int string_length;
+                string_length=1;
+                while ((j<=d) && (string_length !=0))
+                  {
+                        string rat;
+                           ss>>rat;
+                        ptr1=rat.c_str();
+                        string_length=strlen(ptr1);
+                        if (string_length!=0)
+                        { 
+                          if (plrs_readrat (A[i][j], A[0][j], ptr1))
+                                lcm (Lcm[i], A[0][j]);  /* update lcm of denominators */
+                          copy (Temp, A[i][j]);
+                          gcd (Gcd[i], Temp);     /* update gcd of numerators   */
+                          j++;
+                         }
+                   }
+                 ss.clear();
+                 }
+
+
+
+		if (hull)
+		{
+			itomp (ZERO, A[i][0]);	/*for hull, we have to append an extra column of zeroes */
+			if (!one (A[i][1]) || !one (A[0][1]))		/* all rows must have a one in column one */
+				Q->polytope = FALSE;
+		}
+
+		if (!zero (A[i][hull]))	/* for H-rep, are zero in column 0     */
+			Q->homogeneous = FALSE;	/* for V-rep, all zero in column 1     */
+
+
+		storesign (Gcd[i], POS);
+		storesign (Lcm[i], POS);
+
+		if (mp_greater (Gcd[i], mpone) || mp_greater (Lcm[i], mpone))
+			for (j = 0; j <= d; j++)
+		 	{
+		 		divint (A[i][j], Gcd[i], Temp);	/*reduce numerators by Gcd  */
+		 		mulint (Lcm[i], Temp, Temp);	/*remove denominators */
+		 		divint (Temp, A[0][j], A[i][j]);	/*reduce by former denominator */
+		 	}
+	}
+
+
+
+	/* 2010.4.26 patch */
+	if(Q->nonnegative)    /* set up Gcd and Lcm for nonexistent nongative inequalities */
+		for (i=m+1;i<=m+d;i++)
+		{ 
+			itomp (ONE, Lcm[i]);
+			itomp (ONE, Gcd[i]);
+		}
+
+
+	//Make new output node for nonfatal option errors
+	//Make stream to collect prat / pmp data
+	stringstream out_stream;
+
+	if (Q->homogeneous && Q->verbose)
+	{
+		out_stream<<"*Input is homogeneous, column 1 not treated as redundant"<<endl;
+	}
+
+
+	while(input_file.good()){
+			getline(input_file, line);
+			if(line.find("*") == 0){
+				//Ignore lines starting with *
+
+			}else if(line.find("startingcobasis") != string::npos){
+				if(Q->nonnegative){
+					out_stream<<"*Starting cobasis incompatible with nonegative option:skipped"<<endl;
+				}else{
+					
+					Q->givenstart = TRUE;
+					istringstream ss(line);
+					//Trim first word
+					string str;
+					ss >>str;
+					//make string out of facts
+					stringstream facets;
+					facets << ss.rdbuf();
+
+					//Readfacets
+					plrs_readfacets(Q, Q->inequality,facets.str());
+				
+				}	
+		
+			}else if(line.find("restart") != string::npos){
+				
+				Q->restart = TRUE;	
+				istringstream ss(line);
+				//Trim first word
+				string str;
+				ss >>str;
+				//Pipe restart data from string stream
+				if(Q->voronoi){
+					if(!(ss>>Q->count[1]>>Q->count[0]>>Q->count[2]>>P->depth)){
+						printf("\nError reading restart data!\n");
+						exit(1);
+					}
+
+				}else if(hull){
+					if(!(ss>>Q->count[0]>>Q->count[2]>>P->depth)){
+						printf("\nError reading restart data!\n");
+						exit(1);
+					}
+				}else{
+					if(!(ss>>Q->count[1]>>Q->count[0]>>Q->count[2]>>P->depth)){
+						printf("\nError reading restart data!\n");
+						exit(1);					
+					}
+				}
+				//Store starting counts to calculate totals
+				for (int i = 0; i<5; i++){
+					Q->startcount[i] = Q->count[i];
+				}
+	
+				//Make string out of facets
+				stringstream facets;
+				facets << ss.rdbuf();
+				plrs_readfacets(Q,Q->facet,facets.str());
+
+			}else if(line.find("geometric") != string::npos){
+				if(hull && !Q->voronoi)
+					out_stream<<"*Geometric option for H-representation or voronoi only, skipped"<<endl;
+				else
+					Q->geometric = TRUE;
+
+			}else if(line.find("allbases") != string::npos){
+				Q->allbases = TRUE;
+
+			}else if(line.find("countonly") != string::npos){
+				Q->countonly = TRUE;
+
+			}else if(line.find("incidence") != string::npos){
+				Q->incidence = TRUE;
+
+			}else if(line.find("#incidence") != string::npos){
+				Q->printcobasis = TRUE;
+
+			}else if(line.find("printcobasis") != string::npos){
+				istringstream ss(line);
+				//Trim first word
+				string str;
+				ss >>str;
+				if(!(ss>>Q->frequency))
+					Q->frequency = 0;
+				Q->printcobasis = TRUE;
+
+			}else if(line.find("printslack") != string::npos){
+				Q->printslack = TRUE;
+
+			}else if(line.find("maxdepth") != string::npos){
+				istringstream ss(line);
+				//Trim first word
+				string str;
+				ss >>str;
+				if(!(ss>>Q->maxdepth)){
+					Q->maxdepth = 1;
+				}
+				
+
+			}else if(line.find("maxoutput") != string::npos){
+				istringstream ss(line);
+				//Trim first word
+				string str;
+				ss >>str;
+				if(!(ss>>Q->maxoutput)){
+					Q->maxoutput = 100;
+				}
+				
+			}else if(line.find("maxcobases") != string::npos){
+				istringstream ss(line);
+				//Trim first word
+				string str;
+				ss >>str;
+				if(!(ss>>Q->maxcobases)){
+					Q->maxcobases = 1000;
+				}
+				
+                        }else if(line.find("lponly")!= string::npos){
+                                  printf("\nError: lponly option not supported - use lrs!\n");
+                                  exit(1);
+
+			}else if(line.find("mindepth") != string::npos){
+				istringstream ss(line);
+				//Trim first word
+				string str;
+				ss >>str;
+				if(!(ss>>Q->mindepth)){
+					Q->mindepth = 0;
+				}
+			
+			}else if(line.find("estimates") != string::npos){
+				istringstream ss(line);
+				//Trim first word
+				string str;
+				ss >>str;
+				if (!(ss>>Q->runs)){
+					Q->runs=1;
+				}	
+
+			}else if(line.find("subtreesize") != string::npos){
+				istringstream ss(line);
+				//Trim first word
+				string str;
+				ss >>str;
+				if (!(ss>>Q->subtreesize)){
+					Q->subtreesize=MAXD;
+				}	
+
+
+			}else if(line.find("truncate") != string::npos){
+				if(!hull)
+					Q->truncate = TRUE;
+				else
+					out_stream<<"*Truncate option for H-representation only, skipped"<<endl;
+
+			}else if(line.find("verbose") != string::npos){
+				Q->verbose = TRUE;
+
+			}else if(line.find("bound") != string::npos){
+				istringstream ss(line);
+				//Trim first word
+				string str;
+				ss>>str;
+				//get rational number
+				ss>>str;
+				plrs_readrat(Q->boundn, Q->boundd, str.c_str());
+				Q->bound = TRUE;
+
+			}else if(line.find("nonnegative") != string::npos){
+				out_stream<<"*Nonnegative option must come before begin line, skipped"<<endl;
+			}else if(line.find("seed") != string::npos){
+				istringstream ss(line);
+				if(!(ss>>Q->seed)){
+					Q->seed = 3142;
+				}
+				
+
+			}else if(line.find("voronoi") != string::npos || line.find("Voronoi") != string::npos){
+				if(!hull)
+					out_stream<<"*voronoi requires V-representation - option skipped"<<endl;
+				else{
+					Q->voronoi = TRUE;
+					Q->polytope = FALSE;
+				}
+			}
+	}
+
+        if (Q->restart && Q->maxcobases > 0) //2015.4.3 adjust for restart
+               Q->maxcobases = Q->maxcobases + Q->count[2];
+          
+	if (Q->incidence)
+	{
+		Q->printcobasis = TRUE;
+		/* 2010.5.7    No need to reset this, as it may have been set by printcobasis */
+		/*    Q->frequency    = ZERO;                     */                                                      
+	}
+
+	lrs_clear_mp(Temp); lrs_clear_mp(mpone);
+	lrs_clear_mp_vector (oD,d);
+
+	//post output in a nonblocking manner (a consumer thread will manage output)
+	post_output("options warning", out_stream.str().c_str());
+
+}
+
+
+/* read and check facet list for obvious errors during start/restart */
+/* this must be done after linearity option is processed!! */
+void plrs_readfacets (lrs_dat * Q, long facet[], string facets)
+{
+	long i, j;
+	/* assign local variables to structures */
+	long m, d;
+	long *linearity = Q->linearity;
+	m = Q->m;
+	d = Q->inputd;
+
+	istringstream ss(facets);
+	for (j = Q->nlinearity; j < d; j++)	/* note we place these after the linearity indices */
+	{
+		if(!(ss>>facet[j]))
+		{
+			return;
+		}
+
+
+		//fprintf (lrs_ofp, " %ld", facet[j]);
+		/* 2010.4.26 nonnegative option needs larger range of indices */
+		if(Q->nonnegative)
+			if (facet[j] < 1 || facet[j] > m+d)
+			{
+				printf("\nStart/Restart cobasic indices must be in range 1 .. %ld \n",m+d);
+				exit(1);
+			}
+			if(!Q->nonnegative)
+		 		if (facet[j] < 1 || facet[j] > m)
+		 		{
+		  			printf("\nStart/Restart cobasic indices must be in range 1 .. %ld \n",m);
+		  			exit(1);
+		  		}
+			for (i = 0; i < Q->nlinearity; i++)
+				if (linearity[i] == facet[j])
+		  		{
+		    			 printf("\nStart/Restart cobasic indices should not include linearities\n");;
+		    			exit(1);
+		  		}
+				/* bug fix 2011.8.1  reported by Steven Wu*/
+				for (i = Q->nlinearity; i < j; i++)
+				/* end bug fix 2011.8.1 */
+
+			if (facet[i] == facet[j])
+		  	{
+		   		 printf("\nStart/Restart cobasic indices must be distinct\n");
+		    		exit(1);
+		  	}
+	}
+}				/* end of readfacets */
+
+extern int PLRS_DEBUG;
+#endif
+
+
+
+
+
+/*******************************************************/
+/* redund_main is driver for redund.c, removes all     */
+/* redundant rows from an H or V-representation        */
+/* showing function calls intended for public use      */
+/*******************************************************/
+long
+redund_main (int argc, char *argv[])
+
+{
+  lrs_mp_matrix Ain;		/* holds a copy of the input matrix to output at the end */
+
+  long *redineq;		/* redineq[i]=0 if ineq i non-red,1 if red,2 linearity  */
+  long ineq;			/* input inequality number of current index             */
+
+  lrs_dic *P;			/* structure for holding current dictionary and indices */
+  lrs_dat *Q;			/* structure for holding static problem data            */
+
+  lrs_mp_matrix Lin;		/* holds input linearities if any are found             */
+
+  long i, j, d, m;
+  long nlinearity;		/* number of linearities in input file                  */
+  long nredund;			/* number of redundant rows in input file               */
+  long lastdv;
+  long debug;
+  long index;			/* basic index for redundancy test */
+
+/* global variables lrs_ifp and lrs_ofp are file pointers for input and output   */
+/* they default to stdin and stdout, but may be overidden by command line parms. */
+/* Lin is global 2-d array for linearity space if it is found (redund columns)   */
+
+  lrs_ifp = stdin;
+  lrs_ofp = stdout;
+/***************************************************
+ Step 0: 
+  Do some global initialization that should only be done once,
+  no matter how many lrs_dat records are allocated. db
+
+***************************************************/
+
+  if ( !lrs_init ("\n*redund:"))
+    return 1;
+
+  printf ("\n");
+  printf (AUTHOR);
+
+/*********************************************************************************/
+/* Step 1: Allocate lrs_dat, lrs_dic and set up the problem                      */
+/*********************************************************************************/
+
+  Q = lrs_alloc_dat ("LRS globals");	/* allocate and init structure for static problem data */
+
+  if (Q == NULL)
+    return 1;
+
+  if (!lrs_read_dat (Q, argc, argv))	/* read first part of problem data to get dimensions   */
+    return 1;                         	/* and problem type: H- or V- input representation     */
+
+  P = lrs_alloc_dic (Q);	/* allocate and initialize lrs_dic                     */
+  if (P == NULL)
+    return 1;
+
+  if (!lrs_read_dic (P, Q))	/* read remainder of input to setup P and Q            */
+    return 1;
+
+/* if non-negative flag is set, non-negative constraints are not input */
+/* explicitly, and are not checked for redundancy                      */
+
+  m = P->m_A;              /* number of rows of A matrix */   
+  d = P->d;
+  debug = Q->debug;
+
+  redineq = (long *) calloc ((m + 1), sizeof (long));
+  Ain = lrs_alloc_mp_matrix (m, d);	/* make a copy of A matrix for output later            */
+
+  for (i = 1; i <= m; i++)
+    {
+      for (j = 0; j <= d; j++)
+	copy (Ain[i][j], P->A[i][j]);
+
+      if (debug)
+	lrs_printrow ("*", Q, Ain[i], d);
+    }
+
+/*********************************************************************************/
+/* Step 2: Find a starting cobasis from default of specified order               */
+/*         Lin is created if necessary to hold linearity space                   */
+/*********************************************************************************/
+
+  if (!lrs_getfirstbasis (&P, Q, &Lin, TRUE))
+    return 1;
+
+  /* Pivot to a starting dictionary                      */
+  /* There may have been column redundancy               */
+  /* If so the linearity space is obtained and redundant */
+  /* columns are removed. User can access linearity space */
+  /* from lrs_mp_matrix Lin dimensions nredundcol x d+1  */
+
+
+/*********************************************************************************/
+/* Step 3: Test each row of the dictionary to see if it is redundant             */
+/*********************************************************************************/
+
+/* note some of these may have been changed in getting initial dictionary        */
+  m = P->m_A;
+  d = P->d;
+  nlinearity = Q->nlinearity;
+  lastdv = Q->lastdv;
+      if (debug)
+	fprintf (lrs_ofp, "\ncheckindex m=%ld, n=%ld, nlinearity=%ld lastdv=%ld", m,d,nlinearity,lastdv);
+
+/* linearities are not considered for redundancy */
+
+  for (i = 0; i < nlinearity; i++)
+    redineq[Q->linearity[i]] = 2L;
+
+/* rows 0..lastdv are cost, decision variables, or linearities  */
+/* other rows need to be tested                                */
+
+  for (index = lastdv + 1; index <= m + d; index++)
+    {
+      ineq = Q->inequality[index - lastdv];	/* the input inequality number corr. to this index */
+
+      redineq[ineq] = checkindex (P, Q, index);
+      if (debug)
+	fprintf (lrs_ofp, "\ncheck index=%ld, inequality=%ld, redineq=%ld", index, ineq, redineq[ineq]);
+      if (redineq[ineq] == ONE)
+        {
+	fprintf (lrs_ofp, "\n*row %ld was redundant and removed", ineq);
+        fflush  (lrs_ofp);
+        }
+
+    }				/* end for index ..... */
+
+  if (debug)
+    {
+      fprintf (lrs_ofp, "\n*redineq:");
+      for (i = 1; i <= m; i++)
+	fprintf (lrs_ofp, " %ld", redineq[i]);
+    }
+
+  if (!Q->hull)
+    fprintf (lrs_ofp, "\nH-representation");
+  else
+    fprintf (lrs_ofp, "\nV-representation");
+
+/* linearities will be printed first in output */
+
+  if (nlinearity > 0)
+    {
+      fprintf (lrs_ofp, "\nlinearity %ld", nlinearity);
+      for (i = 1; i <= nlinearity; i++)
+	fprintf (lrs_ofp, " %ld", i);
+
+    }
+  nredund = nlinearity;		/* count number of non-redundant inequalities */
+  for (i = 1; i <= m; i++)
+    if (redineq[i] == 0)
+      nredund++;
+  fprintf (lrs_ofp, "\nbegin");
+  fprintf (lrs_ofp, "\n%ld %ld rational", nredund, Q->n);
+
+/* print the linearities first */
+
+  for (i = 0; i < nlinearity; i++)
+    lrs_printrow ("", Q, Ain[Q->linearity[i]], Q->inputd);
+
+  for (i = 1; i <= m; i++)
+    if (redineq[i] == 0)
+      lrs_printrow ("", Q, Ain[i], Q->inputd);
+  fprintf (lrs_ofp, "\nend");
+  fprintf (lrs_ofp, "\n*Input had %ld rows and %ld columns", m, Q->n);
+  fprintf (lrs_ofp, ": %ld row(s) redundant", m - nredund);
+
+/* 2015.9.9  fix memory leak on Gcd Lcm */
+  long savem=P->m;              /* need this to clear Q*/
+  lrs_free_dic (P,Q);           /* deallocate lrs_dic */
+  Q->m=savem;
+
+  lrs_free_dat (Q);             /* deallocate lrs_dat */
+
+
+  lrs_close ("redund:");
+
+  return 0;
+}
+/*********************************************/
+/* end of redund.c                           */
+/*********************************************/
+/*******************/
+/* lrs_printoutput */
+/*******************/
+void 
+lrs_printoutput (lrs_dat * Q, lrs_mp_vector output)
+{
+
+if (Q->countonly)
+   return;
+
+#ifdef PLRS
+	//Make new output node
+	char *type=NULL;
+	
+	//Make stream to collect prat / pmp data
+	stringstream ss;
+
+
+	if (Q->hull || zero (output[0])){
+		/*non vertex */
+		type = "ray";
+		for (int i = 0; i < Q->n; i++)
+			ss<<pmp ("", output[i]);
+	}else{
+		type = "vertex";
+		/* vertex   */
+		ss<<" 1 ";
+		for (int i = 1; i < Q->n; i++)
+			ss<<prat ("", output[i], output[0]);
+	}
+	//post output in a nonblocking manner (a consumer thread will manage output)
+	post_output(type, ss.str().c_str());
+#else
+  long i;
+
+  fprintf (lrs_ofp, "\n");
+
+  if (Q->hull || zero (output[0]))	/*non vertex */
+    {
+      for (i = 0; i < Q->n; i++)
+	pmp ("", output[i]);
+
+    }
+  else
+    {				/* vertex   */
+      fprintf (lrs_ofp, " 1 ");
+      for (i = 1; i < Q->n; i++)
+	prat ("", output[i], output[0]);
+    }
+  fflush(lrs_ofp);
+#endif
+
+}
+/**************************/
+/* end of lrs_printoutput */
+/**************************/
+
+/****************/
+/* lrs_lpoutput */
+/****************/
+void lrs_lpoutput(lrs_dic * P,lrs_dat * Q, lrs_mp_vector output)
+{
+	
+
+#ifndef LRS_QUIET
+  lrs_mp Temp1, Temp2;
+  long i;
+
+  lrs_alloc_mp (Temp1);
+  lrs_alloc_mp (Temp2);
+
+  fprintf (lrs_ofp, "\n*LP solution only requested");
+  prat ("\n\n*Objective function has value ", P->objnum, P->objden);
+
+  fprintf (lrs_ofp, "\n\n*Primal: ");
+  for (i = 1; i < Q->n; i++)
+      {
+        fprintf(lrs_ofp,"x_%ld=",i);
+        prat ("", output[i], output[0]);
+       }
+
+  if(Q->nlinearity > 0)
+      fprintf (lrs_ofp, "\n\n*Linearities in input file - partial dual solution only");
+  fprintf (lrs_ofp, "\n\n*Dual: ");
+
+  for (i = 0; i < P->d; i++)
+	  {
+	        fprintf(lrs_ofp,"y_%ld=",Q->inequality[P->C[i]-Q->lastdv]);
+	        changesign(P->A[0][P->Col[i]]);
+                mulint(Q->Lcm[P->Col[i]],P->A[0][P->Col[i]],Temp1);
+                mulint(Q->Gcd[P->Col[i]],P->det,Temp2);
+	        prat("",Temp1,Temp2);
+	        changesign(P->A[0][P->Col[i]]);
+          }
+  fprintf (lrs_ofp, "\n");
+  lrs_clear_mp (Temp1);
+  lrs_clear_mp (Temp2);
+#endif
+ }
+/***********************/
+/* end of lrs_lpoutput */
+/***********************/
+void 
+lrs_printrow (char name[], lrs_dat * Q, lrs_mp_vector output, long rowd)
+/* print a row of A matrix in output in "original" form  */
+/* rowd+1 is the dimension of output vector                */
+/* if input is H-rep. output[0] contains the RHS      */
+/* if input is V-rep. vertices are scaled by 1/output[1] */
+{
+  long i;
+
+  fprintf (lrs_ofp, "\n%s", name);
+  if (!Q->hull)			/* input was inequalities, print directly */
+
+    {
+
+      for (i = 0; i <= rowd; i++)
+	pmp ("", output[i]);
+      return;
+    }
+
+/* input was vertex/ray */
+
+  if (zero (output[1]))		/*non-vertex */
+    {
+      for (i = 1; i <= rowd; i++)
+	pmp ("", output[i]);
+
+    }
+  else
+    {				/* vertex */
+      fprintf (lrs_ofp, " 1 ");
+      for (i = 2; i <= rowd; i++)
+	prat ("", output[i], output[1]);
+    }
+
+  return;
+
+}				/* end of lrs_printrow */
+
+long 
+lrs_getsolution (lrs_dic * P, lrs_dat * Q, lrs_mp_vector output, long col)
+   /* check if column indexed by col in this dictionary */
+   /* contains output                                   */
+   /* col=0 for vertex 1....d for ray/facet             */
+{
+
+	
+  long j;			/* cobasic index     */
+
+  lrs_mp_matrix A = P->A;
+  long *Row = P->Row;
+
+  if (col == ZERO)		/* check for lexmin vertex */
+    return lrs_getvertex (P, Q, output);
+
+/*  check for rays: negative in row 0 , positive if lponly */
+
+  if (Q->lponly)
+    {
+      if (!positive (A[0][col]))
+	return FALSE;
+    }
+
+  else if (!negative (A[0][col]))
+    return FALSE;
+
+/*  and non-negative for all basic non decision variables */
+
+  j = Q->lastdv + 1;
+  while (j <= P->m && !negative (A[Row[j]][col]))
+    j++;
+
+  if (j <= P->m)
+    return FALSE;
+
+  if (Q->geometric || Q->allbases || lexmin (P, Q, col) || Q->lponly)
+
+    return lrs_getray (P, Q, col, Q->n, output);
+
+  return FALSE;			/* no more output in this dictionary */
+
+}				/* end of lrs_getsolution */
+
+
+long
+lrs_init (char *name)       /* returns TRUE if successful, else FALSE */
+{
+
+  printf ("%s", name);
+  printf (TITLE);
+  printf (VERSION);
+  printf ("(");
+  printf (BIT); 
+  printf (","); 
+  printf (ARITH);
+  if (!lrs_mp_init (ZERO, stdin, stdout))  /* initialize arithmetic */
+    return FALSE;
+  printf (")");
+
+
+  lrs_global_count = 0;
+  lrs_checkpoint_seconds = 0;
+#ifdef SIGNALS
+  setup_signals ();
+#endif
+  return TRUE;
+}
+
+void 
+lrs_close (char *name)
+{
+
+  fprintf (lrs_ofp, "\n*%s", name);
+  fprintf (lrs_ofp, TITLE);
+  fprintf (lrs_ofp, VERSION);
+  fprintf (lrs_ofp, "(");
+  fprintf (lrs_ofp, BIT);
+  fprintf (lrs_ofp, ",");
+  fprintf (lrs_ofp, ARITH);
+  fprintf (lrs_ofp, ")");
+
+#ifdef MP   
+  fprintf (lrs_ofp, " max digits=%ld/%ld", DIG2DEC (lrs_record_digits), DIG2DEC (lrs_digits));
+#endif
+
+#ifdef TIMES
+  ptimes ();
+#endif
+
+  fprintf (lrs_ofp, "\n");
+  fclose (lrs_ifp);
+  if (lrs_ofp != stdout)
+    fclose (lrs_ofp);
+}
+
+/***********************************/
+/* allocate and initialize lrs_dat */
+/***********************************/
+lrs_dat *
+lrs_alloc_dat (const char *name)
+{
+  lrs_dat *Q;
+  long i;
+
+
+  if (lrs_global_count >= MAX_LRS_GLOBALS)
+    {
+      fprintf (stderr,
+	       "Fatal: Attempt to allocate more than %ld global data blocks\n", MAX_LRS_GLOBALS);
+      exit (1);
+
+    }
+
+  Q = (lrs_dat *) malloc (sizeof (lrs_dat));
+  if (Q == NULL)
+    return Q;			/* failure to allocate */
+
+  lrs_global_list[lrs_global_count] = Q;
+  Q->id = lrs_global_count;
+  lrs_global_count++;
+  Q->name=(char *) CALLOC ((unsigned) strlen(name)+1, sizeof (char));
+  strcpy(Q->name,name); 
+
+/* initialize variables */
+  Q->m = 0L;
+  Q->n = 0L;
+  Q->inputd = 0L;
+  Q->deepest = 0L;
+  Q->nlinearity = 0L;
+  Q->nredundcol = 0L;
+  Q->runs = 0L;
+  Q->subtreesize=MAXD;
+  Q->seed = 1234L;
+  Q->totalnodes = 0L;
+  for (i = 0; i < 10; i++)
+    {
+      Q->count[i] = 0L;
+      Q->cest[i] = 0.0;
+      if(i < 5)
+	Q->startcount[i] = 0L;
+    }
+  Q->count[2] = 1L;		/* basis counter */
+  Q->startcount[2] = 0L;	/* starting basis counter */
+/* initialize flags */
+  Q->allbases = FALSE;
+  Q->bound = FALSE;            /* upper/lower bound on objective function given */
+  Q->countonly = FALSE;        /* produce the usual output */
+  Q->debug = FALSE;
+  Q->frequency = 0L;
+  Q->dualdeg = FALSE;          /* TRUE if dual degenerate starting dictionary */
+  Q->geometric = FALSE;
+  Q->getvolume = FALSE;
+  Q->homogeneous = TRUE;
+  Q->polytope = FALSE;
+  Q->hull = FALSE;
+  Q->incidence = FALSE;
+  Q->lponly = FALSE;
+  Q->maxdepth = MAXD;
+  Q->mindepth = -MAXD;
+  Q->maxoutput = 0L;
+  Q->maxcobases = 0L;     /* after maxcobases have been found unexplored subtrees reported */
+  Q->nash = FALSE;
+  Q->nonnegative = FALSE;
+  Q->printcobasis = FALSE;
+  Q->printslack = FALSE;
+  Q->truncate = FALSE;          /* truncate tree when moving from opt vertex        */
+  Q->verbose=FALSE;
+  Q->voronoi = FALSE;
+  Q->maximize = FALSE;		/*flag for LP maximization                          */
+  Q->minimize = FALSE;		/*flag for LP minimization                          */
+  Q->restart = FALSE;		/* TRUE if restarting from some cobasis             */
+  Q->givenstart = FALSE;	/* TRUE if a starting cobasis is given              */
+  Q->strace = -1L;		/* turn on  debug at basis # strace */
+  Q->etrace = -1L;		/* turn off debug at basis # etrace */
+
+  Q->saved_flag = 0;		/* no cobasis saved initially, db */
+  lrs_alloc_mp (Q->Nvolume);
+  lrs_alloc_mp (Q->Dvolume);
+  lrs_alloc_mp (Q->sumdet);
+  lrs_alloc_mp (Q->saved_det);
+  lrs_alloc_mp (Q->boundn);
+  lrs_alloc_mp (Q->boundd);
+  itomp (ZERO, Q->Nvolume);
+  itomp (ONE, Q->Dvolume);
+  itomp (ZERO, Q->sumdet);
+/* 2012.6.1 */
+  Q->unbounded = FALSE;
+
+  return Q;
+}				/* end of allocate and initialize lrs_dat */
+
+/*******************************/
+/*  lrs_read_dat               */
+/*******************************/
+long 
+lrs_read_dat (lrs_dat * Q, int argc, char *argv[])
+{
+  char name[100];
+  long dec_digits = 0;
+  long infile=0;                /*input file number to open if any        */
+  long firstline = TRUE;	/*flag for picking off name at line 1     */
+
+  int c;			/* for fgetc */
+
+
+  if(argc > 1 )
+	  infile=1;
+  if(Q->nash && argc == 2)        /* open second nash input file */
+	  infile=2;
+
+  if (infile > 0)			/* command line argument overides stdin   */
+    {
+      if ((lrs_ifp = fopen (argv[infile], "r")) == NULL)
+	{
+	  printf ("\nBad input file name\n");
+	  return (FALSE);
+	}
+      else
+      { if (infile==1)
+	printf ("\n*Input taken from file %s", argv[infile]);
+      }
+    }
+
+       	  /* command line argument overides stdout   */
+  if ((!Q->nash && argc == 3) || (Q->nash && argc == 4)) 
+    {
+      if ((lrs_ofp = fopen (argv[argc-1], "w")) == NULL)
+	{
+	  printf ("\nBad output file name\n");
+	  return (FALSE);
+	}
+      else
+	printf ("\n*Output sent to file %s\n", argv[argc-1]);
+    }
+
+
+/* process input file */
+  if( fscanf (lrs_ifp, "%s", name) == EOF)
+	    {
+	      fprintf (lrs_ofp, "\nNo begin line");
+	      return (FALSE);
+	    }
+
+  while (strcmp (name, "begin") != 0)	/*skip until "begin" found processing options */
+    {
+      if (strncmp (name, "*", 1) == 0)	/* skip any line beginning with * */
+	{
+	  c = name[0];
+	  while (c != EOF && c != '\n')
+	    c = fgetc (lrs_ifp);
+	}
+
+      else if (strcmp (name, "H-representation") == 0)
+	Q->hull = FALSE;
+      else if ((strcmp (name, "hull") == 0) || (strcmp (name, "V-representation") == 0))
+        {
+	   Q->hull = TRUE;
+           Q->polytope = TRUE;		/* will be updated as input read */
+        }
+      else if (strcmp (name, "digits") == 0)
+	{
+	  if (fscanf (lrs_ifp, "%ld", &dec_digits) == EOF)
+	    {
+	      fprintf (lrs_ofp, "\nNo begin line");
+	      return (FALSE);
+	    }
+          if  (!lrs_set_digits(dec_digits))
+             return (FALSE);
+	}
+      else if (strcmp (name, "linearity") == 0)
+	{
+	  if (!readlinearity (Q))
+	    return FALSE;
+	}
+      else if (strcmp (name, "nonnegative") == 0)
+	{
+         if(Q->nash)
+	  fprintf (lrs_ofp, "\nNash incompatibile with nonnegative option - skipped");
+         else
+	  Q->nonnegative = TRUE;
+	}
+      else if (firstline)
+	{
+	  stringcpy (Q->fname, name);
+	  fprintf (lrs_ofp, "\n%s", Q->fname);
+	  firstline = FALSE;
+	}
+
+      if (fscanf (lrs_ifp, "%s", name) == EOF)
+	{
+	  fprintf (lrs_ofp, "\nNo begin line");
+	  return (FALSE);
+	}
+
+    }				/* end of while */
+
+
+  if (fscanf (lrs_ifp, "%ld %ld %s", &Q->m, &Q->n, name) == EOF)
+    {
+      fprintf (lrs_ofp, "\nNo data in file");
+      return (FALSE);
+    }
+  if (strcmp (name, "integer") != 0 && strcmp (name, "rational") != 0)
+    {
+      fprintf (lrs_ofp, "\nData type must be integer of rational");
+      return (FALSE);
+    }
+
+
+  if (Q->m == 0)
+    {
+      fprintf (lrs_ofp, "\nNo input given");	/* program dies ungracefully */
+      return (FALSE);
+    }
+
+
+
+  /* inputd may be reduced in preprocessing of linearities and redund cols */
+
+  return TRUE;
+}				/* end of lrs_read_dat */
+
+/****************************/
+/* set up lrs_dic structure */
+/****************************/
+long 
+lrs_read_dic (lrs_dic * P, lrs_dat * Q)
+/* read constraint matrix and set up problem and dictionary  */
+
+{
+  lrs_mp Temp,Tempn,Tempd, mpone, mpten;
+  lrs_mp_vector oD;		/* Denom for objective function */
+
+  long i, j;
+  char name[100];
+  int c; /* fgetc actually returns an int. db */
+
+/* assign local variables to structures */
+
+  lrs_mp_matrix A;
+  lrs_mp_vector Gcd, Lcm;
+  long hull = Q->hull;
+  long m, d;
+  long dualperturb=FALSE;    /* dualperturb=TRUE: objective function perturbed */
+
+  lrs_alloc_mp(Temp); lrs_alloc_mp(mpone);
+  lrs_alloc_mp(Tempn); lrs_alloc_mp(Tempd); lrs_alloc_mp(mpten);
+  A = P->A;
+  m = Q->m;
+  d = Q->inputd;
+  Gcd = Q->Gcd;
+  Lcm = Q->Lcm;
+
+  oD = lrs_alloc_mp_vector (d);
+
+
+
+
+  itomp (ONE, mpone);
+  itomp(10L,mpten);
+  itomp (ONE, A[0][0]);
+  itomp (ONE, Lcm[0]);
+  itomp (ONE, Gcd[0]);
+
+  for (i = 1; i <= m; i++)	/* read in input matrix row by row                 */
+    {
+      itomp (ONE, Lcm[i]);	/* Lcm of denominators */
+      itomp (ZERO, Gcd[i]);	/* Gcd of numerators */
+      for (j = hull; j <= d; j++)	/* hull data copied to cols 1..d */
+	{
+	  if (readrat (A[i][j], A[0][j]))
+	    lcm (Lcm[i], A[0][j]);	/* update lcm of denominators */
+	  copy (Temp, A[i][j]);
+	  gcd (Gcd[i], Temp);	/* update gcd of numerators   */
+	}
+
+      if (hull)
+	{
+	  itomp (ZERO, A[i][0]);	/*for hull, we have to append an extra column of zeroes */
+	  if (!one (A[i][1]) || !one (A[0][1]))		/* all rows must have a one in column one */
+	    Q->polytope = FALSE;
+	}
+      if (!zero (A[i][hull]))	/* for H-rep, are zero in column 0     */
+	Q->homogeneous = FALSE;	/* for V-rep, all zero in column 1     */
+
+
+      storesign (Gcd[i], POS);
+      storesign (Lcm[i], POS);
+      if (mp_greater (Gcd[i], mpone) || mp_greater (Lcm[i], mpone))
+	for (j = 0; j <= d; j++)
+	  {
+	    exactdivint (A[i][j], Gcd[i], Temp);	/*reduce numerators by Gcd  */
+	    mulint (Lcm[i], Temp, Temp);	/*remove denominators */
+	    exactdivint (Temp, A[0][j], A[i][j]);	/*reduce by former denominator */
+	  }
+
+    }				/* end of for i=       */
+
+/* 2010.4.26 patch */
+  if(Q->nonnegative)    /* set up Gcd and Lcm for nonexistent nongative inequalities */
+      for (i=m+1;i<=m+d;i++)
+          { itomp (ONE, Lcm[i]);
+            itomp (ONE, Gcd[i]);
+          }
+  
+  if (Q->homogeneous && Q->verbose)
+    {
+      fprintf (lrs_ofp, "\n*Input is homogeneous, column 1 not treated as redundant");
+    }
+
+
+/* read in flags */
+  while (fscanf (lrs_ifp, "%s", name) != EOF)
+    {
+      if (strncmp (name, "*", 1) == 0)	/* skip any line beginning with * */
+	{
+	  c = name[0];
+	  while (c != EOF && c != '\n')
+	    c = fgetc (lrs_ifp);
+	}
+
+
+      if (strcmp (name, "checkpoint") == 0)
+	{
+	  long seconds;
+
+	  if(fscanf (lrs_ifp, "%ld", &seconds) == EOF)
+            {
+              fprintf (lrs_ofp, "\nInvalid checkpoint option");
+              return (FALSE);
+            }
+
+#ifdef SIGNALS
+	  if (seconds > 0)
+	    {
+	      lrs_checkpoint_seconds = seconds;
+	      errcheck ("signal", signal (SIGALRM, timecheck));
+	      alarm (lrs_checkpoint_seconds);
+	    }
+#endif
+	}
+
+      if (strcmp (name, "debug") == 0)
+	{
+          Q->etrace =0;
+	  if(fscanf (lrs_ifp, "%ld %ld", &Q->strace, &Q->etrace)==EOF)
+             Q->strace =0;
+	  fprintf (lrs_ofp, "\n*%s from B#%ld to B#%ld", name, Q->strace, Q->etrace);
+          Q->verbose=TRUE;
+	  if (Q->strace <= 1)
+	    Q->debug = TRUE;
+	}
+      if (strcmp (name, "startingcobasis") == 0)
+	{
+          if(Q->nonnegative)
+	      fprintf (lrs_ofp, "\n*startingcobasis incompatible with nonnegative option:skipped");
+          else
+            {    
+	      fprintf (lrs_ofp, "\n*startingcobasis");
+	      Q->givenstart = TRUE;
+	      if (!readfacets (Q, Q->inequality))
+	          return FALSE;
+	    }
+        }
+
+      if (strcmp (name, "restart") == 0)
+	{
+	  Q->restart = TRUE;
+          if(Q->voronoi)
+           {
+             if(fscanf (lrs_ifp, "%ld %ld %ld %ld", &Q->count[1], &Q->count[0], &Q->count[2], &P->depth)==EOF)
+               return FALSE;
+             fprintf (lrs_ofp, "\n*%s V#%ld R#%ld B#%ld h=%ld data points", name, Q->count[1], Q->count[0], Q->count[2], P->depth);
+            }
+          else if(hull)
+            {
+	    if( fscanf (lrs_ifp, "%ld %ld %ld", &Q->count[0], &Q->count[2], &P->depth)==EOF)
+	     fprintf (lrs_ofp, "\n*%s F#%ld B#%ld h=%ld vertices/rays", name, Q->count[0], Q->count[2], P->depth);
+            }
+          else
+            {
+	     if(fscanf (lrs_ifp, "%ld %ld %ld %ld", &Q->count[1], &Q->count[0], &Q->count[2], &P->depth)==EOF)
+               return FALSE;
+	     fprintf (lrs_ofp, "\n*%s V#%ld R#%ld B#%ld h=%ld facets", name, Q->count[1], Q->count[0], Q->count[2], P->depth);
+            }
+	  if (!readfacets (Q, Q->facet))
+	    return FALSE;
+	}			/* end of restart */
+
+/* The next flag request a LP solution only */
+      if (strcmp (name, "lponly") == 0)
+	{
+	  if (Q->hull)
+	    fprintf (lrs_ofp, "\n*lponly  option not valid for V-representation-skipped");
+	  else
+	    Q->lponly = TRUE;
+	}
+
+
+/* The LP will be solved after initialization to get starting vertex   */
+/* Used also with lponly flag                                          */
+      if (strcmp (name, "maximize") == 0 || strcmp (name, "minimize") == 0)
+	{
+	  if (Q->hull)
+	    fprintf (lrs_ofp, "\n*%s option not valid for V-representation-skipped", name);
+	  else
+	    {
+	      {
+		if (strcmp (name, "maximize") == 0)
+		  Q->maximize = TRUE;
+		else
+		  Q->minimize = TRUE;
+	      }
+	      fprintf (lrs_ofp,"\n*%s", name);
+
+              if(dualperturb)   /* apply a perturbation to objective function */
+                {
+	          fprintf (lrs_ofp, " - Objective function perturbed");
+                  copy(Temp,mpten);
+                  for (j = 0; j <= 10; j++)
+                      mulint(mpten,Temp,Temp);
+                }
+
+              fprintf (lrs_ofp, ":  ");
+
+	      for (j = 0; j <= d; j++)
+		{
+		  if (readrat (A[0][j], oD[j]) || dualperturb )
+		    {
+                      if(dualperturb && j > 0 && j < d )
+                        {
+                         if (Q->maximize)
+                          linrat(A[0][j], oD[j],ONE,mpone,Temp,ONE,Tempn,Tempd);
+                         else
+                          linrat(A[0][j], oD[j],ONE,mpone,Temp,-1L,Tempn,Tempd);
+
+                         copy(A[0][j],Tempn);
+                         copy(oD[j],Tempd);
+                         mulint(mpten,Temp,Temp);
+                        }
+
+		      reduce (A[0][j], oD[j]);
+		      lcm (Q->Lcm[0], oD[j]);	/* update lcm of denominators */
+		    }
+		  prat ("", A[0][j], oD[j]);
+		  if (!Q->maximize)
+		    changesign (A[0][j]);
+		}
+	      storesign (Q->Lcm[0], POS);
+	      if (mp_greater (Q->Lcm[0], mpone))
+		for (j = 0; j <= d; j++)
+		  {
+		    mulint (Q->Lcm[0], A[0][j], A[0][j]);	/*remove denominators */
+		    copy (Temp, A[0][j]);
+		    exactdivint (Temp, oD[j], A[0][j]);
+		  }
+	      if (Q->debug)
+		printA (P, Q);
+	    }
+	}			/* end of LP setup */
+      if (strcmp (name, "volume") == 0)
+	{
+	  fprintf (lrs_ofp, "\n*%s", name);
+	  Q->getvolume = TRUE;
+	}
+      if (strcmp (name, "geometric") == 0)
+	{
+	  fprintf (lrs_ofp, "\n*%s", name);
+	  if (hull & !Q->voronoi)
+	    fprintf (lrs_ofp, " - option for H-representation or voronoi only, skipped");
+	  else
+	    Q->geometric = TRUE;
+	}
+      if (strcmp (name, "allbases") == 0)
+	{
+	  fprintf (lrs_ofp, "\n*%s", name);
+	  Q->allbases = TRUE;
+        }
+
+      if (strcmp (name, "countonly") == 0)
+	{
+	  fprintf (lrs_ofp, "\n*%s", name);
+	  Q->countonly = TRUE;
+	}
+      if (strcmp (name, "dualperturb") == 0)
+	{
+	  dualperturb = TRUE;
+	}
+
+      if (strcmp (name, "incidence") == 0)
+	{
+	  fprintf (lrs_ofp, "\n*%s", name);
+	  Q->incidence = TRUE;
+	}
+
+      if (strcmp (name, "#incidence") == 0) /* number of incident inequalities only */
+	{
+	  Q->printcobasis = TRUE;
+	}
+
+      if (strcmp (name, "printcobasis") == 0)
+	{
+	  if(fscanf (lrs_ifp, "%ld", &Q->frequency)==EOF)
+/*2010.7.7  set default to zero = print only when outputting vertex/ray/facet */
+             Q->frequency=0;
+	  fprintf (lrs_ofp, "\n*%s", name);
+          if (Q->frequency > 0)
+            fprintf(lrs_ofp," %ld", Q->frequency);
+	  Q->printcobasis = TRUE;
+	}
+
+      if (strcmp (name, "printslack") == 0)
+	{
+	  Q->printslack = TRUE;
+	}
+
+      if (strcmp (name, "cache") == 0)
+	{
+	  if(fscanf (lrs_ifp, "%ld", &dict_limit)==EOF)
+              dict_limit=1;
+	  fprintf (lrs_ofp, "\n*cache %ld", dict_limit);
+	  if (dict_limit < 1)
+	    dict_limit = 1;
+	}
+      if (strcmp (name, "linearity") == 0)
+	{
+	  if (!readlinearity (Q))
+	    return FALSE;
+	}
+
+      if (strcmp (name, "maxdepth") == 0)
+	{
+	  if(fscanf (lrs_ifp, "%ld", &Q->maxdepth)==EOF)
+                    Q->maxdepth=MAXD;
+	  fprintf (lrs_ofp, "\n*%s  %ld", name, Q->maxdepth);
+	}
+
+      if (strcmp (name, "maxoutput") == 0)
+	{
+	  if(fscanf (lrs_ifp, "%ld", &Q->maxoutput)==EOF)
+             Q->maxoutput = 100;
+	  fprintf (lrs_ofp, "\n*%s  %ld", name, Q->maxoutput);
+	}
+
+      if (strcmp (name, "maxcobases") == 0)
+	{
+	  if(fscanf (lrs_ifp, "%ld", &Q->maxcobases)==EOF)
+             Q->maxcobases = 1000;
+	  fprintf (lrs_ofp, "\n*%s  %ld", name, Q->maxcobases);
+	}
+
+      if (strcmp (name, "mindepth") == 0)
+	{
+	if( fscanf (lrs_ifp, "%ld", &Q->mindepth)==EOF)
+           Q->mindepth = 0;
+	  fprintf (lrs_ofp, "\n*%s  %ld", name, Q->mindepth);
+	}
+
+      if (strcmp (name, "truncate") == 0)
+        {
+          fprintf (lrs_ofp, "\n*%s", name);
+          if (!hull)
+            Q->truncate = TRUE;
+          else
+            fprintf (lrs_ofp, " - option for H-representation only, skipped");
+        }
+
+
+      if (strcmp (name, "verbose") == 0)
+          Q->verbose = TRUE;
+
+      if (strcmp (name, "bound") == 0)
+         {
+          readrat(Q->boundn,Q->boundd);
+          Q->bound = TRUE;
+         }
+
+      if (strcmp (name, "nonnegative") == 0)
+	{
+	  fprintf (lrs_ofp, "\n*%s", name);
+	  fprintf (lrs_ofp, " - option must come before begin line - skipped");
+	}
+
+      if (strcmp (name, "seed") == 0)
+	{
+	  if(fscanf (lrs_ifp, "%ld", &Q->seed)==EOF)
+               Q->seed = 3142;
+	  fprintf (lrs_ofp, "\n*seed= %ld ", Q->seed);
+	}
+
+      if (strcmp (name, "estimates") == 0)
+	{
+	  if(fscanf (lrs_ifp, "%ld", &Q->runs)==EOF)
+             Q->runs=1;
+	  fprintf (lrs_ofp, "\n*%ld %s", Q->runs, name);
+	}
+
+// 2015.2.9   Estimates will continue until estimate is less than subtree size
+      if (strcmp (name, "subtreesize") == 0)
+        {
+          if(fscanf (lrs_ifp, "%ld", &Q->subtreesize)==EOF)
+             Q->subtreesize=MAXD;
+          fprintf (lrs_ofp, "\n*%s %ld", name, Q->subtreesize);
+        }
+
+
+
+      if ((strcmp (name, "voronoi") == 0) || (strcmp (name, "Voronoi") == 0))
+	{
+	  if (!hull)
+	    fprintf (lrs_ofp, "\n*voronoi requires V-representation - option skipped");
+	  else
+	    {
+	      Q->voronoi = TRUE;
+	      Q->polytope = FALSE;
+	    }
+	}
+
+    }				/* end of while for reading flags */
+
+  if (Q->polytope)
+    Q->getvolume = TRUE; 	/* might as well get volume, it doesn't cost much */
+
+  if (Q->bound && Q->maximize)
+    prat("\n*Lower bound on objective function:",Q->boundn,Q->boundd);
+
+  if (Q->bound && Q->minimize)
+    prat("\n*Upper bound on objective function:",Q->boundn,Q->boundd);
+
+/* Certain options are incompatible, this is fixed here */
+
+  if (Q->restart)
+    Q->getvolume = FALSE;       /* otherwise incorrect volume reported            */
+
+    if (Q->restart && Q->maxcobases > 0) //2015.4.3 adjust for restart
+               Q->maxcobases = Q->maxcobases + Q->count[2];
+
+  if (Q->incidence)
+    {
+      Q->printcobasis = TRUE;
+/* 2010.5.7    No need to reset this, as it may have been set by printcobasis */
+/*    Q->frequency    = ZERO;                     */                                                      
+    }
+
+  if (Q->debug)
+    {
+      printA (P, Q);
+      fprintf (lrs_ofp, "\nexiting lrs_read_dic");
+    }
+  lrs_clear_mp(Temp); lrs_clear_mp(mpone);
+  lrs_clear_mp(Tempn); lrs_clear_mp(Tempd); lrs_clear_mp(mpten);
+  lrs_clear_mp_vector (oD,d);
+  return TRUE;
+
+}
+
+ 		/* end of if(voronoi)     */
+
+/* In lrs_getfirstbasis and lrs_getnextbasis we use D instead of P */
+/* since the dictionary P may change, ie. &P in calling routine    */
+
+#define D (*D_p)
+
+long 
+lrs_getfirstbasis (lrs_dic ** D_p, lrs_dat * Q, lrs_mp_matrix * Lin, long no_output)
+/* gets first basis, FALSE if none              */
+/* P may get changed if lin. space Lin found    */
+/* no_output is TRUE supresses output headers   */
+{
+  lrs_mp scale, Temp;
+
+  long i, j, k;
+
+/* assign local variables to structures */
+
+  lrs_mp_matrix A;
+  long *B, *C, *Col;
+  long *inequality;
+  long *linearity;
+  long hull = Q->hull;
+  long m, d, lastdv, nlinearity, nredundcol;
+
+  lrs_alloc_mp(Temp); lrs_alloc_mp(scale);
+
+  if (Q->lponly)
+    no_output = TRUE;
+  m = D->m;
+  d = D->d;
+
+	
+
+  lastdv = Q->lastdv;
+
+  nredundcol = 0L;		/* will be set after getabasis        */
+  nlinearity = Q->nlinearity;	/* may be reset if new linearity read or in getabasis*/
+  linearity = Q->linearity;
+
+	
+
+  A = D->A;
+  B = D->B;
+  C = D->C;
+  Col = D->Col;
+  inequality = Q->inequality;
+
+
+  if (Q->nlinearity > 0 && Q->nonnegative)
+   {
+	    fprintf (lrs_ofp, "\n*linearity and nonnegative options incompatible");
+	    fprintf (lrs_ofp, " - all linearities are skipped");
+	    fprintf (lrs_ofp, "\n*add nonnegative constraints explicitly and ");
+	    fprintf (lrs_ofp, " remove nonnegative option");
+   }
+
+  if (Q->nlinearity && Q->voronoi){
+    	fprintf (lrs_ofp, "\n*linearity and Voronoi options set - results unpredictable");
+  }
+
+  if (Q->lponly && !Q->maximize && !Q->minimize)
+    	    fprintf (lrs_ofp, "\n*LP has no objective function given - assuming all zero");
+
+
+  if (Q->runs > 0)		/* arrays for estimator */
+    {
+      Q->isave = (long *) CALLOC ((unsigned) (m * d), sizeof (long));
+      Q->jsave = (long *) CALLOC ((unsigned) (m * d), sizeof (long));
+    }
+/* default is to look for starting cobasis using linearies first, then     */
+/* filling in from last rows of input as necessary                         */
+/* linearity array is assumed sorted here                                  */
+/* note if restart/given start inequality indices already in place         */
+/* from nlinearity..d-1                                                    */
+  for (i = 0; i < nlinearity; i++)	/* put linearities first in the order */
+    inequality[i] = linearity[i];
+
+  k = 0;			/* index for linearity array   */
+
+  if (Q->givenstart)
+    k = d;
+  else
+    k = nlinearity;
+  for (i = m; i >= 1; i--)
+    {
+      j = 0;
+      while (j < k && inequality[j] != i)
+	j++;			/* see if i is in inequality  */
+      if (j == k)
+	inequality[k++] = i;
+    }
+	#ifndef PLRS
+  if (Q->debug)
+    {
+	      fprintf (lrs_ofp, "\n*Starting cobasis uses input row order");
+	      for (i = 0; i < m; i++)
+		fprintf (lrs_ofp, " %ld", inequality[i]);
+	
+    }
+	#endif
+/* for voronoi convert to h-description using the transform                  */
+/* a_0 .. a_d-1 -> (a_0^2 + ... a_d-1 ^2)-2a_0x_0-...-2a_d-1x_d-1 + x_d >= 0 */
+/* note constant term is stored in column d, and column d-1 is all ones      */
+/* the other coefficients are multiplied by -2 and shifted one to the right  */
+  if (Q->debug)
+    printA (D, Q);
+  if (Q->voronoi)
+    {
+      Q->hull = FALSE;
+      hull = FALSE;
+      for (i = 1; i <= m; i++)
+	{
+	  if (zero (A[i][1]))
+	    {
+              printf("\nWith voronoi option column one must be all one\n");
+	      return (FALSE);
+	    }
+	  copy (scale, A[i][1]);	/*adjust for scaling to integers of rationals */
+	  itomp (ZERO, A[i][0]);
+	  for (j = 2; j <= d; j++)	/* transform each input row */
+	    {
+	      copy (Temp, A[i][j]);
+	      mulint (A[i][j], Temp, Temp);
+	      linint (A[i][0], ONE, Temp, ONE);
+	      linint (A[i][j - 1], ZERO, A[i][j], -TWO);
+	      mulint (scale, A[i][j - 1], A[i][j - 1]);
+	    }			/* end of for (j=1;..) */
+	  copy (A[i][d], scale);
+	  mulint (scale, A[i][d], A[i][d]);
+	}/* end of for (i=1;..) */
+	#ifndef PLRS			
+	      if (Q->debug)
+		printA (D, Q);
+	#endif
+    }				/* end of if(voronoi)     */
+  if (!Q->maximize && !Q->minimize)
+    for (j = 0; j <= d; j++)
+      itomp (ZERO, A[0][j]);
+
+/* Now we pivot to standard form, and then find a primal feasible basis       */
+/* Note these steps MUST be done, even if restarting, in order to get         */
+/* the same index/inequality correspondance we had for the original prob.     */
+/* The inequality array is used to give the insertion order                   */
+/* and is defaulted to the last d rows when givenstart=FALSE                  */
+
+  if(Q->nonnegative) 
+   {
+/* no need for initial pivots here, labelling already done */
+     Q->lastdv = d;
+     Q->nredundcol = 0;
+   }
+  else
+  {
+     if (!getabasis (D, Q, inequality))
+          return FALSE;
+/* bug fix 2009.12.2 */
+     nlinearity=Q->nlinearity;   /*may have been reset if some lins are redundant*/
+  }
+
+
+
+#ifndef PLRS
+  if(Q->debug)
+  {
+	
+    	fprintf(lrs_ofp,"\nafter getabasis");
+    	printA(D, Q);
+	
+  }
+#endif
+  nredundcol = Q->nredundcol;
+  lastdv = Q->lastdv;
+  d = D->d;
+
+
+
+/********************************************************************/
+/* now we start printing the output file  unless no output requested */
+/********************************************************************/
+
+  if (!no_output || Q->debug)
+    {
+	
+
+	
+      if (Q->voronoi){
+	#ifndef PLRS
+	fprintf (lrs_ofp, "\n*Voronoi Diagram: Voronoi vertices and rays are output");
+	#else
+	char *type = "header";
+	char *data = "*Voronoi Diagram: Voronoi vertices and rays are output";
+	//post output in a nonblocking manner (a consumer thread will manage output)
+	post_output(type,data);
+	#endif
+	}
+      if (hull){
+	#ifndef PLRS
+	fprintf (lrs_ofp, "\nH-representation");
+	#else
+	char *type = "header";
+	char *data = "H-representation";
+	//post output in a nonblocking manner (a consumer thread will manage output)
+	post_output(type, data);
+	#endif
+	}
+      else{
+	#ifndef PLRS
+	fprintf (lrs_ofp, "\nV-representation");
+	#else
+	char *type = "header";
+	char *data = "V-representation";
+	//post output in a nonblocking manner (a consumer thread will manage output)
+	post_output(type,data);
+	#endif
+	}	
+	
+
+/* Print linearity space                 */
+/* Don't print linearity if first column zero in hull computation */
+
+      if (hull && Q->homogeneous)
+	k = 1;			/* 0 normally, 1 for homogeneous case     */
+      else
+	k = 0;
+
+      if (nredundcol > k)
+	{
+		#ifndef PLRS
+	  	fprintf (lrs_ofp, "\nlinearity %ld ", nredundcol - k);	/*adjust nredundcol for homog. */
+		#else
+		stringstream ss;
+		char *type = "header";
+		ss<<"linearity "<<(nredundcol -k);
+		#endif
+	  	for (i = 1; i <= nredundcol - k; i++){
+			#ifndef PLRS
+	    		fprintf (lrs_ofp, " %ld", i);
+			#else
+			ss<<" "<<i;
+			#endif
+		}
+		#ifdef PLRS
+		//post output in a nonblocking manner (a consumer thread will manage output)
+		post_output(type, ss.str().c_str());
+		#endif
+	}			/* end print of linearity space */
+
+	#ifndef PLRS
+      	fprintf (lrs_ofp, "\nbegin");
+      	fprintf (lrs_ofp, "\n***** %ld rational", Q->n);
+	#else
+	char *type = "header";
+	stringstream ss;
+	ss<<"begin"<<endl<<"***** "<<Q->n<<" rational";
+	post_output(type, ss.str().c_str());
+	#endif
+    }
+
+			/* end of if !no_output .......   */
+
+
+/* Reset up the inequality array to remember which index is which input inequality */
+/* inequality[B[i]-lastdv] is row number of the inequality with index B[i]              */
+/* inequality[C[i]-lastdv] is row number of the inequality with index C[i]              */
+
+  for (i = 1; i <= m; i++)
+    inequality[i] = i;
+  if (nlinearity > 0)		/* some cobasic indices will be removed */
+    {
+      for (i = 0; i < nlinearity; i++)	/* remove input linearity indices */
+	inequality[linearity[i]] = 0;
+      k = 1;			/* counter for linearities         */
+      for (i = 1; i <= m - nlinearity; i++)
+	{
+	  while (k <= m && inequality[k] == 0)
+	    k++;		/* skip zeroes in corr. to linearity */
+	  inequality[i] = inequality[k++];
+	}
+    }
+				/* end if linearity */
+
+ #ifndef PLRS
+  if (Q->debug)
+    {
+      fprintf (lrs_ofp, "\ninequality array initialization:");
+      for (i = 1; i <= m - nlinearity; i++)
+	fprintf (lrs_ofp, " %ld", inequality[i]);
+	
+    }
+ #endif
+
+
+
+
+  if (nredundcol > 0)
+    {
+      const unsigned int Qn = Q->n;
+      *Lin = lrs_alloc_mp_matrix (nredundcol, Qn);
+
+      for (i = 0; i < nredundcol; i++)
+	{
+		
+
+	  if (!(Q->homogeneous && Q->hull && i == 0))	/* skip redund col 1 for homog. hull */
+	    {
+		
+	      lrs_getray (D, Q, Col[0], D->C[0] + i - hull, (*Lin)[i]);		/* adjust index for deletions */
+	    }
+
+	  if (!removecobasicindex (D, Q, 0L))
+	    {
+	      lrs_clear_mp_matrix (*Lin, nredundcol, Qn);
+	      return FALSE;
+	    }
+	}
+	
+
+    }				/* end if nredundcol > 0 */
+
+
+
+  	if (Q->lponly || Q->nash ){
+		if (Q->verbose)
+		{
+			fprintf (lrs_ofp, "\nNumber of pivots for starting dictionary: %ld",Q->count[3]);
+			if(Q->lponly)
+			     printA (D, Q);
+		}
+        }
+
+/* Do dual pivots to get primal feasibility */
+  if (!primalfeasible (D, Q))
+    {
+#ifndef LRS_QUIET
+          fprintf (lrs_ofp, "\nNo feasible solution\n");
+#endif
+     if (Q->nash && Q->verbose )
+      {
+          fprintf (lrs_ofp, "\nNumber of pivots for feasible solution: %ld",Q->count[3]);
+          fprintf (lrs_ofp, " - No feasible solution");
+      }
+      return FALSE;
+    }
+
+  if (Q->lponly || Q->nash )
+      if (Q->verbose)
+     {
+      fprintf (lrs_ofp, "\nNumber of pivots for feasible solution: %ld",Q->count[3]);
+      if(Q->lponly)
+	      printA (D, Q);
+     }
+
+
+
+/* Now solve LP if objective function was given */
+  if (Q->maximize || Q->minimize)
+    {
+      Q->unbounded = !lrs_solvelp (D, Q, Q->maximize);
+      if (Q->lponly)		
+        {
+
+	#ifndef PLRS
+         if (Q->verbose)
+         {
+           fprintf (lrs_ofp, "\nNumber of pivots for optimum solution: %ld",Q->count[3]);
+           printA (D, Q);
+          }
+	#endif
+          lrs_clear_mp(Temp); lrs_clear_mp(scale);
+          return TRUE;
+        }
+
+      else                         /* check to see if objective is dual degenerate */
+       {
+	  j = 1;
+	  while (j <= d && !zero (A[0][j]))
+	    j++;
+	  if (j <= d)
+	    Q->dualdeg = TRUE;
+	}
+    }
+  else
+/* re-initialize cost row to -det */
+    {
+      for (j = 1; j <= d; j++)
+	{
+	  copy (A[0][j], D->det);
+	  storesign (A[0][j], NEG);
+	}
+
+      itomp (ZERO, A[0][0]);	/* zero optimum objective value */
+    }
+
+
+/* reindex basis to 0..m if necessary */
+/* we use the fact that cobases are sorted by index value */
+#ifndef PLRS
+  if (Q->debug)
+    printA (D, Q);
+#endif
+
+
+
+  while (C[0] <= m)
+    {
+      i = C[0];
+      j = inequality[B[i] - lastdv];
+      inequality[B[i] - lastdv] = inequality[C[0] - lastdv];
+      inequality[C[0] - lastdv] = j;
+      C[0] = B[i];
+      B[i] = i;
+      reorder1 (C, Col, ZERO, d);
+    }
+
+
+
+#ifndef PLRS
+  if (Q->debug)
+    {
+      fprintf (lrs_ofp, "\n*Inequality numbers for indices %ld .. %ld : ", lastdv + 1, m + d);
+      for (i = 1; i <= m - nlinearity; i++)
+	fprintf (lrs_ofp, " %ld ", inequality[i]);
+      printA (D, Q);
+    }
+#endif
+
+
+
+  if (Q->restart)
+    {
+	#ifndef PLRS
+      if (Q->debug)
+	fprintf (lrs_ofp, "\nPivoting to restart co-basis");
+	#endif
+      if (!restartpivots (D, Q))
+	return FALSE;
+      D->lexflag = lexmin (D, Q, ZERO);		/* see if lexmin basis */
+	#ifndef PLRS
+      if (Q->debug)
+	printA (D, Q);
+	#endif
+    }
+
+/* Check to see if necessary to resize */
+  if (Q->inputd > d)
+    *D_p = resize (D, Q);
+
+  lrs_clear_mp(Temp); lrs_clear_mp(scale);
+  return TRUE;
+}
+/********* end of lrs_getfirstbasis  ***************/
+
+
+/*****************************************/
+/* getnextbasis in reverse search order  */
+/*****************************************/
+
+
+long 
+lrs_getnextbasis (lrs_dic ** D_p, lrs_dat * Q, long backtrack)
+	 /* gets next reverse search tree basis, FALSE if none  */
+	 /* switches to estimator if maxdepth set               */
+	 /* backtrack TRUE means backtrack from here            */
+
+{
+  /* assign local variables to structures */
+  long i = 0L, j = 0L;
+  long m = D->m;
+  long d = D->d;
+  long saveflag;
+  long cob_est=0;     /* estimated number of cobases in subtree from current node */
+
+
+  if (backtrack && D->depth == 0)
+    return FALSE;                       /* cannot backtrack from root      */
+
+  if (Q->maxoutput > 0 && Q->count[0]+Q->count[1]-Q->hull >= Q->maxoutput)
+     return FALSE;                      /* output limit reached            */
+
+  while ((j < d) || (D->B[m] != m))	/*main while loop for getnextbasis */
+    {
+      if (D->depth >= Q->maxdepth)
+        {
+          if (Q->runs > 0 && !backtrack )     /*get an estimate of remaining tree */
+            {
+
+//2015.2.9 do iterative estimation backtracking when estimate is small
+
+                 saveflag=Q->printcobasis;
+                 Q->printcobasis=FALSE;
+                 cob_est=lrs_estimate (D, Q);
+                 Q->printcobasis=saveflag;
+                 if(cob_est <= Q->subtreesize) /* stop iterative estimation */
+                  {
+                    if(cob_est > 0)   /* when zero we are at a leaf */
+                       {  lrs_printcobasis(D,Q,ZERO);
+#ifndef PLRS
+                        fprintf(lrs_ofp," cob_est= %ld *subtree",cob_est);
+#else
+                        if (PLRS_DEBUG)
+			{
+				stringstream ss;
+				ss<< "cob_est= " << cob_est << " *subtree" << endl;
+				post_output("debug", ss.str().c_str());
+			}
+#endif
+                       }
+                    backtrack=TRUE;
+                  }
+
+            }
+            else    // either not estimating or we are backtracking
+
+              if (!backtrack && !Q->printcobasis) 
+                 if(!lrs_leaf(D,Q))    /* 2015.6.5 cobasis returned if not a leaf */
+                      lrs_printcobasis(D,Q,ZERO);
+
+            backtrack = TRUE;
+
+ 
+	     if (Q->maxdepth == 0 && cob_est <= Q->subtreesize)	/* root estimate only */
+	       return FALSE;	/* no nextbasis  */
+       }     // if (D->depth >= Q->maxdepth)
+
+
+/*      if ( Q->truncate && negative(D->A[0][0]))*/   /* truncate when moving from opt. vertex */
+/*          backtrack = TRUE;    2011.7.14 */
+
+      if (backtrack)		/* go back to prev. dictionary, restore i,j */
+	{
+	  backtrack = FALSE;
+
+	  if (check_cache (D_p, Q, &i, &j))
+	    {
+	      if (Q->debug)
+		fprintf (lrs_ofp, "\n Cached Dict. restored to depth %ld\n", D->depth);
+	    }
+	  else
+	    {
+	      D->depth--;
+	      selectpivot (D, Q, &i, &j);
+	      pivot (D, Q, i, j);
+	      update (D, Q, &i, &j);	/*Update B,C,i,j */
+	    }
+
+	  if (Q->debug)
+	    {
+	      fprintf (lrs_ofp, "\n Backtrack Pivot: indices i=%ld j=%ld depth=%ld", i, j, D->depth);
+	      printA (D, Q);
+	    };
+
+	  j++;			/* go to next column */
+	}			/* end of if backtrack  */
+
+      if (D->depth < Q->mindepth)
+	break;
+
+      /* try to go down tree */
+
+/* 2011.7.14 patch */
+      while ((j < d) && (!reverse (D, Q, &i, j) || (Q->truncate && Q->minratio[D->m]==1)))
+	j++;
+      if (j == d )
+	backtrack = TRUE;
+      else
+	/*reverse pivot found */
+	{
+	  cache_dict (D_p, Q, i, j);
+	  /* Note that the next two lines must come _after_ the 
+	     call to cache_dict */
+
+	  D->depth++;
+	  if (D->depth > Q->deepest)
+	    Q->deepest++;
+
+	  pivot (D, Q, i, j);
+	  update (D, Q, &i, &j);	/*Update B,C,i,j */
+
+	  D->lexflag = lexmin (D, Q, ZERO);	/* see if lexmin basis */
+	  Q->count[2]++;
+	  Q->totalnodes++;
+
+	  save_basis (*D_p, Q);
+	  if (Q->strace == Q->count[2])
+	    Q->debug = TRUE;
+	  if (Q->etrace == Q->count[2])
+	    Q->debug = FALSE;
+	  return TRUE;		/*return new dictionary */
+	}
+
+
+    }				/* end of main while loop for getnextbasis */
+  return FALSE;			/* done, no more bases */
+}				/*end of lrs_getnextbasis */
+
+/*************************************/
+/* print out one line of output file */
+/*************************************/
+long 
+lrs_getvertex (lrs_dic * P, lrs_dat * Q, lrs_mp_vector output)
+/*Print out current vertex if it is lexmin and return it in output */
+/* return FALSE if no output generated  */
+{
+  lrs_mp_matrix A = P->A;
+
+  long i;
+  long ind;			/* output index                                  */
+  long ired;			/* counts number of redundant columns            */
+/* assign local variables to structures */
+  long *redundcol = Q->redundcol;
+  long *count = Q->count;
+  long *B = P->B;
+  long *Row = P->Row;
+
+  long lastdv = Q->lastdv;
+
+  long hull;
+  long lexflag;
+
+	
+
+  hull = Q->hull;
+  lexflag = P->lexflag;
+  if (lexflag || Q->allbases)
+    ++(Q->count[1]);
+#ifdef PLRS
+        // do not print vertex again in PLRS at root
+	if(P->depth == Q->mindepth ){
+		return FALSE;
+	}
+
+#else
+	//If we are at minimum depth and not at root do not print vertex
+	if(P->depth == Q->mindepth && Q->mindepth != 0){
+		return FALSE;
+	}
+#endif
+
+  if (Q->debug)
+    printA (P, Q);
+
+  linint (Q->sumdet, 1, P->det, 1);
+  if (Q->getvolume)
+   {
+    updatevolume (P, Q);
+    if(Q->verbose)   /* this will print out a triangulation */
+	lrs_printcobasis(P,Q,ZERO);
+   }
+
+
+  /*print cobasis if printcobasis=TRUE and count[2] a multiple of frequency */
+  /* or for lexmin basis, except origin for hull computation - ugly!        */
+
+  if (Q->printcobasis)
+    if ((lexflag && !hull)  || ((Q->frequency > 0) && (count[2] == (count[2] / Q->frequency) * Q->frequency)))
+	if(P->depth != Q->mindepth || Q->mindepth == 0) //Don't print cobasis if this is a restart cobasis
+		lrs_printcobasis(P,Q,ZERO);
+
+  if (hull)
+    return FALSE;		/* skip printing the origin */
+
+  if (!lexflag && !Q->allbases && !Q->lponly)	/* not lexmin, and not printing forced */
+    return FALSE;
+
+
+  /* copy column 0 to output */
+
+  i = 1;
+  ired = 0;
+  copy (output[0], P->det);
+
+  for (ind = 1; ind < Q->n; ind++)	/* extract solution */
+
+    if ((ired < Q->nredundcol) && (redundcol[ired] == ind))
+      /* column was deleted as redundant */
+      {
+	itomp (ZERO, output[ind]);
+	ired++;
+      }
+    else
+      /* column not deleted as redundant */
+      {
+	getnextoutput (P, Q, i, ZERO, output[ind]);
+	i++;
+      }
+
+  reducearray (output, Q->n);
+  if (lexflag && one(output[0]))
+      ++Q->count[4];               /* integer vertex */
+
+
+/* uncomment to print nonzero basic variables 
+
+ printf("\n nonzero basis: vars");
+  for(i=1;i<=lastdv; i++)
+   {
+    if ( !zero(A[Row[i]][0]) )
+         printf(" %ld ",B[i]);
+   }
+*/
+
+/* printslack inequality indices  */
+
+   if (Q->printslack)
+    {
+       fprintf(lrs_ofp,"\nslack ineq:");
+       for(i=lastdv+1;i<=P->m; i++)
+         {
+           if (!zero(A[Row[i]][0]))
+                 fprintf(lrs_ofp," %ld ", Q->inequality[B[i]-lastdv]);
+         }
+    }
+
+
+
+  return TRUE;
+}				/* end of lrs_getvertex */
+
+long 
+lrs_getray (lrs_dic * P, lrs_dat * Q, long col, long redcol, lrs_mp_vector output)
+/*Print out solution in col and return it in output   */
+/*redcol =n for ray/facet 0..n-1 for linearity column */
+/*hull=1 implies facets will be recovered             */
+/* return FALSE if no output generated in column col  */
+{
+  long i;
+  long ind;			/* output index                                  */
+  long ired;			/* counts number of redundant columns            */
+/* assign local variables to structures */
+  long *redundcol = Q->redundcol;
+  long *count = Q->count;
+  long hull = Q->hull;
+  long n = Q->n;
+
+  long *B = P->B;
+  long *Row = P->Row;
+  long lastdv = Q->lastdv;
+
+#ifdef PLRS
+        // do not print vertex again in PLRS at root
+	if(P->depth == Q->mindepth ){
+		return FALSE;
+	}
+
+#else
+	//If we are at minimum depth and not at origin do not print ray
+	if(P->depth == Q->mindepth && Q->mindepth != 0){
+		return FALSE;
+	}
+#endif
+
+
+  if (Q->debug)
+    {
+      printA (P, Q);
+      for (i = 0; i < Q->nredundcol; i++)
+	fprintf (lrs_ofp, " %ld", redundcol[i]);
+      fflush(lrs_ofp);
+    }
+
+  if (redcol == n)
+    {
+      ++count[0];
+      if (Q->printcobasis)
+	if(P->depth != Q->mindepth || Q->mindepth == 0) //Don't print cobasis if this is a restart cobasis
+		lrs_printcobasis(P,Q,col);
+
+    }
+
+  i = 1;
+  ired = 0;
+
+  for (ind = 0; ind < n; ind++)	/* print solution */
+    {
+      if (ind == 0 && !hull)	/* must have a ray, set first column to zero */
+	itomp (ZERO, output[0]);
+
+      else if ((ired < Q->nredundcol) && (redundcol[ired] == ind))
+	/* column was deleted as redundant */
+	{
+	  if (redcol == ind)	/* true for linearity on this cobasic index */
+	    /* we print reduced determinant instead of zero */
+	    copy (output[ind], P->det);
+	  else
+	    itomp (ZERO, output[ind]);
+	  ired++;
+	}
+      else
+	/* column not deleted as redundant */
+	{
+	  getnextoutput (P, Q, i, col, output[ind]);
+	  i++;
+	}
+
+    }
+  reducearray (output, n);
+/* printslack for rays: 2006.10.10 */
+/* printslack inequality indices  */
+
+   if (Q->printslack)
+    {
+       fprintf(lrs_ofp,"\nslack ineq:");
+       for(i=lastdv+1;i<=P->m; i++)
+         {
+           if (!zero(P->A[Row[i]][col]))
+                 fprintf(lrs_ofp," %ld ", Q->inequality[B[i]-lastdv]);
+         }
+    }
+
+	
+
+  return TRUE;
+}				/* end of lrs_getray */
+
+
+void 
+getnextoutput (lrs_dic * P, lrs_dat * Q, long i, long col, lrs_mp out)
+      /* get A[B[i][col] and copy to out */
+{
+  long row;
+  long m = P->m;
+  long d = P->d;
+  long lastdv = Q->lastdv;
+  lrs_mp_matrix A = P->A;
+  long *B = P->B;
+  long *Row = P->Row;
+  long j;
+
+  if (i == d && Q->voronoi)
+    return;			/* skip last column if voronoi set */
+
+  row = Row[i];
+
+  if (Q->nonnegative) 	/* if m+i basic get correct value from dictionary          */
+                        /* the slack for the inequality m-d+i contains decision    */
+                        /* variable x_i. We first see if this is in the basis      */
+                        /* otherwise the value of x_i is zero, except for a ray    */
+                        /* when it is one (det/det) for the actual column it is in */
+    {
+	for (j = lastdv+ 1; j <= m; j++)
+           {
+             if ( Q->inequality[B[j]-lastdv] == m-d+i )
+	          {
+                    copy (out, A[Row[j]][col]);
+                    return;
+                  }
+           }
+/* did not find inequality m-d+i in basis */
+      if ( i == col )
+            copy(out,P->det);
+      else
+            itomp(ZERO,out);
+	
+    }
+  else
+    copy (out, A[row][col]);
+
+}				/* end of getnextoutput */
+
+void 
+lrs_printcobasis (lrs_dic * P, lrs_dat * Q, long col)
+/* col is output column being printed */
+{
+
+	#ifdef PLRS
+	long i;
+	long rflag;/* used to find inequality number for ray column */
+	/* assign local variables to structures */
+	lrs_mp_matrix A = P->A;
+	lrs_mp Nvol, Dvol;	/* hold rescaled det of current basis */
+	long *B = P->B;
+	long *C = P->C;
+	long *Col = P->Col;
+	long *Row = P->Row;
+	long *inequality = Q->inequality;
+	long *temparray = Q->temparray;
+	long *count = Q->count;
+	long hull = Q->hull;
+	long d = P->d;
+	long lastdv = Q->lastdv;
+	long m=P->m;
+	long firstime=TRUE;
+	long nincidence;	/* count number of tight inequalities */
+
+	//Make new output node
+	char *type = "cobasis";
+	//Make stream to collect prat / pmp data
+	stringstream ss;
+
+
+	lrs_alloc_mp(Nvol); lrs_alloc_mp(Dvol);
+
+	if (hull)
+		ss<<"F#"<<count[0]<<" B#"<<count[2]<<" h="<<P->depth<<" vertices/rays ";
+	else if (Q->voronoi)
+		ss<<"V#"<<count[1]<<" R#"<<count[0]<<" B#"<<count[2]<<" h="<<P->depth<<" data points ";
+	else
+		ss<<"V#"<<count[1]<<" R#"<<count[0]<<" B#"<<count[2]<<" h="<<P->depth<<" facets ";
+
+	rflag = (-1);
+	for (i = 0; i < d; i++)
+	{
+		temparray[i] = inequality[C[i] - lastdv];
+	if (Col[i] == col)
+		rflag = temparray[i];	/* look for ray index */
+	}
+	for (i = 0; i < d; i++)
+		reorder (temparray, d);
+	for (i = 0; i < d; i++)
+	{
+		ss<<" "<<temparray[i];
+
+		/* missing cobasis element for ray */
+		if (!(col == ZERO) && (rflag == temparray[i])){ 
+		  	ss<<"*";
+			type = "V cobasis";
+		}
+
+	}
+
+	/* get and print incidence information */
+	if ( col == 0 )
+		nincidence = d;
+	else
+		nincidence = d-1;
+
+	for(i=lastdv+1;i<=m;i++)
+	if ( zero (A[Row[i]][0] ))
+	if( ( col == ZERO ) || zero (A[Row[i]] [col]) ){
+		nincidence++;
+		if( Q->incidence ){
+			if (firstime){
+		    		ss<<" :";
+		    		firstime = FALSE;
+		   	}
+			ss<<inequality[B[i] - lastdv ];
+		}
+	}
+	 
+	ss<<" I#"<<nincidence;
+
+	ss<<pmp (" det=", P->det);
+	//fflush (lrs_ofp);
+	rescaledet (P, Q, Nvol, Dvol); 	/* scales determinant in case input rational */
+
+	ss<<prat(" in_det=",Nvol,Dvol);
+
+	//pipe stream into output node
+	//post output in a nonblocking manner (a consumer thread will manage output)
+	post_output(type, ss.str().c_str());
+
+	lrs_clear_mp(Nvol); lrs_clear_mp(Dvol);
+
+
+	#else
+	long i;
+	long rflag;			/* used to find inequality number for ray column */
+	/* assign local variables to structures */
+	lrs_mp_matrix A = P->A;
+	lrs_mp Nvol, Dvol;		/* hold rescaled det of current basis */
+	long *B = P->B;
+	long *C = P->C;
+	long *Col = P->Col;
+	long *Row = P->Row;
+	long *inequality = Q->inequality;
+	long *temparray = Q->temparray;
+	long *count = Q->count;
+	long hull = Q->hull;
+	long d = P->d;
+	long lastdv = Q->lastdv;
+	long m=P->m;
+	long firstime=TRUE;
+	long nincidence;       /* count number of tight inequalities */
+
+	lrs_alloc_mp(Nvol); lrs_alloc_mp(Dvol);
+
+	if (hull)
+	fprintf (lrs_ofp, "\nF#%ld B#%ld h=%ld vertices/rays ", count[0], count[2], P->depth);
+	else if (Q->voronoi)
+	fprintf (lrs_ofp, "\nV#%ld R#%ld B#%ld h=%ld data points ", count[1], count[0], count[2], P->depth);
+	else
+	fprintf (lrs_ofp, "\nV#%ld R#%ld B#%ld h=%ld facets ", count[1], count[0], count[2], P->depth);
+
+	rflag = (-1);
+	for (i = 0; i < d; i++)
+	{
+	temparray[i] = inequality[C[i] - lastdv];
+	if (Col[i] == col)
+	rflag = temparray[i];	/* look for ray index */
+	}
+	for (i = 0; i < d; i++)
+	reorder (temparray, d);
+	for (i = 0; i < d; i++)
+	{
+	fprintf (lrs_ofp, " %ld", temparray[i]);
+
+	if (!(col == ZERO) && (rflag == temparray[i])) /* missing cobasis element for ray */
+	   fprintf (lrs_ofp, "*");
+
+	}
+
+	/* get and print incidence information */
+	if ( col == 0 )
+	nincidence = d;
+	else
+	nincidence = d-1;
+
+	for(i=lastdv+1;i<=m;i++)
+	if ( zero (A[Row[i]][0] ))
+	if( ( col == ZERO ) || zero (A[Row[i]] [col]) )
+	  { 
+	    nincidence++;
+	    if( Q->incidence )
+	      {
+		if (firstime)
+		  {
+		    fprintf (lrs_ofp," :");
+		    firstime = FALSE;
+		   }
+		fprintf(lrs_ofp," %ld",inequality[B[i] - lastdv ] );
+	      }
+	   }
+	 
+	fprintf(lrs_ofp," I#%ld",nincidence);
+
+	pmp (" det=", P->det);
+	fflush (lrs_ofp);
+	rescaledet (P, Q, Nvol, Dvol); 	/* scales determinant in case input rational */
+	prat(" in_det=",Nvol,Dvol);
+        prat (" z=", P->objnum, P->objden);
+	lrs_clear_mp(Nvol); lrs_clear_mp(Dvol);
+	#endif
+  
+
+}				/* end of lrs_printcobasis */
+
+/*********************/
+/* print final totals */
+/*********************/
+void 
+lrs_printtotals (lrs_dic * P, lrs_dat * Q)
+{
+
+
+#ifdef PLRS
+
+	long *count = Q->count;
+	long *startcount = Q->startcount;
+	std::stringstream ss;
+
+	//output node number of basis
+	ss.str("");
+	ss<<count[2] - startcount[2];
+	post_output("basis count", ss.str().c_str());
+
+	if(Q->hull){
+		//output node for number of facets
+		ss.str("");
+		ss<<count[0] - startcount[0];
+		post_output("facet count", ss.str().c_str());
+
+
+      		rescalevolume (P, Q, Q->Nvolume, Q->Dvolume);
+
+		ss.str("");
+		string str1 = prat("",Q->Nvolume,Q->Dvolume);
+//strip trailing blank introduced by prat
+//for some reason next line fails for mp library !   2014.12.3 so no volume is reported!
+#if (defined(LRSLONG) || defined(GMP))
+                ss << str1.substr (0,str1.length()-1);
+#endif
+		post_output("volume", ss.str().c_str());
+
+
+
+	}else{
+		//output node for number of vertices
+		ss.str("");
+		ss<<count[1] - startcount[1];
+		post_output("vertex count", ss.str().c_str());
+
+		//output node for number of rays  
+		ss.str("");
+		ss<<count[0] - startcount[0];
+		post_output("ray count", ss.str().c_str());
+
+		//output node for number of integer vertices
+/* inaccurate for plrs as restart command does not contain number of integer vertices : an overcount is produced */ 
+
+		ss.str("");
+		ss<<count[4] - startcount[4];
+		post_output("integer vertex count", ss.str().c_str());
+	}
+	
+#else
+  long i;
+  double x;
+/* local assignments */
+  double *cest = Q->cest;
+  long *count = Q->count;
+  long *inequality = Q->inequality;
+  long *linearity = Q->linearity;
+  long *temparray = Q->temparray;
+
+  long *C = P->C;
+
+  long hull = Q->hull;
+  long homogeneous = Q->homogeneous;
+  long nlinearity = Q->nlinearity;
+  long nredundcol = Q->nredundcol;
+  long d, lastdv;
+  d = P->d;
+  lastdv = Q->lastdv;
+
+  fprintf (lrs_ofp, "\nend");
+  if (Q->dualdeg)
+    {
+      fprintf (lrs_ofp, "\n*Warning: Starting dictionary is dual degenerate");
+      fprintf (lrs_ofp, "\n*Complete enumeration may not have been produced");
+      if (Q->maximize)
+         fprintf(lrs_ofp,"\n*Recommendation: Add dualperturb option before maximize in input file\n");
+      else
+         fprintf(lrs_ofp,"\n*Recommendation: Add dualperturb option before minimize in input file\n");
+    }
+
+  if (Q->unbounded)
+    {
+      fprintf (lrs_ofp, "\n*Warning: Starting dictionary contains rays");
+      fprintf (lrs_ofp, "\n*Complete enumeration may not have been produced");
+      if (Q->maximize)
+        fprintf(lrs_ofp,"\n*Recommendation: Change or remove maximize option or add bounds\n");
+      else
+        fprintf(lrs_ofp,"\n*Recommendation: Change or remove minimize option or add bounds\n");
+    }
+  if (Q->truncate)
+    fprintf(lrs_ofp,"\n*Tree truncated at each new vertex");
+  if (Q->maxdepth < MAXD)
+    fprintf (lrs_ofp, "\n*Tree truncated at depth %ld", Q->maxdepth);
+  if (Q->maxoutput > 0L)
+    fprintf (lrs_ofp, "\n*Maximum number of output lines = %ld", Q->maxoutput);
+
+
+#ifdef LRSLONG
+  fprintf (lrs_ofp, "\n*Caution: no overflow checking with long integer arithemtic");
+#else
+  if( Q->verbose)
+    {
+      fprintf (lrs_ofp, "\n*Sum of det(B)=");
+      pmp ("", Q->sumdet);
+    }
+#endif
+
+/* next block with volume rescaling must come before estimates are printed */
+
+  if (Q->getvolume)
+    {
+      rescalevolume (P, Q, Q->Nvolume, Q->Dvolume);
+
+      if (Q->polytope)
+	prat ("\n*Volume=", Q->Nvolume, Q->Dvolume);
+      else
+	prat ("\n*Pseudovolume=", Q->Nvolume, Q->Dvolume);
+    }
+
+  if (hull)     /* output things that are specific to hull computation */
+    {
+      fprintf (lrs_ofp, "\n*Totals: facets=%ld bases=%ld", count[0], count[2]);
+
+      if (nredundcol > homogeneous)	/* don't count column 1 as redundant if homogeneous */
+       {
+        fprintf (lrs_ofp, " linearities=%ld", nredundcol - homogeneous);
+        fprintf (lrs_ofp, " facets+linearities=%ld",nredundcol-homogeneous+count[0]);
+       }
+     if(lrs_ofp != stdout)
+       {
+           printf ("\n*Totals: facets=%ld bases=%ld", count[0], count[2]);
+
+           if (nredundcol > homogeneous)	/* don't count column 1 as redundant if homogeneous */
+           {
+            printf (" linearities=%ld", nredundcol - homogeneous);
+            printf (" facets+linearities=%ld",nredundcol-homogeneous+count[0]);
+           }
+       }
+
+
+      if ((cest[2] > 0) || (cest[0] > 0))
+      {
+	fprintf (lrs_ofp, "\n*Estimates: facets=%.0f bases=%.0f", count[0] + cest[0], count[2] + cest[2]);
+        if (Q->getvolume)
+	  {
+	    rattodouble (Q->Nvolume, Q->Dvolume, &x);
+	    for (i = 2; i < d; i++)
+	      cest[3] = cest[3] / i;	/*adjust for dimension */
+	    fprintf (lrs_ofp, " volume=%g", cest[3] + x);
+	  }
+
+      fprintf (lrs_ofp, "\n*Total number of tree nodes evaluated: %ld", Q->totalnodes);
+#ifdef TIMES
+      fprintf (lrs_ofp, "\n*Estimated total running time=%.1f secs ",(count[2]+cest[2])/Q->totalnodes*get_time () );
+#endif
+
+      }
+/*    Should not happen since we homogenize    */
+/*
+      if ( Q-> restart || Q->allbases || (count[0] > 1 && !Q->homogeneous && !Q->polytope))
+	    fprintf (lrs_ofp, "\n*Note! Duplicate facets may be present");
+*/
+
+    }
+  else         /* output things specific to vertex/ray computation */
+    {
+      fprintf (lrs_ofp, "\n*Totals: vertices=%ld rays=%ld bases=%ld", count[1], count[0], count[2]);
+
+      fprintf (lrs_ofp, " integer_vertices=%ld ",count[4]);
+
+      if (nredundcol > 0)
+        fprintf (lrs_ofp, " linearities=%ld", nredundcol);
+      if ( count[0] + nredundcol > 0 )
+         {
+           fprintf (lrs_ofp, " vertices+rays");
+           if ( nredundcol > 0 )
+              fprintf (lrs_ofp, "+linearities");
+           fprintf (lrs_ofp, "=%ld",nredundcol+count[0]+count[1]);
+         }
+
+      if(lrs_ofp != stdout)
+       { 
+	printf ("\n*Totals: vertices=%ld rays=%ld bases=%ld", count[1], count[0], count[2]);
+
+        printf (" integer_vertices=%ld ",count[4]);
+
+      if (nredundcol > 0)
+           printf (" linearities=%ld", nredundcol);
+      if ( count[0] + nredundcol > 0 )
+         {
+           printf (" vertices+rays");
+           if ( nredundcol > 0 )
+              printf ("+linearities");
+              printf ("=%ld",nredundcol+count[0]+count[1]);
+         }
+        } /* end lrs_ofp != stdout */
+
+
+      if ((cest[2] > 0) || (cest[0] > 0))
+        {
+	fprintf (lrs_ofp, "\n*Estimates: vertices=%.0f rays=%.0f", count[1]+cest[1], count[0]+cest[0]);
+	fprintf (lrs_ofp, " bases=%.0f integer_vertices=%.0f ",count[2]+cest[2], count[4]+cest[4]);
+
+         if (Q->getvolume)
+	   {
+	     rattodouble (Q->Nvolume, Q->Dvolume, &x);
+	     for (i = 2; i <= d-homogeneous; i++)
+	       cest[3] = cest[3] / i;	/*adjust for dimension */
+	     fprintf (lrs_ofp, " pseudovolume=%g", cest[3] + x);
+	   }
+         fprintf (lrs_ofp, "\n*Total number of tree nodes evaluated: %ld", Q->totalnodes);
+#ifdef TIMES
+         fprintf (lrs_ofp, "\n*Estimated total running time=%.1f secs ",(count[2]+cest[2])/Q->totalnodes*get_time () );
+#endif
+        }
+
+      if (Q->restart || Q->allbases)        /* print warning  */
+          fprintf (lrs_ofp, "\n*Note! Duplicate vertices/rays may be present");
+
+      else if ( (count[0] > 1 && !Q->homogeneous))
+          fprintf (lrs_ofp, "\n*Note! Duplicate rays may be present");
+
+    }				/* end of output for vertices/rays */
+
+  fprintf (lrs_ofp, "\n*Dictionary Cache: max size= %ld misses= %ld/%ld   Tree Depth= %ld", dict_count, cache_misses, cache_tries, Q->deepest);
+  if(lrs_ofp != stdout)
+      printf ("\n*Dictionary Cache: max size= %ld misses= %ld/%ld   Tree Depth= %ld", dict_count, cache_misses, cache_tries, Q->deepest);
+
+  if(!Q->verbose)
+     return;
+
+  fprintf (lrs_ofp, "\n*Input size m=%ld rows n=%ld columns", P->m, Q->n);
+  if (hull)
+    fprintf (lrs_ofp, " working dimension=%ld", d - 1 + homogeneous);
+  else
+    fprintf (lrs_ofp, " working dimension=%ld", d);
+
+  fprintf (lrs_ofp, "\n*Starting cobasis defined by input rows");
+  for (i = 0; i < nlinearity; i++)
+    temparray[i] = linearity[i];
+  for (i = nlinearity; i < lastdv; i++)
+    temparray[i] = inequality[C[i - nlinearity] - lastdv];
+  for (i = 0; i < lastdv; i++)
+    reorder (temparray, lastdv);
+  for (i = 0; i < lastdv; i++)
+    fprintf (lrs_ofp, " %ld", temparray[i]);
+
+
+#endif
+
+
+}				/* end of lrs_printtotals */
+/************************/
+/*  Estimation function */
+/************************/
+long  
+lrs_estimate (lrs_dic * P, lrs_dat * Q)
+		   /*returns estimate of subtree size (no. cobases) from current node    */
+		   /*current node is not counted.                   */
+		   /*cest[0]rays [1]vertices [2]bases [3]volume     */
+                   /*    [4] integer vertices                       */
+{
+
+  lrs_mp_vector output;		/* holds one line of output; ray,vertex,facet,linearity */
+  lrs_mp Nvol, Dvol;		/* hold volume of current basis */
+  long estdepth = 0;		/* depth of basis/vertex in subtree for estimate */
+  long i = 0, j = 0, k, nchild, runcount, col;
+  double prod = 0.0;
+  double cave[] =
+  {0.0, 0.0, 0.0, 0.0, 0.0};
+  double nvertices, nbases, nrays, nvol, nivertices;
+  long rays = 0;
+  double newvol = 0.0;
+/* assign local variables to structures */
+  lrs_mp_matrix A = P->A;
+  long *isave = Q->isave;
+  long *jsave = Q->jsave;
+  double *cest = Q->cest;
+  long d = P->d;
+  lrs_alloc_mp(Nvol); lrs_alloc_mp(Dvol);
+/* Main Loop of Estimator */
+
+  output = lrs_alloc_mp_vector (Q->n);	/* output holds one line of output from dictionary     */
+
+  for (runcount = 1; runcount <= Q->runs; runcount++)
+    {				/* runcount counts number of random probes */
+      j = 0;
+      nchild = 1;
+      prod = 1;
+      nvertices = 0.0;
+      nbases = 0.0;
+      nrays = 0.0;
+      nvol = 0.0;
+      nivertices =0.0;
+
+      while (nchild != 0)	/* while not finished yet */
+	{
+
+	  nchild = 0;
+	  while (j < d)
+	    {
+	      if (reverse (P, Q, &i, j))
+		{
+		  isave[nchild] = i;
+		  jsave[nchild] = j;
+		  nchild++;
+		}
+	      j++;
+	    }
+
+	  if (estdepth == 0 && nchild == 0)
+	    {
+	      cest[0] = cest[0] + rays;		/* may be some rays here */
+              lrs_clear_mp(Nvol); lrs_clear_mp(Dvol);
+              lrs_clear_mp_vector(output, Q->n);
+	      return(0L);		/*subtree is a leaf */
+	    }
+
+	  prod = prod * nchild;
+	  nbases = nbases + prod;
+	  if (Q->debug)
+	    {
+	      fprintf (lrs_ofp, "   degree= %ld ", nchild);
+	      fprintf (lrs_ofp, "\nPossible reverse pivots: i,j=");
+	      for (k = 0; k < nchild; k++)
+		fprintf (lrs_ofp, "%ld,%ld ", isave[k], jsave[k]);
+	    }
+
+	  if (nchild > 0)	/*reverse pivot found choose random child */
+	    {
+	      k = myrandom (Q->seed, nchild);
+	      Q->seed = myrandom (Q->seed, 977L);
+	      i = isave[k];
+	      j = jsave[k];
+	      if (Q->debug)
+		fprintf (lrs_ofp, "  selected pivot k=%ld seed=%ld ", k, Q->seed);
+	      estdepth++;
+	      Q->totalnodes++;	/* calculate total number of nodes evaluated */
+	      pivot (P, Q, i, j);
+	      update (P, Q, &i, &j);	/*Update B,C,i,j */
+	      if (lexmin (P, Q, ZERO))	/* see if lexmin basis for vertex */
+                {
+		  nvertices = nvertices + prod;
+                                                    /* integer vertex estimate */
+                  if( lrs_getvertex(P,Q,output))
+                      { --Q->count[1];
+                       if  (one(output[0] ))
+                         { --Q->count[4];
+                           nivertices = nivertices + prod;
+                         }
+                      }
+                }
+
+	      rays = 0;
+	      for (col = 1; col <= d; col++)
+		if (negative (A[0][col]) && (lrs_ratio (P, Q, col) == 0) && lexmin (P, Q, col))
+		  rays++;
+	      nrays = nrays + prod * rays;	/* update ray info */
+
+	      if (Q->getvolume)
+		{
+		  rescaledet (P, Q, Nvol, Dvol);	/* scales determinant in case input rational */
+		  rattodouble (Nvol, Dvol, &newvol);
+		  nvol = nvol + newvol * prod;	/* adjusts volume for degree              */
+		}
+	      j = 0;
+	    }
+	}
+      cave[0] = cave[0] + nrays;
+      cave[1] = cave[1] + nvertices;
+      cave[2] = cave[2] + nbases;
+      cave[3] = cave[3] + nvol;
+      cave[4] = cave[4] + nivertices;
+
+/*  backtrack to root and do it again */
+
+      while (estdepth > 0)
+	{
+	  estdepth = estdepth - 1;
+	  selectpivot (P, Q, &i, &j);
+	  pivot (P, Q, i, j);
+	  update (P, Q, &i, &j);	/*Update B,C,i,j */
+	  /*fprintf(lrs_ofp,"\n0  +++"); */
+	  if (Q->debug)
+	    {
+	      fprintf (lrs_ofp, "\n Backtrack Pivot: indices i,j %ld %ld ", i, j);
+	      printA (P, Q);
+	    }
+	  j++;
+	}
+
+    }				/* end of for loop on runcount */
+
+  j=(long) cave[2]/Q->runs;
+
+//2015.2.9   Do not update totals if we do iterative estimating and subtree is too big
+  if(Q->subtreesize == 0  || j <= Q->subtreesize )
+       for (i = 0; i < 5; i++)
+           cest[i] = cave[i] / Q->runs + cest[i];
+
+  
+  lrs_clear_mp(Nvol); lrs_clear_mp(Dvol);
+  lrs_clear_mp_vector(output, Q->n);
+  return(j);
+}				/* end of lrs_estimate  */
+
+
+/*********************************/
+/* Internal functions            */
+/*********************************/
+/* Basic Dictionary functions    */
+/******************************* */
+
+long 
+reverse (lrs_dic * P, lrs_dat * Q, long *r, long s)
+/*  find reverse indices  */
+/* TRUE if B[*r] C[s] is a reverse lexicographic pivot */
+{
+  long i, j, enter, row, col;
+
+/* assign local variables to structures */
+  lrs_mp_matrix A = P->A;
+  long *B = P->B;
+  long *C = P->C;
+  long *Row = P->Row;
+  long *Col = P->Col;
+  long d = P->d;
+
+  enter = C[s];
+  col = Col[s];
+  if (Q->debug)
+    {
+      fprintf (lrs_ofp, "\n+reverse: col index %ld C %ld Col %ld ", s, enter, col);
+      fflush (lrs_ofp);
+    }
+  if (!negative (A[0][col]))
+    {
+      if (Q->debug)
+	fprintf (lrs_ofp, " Pos/Zero Cost Coeff");
+      Q->minratio[P->m]=0;  /* 2011.7.14 */
+      return (FALSE);
+    }
+
+  *r = lrs_ratio (P, Q, col);
+  if (*r == 0)			/* we have a ray */
+    {
+      if (Q->debug)
+	fprintf (lrs_ofp, " Pivot col non-negative:  ray found");
+      Q->minratio[P->m]=0;  /* 2011.7.14 */
+      return (FALSE);
+    }
+
+  row = Row[*r];
+
+/* check cost row after "pivot" for smaller leaving index    */
+/* ie. j s.t.  A[0][j]*A[row][col] < A[0][col]*A[row][j]     */
+/* note both A[row][col] and A[0][col] are negative          */
+
+  for (i = 0; i < d && C[i] < B[*r]; i++)
+    if (i != s)
+      {
+	j = Col[i];
+	if (positive (A[0][j]) || negative (A[row][j]))		/*or else sign test fails trivially */
+	  if ((!negative (A[0][j]) && !positive (A[row][j])) ||
+	      comprod (A[0][j], A[row][col], A[0][col], A[row][j]) == -1)
+	    {			/*+ve cost found */
+	      if (Q->debug)
+               {
+		fprintf (lrs_ofp, "\nPositive cost found: index %ld C %ld Col %ld", i, C[i], j);
+                fflush(lrs_ofp);
+               }
+              Q->minratio[P->m]=0;  /* 2011.7.14 */
+
+	      return (FALSE);
+	    }
+      }
+  if (Q->debug)
+    {
+      fprintf (lrs_ofp, "\n+end of reverse : indices r %ld s %ld \n", *r, s);
+      fflush (stdout);
+    }
+  return (TRUE);
+}				/* end of reverse */
+
+long 
+selectpivot (lrs_dic * P, lrs_dat * Q, long *r, long *s)
+/* select pivot indices using lexicographic rule   */
+/* returns TRUE if pivot found else FALSE          */
+/* pivot variables are B[*r] C[*s] in locations Row[*r] Col[*s] */
+{
+  long j, col;
+/* assign local variables to structures */
+  lrs_mp_matrix A = P->A;
+  long *Col = P->Col;
+  long d = P->d;
+
+  *r = 0;
+  *s = d;
+  j = 0;
+
+/*find positive cost coef */
+  while ((j < d) && (!positive (A[0][Col[j]])))
+    j++;
+
+  if (j < d)			/* pivot column found! */
+    {
+      *s = j;
+      col = Col[j];
+
+      /*find min index ratio */
+      *r = lrs_ratio (P, Q, col);
+      if (*r != 0)
+	return (TRUE);		/* unbounded */
+    }
+  return (FALSE);
+}				/* end of selectpivot        */
+/******************************************************* */
+
+void 
+pivot (lrs_dic * P, lrs_dat * Q, long bas, long cob)
+		     /* Qpivot routine for array A              */
+		     /* indices bas, cob are for Basis B and CoBasis C    */
+		     /* corresponding to row Row[bas] and column       */
+		     /* Col[cob]   respectively                       */
+{
+  long r, s;
+  long i, j;
+  lrs_mp Ns, Nt, Ars;
+/* assign local variables to structures */
+
+  lrs_mp_matrix A = P->A;
+  long *B = P->B;
+  long *C = P->C;
+  long *Row = P->Row;
+  long *Col = P->Col;
+  long d, m_A;
+
+  lrs_alloc_mp(Ns); lrs_alloc_mp(Nt); lrs_alloc_mp(Ars);
+
+  d = P->d;
+  m_A = P->m_A;
+  Q->count[3]++;    /* count the pivot */
+
+  r = Row[bas];
+  s = Col[cob];
+
+/* Ars=A[r][s]    */
+  if (Q->debug)
+    {
+      fprintf (lrs_ofp, "\n pivot  B[%ld]=%ld  C[%ld]=%ld ", bas, B[bas], cob, C[cob]);
+      printA(P,Q);
+      fflush (stdout);
+    }
+  copy (Ars, A[r][s]);
+  storesign (P->det, sign (Ars));	/*adjust determinant to new sign */
+
+
+  for (i = 0; i <= m_A; i++)
+    if (i != r)
+      for (j = 0; j <= d; j++)
+	if (j != s)
+	  {
+/*          A[i][j]=(A[i][j]*Ars-A[i][s]*A[r][j])/P->det; */
+
+	    mulint (A[i][j], Ars, Nt);
+	    mulint (A[i][s], A[r][j], Ns);
+	    decint (Nt, Ns);
+	    exactdivint (Nt, P->det, A[i][j]);
+	  }			/* end if j ....  */
+
+  if (sign (Ars) == POS)
+    {
+      for (j = 0; j <= d; j++)	/* no need to change sign if Ars neg */
+	/*   A[r][j]=-A[r][j];              */
+	if (!zero (A[r][j]))
+	  changesign (A[r][j]);
+    }				/* watch out for above "if" when removing this "}" ! */
+  else
+    for (i = 0; i <= m_A; i++)
+      if (!zero (A[i][s]))
+	changesign (A[i][s]);
+
+  /*  A[r][s]=P->det;                  */
+
+  copy (A[r][s], P->det);		/* restore old determinant */
+  copy (P->det, Ars);
+  storesign (P->det, POS);		/* always keep positive determinant */
+
+
+  if (Q->debug)
+    {
+      fprintf (lrs_ofp, " depth=%ld ", P->depth);
+      pmp ("det=", P->det);
+      fflush(stdout);
+    }
+/* set the new rescaled objective function value */
+
+  mulint (P->det, Q->Lcm[0], P->objden);
+  mulint (Q->Gcd[0], A[0][0], P->objnum);
+
+  if (!Q->maximize)
+        changesign (P->objnum);
+  if (zero (P->objnum))
+        storesign (P->objnum, POS);
+  else
+	reduce (P->objnum,P->objden);
+
+
+  lrs_clear_mp(Ns); lrs_clear_mp(Nt); lrs_clear_mp(Ars);
+}				/* end of pivot */
+
+long 
+primalfeasible (lrs_dic * P, lrs_dat * Q)
+/* Do dual pivots to get primal feasibility */
+/* Note that cost row is all zero, so no ratio test needed for Dual Bland's rule */
+{
+  long primalinfeasible = TRUE;
+  long i, j;
+/* assign local variables to structures */
+  lrs_mp_matrix A = P->A;
+  long *Row = P->Row;
+  long *Col = P->Col;
+  long m, d, lastdv;
+  m = P->m;
+  d = P->d;
+  lastdv = Q->lastdv;
+
+/*temporary: try to get new start after linearity */
+
+  while (primalinfeasible)
+    {
+      i=lastdv+1;
+      while (i <= m && !negative (A[Row[i]][0]) )
+        i++;
+      if (i <= m )
+	{
+	      j = 0;		/*find a positive entry for in row */
+	      while (j < d && !positive (A[Row[i]][Col[j]]))
+		j++;
+	      if (j >= d)
+		return (FALSE);	/* no positive entry */
+	      pivot (P, Q, i, j);
+	      update (P, Q, &i, &j);
+        }
+      else
+         primalinfeasible = FALSE;
+    }				/* end of while primalinfeasibile */
+  return (TRUE);
+}				/* end of primalfeasible */
+
+
+long 
+lrs_solvelp (lrs_dic * P, lrs_dat * Q, long maximize)
+/* Solve primal feasible lp by Dantzig`s rule and lexicographic ratio test */
+/* return TRUE if bounded, FALSE if unbounded                              */
+{
+  long i, j;
+/* assign local variables to structures */
+  long d = P->d;
+
+  while (dan_selectpivot (P, Q, &i, &j))
+    {
+      Q->count[3]++;
+      pivot (P, Q, i, j);
+      update (P, Q, &i, &j);	/*Update B,C,i,j */
+    }
+  if (Q->debug)
+    printA (P, Q);
+
+  if (j < d && i == 0)		/* selectpivot gives information on unbounded solution */
+    {
+#ifndef LRS_QUIET
+      if (Q->lponly)
+	fprintf (lrs_ofp, "\n*Unbounded solution");
+#endif
+      return FALSE;
+    }
+  return TRUE;
+}				/* end of lrs_solvelp  */
+
+long 
+getabasis (lrs_dic * P, lrs_dat * Q, long order[])
+
+/* Pivot Ax<=b to standard form */
+/*Try to find a starting basis by pivoting in the variables x[1]..x[d]        */
+/*If there are any input linearities, these appear first in order[]           */
+/* Steps: (a) Try to pivot out basic variables using order                    */
+/*            Stop if some linearity cannot be made to leave basis            */
+/*        (b) Permanently remove the cobasic indices of linearities           */
+/*        (c) If some decision variable cobasic, it is a linearity,           */
+/*            and will be removed.                                            */
+
+{
+  long i, j, k;
+/* assign local variables to structures */
+  lrs_mp_matrix A = P->A;
+  long *B = P->B;
+  long *C = P->C;
+  long *Row = P->Row;
+  long *Col = P->Col;
+  long *linearity = Q->linearity;
+  long *redundcol = Q->redundcol;
+  long m, d, nlinearity;
+  long nredundcol = 0L;		/* will be calculated here */
+  m = P->m;
+  d = P->d;
+  nlinearity = Q->nlinearity;
+
+
+
+  if (Q->debug)
+    {
+      fprintf (lrs_ofp, "\ngetabasis from inequalities given in order");
+      for (i = 0l; i < m; i++)
+	fprintf (lrs_ofp, " %ld", order[i]);
+    }
+	for (j = 0l; j < m; j++)
+	{
+		i = 0l;
+		while (i <= m && B[i] != d + order[j])
+			i++;			/* find leaving basis index i */
+		if (j < nlinearity && i > m)	/* cannot pivot linearity to cobasis */
+		{
+	  		if (Q->debug)
+	    			printA (P, Q);
+			#ifndef LRS_QUIET
+	  		fprintf (lrs_ofp, "\nCannot find linearity in the basis");
+			#endif
+	  		return FALSE;
+		}
+		if (i <= m)
+		{			/* try to do a pivot */
+	  		k = 0l;
+	  		while (C[k] <= d && zero (A[Row[i]][Col[k]])){
+	    			k++;
+			}
+	  		if (C[k] <= d)
+	    		{
+				
+	      			pivot (P, Q, i, k);
+	      			update (P, Q, &i, &k);
+    			}
+	  		else if (j < nlinearity)
+	    		{			/* cannot pivot linearity to cobasis */
+				if (zero (A[Row[i]][0]))
+				{
+					#ifndef LRS_QUIET
+			  		fprintf (lrs_ofp, "\n*Input linearity in row %ld is redundant--converted to inequality", order[j]);
+					#endif
+			  		linearity[j]=0l;
+				}
+		      		else
+				{
+			  		if (Q->debug)
+			    			printA (P, Q);
+					#ifndef LRS_QUIET
+			  		fprintf (lrs_ofp, "\n*Input linearity in row %ld is inconsistent with earlier linearities", order[j]);
+			  		fprintf (lrs_ofp, "\n*No feasible solution");
+					#endif
+			  		return FALSE;
+				}
+	    		}
+			
+
+		}
+	}
+
+
+/* update linearity array to get rid of redundancies */
+  i = 0;
+  k = 0;			/* counters for linearities         */
+  while (k < nlinearity)
+    {
+      while (k < nlinearity && linearity[k] == 0)
+	k++;
+      if (k < nlinearity)
+	linearity[i++] = linearity[k++];
+    }
+
+  nlinearity = i;
+/* bug fix, 2009.6.27 */    Q->nlinearity = i;
+
+/* column dependencies now can be recorded  */
+/* redundcol contains input column number 0..n-1 where redundancy is */
+  k = 0;
+  while (k < d && C[k] <= d)
+    {
+      if (C[k] <= d){		/* decision variable still in cobasis */
+	redundcol[nredundcol++] = C[k] - Q->hull;	/* adjust for hull indices */
+		
+	}
+      k++;
+    }
+
+/* now we know how many decision variables remain in problem */
+  Q->nredundcol = nredundcol;
+  Q->lastdv = d - nredundcol;
+  if (Q->debug)
+    {
+      fprintf (lrs_ofp, "\nend of first phase of getabasis: ");
+      fprintf (lrs_ofp, "lastdv=%ld nredundcol=%ld", Q->lastdv, Q->nredundcol);
+      fprintf (lrs_ofp, "\nredundant cobases:");
+      for (i = 0; i < nredundcol; i++)
+	fprintf (lrs_ofp, " %ld", redundcol[i]);
+      printA (P, Q);
+    }
+
+/* Remove linearities from cobasis for rest of computation */
+/* This is done in order so indexing is not screwed up */
+
+  for (i = 0; i < nlinearity; i++)
+    {				/* find cobasic index */
+      k = 0;
+      while (k < d && C[k] != linearity[i] + d)
+	k++;
+      if (k >= d)
+	{
+	  fprintf (lrs_ofp, "\nError removing linearity");
+	  return FALSE;
+	}
+      if (!removecobasicindex (P, Q, k))
+	return FALSE;
+      d = P->d;
+    }
+  if (Q->debug && nlinearity > 0)
+    printA (P, Q);
+/* set index value for first slack variable */
+
+/* Check feasability */
+  if (Q->givenstart)
+    {
+      i = Q->lastdv + 1;
+      while (i <= m && !negative (A[Row[i]][0]))
+	i++;
+      if (i <= m)
+	fprintf (lrs_ofp, "\n*Infeasible startingcobasis - will be modified");
+    }
+  return TRUE;
+}				/*  end of getabasis */
+
+long 
+removecobasicindex (lrs_dic * P, lrs_dat * Q, long k)
+/* remove the variable C[k] from the problem */
+/* used after detecting column dependency    */
+{
+  long i, j, cindex, deloc;
+/* assign local variables to structures */
+  lrs_mp_matrix A = P->A;
+  long *B = P->B;
+  long *C = P->C;
+  long *Col = P->Col;
+  long m, d;
+  m = P->m;
+  d = P->d;
+
+  if (Q->debug)
+    fprintf (lrs_ofp, "\nremoving cobasic index k=%ld C[k]=%ld", k, C[k]);
+  cindex = C[k];		/* cobasic index to remove              */
+  deloc = Col[k];		/* matrix column location to remove     */
+
+  for (i = 1; i <= m; i++)	/* reduce basic indices by 1 after index */
+    if (B[i] > cindex)
+      B[i]--;
+  
+  for (j = k; j < d; j++)	/* move down other cobasic variables    */
+    {
+      C[j] = C[j + 1] - 1;	/* cobasic index reduced by 1           */
+      Col[j] = Col[j + 1];
+    }
+
+  if (deloc != d)               
+    {
+  /* copy col d to deloc */
+      for (i = 0; i <= m; i++)
+        copy (A[i][deloc], A[i][d]);
+
+  /* reassign location for moved column */
+      j = 0;
+      while (Col[j] != d)
+        j++;
+      
+    Col[j] = deloc;
+  }
+
+  P->d--;
+  if (Q->debug)
+    printA (P, Q);
+  return TRUE;
+}				/* end of removecobasicindex */
+
+lrs_dic *
+resize (lrs_dic * P, lrs_dat * Q)
+	/* resize the dictionary after some columns are deleted, ie. inputd>d */
+	/* a new lrs_dic record is created with reduced size, and items copied over */
+{
+  lrs_dic *P1;			/* to hold new dictionary in case of resizing */
+
+  long i, j;
+  long m, d, m_A;
+
+
+  m = P->m;
+  d = P->d;
+  m_A = P->m_A;
+
+/* get new dictionary record */
+
+  P1 = new_lrs_dic (m, d, m_A);
+
+/* copy data from P to P1    */
+  P1->i = P->i;
+  P1->j = P->j;
+  P1->depth = P->depth;
+  P1->m = P->m;
+  P1->d = P1->d_orig = d;
+  P1->lexflag = P->lexflag;
+  P1->m_A = P->m_A;
+  copy (P1->det, P->det);
+  copy (P1->objnum, P->objnum);
+  copy (P1->objden, P->objden);
+
+  for (i = 0; i <= m; i++)
+    {
+      P1->B[i] = P->B[i];
+      P1->Row[i] = P->Row[i];
+    }
+  for (i = 0; i <= m_A; i++)
+    {
+      for (j = 0; j <= d; j++)
+	copy (P1->A[i][j], P->A[i][j]);
+    }
+
+
+  for (j = 0; j <= d; j++)
+    {
+      P1->Col[j] = P->Col[j];
+      P1->C[j] = P->C[j];
+    }
+
+  if (Q->debug)
+    {
+      fprintf (lrs_ofp, "\nDictionary resized from d=%ld to d=%ld", Q->inputd, P->d);
+      printA (P1, Q);
+    }
+
+  lrs_free_dic (P,Q);
+
+/* Reassign cache pointers */
+
+  Q->Qhead = P1;
+  Q->Qtail = P1;
+  P1->next = P1;
+  P1->prev = P1;
+
+  return P1;
+
+}
+/********* resize                    ***************/
+
+
+long 
+restartpivots (lrs_dic * P, lrs_dat * Q)
+/* facet contains a list of the inequalities in the cobasis for the restart */
+/* inequality contains the relabelled inequalities after initialization     */
+{
+  long i, j, k;
+  long *Cobasic;		/* when restarting, Cobasic[j]=1 if j is in cobasis */
+/* assign local variables to structures */
+  lrs_mp_matrix A = P->A;
+  long *B = P->B;
+  long *C = P->C;
+  long *Row = P->Row;
+  long *Col = P->Col;
+  long *inequality = Q->inequality;
+  long *facet = Q->facet;
+  long nlinearity = Q->nlinearity;
+  long m, d, lastdv;
+  m = P->m;
+  d = P->d;
+  lastdv = Q->lastdv;
+
+  Cobasic = (long *) CALLOC ((unsigned) m + d + 2, sizeof (long));
+
+  if (Q->debug)
+    fprintf(lrs_ofp,"\nCobasic flags in restartpivots");
+  /* set Cobasic flags */
+  for (i = 0; i < m + d + 1; i++)
+    Cobasic[i] = 0;
+  for (i = 0; i < d; i++)	/* find index corresponding to facet[i] */
+    {
+      j = 1;
+      while (facet[i + nlinearity] != inequality[j])
+	j++;
+      Cobasic[j + lastdv] = 1;
+      if (Q->debug)
+        fprintf(lrs_ofp," %ld %ld;",facet[i+nlinearity],j+lastdv);
+    }
+
+  /* Note that the order of doing the pivots is important, as */
+  /* the B and C vectors are reordered after each pivot       */
+
+/* code below replaced 2006.10.30 */
+/*
+
+  for (i = m; i >= d + 1; i--)	
+    if (Cobasic[B[i]])	
+      {
+	k = d - 1;
+	while ((k >= 0) &&
+	       (zero (A[Row[i]][Col[k]]) || Cobasic[C[k]]))
+	  k--;
+	if (k >= 0)
+	  {
+	    pivot (P, Q, i, k);
+	    update (P, Q, &i, &k);
+	  }
+	else
+	  {
+	    fprintf (lrs_ofp, "\nInvalid Co-basis - does not have correct rank");
+            free(Cobasic);
+	    return FALSE;
+	  }
+      }
+*/     
+/*end of code that was replaced */	
+
+/* Suggested new code from db starts */
+  i=m;
+  while (i>d){
+    while(Cobasic[B[i]]){
+      k = d - 1;
+      while ((k >= 0) && (zero (A[Row[i]][Col[k]]) || Cobasic[C[k]])){
+       k--;
+      }
+      if (k >= 0)  {
+           /*db asks: should i really be modified here? (see old code) */
+           /*da replies: modifying i only makes is larger, and so      */
+           /*the second while loop will put it back where it was       */
+           /*faster (and safer) as done below                          */
+       long  ii=i;
+       pivot (P, Q, ii, k);
+       update (P, Q, &ii, &k);
+      } else {
+       fprintf (lrs_ofp, "\nInvalid Co-basis - does not have correct rank");
+       free(Cobasic);
+       return FALSE;
+      }
+    }
+    i--;
+  }
+/* Suggested new code from db ends */
+
+  if (lexmin (P, Q, ZERO))
+    --Q->count[1];		/* decrement vertex count if lexmin */
+/* check restarting from a primal feasible dictionary               */
+  for (i = lastdv + 1; i <= m; i++)
+    if (negative (A[Row[i]][0]))
+      {
+	fprintf (lrs_ofp, "\nTrying to restart from infeasible dictionary");
+        free(Cobasic);
+	return FALSE;
+      }
+  free(Cobasic);
+  return TRUE;
+
+}				/* end of restartpivots */
+
+
+long 
+lrs_ratio (lrs_dic * P, lrs_dat * Q, long col)	/*find lex min. ratio */
+		  /* find min index ratio -aig/ais, ais<0 */
+		  /* if multiple, checks successive basis columns */
+		  /* recoded Dec 1997                     */
+{
+  long i, j, comp, ratiocol, basicindex, start, nstart, cindex, bindex;
+  long firstime;		/*For ratio test, true on first pass,else false */
+  lrs_mp Nmin, Dmin;
+  long degencount, ndegencount;
+/* assign local variables to structures */
+  lrs_mp_matrix A = P->A;
+  long *B = P->B;
+  long *Row = P->Row;
+  long *Col = P->Col;
+  long *minratio = Q->minratio;
+  long m, d, lastdv;
+
+  m = P->m;
+  d = P->d;
+  lastdv = Q->lastdv;
+
+
+  nstart=0;
+  ndegencount=0;
+  degencount = 0;
+  minratio[P->m]=1;   /*2011.7.14 non-degenerate pivot flag */
+
+  for (j = lastdv + 1; j <= m; j++)
+    {
+      /* search rows with negative coefficient in dictionary */
+      /*  minratio contains indices of min ratio cols        */
+      if (negative (A[Row[j]][col]))
+       {
+	minratio[degencount++] = j;
+        if(zero (A[Row[j]][0]))
+          minratio[P->m]=0;   /*2011.7.14 degenerate pivot flag */
+       }
+    }				/* end of for loop */
+  if (Q->debug)
+    {
+      fprintf (lrs_ofp, "  Min ratios: ");
+      for (i = 0; i < degencount; i++)
+	fprintf (lrs_ofp, " %ld ", B[minratio[i]]);
+    }
+  if (degencount == 0)
+    return (degencount);	/* non-negative pivot column */
+
+  lrs_alloc_mp(Nmin); lrs_alloc_mp(Dmin);
+  ratiocol = 0;			/* column being checked, initially rhs */
+  start = 0;			/* starting location in minratio array */
+  bindex = d + 1;		/* index of next basic variable to consider */
+  cindex = 0;			/* index of next cobasic variable to consider */
+  basicindex = d;		/* index of basis inverse for current ratio test, except d=rhs test */
+  while (degencount > 1)	/*keep going until unique min ratio found */
+    {
+      if (B[bindex] == basicindex)	/* identity col in basis inverse */
+	{
+	  if (minratio[start] == bindex)
+	    /* remove this index, all others stay */
+	    {
+	      start++;
+	      degencount--;
+	    }
+	  bindex++;
+	}
+      else
+	/* perform ratio test on rhs or column of basis inverse */
+	{
+	  firstime = TRUE;
+	  /*get next ratio column and increment cindex */
+	  if (basicindex != d)
+	    ratiocol = Col[cindex++];
+	  for (j = start; j < start + degencount; j++)
+	    {
+	      i = Row[minratio[j]];	/* i is the row location of the next basic variable */
+	      comp = 1;		/* 1:  lhs>rhs;  0:lhs=rhs; -1: lhs<rhs */
+	      if (firstime)
+		firstime = FALSE;	/*force new min ratio on first time */
+	      else
+		{
+		  if (positive (Nmin) || negative (A[i][ratiocol]))
+		    {
+		      if (negative (Nmin) || positive (A[i][ratiocol]))
+			comp = comprod (Nmin, A[i][col], A[i][ratiocol], Dmin);
+		      else
+			comp = -1;
+		    }
+
+		  else if (zero (Nmin) && zero (A[i][ratiocol]))
+		    comp = 0;
+
+		  if (ratiocol == ZERO)
+		    comp = -comp;	/* all signs reversed for rhs */
+		}
+	      if (comp == 1)
+		{		/*new minimum ratio */
+		  nstart = j;
+		  copy (Nmin, A[i][ratiocol]);
+		  copy (Dmin, A[i][col]);
+		  ndegencount = 1;
+		}
+	      else if (comp == 0)	/* repeated minimum */
+		minratio[nstart + ndegencount++] = minratio[j];
+
+	    }			/* end of  for (j=start.... */
+	  degencount = ndegencount;
+	  start = nstart;
+	}			/* end of else perform ratio test statement */
+      basicindex++;		/* increment column of basis inverse to check next */
+      if (Q->debug)
+	{
+	  fprintf (lrs_ofp, " ratiocol=%ld degencount=%ld ", ratiocol, degencount);
+	  fprintf (lrs_ofp, "  Min ratios: ");
+	  for (i = start; i < start + degencount; i++)
+	    fprintf (lrs_ofp, " %ld ", B[minratio[i]]);
+	}
+    }				/*end of while loop */
+  lrs_clear_mp(Nmin); lrs_clear_mp(Dmin);
+  return (minratio[start]);
+}				/* end of ratio */
+
+
+
+long 
+lexmin (lrs_dic * P, lrs_dat * Q, long col)
+  /*test if basis is lex-min for vertex or ray, if so TRUE */
+  /* FALSE if a_r,g=0, a_rs !=0, r > s          */
+{
+/*do lexmin test for vertex if col=0, otherwise for ray */
+  long r, s, i, j;
+/* assign local variables to structures */
+  lrs_mp_matrix A = P->A;
+  long *B = P->B;
+  long *C = P->C;
+  long *Row = P->Row;
+  long *Col = P->Col;
+  long m = P->m;
+  long d = P->d;
+
+  for (i = Q->lastdv + 1; i <= m; i++)
+    {
+      r = Row[i];
+      if (zero (A[r][col]))	/* necessary for lexmin to fail */
+	for (j = 0; j < d; j++)
+	  {
+	    s = Col[j];
+	    if (B[i] > C[j])	/* possible pivot to reduce basis */
+	      {
+		if (zero (A[r][0]))	/* no need for ratio test, any pivot feasible */
+		  {
+		    if (!zero (A[r][s]))
+		      return (FALSE);
+		  }
+		else if (negative (A[r][s]) && ismin (P, Q, r, s))
+		  {
+		    return (FALSE);
+		  }
+	      }			/* end of if B[i] ... */
+	  }
+    }
+  if ((col != ZERO) && Q->debug)
+    {
+      fprintf (lrs_ofp, "\n lexmin ray in col=%ld ", col);
+      printA (P, Q);
+    }
+  return (TRUE);
+}				/* end of lexmin */
+
+long 
+ismin (lrs_dic * P, lrs_dat * Q, long r, long s)
+/*test if A[r][s] is a min ratio for col s */
+{
+  long i;
+/* assign local variables to structures */
+  lrs_mp_matrix A = P->A;
+  long m_A = P->m_A;
+
+  for (i = 1; i <= m_A; i++)
+    if ((i != r) &&
+	negative (A[i][s]) && comprod (A[i][0], A[r][s], A[i][s], A[r][0]))
+      {
+	return (FALSE);
+      }
+
+  return (TRUE);
+}
+
+void 
+update (lrs_dic * P, lrs_dat * Q, long *i, long *j)
+ /*update the B,C arrays after a pivot */
+ /*   involving B[bas] and C[cob]           */
+{
+
+  long leave, enter;
+/* assign local variables to structures */
+  long *B = P->B;
+  long *C = P->C;
+  long *Row = P->Row;
+  long *Col = P->Col;
+  long m = P->m;
+  long d = P->d;
+
+  leave = B[*i];
+  enter = C[*j];
+  B[*i] = enter;
+  reorder1 (B, Row, *i, m + 1);
+  C[*j] = leave;
+  reorder1 (C, Col, *j, d);
+/* restore i and j to new positions in basis */
+  for (*i = 1; B[*i] != enter; (*i)++);		/*Find basis index */
+  for (*j = 0; C[*j] != leave; (*j)++);		/*Find co-basis index */
+}				/* end of update */
+
+long 
+lrs_degenerate (lrs_dic * P, lrs_dat * Q)
+/* TRUE if the current dictionary is primal degenerate */
+/* not thoroughly tested   2000/02/15                  */
+{
+  long i;
+  long *Row;
+
+  lrs_mp_matrix A = P->A;
+  long d = P->d;
+  long m = P->m;
+
+  Row = P->Row;
+
+  for (i = d + 1; i <= m; i++)
+    if (zero (A[Row[i]][0]))
+      return TRUE;
+
+  return FALSE;
+}
+
+
+/*********************************************************/
+/*                 Miscellaneous                         */
+/******************************************************* */
+
+void 
+reorder (long a[], long range)
+/*reorder array in increasing order with one misplaced element */
+{
+  long i, temp;
+  for (i = 0; i < range - 1; i++)
+    if (a[i] > a[i + 1])
+      {
+	temp = a[i];
+	a[i] = a[i + 1];
+	a[i + 1] = temp;
+      }
+  for (i = range - 2; i >= 0; i--)
+    if (a[i] > a[i + 1])
+      {
+	temp = a[i];
+	a[i] = a[i + 1];
+	a[i + 1] = temp;
+      }
+
+}				/* end of reorder */
+
+void 
+reorder1 (long a[], long b[], long newone, long range)
+/*reorder array a in increasing order with one misplaced element at index newone */
+/*elements of array b are updated to stay aligned with a */
+{
+  long temp;
+  while (newone > 0 && a[newone] < a[newone - 1])
+    {
+      temp = a[newone];
+      a[newone] = a[newone - 1];
+      a[newone - 1] = temp;
+      temp = b[newone];
+      b[newone] = b[newone - 1];
+      b[--newone] = temp;
+    }
+  while (newone < range - 1 && a[newone] > a[newone + 1])
+    {
+      temp = a[newone];
+      a[newone] = a[newone + 1];
+      a[newone + 1] = temp;
+      temp = b[newone];
+      b[newone] = b[newone + 1];
+      b[++newone] = temp;
+    }
+}				/* end of reorder1 */
+
+
+void 
+rescaledet (lrs_dic * P, lrs_dat * Q, lrs_mp Vnum, lrs_mp Vden)
+  /* rescale determinant to get its volume */
+  /* Vnum/Vden is volume of current basis  */
+{
+  lrs_mp gcdprod;		/* to hold scale factors */
+  long i;
+/* assign local variables to structures */
+  long *B = P->B;
+  long *C = P->C;
+  long m, d, lastdv;
+
+  lrs_alloc_mp(gcdprod);
+  m = P->m;
+  d = P->d;
+  lastdv = Q->lastdv;
+
+  itomp (ONE, gcdprod);
+  itomp (ONE, Vden);
+  for (i = 0; i < d; i++)
+    if (B[i] <= m)
+      {
+	mulint (Q->Gcd[Q->inequality[C[i] - lastdv]], gcdprod, gcdprod);
+	mulint (Q->Lcm[Q->inequality[C[i] - lastdv]], Vden, Vden);
+      }
+  mulint (P->det, gcdprod, Vnum);
+  reduce (Vnum, Vden);
+  lrs_clear_mp(gcdprod);
+}				/* end rescaledet */
+
+void 
+rescalevolume (lrs_dic * P, lrs_dat * Q, lrs_mp Vnum, lrs_mp Vden)
+/* adjust volume for dimension */
+{
+  lrs_mp temp, dfactorial;
+/* assign local variables to structures */
+  long lastdv = Q->lastdv;
+
+  lrs_alloc_mp(temp); lrs_alloc_mp(dfactorial);
+
+  /*reduce Vnum by d factorial  */
+  getfactorial (dfactorial, lastdv);
+  mulint (dfactorial, Vden, Vden);
+  if (Q->hull && !Q->homogeneous)
+    {				/* For hull option multiply by d to correct for lifting */
+      itomp (lastdv, temp);
+      mulint (temp, Vnum, Vnum);
+    }
+
+  reduce (Vnum, Vden);
+  lrs_clear_mp(temp); lrs_clear_mp(dfactorial);
+}
+
+
+void 
+updatevolume (lrs_dic * P, lrs_dat * Q)		/* rescale determinant and update the volume */
+{
+  lrs_mp tN, tD, Vnum, Vden;
+  lrs_alloc_mp(tN); lrs_alloc_mp(tD); lrs_alloc_mp(Vnum); lrs_alloc_mp(Vden);
+  rescaledet (P, Q, Vnum, Vden);
+  copy (tN, Q->Nvolume);
+  copy (tD, Q->Dvolume);
+  linrat (tN, tD, ONE, Vnum, Vden, ONE, Q->Nvolume, Q->Dvolume);
+  if (Q->debug)
+    {
+      prat ("\n*Volume=", Q->Nvolume, Q->Dvolume);
+      pmp (" Vnum=", Vnum);
+      pmp (" Vden=", Vden);
+    }
+  lrs_clear_mp(tN); lrs_clear_mp(tD); lrs_clear_mp(Vnum); lrs_clear_mp(Vden);
+
+}				/* end of updatevolume */
+
+
+
+/***************************************************/
+/* Routines for redundancy checking                */
+/***************************************************/
+
+long 
+checkredund (lrs_dic * P, lrs_dat * Q)
+/* Solve primal feasible lp by least subscript and lex min basis method */
+/* to check redundancy of a row in objective function                   */
+/* returns TRUE if redundant, else FALSE                                */
+{
+  lrs_mp Ns, Nt;
+  long i, j;
+  long r, s;
+
+/* assign local variables to structures */
+  lrs_mp_matrix A = P->A;
+  long *Row, *Col;
+  long d = P->d;
+
+  lrs_alloc_mp(Ns); lrs_alloc_mp(Nt);
+  Row = P->Row;
+  Col = P->Col;
+
+  while (selectpivot (P, Q, &i, &j))
+    {
+      Q->count[2]++;
+
+/* sign of new value of A[0][0]            */
+/* is      A[0][s]*A[r][0]-A[0][0]*A[r][s] */
+
+      r = Row[i];
+      s = Col[j];
+
+      mulint (A[0][s], A[r][0], Ns);
+      mulint (A[0][0], A[r][s], Nt);
+
+      if (mp_greater (Ns, Nt))
+        {
+          lrs_clear_mp(Ns); lrs_clear_mp(Nt);
+	  return FALSE;		/* non-redundant */
+        }
+
+      pivot (P, Q, i, j);
+      update (P, Q, &i, &j);	/*Update B,C,i,j */
+
+    }
+
+  lrs_clear_mp(Ns); lrs_clear_mp(Nt);
+
+  return !(j < d && i == 0);	/* unbounded is also non-redundant */
+
+
+}				/* end of checkredund  */
+
+long 
+checkcobasic (lrs_dic * P, lrs_dat * Q, long index)
+/* TRUE if index is cobasic and nonredundant                         */
+/* FALSE if basic, or degen. cobasic, where it will get pivoted out  */
+
+{
+
+/* assign local variables to structures */
+
+  lrs_mp_matrix A = P->A;
+  long *B, *C, *Row, *Col;
+  long d = P->d;
+  long m = P->m;
+  long debug = Q->debug;
+  long i = 0;
+  long j = 0;
+  long s;
+
+  B = P->B;
+  C = P->C;
+  Row = P->Row;
+  Col = P->Col;
+
+
+  while ((j < d) && C[j] != index)
+    j++;
+
+  if (j == d)
+    return FALSE;		/* not cobasic index */
+
+
+/* index is cobasic */
+
+  if (debug)
+    fprintf (lrs_ofp, "\nindex=%ld cobasic", index);
+/* not debugged for new LOC 
+   s=LOC[index];
+ */
+  s = Col[j];
+  i = Q->lastdv + 1;
+
+  while ((i <= m) &&
+	 (zero (A[Row[i]][s]) || !zero (A[Row[i]][0])))
+    i++;
+
+  if (i > m)
+    {
+      if (debug)
+	fprintf (lrs_ofp, " is non-redundant");
+      return TRUE;
+    }
+  if (debug)
+    fprintf (lrs_ofp, " is degenerate B[i]=%ld", B[i]);
+
+  pivot (P, Q, i, j);
+  update (P, Q, &i, &j);	/*Update B,C,i,j */
+
+  return FALSE;			/*index is no longer cobasic */
+
+}				/* end of checkcobasic */
+
+long 
+checkindex (lrs_dic * P, lrs_dat * Q, long index)
+/* 0 if index is non-redundant inequality */
+/* 1 if index is redundant     inequality */
+/* 2 if index is input linearity          */
+/*NOTE: row is returned all zero if redundant!! */
+{
+  long i, j;
+
+  lrs_mp_matrix A = P->A;
+  long *Row = P->Row;
+  long *B = P->B;
+  long d = P->d;
+  long m = P->m;
+
+  if (Q->debug)
+    printA (P, Q);
+
+/* each slack index must be checked for redundancy */
+/* if in cobasis, it is pivoted out if degenerate */
+/* else it is non-redundant                       */
+
+  if (checkcobasic (P, Q, index))
+    return ZERO;
+
+/* index is basic   */
+/* not debugged for new LOC 
+   i=LOC[index];
+ */
+  j = 1;
+  while ((j <= m) && (B[j] != index))
+    j++;
+
+  i = Row[j];
+
+  /* copy row i to cost row, and set it to zero */
+
+  for (j = 0; j <= d; j++)
+    {
+      copy (A[0][j], A[i][j]);
+      changesign (A[0][j]);
+      itomp (ZERO, A[i][j]);
+    }
+
+
+  if (checkredund (P, Q))
+    return ONE;
+
+/* non-redundant, copy back and change sign */
+
+  for (j = 0; j <= d; j++)
+    {
+      copy (A[i][j], A[0][j]);
+      changesign (A[i][j]);
+    }
+
+  return ZERO;
+
+}				/* end of checkindex */
+
+
+/***************************************************************/
+/*                                                             */
+/*            Package of I/O routines                          */
+/*                                                             */
+/***************************************************************/
+
+void
+lprat (const char *name, long Nt, long Dt)
+/*print the long precision rational Nt/Dt without reducing  */
+{
+  if ( Nt > 0 ) 
+    fprintf (lrs_ofp, " ");
+  fprintf (lrs_ofp, "%s%ld", name, Nt);
+  if (Dt != 1)   
+    fprintf (lrs_ofp, "/%ld", Dt);
+  fprintf (lrs_ofp, " ");
+}                               /* lprat */
+
+long
+lreadrat (long *Num, long *Den)
+ /* read a rational string and convert to long    */
+ /* returns true if denominator is not one        */
+{
+  char in[MAXINPUT], num[MAXINPUT], den[MAXINPUT];
+  if(fscanf (lrs_ifp, "%s", in) == EOF)
+     return(FALSE);
+  atoaa (in, num, den);         /*convert rational to num/dem strings */
+  *Num = atol (num);
+  if (den[0] == '\0')
+    {
+      *Den = 1L;
+      return (FALSE);
+    }
+  *Den = atol (den);
+  return (TRUE);
+}
+
+void
+lrs_getinput(lrs_dic *P,lrs_dat *Q,long *num,long *den, long m, long d)
+/* code for reading data matrix in lrs/cdd format */
+{
+  long j,row;
+
+  printf("\nEnter each row: b_i  a_ij j=1..%ld",d);
+  for (row=1;row<=m;row++)
+    {
+      printf("\nEnter row %ld: ",row );
+      for(j=0;j<=d;j++)
+       {
+         lreadrat(&num[j],&den[j]);
+         lprat(" ",num[j],den[j]);
+       }
+
+       lrs_set_row(P,Q,row,num,den,GE);
+     }
+
+ printf("\nEnter objective row c_j j=1..%ld: ",d);
+ num[0]=0; den[0]=1;
+ for(j=1;j<=d;j++)
+       {
+         lreadrat(&num[j],&den[j]);
+         lprat(" ",num[j],den[j]);
+       }
+
+ lrs_set_obj(P,Q,num,den,MAXIMIZE);
+}
+
+
+long 
+readlinearity (lrs_dat * Q)	/* read in and check linearity list */
+{
+  long i, j;
+  long nlinearity;
+  if(fscanf (lrs_ifp, "%ld", &nlinearity)==EOF )
+    {
+      fprintf (lrs_ofp, "\nLinearity option invalid, no indices ");
+      return (FALSE);
+    } 
+  if (nlinearity < 1)
+    {
+      fprintf (lrs_ofp, "\nLinearity option invalid, indices must be positive");
+      return (FALSE);
+    } 
+
+  Q->linearity  = (long int*) CALLOC ((nlinearity + 1), sizeof (long));
+
+  for (i = 0; i < nlinearity; i++)
+    {
+      if(fscanf (lrs_ifp, "%ld", &j)==EOF)
+      {
+      fprintf (lrs_ofp, "\nLinearity option invalid, missing indices");
+      return (FALSE);
+      } 
+      Q->linearity[i] = j;	
+
+    }
+  for (i = 1; i < nlinearity; i++)	/*sort in order */
+    reorder (Q->linearity, nlinearity);
+
+  Q->nlinearity = nlinearity;
+  Q->polytope = FALSE;
+  return TRUE;
+}				/* end readlinearity */
+
+#ifdef PLRS
+void plrs_readlinearity(lrs_dat *Q, string line){
+	istringstream ss(line);
+	long nlinearity;
+	if(!(ss>>nlinearity)){
+		printf("\nLinearity option invalid, no indices\n");
+		exit(1);
+	}
+	if(nlinearity < 1)
+	{
+		printf("\nLinearity option invalid, indices must be positive\n");
+		exit(1);
+	}
+
+	Q->linearity = (long int*) CALLOC ((nlinearity + 1), sizeof (long));
+
+	for (int i = 0; i < nlinearity; i++)
+	{
+		if(!(ss>>Q->linearity[i])){
+			printf("\nLinearity option invalid, missing indices\n");
+			exit(1);
+		}
+	}
+
+	for(int i = 1; i < nlinearity; i++)
+		 reorder (Q->linearity, nlinearity);
+
+ 	Q->nlinearity = nlinearity;
+	Q->polytope = FALSE;
+
+}
+#endif
+
+long 
+readfacets (lrs_dat * Q, long facet[])
+/* read and check facet list for obvious errors during start/restart */
+/* this must be done after linearity option is processed!! */
+{
+  long i, j;
+/* assign local variables to structures */
+  long m, d;
+  long *linearity = Q->linearity;
+  m = Q->m;
+  d = Q->inputd;
+
+  for (j = Q->nlinearity; j < d; j++)	/* note we place these after the linearity indices */
+    {
+      if(fscanf (lrs_ifp, "%ld", &facet[j])==EOF)
+        {
+      fprintf (lrs_ofp, "\nrestart: facet list missing indices");                 
+      return (FALSE);
+      }
+
+
+      fprintf (lrs_ofp, " %ld", facet[j]);
+/* 2010.4.26 nonnegative option needs larger range of indices */
+      if(Q->nonnegative)
+         if (facet[j] < 1 || facet[j] > m+d)
+	  {
+	  fprintf (lrs_ofp, "\n Start/Restart cobasic indices must be in range 1 .. %ld ", m+d);
+	  return FALSE;
+	  }
+      if(!Q->nonnegative)
+         if (facet[j] < 1 || facet[j] > m)
+	  {
+	  fprintf (lrs_ofp, "\n Start/Restart cobasic indices must be in range 1 .. %ld ", m);
+	  return FALSE;
+	  }
+      for (i = 0; i < Q->nlinearity; i++)
+	if (linearity[i] == facet[j])
+	  {
+	    fprintf (lrs_ofp, "\n Start/Restart cobasic indices should not include linearities");
+	    return FALSE;
+	  }
+/* bug fix 2011.8.1  reported by Steven Wu*/
+      for (i = Q->nlinearity; i < j; i++)
+/* end bug fix 2011.8.1 */
+
+	if (facet[i] == facet[j])
+	  {
+	    fprintf (lrs_ofp, "\n  Start/Restart cobasic indices must be distinct");
+	    return FALSE;
+	  }
+    }
+  return TRUE;
+}				/* end of readfacets */
+
+void 
+printA (lrs_dic * P, lrs_dat * Q)	/* print the integer m by n array A 
+					   with B,C,Row,Col vectors         */
+{
+  long i, j;
+/* assign local variables to structures */
+  lrs_mp_matrix A = P->A;
+  long *B = P->B;
+  long *C = P->C;
+  long *Row = P->Row;
+  long *Col = P->Col;
+  long m, d;
+  m = P->m;
+  d = P->d;
+
+  fprintf (lrs_ofp, "\n Basis    ");
+  for (i = 0; i <= m; i++)
+    fprintf (lrs_ofp, "%ld ", B[i]);
+  fprintf (lrs_ofp, " Row ");
+  for (i = 0; i <= m; i++)
+    fprintf (lrs_ofp, "%ld ", Row[i]);
+  fprintf (lrs_ofp, "\n Co-Basis ");
+  for (i = 0; i <= d; i++)
+    fprintf (lrs_ofp, "%ld ", C[i]);
+  fprintf (lrs_ofp, " Column ");
+  for (i = 0; i <= d; i++)
+    fprintf (lrs_ofp, "%ld ", Col[i]);
+  pmp (" det=", P->det);
+  fprintf (lrs_ofp, "\n");
+  i=0;
+  while ( i<= m )
+    {
+      for (j = 0; j <= d; j++)
+	pimat (P, i, j, A[Row[i]][Col[j]], "A");
+      fprintf (lrs_ofp, "\n");
+      if (i==0 && Q->nonnegative)  /* skip basic rows - don't exist! */
+          i=d;
+      i++;
+  fflush (stdout);
+    }
+  fflush (stdout);
+}
+
+
+void 
+pimat (lrs_dic * P, long r, long s, lrs_mp Nt, char name[])
+ /*print the long precision integer in row r col s of matrix A */
+{
+  long *B = P->B;
+  long *C = P->C;
+  if (s == 0)
+    fprintf (lrs_ofp, "%s[%ld][%ld]=", name, B[r], C[s]);
+  else
+    fprintf (lrs_ofp, "[%ld]=", C[s]);
+  pmp ("", Nt);
+
+}
+
+/***************************************************************/
+/*                                                             */
+/*     Routines for caching, allocating etc.                   */
+/*                                                             */
+/***************************************************************/
+
+/* From here mostly Bremner's handiwork */
+
+static void
+cache_dict (lrs_dic ** D_p, lrs_dat * global, long i, long j)
+{
+
+  if (dict_limit > 1)
+    {
+      /* save row, column indicies */
+      (*D_p)->i = i;
+      (*D_p)->j = j;
+
+/* Make a new, blank spot at the end of the queue to copy into */ 
+
+      pushQ (global, (*D_p)->m, (*D_p)->d, (*D_p)->m_A);
+
+      copy_dict (global, global->Qtail, *D_p);	/* Copy current dictionary */
+    }
+  *D_p = global->Qtail;
+
+}
+
+void 
+copy_dict (lrs_dat * global, lrs_dic * dest, lrs_dic * src)
+{
+  long m = src->m;
+  long m_A = src->m_A;        /* number of rows in A */
+  long d = src->d;
+  long r,s;
+
+#ifdef GMP 
+  for ( r=0;r<=m_A;r++)
+    for( s=0;s<=d;s++)
+       copy(dest->A[r][s],src->A[r][s]);
+
+#else            /* fast copy for MP and LRSLONG arithmetic */
+  /* Note that the "A" pointer trees need not be copied, since they
+     always point to the same places within the corresponding space
+*/
+/* I wish I understood the above remark. For the time being, do it the easy way for Nash */
+  if(global->nash)
+  {
+  for ( r=0;r<=m_A;r++)
+    for( s=0;s<=d;s++)
+       copy(dest->A[r][s],src->A[r][s]);
+  }
+  else
+  memcpy (dest->A[0][0], (global->Qtail->prev)->A[0][0],
+          (d + 1) * (lrs_digits + 1) * (m_A + 1) * sizeof (long));
+
+#endif
+
+  dest->i = src->i;
+  dest->j = src->j;
+  dest->m = m;
+  dest->d = d;
+  dest->m_A  = src->m_A;
+
+  dest->depth = src->depth;
+  dest->lexflag = src->lexflag;
+
+  copy (dest->det, src->det);
+  copy (dest->objnum, src->objnum);
+  copy (dest->objden, src->objden);
+
+  if (global->debug)
+    fprintf (lrs_ofp, "\nSaving dict at depth %ld\n", src->depth);
+
+  memcpy (dest->B, src->B, (m + 1) * sizeof (long));
+  memcpy (dest->C, src->C, (d + 1) * sizeof (long));
+  memcpy (dest->Row, src->Row, (m + 1) * sizeof (long));
+  memcpy (dest->Col, src->Col, (d + 1) * sizeof (long));
+}
+
+/* 
+ * pushQ(lrs_dat *globals,m,d):
+ * this routine ensures that Qtail points to a record that 
+ * may be copied into.
+ *
+ * It may create a new record, or it may just move the head pointer
+ * forward so that know that the old record has been overwritten.
+ */
+#if 0
+#define TRACE(s) fprintf(stderr,"\n%s %p %p\n",s,global->Qhead,global->Qtail);
+#else
+#define TRACE(s)
+#endif
+
+static void
+pushQ (lrs_dat * global, long m, long d ,long m_A)
+{
+
+  if ((global->Qtail->next) == global->Qhead)
+    {
+      /* the Queue is full */
+      if (dict_count < dict_limit)
+	{
+	  /* but we are allowed to create more */
+	  lrs_dic *p;
+
+	  p = new_lrs_dic (m, d, m_A);
+
+	  if (p)
+	    {
+
+	      /* we successfully created another record */
+
+	      p->next = global->Qtail->next;
+	      (global->Qtail->next)->prev = p;
+	      (global->Qtail->next) = p;
+	      p->prev = global->Qtail;
+
+	      dict_count++;
+	      global->Qtail = p;
+
+	      TRACE ("Added new record to Q");
+
+	    }
+	  else
+	    {
+	      /* virtual memory exhausted. bummer */
+	      global->Qhead = global->Qhead->next;
+	      global->Qtail = global->Qtail->next;
+
+	      TRACE ("VM exhausted");
+	    }
+	}
+      else
+	{
+	  /*
+	   * user defined limit reached. start overwriting the
+	   * beginning of Q 
+	   */
+	  global->Qhead = global->Qhead->next;
+	  global->Qtail = global->Qtail->next;
+	  TRACE ("User  limit");
+
+	}
+    }
+
+  else
+    {
+      global->Qtail = global->Qtail->next;
+      TRACE ("Reusing");
+    }
+}
+
+lrs_dic *
+lrs_getdic(lrs_dat *Q)
+/* create another dictionary for Q without copying any values */
+/* derived from lrs_alloc_dic,  used by nash.c                */
+{
+lrs_dic *p;
+
+  long m;
+
+  m = Q->m;
+
+/* nonnegative flag set means that problem is d rows "bigger"     */
+/* since nonnegative constraints are not kept explicitly          */
+
+
+  if(Q->nonnegative)
+    m = m+Q->inputd;
+
+  p = new_lrs_dic (m, Q->inputd, Q->m);
+  if (!p)
+    return NULL;
+
+  p->next = p;
+  p->prev = p;
+  Q->Qhead = p;
+  Q->Qtail = p;
+
+  return p;
+}
+
+
+#define NULLRETURN(e) if (!(e)) return NULL;
+
+static lrs_dic *
+new_lrs_dic (long m, long d, long m_A)
+{
+  lrs_dic *p;
+
+  NULLRETURN (p = (lrs_dic *) malloc (sizeof (lrs_dic)));
+
+
+  NULLRETURN (p->B = (long int*) calloc ((m + 1), sizeof (long)));
+  NULLRETURN (p->Row = (long int*) calloc ((m + 1), sizeof (long)));
+
+  NULLRETURN (p->C =  (long int*) calloc ((d + 1), sizeof (long)));
+  NULLRETURN (p->Col = (long int*) calloc ((d + 1), sizeof (long)));
+
+#ifdef GMP
+  lrs_alloc_mp(p->det);
+  lrs_alloc_mp(p->objnum);
+  lrs_alloc_mp(p->objden);
+#endif
+
+  p->d_orig=d;
+  p->A=lrs_alloc_mp_matrix(m_A,d);
+
+
+  return p;
+}
+
+void 
+lrs_free_dic (lrs_dic * P, lrs_dat *Q)
+{
+/* do the same steps as for allocation, but backwards */
+/* gmp variables cannot be cleared using free: use lrs_clear_mp* */
+  lrs_dic *P1;
+
+/* repeat until cache is empty */
+
+  do
+  {
+    /* I moved these here because I'm not certain the cached dictionaries
+       need to be the same size. Well, it doesn't cost anything to be safe. db */
+
+  long d = P->d_orig;
+  long m_A = P->m_A;
+
+  lrs_clear_mp_matrix (P->A,m_A,d);
+
+/* "it is a ghastly error to free something not assigned my malloc" KR167 */
+/* so don't try: free (P->det);                                           */
+
+  lrs_clear_mp (P->det);      
+  lrs_clear_mp (P->objnum);      
+  lrs_clear_mp (P->objden);      
+
+  free (P->Row);
+  free (P->Col);
+  free (P->C);
+  free (P->B);
+
+/* go to next record in cache if any */
+  P1 =P->next;
+  free (P);
+  P=P1;
+
+  }  while (Q->Qhead != P );
+
+
+}
+
+void 
+lrs_free_dic2 (lrs_dic * P, lrs_dat *Q)
+{
+/* do the same steps as for allocation, but backwards */
+/* same as lrs_free_dic except no cache for P */
+    /* I moved these here because I'm not certain the cached dictionaries
+       need to be the same size. Well, it doesn't cost anything to be safe. db */
+
+  long d = P->d_orig;
+  long m_A = P->m_A;
+
+
+  lrs_clear_mp_matrix (P->A,m_A,d);
+
+/* "it is a ghastly error to free something not assigned my malloc" KR167 */
+/* so don't try: free (P->det);                                           */
+
+printf("\n hello 2"); fflush(stdout);
+  lrs_clear_mp (P->det);      
+  lrs_clear_mp (P->objnum);      
+  lrs_clear_mp (P->objden);      
+printf("\n hello 2"); fflush(stdout);
+
+  free (P->Row);
+  free (P->Col);
+  free (P->C);
+  free (P->B);
+
+printf("\n hello 2"); fflush(stdout);
+  free (P);
+
+}
+
+void
+lrs_free_dat ( lrs_dat *Q )
+{
+   long m=Q->m;
+
+/* most of these items were allocated in lrs_alloc_dic */
+
+  lrs_clear_mp_vector (Q->Gcd,m);
+  lrs_clear_mp_vector (Q->Lcm,m);
+
+  lrs_clear_mp (Q->sumdet);
+  lrs_clear_mp (Q->Nvolume);
+  lrs_clear_mp (Q->Dvolume);
+  lrs_clear_mp (Q->saved_det);
+  lrs_clear_mp (Q->boundd);
+  lrs_clear_mp (Q->boundn);
+
+  free (Q->inequality);
+  free (Q->facet);
+  free (Q->redundcol);
+  free (Q->linearity);
+  free (Q->minratio);
+  free (Q->temparray);
+
+  free (Q->name);  
+  free (Q->saved_C);
+
+  lrs_global_count--;
+
+  free(Q);
+}
+
+
+static long
+check_cache (lrs_dic ** D_p, lrs_dat * global, long *i_p, long *j_p)
+{
+/* assign local variables to structures */
+
+
+
+  cache_tries++;
+
+  if (global->Qtail == global->Qhead)
+    {
+      TRACE ("cache miss");
+      /* Q has only one element */
+      cache_misses++;
+      return 0;
+
+    }
+  else
+    {
+      global->Qtail = global->Qtail->prev;
+
+      *D_p = global->Qtail;
+
+      *i_p = global->Qtail->i;
+      *j_p = global->Qtail->j;
+
+      TRACE ("restoring dict");
+      return 1;
+    }
+}
+
+
+lrs_dic *
+lrs_alloc_dic (lrs_dat * Q)
+/* allocate and initialize lrs_dic */
+{
+
+  lrs_dic *p;
+  long i, j;
+  long m, d, m_A;
+
+  if (Q->hull)                       /* d=col dimension of A */
+    Q->inputd = Q->n;                /* extra column for hull */
+  else
+    Q->inputd = Q->n - 1;
+
+  m = Q->m;
+  d = Q->inputd;
+  m_A = m;   /* number of rows in A */
+
+/* nonnegative flag set means that problem is d rows "bigger"     */
+/* since nonnegative constraints are not kept explicitly          */
+
+  if(Q->nonnegative)
+    m = m+d;
+
+  p = new_lrs_dic (m, d, m_A);
+  if (!p)
+    return NULL;
+
+  p->next = p;
+  p->prev = p;
+  Q->Qhead = p;
+  Q->Qtail = p;
+
+
+  dict_count = 1;
+  dict_limit = 50; 
+  cache_tries = 0;
+  cache_misses = 0;
+
+/* Initializations */
+
+  p->d = p->d_orig = d;
+  p->m = m;
+  p->m_A  = m_A;
+  p->depth = 0L;
+  p->lexflag = TRUE;
+  itomp (ONE, p->det);
+  itomp (ZERO, p->objnum);
+  itomp (ONE, p->objden);
+
+/*m+d+1 is the number of variables, labelled 0,1,2,...,m+d  */
+/*  initialize array to zero   */
+  for (i = 0; i <= m_A; i++)
+    for (j = 0; j <= d; j++)
+      itomp (ZERO, p->A[i][j]);
+
+  Q->inequality = (long int*) CALLOC ((m + 1), sizeof (long));
+  if (Q->nlinearity == ZERO)   /* linearity may already be allocated */
+      Q->linearity  = (long int*) CALLOC ((m + 1), sizeof (long));
+
+  Q->facet =  (long int*) CALLOC ((unsigned) d + 1, sizeof (long));
+  Q->redundcol = (long int*) CALLOC ((d + 1), sizeof (long));
+  Q->minratio = (long int*) CALLOC ((m + 1), sizeof (long));
+                         /*  2011.7.14  minratio[m]=0 for degen =1 for nondegen pivot*/
+  Q->temparray = (long int*) CALLOC ((unsigned) d + 1, sizeof (long));
+
+  Q->inequality[0] = 2L;
+  Q->Gcd = lrs_alloc_mp_vector(m);
+  Q->Lcm = lrs_alloc_mp_vector(m);
+  Q->saved_C = (long int*) CALLOC (d + 1, sizeof (long));
+
+  Q->lastdv = d;      /* last decision variable may be decreased */
+                      /* if there are redundant columns          */
+
+/*initialize basis and co-basis indices, and row col locations */
+/*if nonnegative, we label differently to avoid initial pivots */
+/* set basic indices and rows */
+ if(Q->nonnegative)
+  for (i = 0; i <= m; i++)
+    {
+      p->B[i] = i;
+      if (i <= d )
+          p->Row[i]=0; /* no row for decision variables */
+      else 
+          p->Row[i]=i-d;
+    }
+ else
+   for (i = 0; i <= m; i++)
+    {
+      if (i == 0 )
+          p->B[0]=0;
+      else
+          p->B[i] = d + i;
+      p->Row[i] = i;
+    }
+
+  for (j = 0; j < d; j++)
+    {
+      if(Q->nonnegative)
+          p->C[j] = m+j+1;
+      else
+          p->C[j] = j + 1;
+      p->Col[j] = j + 1;
+    }
+  p->C[d] = m + d + 1;
+  p->Col[d] = 0;
+  return p;
+}				/* end of lrs_alloc_dic */
+
+
+/* 
+   this routine makes a copy of the information needed to restart, 
+   so that we can guarantee that if a signal is received, we 
+   can guarantee that nobody is messing with it.
+   This as opposed to adding all kinds of critical regions in 
+   the main line code.
+
+   It is also used to make sure that in case of overflow, we
+   have a valid cobasis to restart from.
+ */
+static void
+save_basis (lrs_dic * P, lrs_dat * Q)
+{
+  int i;
+/* assign local variables to structures */
+  long *C = P->C;
+  long d;
+
+#ifdef SIGNALS
+  sigset_t oset, blockset;
+  sigemptyset (&blockset);
+  sigaddset (&blockset, SIGTERM);
+  sigaddset (&blockset, SIGHUP);
+  sigaddset (&blockset, SIGUSR1);
+
+  errcheck ("sigprocmask", sigprocmask (SIG_BLOCK, &blockset, &oset));
+#endif
+  d = P->d;
+
+  Q->saved_flag = 1;
+
+  for (i = 0; i < 3; i++)
+    Q->saved_count[i] = Q->count[i];
+
+  for (i = 0; i < d + 1; i++)
+    Q->saved_C[i] = C[i];
+
+  copy (Q->saved_det, P->det);
+
+  Q->saved_d = P->d;
+  Q->saved_depth = P->depth;
+
+#ifdef SIGNALS
+  errcheck ("sigprocmask", sigprocmask (SIG_SETMASK, &oset, 0));
+#endif
+}
+
+/* digits overflow is a call from lrs_mp package */
+
+void 
+digits_overflow ()
+{
+  fprintf (lrs_ofp, "\nOverflow at digits=%ld", DIG2DEC (lrs_digits));
+  fprintf (lrs_ofp, "\nRerun with option: digits n, where n > %ld\n", DIG2DEC (lrs_digits));
+  lrs_dump_state ();
+
+  notimpl("");
+}
+
+static void 
+lrs_dump_state ()
+{
+  long i;
+
+  fprintf (stderr, "\n\nlrs_lib: checkpointing:\n");
+
+  fprintf (stderr, "lrs_lib: Current digits at %ld out of %ld\n",
+	   DIG2DEC (lrs_record_digits),
+	   DIG2DEC (lrs_digits));
+
+  for (i = 0; i < lrs_global_count; i++)
+    {
+      print_basis (stderr, lrs_global_list[i]);
+    }
+  fprintf (stderr, "lrs_lib: checkpoint finished\n");
+}
+
+
+/* print out the saved copy of the basis */
+void 
+print_basis (FILE * fp, lrs_dat * global)
+{
+  int i;
+/* assign local variables to structures */
+  fprintf (fp, "lrs_lib: State #%ld: (%s)\t", global->id, global->name);
+
+  if (global->saved_flag)
+    {
+
+      fprintf (fp, "V#%ld R#%ld B#%ld h=%ld facets ",
+	       global->saved_count[1],
+	       global->saved_count[0],
+	       global->saved_count[2],
+	       global->saved_depth);
+      for (i = 0; i < global->saved_d; i++)
+	fprintf (fp, "%ld ",
+		 global->inequality[global->saved_C[i] - global->lastdv]);
+      pmp (" det=", global->saved_det);
+      fprintf (fp, "\n");
+
+    }
+  else
+    {
+      fprintf (fp, "lrs_lib: Computing initial basis\n");
+    }
+
+
+  fflush (fp);
+}
+
+#ifdef SIGNALS
+
+/*
+   If given a signal
+   USR1            print current cobasis and continue
+   TERM            print current cobasis and terminate
+   INT (ctrl-C) ditto
+   HUP                     ditto
+ */
+static void
+setup_signals ()
+{
+  errcheck ("signal", signal (SIGTERM, die_gracefully));
+  errcheck ("signal", signal (SIGALRM, timecheck));
+  errcheck ("signal", signal (SIGHUP, die_gracefully));
+  errcheck ("signal", signal (SIGINT, die_gracefully));
+  errcheck ("signal", signal (SIGUSR1, checkpoint));
+}
+
+static void
+timecheck ()
+{
+  lrs_dump_state ();
+  errcheck ("signal", signal (SIGALRM, timecheck));
+  alarm (lrs_checkpoint_seconds);
+}
+
+static void
+checkpoint ()
+{
+  lrs_dump_state ();
+  errcheck ("signal", signal (SIGUSR1, checkpoint));
+}
+
+static void
+die_gracefully ()
+{
+  lrs_dump_state ();
+
+  exit (1);
+}
+
+#endif
+
+#ifdef TIMES
+/* 
+ * Not sure about the portability of this yet, 
+ *              - db
+ */
+#include <sys/resource.h>
+#define double_time(t) ((double)(t.tv_sec)+(double)(t.tv_usec)/1000000)
+
+static void
+ptimes ()
+{
+  struct rusage rusage;
+  getrusage (RUSAGE_SELF, &rusage);
+  fprintf (lrs_ofp, "\n*%0.3fu %0.3fs %ldKb %ld flts %ld swaps %ld blks-in %ld blks-out \n",
+	   double_time (rusage.ru_utime),
+	   double_time (rusage.ru_stime),
+	   rusage.ru_maxrss, rusage.ru_majflt, rusage.ru_nswap,
+	   rusage.ru_inblock, rusage.ru_oublock);
+  if(lrs_ofp != stdout)
+     printf ("\n*%0.3fu %0.3fs %ldKb %ld flts %ld swaps %ld blks-in %ld blks-out \n",
+	   double_time (rusage.ru_utime),
+	   double_time (rusage.ru_stime),
+	   rusage.ru_maxrss, rusage.ru_majflt, rusage.ru_nswap,
+	   rusage.ru_inblock, rusage.ru_oublock);
+}
+
+static double get_time()
+{
+  struct rusage rusage;
+  getrusage (RUSAGE_SELF, &rusage);
+  return   ( double_time (rusage.ru_utime));
+
+}
+
+#endif
+
+/* Routines based on lp_solve */
+
+void
+lrs_set_row(lrs_dic *P, lrs_dat *Q, long row, long num[], long den[], long ineq)
+/* convert to lrs_mp then call lrs_set_row */
+{
+ lrs_mp_vector Num, Den;
+ long d;
+ long j;
+ 
+  d = P->d;
+
+  Num=lrs_alloc_mp_vector(d+1);
+  Den=lrs_alloc_mp_vector(d+1);
+
+  for (j=0;j<=d;j++)
+  {
+  itomp(num[j],Num[j]);
+  itomp(den[j],Den[j]);
+  }
+ 
+  lrs_set_row_mp(P,Q,row,Num,Den,ineq);
+
+  lrs_clear_mp_vector(Num,d+1);
+  lrs_clear_mp_vector(Den,d+1);
+
+}
+
+void
+lrs_set_row_mp(lrs_dic *P, lrs_dat *Q, long row, lrs_mp_vector num, lrs_mp_vector den, long ineq)
+/* set row of dictionary using num and den arrays for rational input */
+/* ineq = 1 (GE)   - ordinary row  */
+/*      = 0 (EQ)   - linearity     */
+{
+  lrs_mp Temp, mpone;
+  lrs_mp_vector oD;             /* denominator for row  */
+
+  long i, j;
+
+/* assign local variables to structures */
+
+  lrs_mp_matrix A;
+  lrs_mp_vector Gcd, Lcm;
+  long hull;
+  long m, d;
+  lrs_alloc_mp(Temp); lrs_alloc_mp(mpone);
+  hull = Q->hull;
+  A = P->A;
+  m = P->m;
+  d = P->d;
+  Gcd = Q->Gcd;
+  Lcm = Q->Lcm;
+
+  oD = lrs_alloc_mp_vector (d);
+  itomp (ONE, mpone);
+  itomp (ONE, oD[0]);
+
+  i=row;
+  itomp (ONE, Lcm[i]);      /* Lcm of denominators */
+  itomp (ZERO, Gcd[i]);     /* Gcd of numerators */
+  for (j = hull; j <= d; j++)       /* hull data copied to cols 1..d */
+        {
+          copy( A[i][j],num[j-hull]);
+          copy(oD[j],den[j-hull]);
+          if (!one(oD[j]))
+            lcm (Lcm[i], oD[j]);      /* update lcm of denominators */
+          copy (Temp, A[i][j]);
+          gcd (Gcd[i], Temp);   /* update gcd of numerators   */
+        }
+
+  if (hull)
+        {
+          itomp (ZERO, A[i][0]);        /*for hull, we have to append an extra column of zeroes */
+          if (!one (A[i][1]) || !one (oD[1]))         /* all rows must have a one in column one */
+            Q->polytope = FALSE;
+        }
+  if (!zero (A[i][hull]))   /* for H-rep, are zero in column 0     */
+        Q->homogeneous = FALSE; /* for V-rep, all zero in column 1     */
+
+  storesign (Gcd[i], POS);
+  storesign (Lcm[i], POS);
+  if (mp_greater (Gcd[i], mpone) || mp_greater (Lcm[i], mpone))
+        for (j = 0; j <= d; j++)
+          {
+            exactdivint (A[i][j], Gcd[i], Temp);        /*reduce numerators by Gcd  */
+            mulint (Lcm[i], Temp, Temp);        /*remove denominators */
+            exactdivint (Temp, oD[j], A[i][j]);       /*reduce by former denominator */
+          }
+
+  if ( ineq == EQ )        /* input is linearity */
+     {
+      Q->linearity[Q->nlinearity]=row;
+      Q->nlinearity++;
+     }
+
+/* 2010.4.26   Set Gcd and Lcm for the non-existant rows when nonnegative set */
+
+
+  if(Q->nonnegative && row==m)
+      for(j=1;j<=d;j++)
+         { itomp (ONE, Lcm[m+j]);
+           itomp (ONE, Gcd[m+j]);
+         }
+
+
+  lrs_clear_mp_vector (oD,d);
+  lrs_clear_mp(Temp); lrs_clear_mp(mpone);
+}          /* end of lrs_set_row_mp */
+
+void 
+lrs_set_obj(lrs_dic *P, lrs_dat *Q, long num[], long den[], long max)
+{
+  long i;
+
+  if (max == MAXIMIZE)
+       Q->maximize=TRUE;
+  else
+       {
+       Q->minimize=TRUE;
+       for(i=0;i<=P->d;i++)
+         num[i]=-num[i];
+       }
+
+  lrs_set_row(P,Q,0L,num,den,GE);
+}
+
+void
+lrs_set_obj_mp(lrs_dic *P, lrs_dat *Q, lrs_mp_vector num, lrs_mp_vector den, long max)
+{
+  long i;
+
+  if (max == MAXIMIZE)
+       Q->maximize=TRUE;
+  else
+       {
+       Q->minimize=TRUE;
+       for(i=0;i<=P->d;i++)
+         changesign(num[i]);
+       }
+
+  lrs_set_row_mp(P,Q,0L,num,den,GE);
+}
+
+
+long
+lrs_solve_lp(lrs_dic *P, lrs_dat *Q)
+/* user callable function to solve lp only */
+{
+  lrs_mp_matrix Lin;		/* holds input linearities if any are found             */
+  long col;
+
+  Q->lponly = TRUE;
+
+  if (!lrs_getfirstbasis (&P, Q, &Lin, FALSE))
+    return FALSE;
+
+/* There may have been column redundancy                */
+/* If so the linearity space is obtained and redundant  */
+/* columns are removed. User can access linearity space */
+/* from lrs_mp_matrix Lin dimensions nredundcol x d+1   */
+
+  for (col = 0; col < Q->nredundcol; col++)	/* print linearity space               */
+    lrs_printoutput (Q, Lin[col]);   	        /* Array Lin[][] holds the coeffs.     */
+
+  return TRUE;
+} /* end of lrs_solve_lp */
+
+
+long 
+dan_selectpivot (lrs_dic * P, lrs_dat * Q, long *r, long *s)
+/* select pivot indices using dantzig simplex method             */
+/* largest coefficient with lexicographic rule to avoid cycling  */
+/* Bohdan Kaluzny's handiwork                                    */
+/* returns TRUE if pivot found else FALSE                        */
+/* pivot variables are B[*r] C[*s] in locations Row[*r] Col[*s]  */
+{
+  long j,k,col;
+  lrs_mp coeff;
+/* assign local variables to structures */
+  lrs_mp_matrix A = P->A;
+  long *Col = P->Col;
+  long d = P->d;
+
+  lrs_alloc_mp (coeff);
+  *r = 0;
+  *s = d;
+  j = 0;
+  k = 0;
+ 
+  itomp(0,coeff);
+/*find positive cost coef */
+  while (k < d) 
+     {
+       if(mp_greater(A[0][Col[k]],coeff))
+        {
+          j = k;
+          copy(coeff,A[0][Col[j]]);
+	}
+      k++;
+     }
+
+  if (positive(coeff))			/* pivot column found! */
+    {
+      *s = j;
+      col = Col[j];
+
+      /*find min index ratio */
+      *r = lrs_ratio (P, Q, col);
+      if (*r != 0)
+        {
+        lrs_clear_mp(coeff);
+	return (TRUE);		/* unbounded */
+        }
+    }
+  lrs_clear_mp(coeff);
+  return (FALSE);
+}				/* end of dan_selectpivot        */
+
+
+long
+phaseone (lrs_dic * P, lrs_dat * Q)
+/* Do a dual pivot to get primal feasibility (pivot in X_0)*/
+/* Bohdan Kaluzny's handiwork                                    */
+{
+  long i, j, k;
+/* assign local variables to structures */
+  lrs_mp_matrix A = P->A;
+  long *Row = P->Row;
+  long *Col = P->Col;
+  long m, d;
+  lrs_mp b_vector;
+  lrs_alloc_mp (b_vector);
+  m = P->m;
+  d = P->d;
+  i = 0;
+  k = d+1;
+
+  itomp(0,b_vector);
+
+  fprintf (lrs_ofp, "\nLP: Phase One: Dual pivot on artificial variable");
+
+/*find most negative b vector */
+  while (k <= m)
+     {
+       if(mp_greater(b_vector,A[Row[k]][0]))
+        {
+          i = k;
+          copy(b_vector,A[Row[i]][0]);
+        }
+      k++;
+     }
+
+  if (negative(b_vector))                       /* pivot row found! */
+    {
+      j = 0;            /*find a positive entry for in row */
+      while (j < d && !positive (A[Row[i]][Col[j]]))
+        j++;
+      if (j >= d)
+        {
+          lrs_clear_mp (b_vector);
+          return (FALSE);       /* no positive entry */
+        }
+      pivot (P, Q, i, j);
+      update (P, Q, &i, &j);
+    }
+  lrs_clear_mp (b_vector);
+  return (TRUE);
+}
+
+
+long
+lrs_set_digits(long dec_digits)
+{
+/* convert user specified decimal digits to mp digits */
+
+  fprintf (lrs_ofp, "\n*digits %ld", dec_digits);
+  if (dec_digits > 0)
+    lrs_digits = DEC2DIG (dec_digits);
+  if (lrs_digits > MAX_DIGITS)
+    {
+      fprintf (lrs_ofp, "\nDigits must be at most %ld\nChange MAX_DIGITS and recompile",
+	       DIG2DEC (MAX_DIGITS));
+      fflush(stdout);
+      return (FALSE);
+    }
+  return (TRUE);
+}
+
+long
+lrs_checkbound(lrs_dic *P, lrs_dat *Q)
+{
+/* check bound on objective and return TRUE if exceeded */
+
+  if(!Q->bound)
+     return FALSE;
+
+  if( Q->maximize && comprod(Q->boundn,P->objden,P->objnum,Q->boundd) == 1 )
+       {
+	#ifndef PLRS
+        if(Q->verbose)
+             {
+              prat(" \nObj value: ",P->objnum,P->objden);
+              fprintf(lrs_ofp," Pruning ");
+              }
+	#endif
+         return TRUE;
+       }
+  if( Q->minimize && comprod(Q->boundn,P->objden,P->objnum,Q->boundd) == -1 )
+       {
+	#ifndef PLRS
+        if(Q->verbose)
+             {
+              prat(" \nObj value: ",P->objnum,P->objden);
+              fprintf(lrs_ofp," Pruning ");
+              }
+	#endif
+         return TRUE;
+       }
+  return FALSE;
+}
+
+
+long
+lrs_leaf(lrs_dic *P, lrs_dat *Q)
+{
+/* check if current dictionary is a leaf of reverse search tree */
+  long    col=0;
+  long    tmp=0;
+
+  while (col < P->d && !reverse(P,Q,&tmp,col))
+                 col++;
+  if(col < P->d) 
+     return 0;            /* dictionary is not a leaf */
+  else
+     return 1;
+}