/* polyv.c
 *
 * Produces logical formulas in SMT-LIB standard 2.6, looking for
 * counter-examples to show two things are distinct.
 *
 * TODO: refactor, merging at least run_hhcheck, run_hvcheck, run_vvcheck
 *       possibly can merge bodies/etc too
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>

#define VERSION "20220316"

#ifdef LRSLONG
#define ARITH "lrslong.h"    /* lrs long integer arithmetic package */
#else
#if defined(GMP) || defined(FLINT)
#define ARITH "lrsgmp.h"     /* lrs wrapper for gmp multiple precision arithmetic    */
#else
#define ARITH "lrsmp.h"      /* lrs multiple precision arithmetic    */
#endif
#endif

#include ARITH

#define printusage() \
  printf("Usage: %s <file.ine>\n" \
	 "To check redundancy of the first inequality in redund_list or redund when eliminating/projecting as specified by the eliminate/project option.\n\n%s -c <1/2> <file.ine>\nTo generate lrs input file <1/2> checking a witness (given on stdin) of non-redundancy generated by z3\n\n%s <file1> <file2>\nTo check whether H/V-representations file1 and file2 define different polyhedra, where H-representations may contain projections.\n\n%s <file1> <file2> <file3>\nTo check whether the intersection of H/V-representations file1 and file2 is different from H/V-representation file3, where H-representations may contain projections.", \
	 argv[0],argv[0],argv[0],argv[0]);

int debug = 0;
FILE *lrs_ifp; /* hook to make readrat work */
FILE *lrs_ofp; /* hook to be able to link to lrsgmp/etc */
int mode = 0; /* 0: normal, 1: generate certificate 1, 2: generate certificate 2 */
char *filename1=NULL;
char *filename2=NULL;
char *filename3=NULL;

int force[3]={0}; /* 1: force the file to V, 2: force the file to H */

struct mat {
	lrs_mp_matrix num, den;
	unsigned int m, n;	/* m inequalities, n variables x1...xn */
	unsigned int nlinearity; /* number of linearities */
	unsigned int *linearities; /* indices here are equations */
	unsigned int *linmask; /* Boolean flags for equations */
	unsigned int *varmask; /* Boolean flags for variables to eliminate */
	unsigned int vrep; /* 0: h-rep, 1: v-rep */
	unsigned int nrays; /* number of rays in v-rep */
	unsigned int *rays; /* indices here are rays in v-rep */
	unsigned int *raymask; /* Boolean flags for rays */
	unsigned long *vars; /* variables projecting out */
	unsigned int nv;     /* number of variables projecting out */
	unsigned int elim;   /* eliminate/project option present? */
	unsigned int redund; /* redund/redund_list option present? */
	};

struct mat *parse_input(FILE *, unsigned long *, int);
void print_formula(struct mat*, unsigned int, unsigned long *, unsigned long, const char *, const int, char **);
void free_mat(struct mat*);
void print_cert(struct mat*, unsigned int, unsigned long *, unsigned long, const char *, const int, char **);
int get_options(int, char **, FILE **, FILE **, FILE **);
int run_checkpred(struct mat *, unsigned long, int, char **);
int run_hvcheck(struct mat *, struct mat *, int, char **);
int run_hhcheck(struct mat *, struct mat *, int, char **);
int run_vvcheck(struct mat *, struct mat *, int, char **);
int run_intcheck(struct mat*, struct mat *, struct mat *, int, char **);
unsigned int *make_newn(struct mat *);
void printsystem(struct mat *, unsigned long, unsigned int *);
int projvar(struct mat *, unsigned int);
void setup_varmask(struct mat *);
int sanity_check(struct mat *, struct mat *, struct mat *);
int setup_mat(FILE *, struct mat **, unsigned long *, int);

int main(int argc, char **argv)
{
	FILE *f1=NULL, *f2=NULL, *f3=NULL;
	struct mat *mat=NULL, *mat2=NULL, *mat3=NULL;
	unsigned long ineq=0; /* ineq ineq when projecting out vars */

	lrs_ofp = stdout;

	if (get_options(argc, argv, &f1, &f2, &f3)==0)
	{
		printusage();
		if (f1!=NULL)
			fclose(f1);
		if (f2!=NULL)
			fclose(f2);
		if (f3!=NULL)
			fclose(f3);
		return 0;
	}

	if (setup_mat(f1, &mat,&ineq,1)==0)
		return 0;
	if (setup_mat(f2, &mat2,&ineq,2)==0)
	{
		free_mat(mat);
		return 0;
	}
	if (setup_mat(f3, &mat3,&ineq,3)==0)
	{
		free_mat(mat);
		free_mat(mat2);
		return 0;
	}

	if (!sanity_check(mat,mat2,mat3))
		return 0;

	if (mode>=0 && mode<=2)
		return run_checkpred(mat, ineq, argc, argv);

	else if (mode==4)
		return run_hvcheck(mat, mat2, argc, argv);
	else if (mode==5)
		return run_hhcheck(mat, mat2, argc, argv);
	else if (mode==6)
		return run_vvcheck(mat, mat2, argc, argv);
	else if (mode==7)
		return run_intcheck(mat, mat2, mat3, argc, argv);
	else
		return 0;
}

/* parse the file into the mat, returning 0 if a problem
 * close file when done
 */
int setup_mat(FILE *f, struct mat **mat, unsigned long *ineq, int num)
{
	if (f!=NULL)
	{
		*mat = parse_input(f, ineq, num);
		fclose(f);
		if (*mat == NULL)
		{
			fprintf(stderr,"Error parsing input\n");
			return 0;
		}
		setup_varmask(*mat);
		return 1;
	}
	return 1;
}

/* We have inputs mat1, mat2 (possibly NULL)
 * check them, set mode:
 *  0   checkpred mode, check single inequality for projected redundancy
 *  1/2 checkpred certificate mode (mode already set in this case
 *  3   fel verification mode (removed, do as special case of mode 5)
 *  4   H/V verification mode
 *  5   H/H verification mode
 *  6   V/V verification mode
 *  7   mat1 intersect mat2 == mat3 mode
 *
 * return 0: abort run, we cleaned up mat1,mat2,mat3
 * return 1: okay to go, mode is set
 */

int sanity_check(struct mat *mat1, struct mat *mat2, struct mat *mat3)
{
	const char *msg;
	const char *emsg=NULL;
	if (mode==1 || mode==2)
	{
		if (mat3!=NULL || mat2!=NULL || mat1==NULL)
		{
			fprintf(stderr, "certificate generation takes one input\n");
			free_mat(mat1);
			free_mat(mat2);
			free_mat(mat3);
			return 0;
		}
	}
	else if (mat1==NULL)
	{
		fprintf(stderr, "no input file given\n");
		free_mat(mat1);
		free_mat(mat2);
		return 0;
	}
	else if (mat1!=NULL && mat2==NULL && mat3==NULL)
	{
		if (mat1->vrep==0) /* one h rep, mode 0 */
			mode=0;
		else /* single v rep, not supported */
		{
			fprintf(stderr, "no operation for single V-representation\n");
			free_mat(mat1);
			return 0;
		}
	}
	else if (mat1!=NULL && mat2!=NULL && mat3==NULL) /* two inputs */
	{
		if (mat1->vrep==0 && mat2->vrep==0) /* HH */
			mode=5;
		else if (mat1->vrep==1 && mat2->vrep==1) /* VV */
			mode=6;
		else /* HV or VH */
			mode=4;
	}
	else /* three inputs */
		mode = 7;

	if (mode == 0)
		msg="*polv: checkpred mode";
	else if (mode==1 || mode==2)
		msg="*polyv: checkpred certificate mode";
	else if (mode == 4)
		msg="*polyv: H/V verification mode";
	else if (mode == 5)
		msg="*polyv: H/H verification mode";
	else if (mode == 6)
		msg="*polyv: V/V verification mode";
	else if (mode == 7)
		msg="*polyv: intersection verification mode";
	else
		msg="*polyv: very confused mode, good luck!";

	fprintf(stderr, "%s\n", msg);

	if (mode<=2 && mat1->redund==0)
		emsg="redund/redund_list line missing";
	if ( mode<=2 && mat1->elim==0 )
		emsg="eliminate/project line missing";

	if (emsg!=NULL)
	{
		fprintf(stderr, "error: %s\n", emsg);
		free_mat(mat1);
		free_mat(mat2);
		free_mat(mat3);
		return 0;
	}

	return 1;

}

void print_header(int argc, char **argv, const char *what)
{
	unsigned int i;
	printf("; polyv %s formula for verifying %s\n",
		VERSION, what);
	printf("; invocation was:");
	for (i=0; i<argc; i++)
		printf(" %s", argv[i]);
	printf("\n");
	printf("(set-logic LRA)\n\n");
	return;
}

void setup_varmask(struct mat *mat)
{
	unsigned int i;
	unsigned long *vars=mat->vars;
	mat->varmask = calloc(mat->n+1, sizeof(unsigned int));
	for (i=0; i<mat->nv; i++)
		mat->varmask[vars[i]] = 1;
}

int run_checkpred(struct mat *mat, unsigned long ineq, int argc, char **argv)
{
	unsigned long *vars = mat->vars;
	unsigned int i, flag=0, nv=mat->nv;
	for (i=0; i<nv; i++)
		if (vars[i]<1 || vars[i] >= mat->n)
			flag = 1;
	if (flag || ineq<1 || ineq>mat->m)
	{
		fprintf(stderr,"Bounds error: variables 1...%d, inequalities 1...%d\n", mat->n, mat->m);
		free_mat(mat);
		return 0;
	}

	fprintf(stderr, "*Mode %d checking non-redundancy of inequality %lu after eliminating variables", mode, ineq);
	for (i=0; i<nv; i++)
		fprintf(stderr, " %lu", vars[i]);
	fprintf(stderr, " of %s\n", filename1);

	if (mode == 0)
		print_formula(mat, nv, vars, ineq, filename1, argc, argv);
	else /* print a certificate */
		print_cert(mat, nv, vars, ineq, filename1, argc, argv);
	free_mat(mat);

	return 0;
}

int get_options(int argc, char **argv, FILE **f1, FILE **f2, FILE **f3)
{
	int i, j, len, num=0;

	if (argc < 2 || argc > 4)
		return 0;

	for (i=1; i<argc; i++)
	{
		if (!strcmp(argv[i], "-c"))
		{
			if (i+1==argc)
				return 0;
			mode = atoi(argv[i+1]);
			if (mode<0 || mode>2)
				return 0;
			i++; /* skip argument too */
		}
		else if (argv[i][0]=='-')
		{
			len = strlen(argv[i]);
			if (len<2 || len>4)
			{
				fprintf(stderr, "unknown option %s\n",
					argv[i]);
				return 0;
			}
			for (j=1; j<len; j++)
			{
				if (argv[i][j]=='v')
					force[j-1]=1;
				else if (argv[i][j]=='h')
					force[j-1]=2;
				else
				{
					fprintf(stderr, "unknown option %s\n",
						argv[i]);
					return 0;
				}
			}	
		}
		/* at end */
		else if (num==2)
		{
			*f3 = fopen(argv[i], "r");
			if (*f3==NULL)
				return 0;
			filename3=argv[i];
			num++;
		}
		else if (num==1)
		{
			*f2 = fopen(argv[i], "r");
			if (*f2==NULL)
				return 0;
			filename2=argv[i];
			num++;
		}
		else if (num==0)
		{
			*f1 = fopen(argv[i], "r");
			if (*f1==NULL)
				return 0;
			filename1=argv[i];
			num++;
		}
		else
			return 0;
	}

	
	return num;
}

void parse_error(struct mat *mat, const char *s)
{
	fprintf(stderr, "Parse error: %s\n", s);
	if (mat == NULL)
		return;
	if (mat->num != NULL)
		lrs_clear_mp_matrix(mat->num,mat->m,mat->n);
	if (mat->den != NULL)
		lrs_clear_mp_matrix(mat->den,mat->m,mat->n);
	free(mat);
	exit(0);
}

void eatspace(FILE *f)
{
	int c;
	while (isspace(c=getc(f)));
	ungetc(c,f);
}

/* just read name (either "redund" or "redund_list"
 * handle the line, put the first ineq in ineq
 */
void handle_redund(FILE *f, char *name, struct mat *mat, unsigned long *ineq)
{
	unsigned int i, j, tmp;
	int ret;
	eatspace(f);
	if (!strcmp("redund", name))
	{
		*ineq = 1;
		return;
	}
	/* redund_list */
	ret = fscanf(f, "%u", &j);
	if (ret!=1 || j==0)
		parse_error(mat, "broken redund_list line");
	for (i=0; i<j; i++)
	{
		eatspace(f);
		ret = fscanf(f, "%u", &tmp);
		if (ret!=1)
			parse_error(mat, "broken redund_list line");
		if (i==0)
			*ineq = tmp;
	}
	eatspace(f);
	return;
}

/* just read name (either "project" or "eliminate"
 * handle the line, set nv, allocate and set vars
 * note: could be multiple lines, last over-writes earlier ones
 * so realloc as needed
 * TODO: redo comment, put vars and nv in mat instead
 */
void handle_project(FILE *f, char *name, struct mat *mat, unsigned int* nv,
		    unsigned long **vars)
{
	unsigned int i, j, num;
	unsigned long *tmp;
	unsigned long *tmp2;
	int ret;

	eatspace(f);
	ret = fscanf(f, "%u", &num);
	if (ret!=1 || num==0)
		parse_error(mat, "broken project/eliminate line");
	tmp = malloc(sizeof(unsigned long) * num);
	for (i=0; i<num; i++)
	{
		eatspace(f);
		ret = fscanf(f, "%lu", tmp+i);
		if (ret!=1 || tmp[i]<1 || tmp[i]>mat->m)
			parse_error(mat, "invalid row index in project/eliminate");
	}
	eatspace(f);
	if (!strcmp(name, "eliminate"))
	{
		*nv = num;
		tmp2 = realloc(*vars, sizeof(unsigned long)*num);
		if (tmp2 == NULL)
			parse_error(mat, "memory allocation error");
		*vars = tmp2;
		for (i=0; i<num; i++)
			tmp2[i] = tmp[i];
		free(tmp);
		return;
	}
	/* project, so flip indices */
	tmp2 = calloc(mat->n, sizeof(unsigned long));
	for (i=0; i<num; i++)
		tmp2[tmp[i]] = 1;
	/* tmp2[i] is 1 if i is to be kept, 0 if to be eliminated */
	free(tmp);
	*nv = mat->n - 1 - num;
	tmp = realloc(*vars, sizeof(unsigned long)* *nv);
	if (tmp == NULL)
		parse_error(mat, "memory allocation error");
	*vars = tmp;
	j=0;
	for (i=1; i<mat->n; i++)
		if (tmp2[i] == 0)
			tmp[j++] = i;
	if (j!=*nv)
		parse_error(mat, "horribly broken");
	free(tmp2);	
	return;
}

struct mat *parse_input(FILE *f, unsigned long *ineq, int num)
{
	struct mat *mat;
	char *p, name[4096]={0}; /* lrslib uses 1000 */
	unsigned int cur, i, j, k, offs=0;
	unsigned int nrays = 0, *nv;
	int ret, end=0;
	unsigned long **vars;

	mat = calloc(1, sizeof(struct mat));	
	vars=&mat->vars;
	nv=&mat->nv;

	p = fgets(name, 4000, f);
	if (p == NULL)
		parse_error(mat, "file error");
	if (debug)
		fprintf(stderr, "%s", name);
	while (strncmp(name,"begin",5))
	{
		if (!strncmp(name, "V-representation", 16))
			mat->vrep=1;
		if (!strncmp(name, "linearity ", 10))
		{
			/* read linearity option */
			p = name + 10;
			ret = sscanf(name+10, "%u%n", &k, &cur);
			if (ret == 0)
				parse_error(mat,"linearity error");
			mat->nlinearity = k;
			mat->linearities = calloc(k, sizeof(unsigned int));
			for (offs=cur, i=0; i<k; i++)
			{
				ret = sscanf(name+10+offs, "%u%n", 
					     &mat->linearities[i], &cur);
				if (ret == 0)
					parse_error(mat,"linearity error");
				offs += cur;
			}
			if (debug)
			{
				fprintf(stderr, "linearities:");
				for (i=0; i<mat->nlinearity; i++)
					fprintf(stderr, " %u", 
						mat->linearities[i]);
				fprintf(stderr, "\n");
			}
		} /* linearity option */
		fgets(name, 4000, f);
	}
	/* got begin */
	ret = fscanf(f, "%u %u", &mat->m, &mat->n);
	if (ret != 2)
		parse_error(mat,"missing or broken dimensions");
	mat->linmask = calloc(mat->m+1, sizeof(unsigned int));

	/* support for -hh, -hvh, etc over-riding input files */
	if (force[num-1] == 1 && !mat->vrep) /* force file to V rep */
	{
		fprintf(stderr, "*polyv: treating H-representation as V-representation\n");
		mat->vrep = 1;
	}
	else if (force[num-2] == 2 && mat->vrep) /* force file to H rep */
	{
		fprintf(stderr, "*polyvv treating V-representation as H-representation\n");
		mat->vrep = 0;
	}

	if (mat->vrep)
		mat->raymask = calloc(mat->m+1, sizeof(unsigned int));
	for (i=0; i<mat->nlinearity; i++)
		mat->linmask[mat->linearities[i]] = 1;
	fgets(name, 4000, f);
	for (i=0; isspace(name[i]); i++);
	if (strncmp(name+i, "integer", 7)!=0 && strncmp(name+i, "rational", 8)!=0)
		parse_error(mat,"data type must be integer or rational");

	mat->num = lrs_alloc_mp_matrix(mat->m, mat->n);
	mat->den = lrs_alloc_mp_matrix(mat->m, mat->n);

	lrs_ifp = f; /* hook to make readrat work */

	/* read the matrix */
	for (i=1; i<=mat->m; i++)
		for (j=0; j<mat->n; j++)
		{
			readrat(mat->num[i][j], mat->den[i][j]);
			if (j==0 && mat->vrep && zero(mat->num[i][j]) && !mat->linmask[i]) /* a ray */
			{
				nrays++;
				mat->raymask[i]=1;
			}
		}

	/* read options - required project or eliminate and redund or
	 * redund_list after end. they could be missing or too many
	 * numbers
	 */
	ret = 0;
	
	while (ret != EOF)
	{
		eatspace(f);
		ret = fscanf(f, "%s", name);
		if (ret != 1)
			break;
		else if (!strcmp("end", name))
		{
			end++;
			continue;
		}
		else if (!strcmp("redund_list", name)|| !strcmp("redund", name))
			mat->redund++,handle_redund(f, name, mat, ineq);
		else if (!strcmp("eliminate", name) || !strcmp("project", name))
			mat->elim++,handle_project(f, name, mat, nv, vars);
	}
	if (end!=1)
		parse_error(mat, "single end line required");
	if (mat->vrep && (mat->redund!=0 || mat->elim!=0))
		parse_error(mat, "eliminate/project line not allowed in V-representation inputs\n");
	if (mat->vrep && nrays>0)
	{
		mat->rays = malloc(sizeof(unsigned int)*nrays);
		j=0;
		for (i=1; i<=mat->m; i++)
			if (mat->raymask[i])
				mat->rays[j++] = i;
		mat->nrays = nrays;
	}

	return mat;
}

/* 2020.8.21 do negation the annoying way to make MathSAT happy */
void print_prat(lrs_mp Nt, lrs_mp Dt)
{
	lrs_mp pNt, pDt;
	char *tmp1 = cpmp("", Nt), *tmp2=cpmp("", Dt);
	int neg1 = negative(Nt), neg2 = negative(Dt);
	int neg = (neg1 != neg2); /* Nt/Dt negative? */

	lrs_alloc_mp(pNt); lrs_alloc_mp(pDt);
	copy(pNt, Nt); copy(pDt, Dt);
	if (neg1)
		changesign(pNt);
	if (neg2)
		changesign(pDt);

	tmp1 = cpmp("", pNt);
	tmp2 = cpmp("", pDt);

	if (neg)
		printf("(- ");

	if (one(Dt))
		printf("%s", tmp1);
	else
		printf("(/ %s %s)", tmp1, tmp2);

	if (neg)
		printf(")"); /* close negation */
	free(tmp1); free(tmp2);
	lrs_clear_mp(pNt); lrs_clear_mp(pDt);
}

void print_orig_header(unsigned long *vars, unsigned int nv, unsigned long ineq, const char *filename, const int argc, char **argv)
{
	unsigned int i;
	printf("; polyv %s formula for checking non-redundancy of inequality %lu after eliminating variables", VERSION, ineq);
	for (i=0; i<nv; i++)
		printf(" %lu", vars[i]);
	printf(", using file %s\n", filename);
	printf("; invocation was:");
	for (i=0; i<argc; i++)
		printf(" %s", argv[i]);
	printf("\n");

	printf("(set-logic LRA)\n\n");
	return;
}

/* bool, is inequality i actually an equality? */
int islin(struct mat *mat, unsigned int i)
{
	return mat->linmask[i];
}

/* print inequality i, ax + by <= c, but done as 0 <= ...
 * invert a flag, print the <= if false, print > if true
 * also handle linearities.
 * if newn is not NULL, we need to use new names and newn[i]==k means
 * we should use x_k instead of x_i
 * so far only used with mode=3, normally newn==NULL
 */
void printineq(struct mat *mat, unsigned int i, int invert, unsigned int *newn)
{
	unsigned int j, n = mat->n;
	int lin = islin(mat, i);
	int one;

	printf("\n         (");
	if (invert == 0)
		printf("%s", (lin ? "=" : "<="));
	else /* inverting, so <= -> >, = -> noteq */
		printf("%s", (lin ? "not (=" : ">"));
	printf(" 0\n          (+ ");
	print_prat(mat->num[i][0], mat->den[i][0]);
	for (j=1; j<n; j++)
	{
		if (zero(mat->num[i][j]))
			continue; /* don't bother: 0 * x */
		one = 0;
		if (!mpz_cmp(mat->num[i][j], mat->den[i][j]))
			one = 1;
		
		if (!one) /* just do x_j instead of (* (/ a a) x_j) */
		{
			printf(" (* ");
			print_prat(mat->num[i][j], mat->den[i][j]);
		}
		printf(" x_%u", (newn==NULL?j:newn[j]));
		if (!one)
			printf(")");
	}
	if (invert && lin) /* printed "not (=, so close the extra paren here */
		printf(")");
	printf("))"); /* close + and <= */
}

unsigned int *make_newn(struct mat *mat)
{
	unsigned int i, j=1, n=mat->n;
	unsigned int *newn = calloc(n+1, sizeof(unsigned int));
	unsigned int *varmask=mat->varmask;

	if (newn==NULL)
		return NULL;

	/* varmask[] gives the variables projected out, we want the
	 * ones kept
	 */
	for (i=1; i<=mat->n; i++)
	{
		if (varmask[i] == 0)
			newn[j++] = i;
	}
	return newn;
}

/* print Ax + By <= C from diagram, but 
 * done as 0 <= ...
 * Omit inequality ineq per diagram
 * if offs is not NULL, we need to use offsets and offs[i]==k means
 * we've projected out k variables at this point and should use x_{i+k}
 * instead of x_i
 * so far only used with mode=3, normally offs==NULL
 */
void printsystem(struct mat *mat, unsigned long ineq, unsigned int *offs)
{
	unsigned int i, m=mat->m;
	printf("      (and"); 
	for (i=1; i<=m; i++)
	{
		if (i == ineq)
			continue;
		printineq(mat, i, 0, offs);
	}
	printf(")"); /* close and */
}

int projvar(struct mat *mat, unsigned int i)
{
	return mat->varmask[i];
}

/* rows give a_0 + a_1x_1 + a_2x_2 ... a_{n-1}x_{n_1} \ge 0 */
void print_formulabody(struct mat *mat, unsigned int nv, unsigned long *vars,
		       unsigned long ineq)
{
	unsigned int i, n=mat->n;

	/* declare constants x_1 ... x_{n-1} so (get-model) gives model */
	for (i=1; i<n; i++)
#ifdef USEEXISTS
		if (!projvar(mat, i)) /* omit variables projecting out */
#endif
			printf("(declare-const x_%u Real)\n", i);

	printf("\n(assert (and\n");

	/* \exists y st Ax + By \le C from diagram */
	/* in loop in case we can generalize to projecting out multiple vars */
#ifdef USEEXISTS
	printf("  (exists (");
	for (i=1; i<n; i++)
		if (projvar(mat, i))
			printf("(x_%u Real) ", i);
	printf(")");
#endif
	printsystem(mat, ineq, NULL);
#ifdef USEEXISTS
	printf("\n  )\n"); /* close exists */
#endif

	/* forall y : Ax+By\le C from diagram */
	printf("  (forall (");
	for (i=1; i<n; i++)
		if (projvar(mat, i))
			printf("(x_%u Real) ", i);
	printf(")");
	printf("\n    (=> ");
	printsystem(mat, ineq, NULL);

	/* implies  ax + by > c from diagram */
	printineq(mat, ineq, 1, NULL);
	printf("    )"); /* close => */
	printf("\n  )\n))"); /* close forall, and,assert */

	return;
}

void print_formula(struct mat *mat, unsigned int nv, unsigned long *vars,
		   unsigned long ineq, const char *filename, const int argc,
		   char **argv)
{
	print_orig_header(vars, nv, ineq, filename, argc, argv);
	print_formulabody(mat, nv, vars, ineq);
	printf("\n(check-sat)\n(get-model)\n");
	return;
}

/* f is pointing at a floating point representation of a number,
 * read it into v[i]
 */
void readmp_float(FILE *f, lrs_mp_vector v, unsigned int i)
{
	char s[1024];
	char c;

	eatspace(f);

	fscanf(f, "%[0-9]", s);
	atomp(s, v[i]);

	eatspace(f);
	
	c = getc(f);
	if (c!='.')
	{
		if (isdigit(c)) /* cvc4 does not always add .0 to ints */
		{
			ungetc(c,f);
			c = ' ';
		}
		else if (c!=')') /* cvc4 does not always add .0 to ints */
		{
			printf("parse error 1 %c\n", c);
			exit(1);
		}
	}
	while (!isspace(c) && c!=')')
	{
		c = getc(f);
		if (isdigit(c) && c!='0')
		{
			printf("unexpected decimal part %c\n", c);
			exit(1);
		}
	}
}

/* f is pointing at a z3 representation of a number, 
 * e.g. 15.0 or (/ 3.0 7.0)
 * or (- a b) etc.
 * read this into num[i] / den[i]
 */
void readmp_z3(FILE *f, lrs_mp_vector num, lrs_mp_vector den, unsigned int i)
{
	char c;

	eatspace(f);
	c=getc(f);
	if (isdigit(c))
	{
		ungetc(c, f);
		/* now pointing at a floating point representation */
		readmp_float(f, num, i);
		itomp(ONE, den[i]);
		return;
	}
	/* now pointing at "/ a b)" */
	c = getc(f);
	if (c == '/')
	{
		readmp_float(f, num, i);
		readmp_float(f, den, i);
		c = getc(f);
		if (c != ')')
		{
			printf("parse error 2 %c\n", c);
			exit(1);
		}
	}
	else if (c == '-')
	{
		readmp_z3(f, num, den, i);
		changesign(num[i]);
		c = getc(f);
		if (c != ')')
		{
			printf("parse error 3 %c\n", c);
			exit(1);
		}
	}
	else
	{
		printf("unknown character %c\n", c);
		exit(1);
	}
	
}

/* read z3 output from stdin, put it in num/den
 * notes: x_i read to index i, index 0 unused.
 *        z3 outputs reals, we read them in as naturals or rationals
 */
void parse_witness(lrs_mp_vector num, lrs_mp_vector den, unsigned int n)
{
	FILE *f = stdin;
	unsigned int i;
	unsigned int index=0;

	lrs_ifp = f; /* hook to make readmp work */

	eatspace(f); fscanf(f, "sat"); eatspace(f); fscanf(f, "(model");
	for (i=1; i<n; i++)
	{
		eatspace(f);
		fscanf(f, "(define-fun x_%u () Real", &index);
		readmp_z3(f, num, den, index);
	}
	
}

/* print the linearity line if needed */
void cert_print_linearity(struct mat *mat, unsigned int nv, unsigned long ineq)
{
	unsigned int i;
	unsigned int bound=0; /* set for warning removal */
	unsigned int nlin;
	unsigned int lin;

	if (mode == 1)
		bound = mat->n-1;
	if (mode == 2)
		bound = mat->n-1 - nv;

	nlin = mat->nlinearity + bound;

	/* if we're checking a linearity, it gets omitted in certificate 1 */
	if (mode == 1 && islin(mat, ineq))
		nlin--;

	printf("linearity %u", nlin);

	for (i=0; i<mat->nlinearity; i++)
	{
		if (mode==1 && mat->linearities[i] == ineq)
			continue;
		lin = mat->linearities[i];
		/* mode 1 removes ineq, so later linearities shift down */
		if (mode==1 && lin>ineq)
			lin--;
		printf(" %u", lin);
	}

	for (i=1; i<=bound; i++)
		printf(" %u", mat->m + (mode==1?i-1:i));

	printf("\n");
}

/* print the header for the lrs certificate, which may have linearities */
void cert_print_header(struct mat *mat, unsigned int nv, const char *filename,
		       const int argc, char **argv, unsigned long ineq)
{
	printf("certificate %d for %s\nH-representation\n", mode, filename);
	cert_print_linearity(mat, nv, ineq);
	printf("begin\n");
	if (mode == 1)
		printf("%u", mat->m -1 + mat->n-1);
	else /* mode == 2 */
		printf("%u", mat->m + mat->n-1 - nv);
	printf(" %u rational\n", mat->n);
}

/* print the original system: note omit the inequality if mode==1 */
void lrs_print_system(struct mat *mat, unsigned long ineq)
{
	unsigned int i, j;

	for (i=1; i<=mat->m; i++)
	{
		if (i == ineq && mode == 1) /* omit this one for mode==1 */
			continue;
		for (j=0; j<mat->n; j++)
			prat("", mat->num[i][j], mat->den[i][j]);
		printf("\n");
	}
}

/* print the witness part */
void lrs_print_witness(struct mat *mat, lrs_mp_vector wn, lrs_mp_vector wd)
{
	unsigned int i, j;

	for (i=1; i<mat->n; i++)
	{
		if (mode == 2 && mat->varmask[i]) /* don't include y */
			continue;
		prat("", wn[i], wd[i]);
		for (j=1; j<mat->n; j++)
			printf(" %s", (j==i? "-1" : "0"));
		printf("\n");
	}
}

void print_cert(struct mat *mat, unsigned int nv, unsigned long *vars,
		   unsigned long ineq, const char *filename, const int argc,
		   char **argv)
{
	lrs_mp_vector witness_n = lrs_alloc_mp_vector(mat->n);
	lrs_mp_vector witness_d = lrs_alloc_mp_vector(mat->n);
	parse_witness(witness_n, witness_d, mat->n);

	/* print new header */
	cert_print_header(mat, nv, filename, argc, argv, ineq);
	/* print original input, omitting one inequality if mode==1 */
	lrs_print_system(mat, ineq);
	/* print new part */
	lrs_print_witness(mat, witness_n, witness_d);

	printf("end\nlponly\n");

	lrs_clear_mp_vector(witness_n,mat->n);
	lrs_clear_mp_vector(witness_d,mat->n);
	return;
}

void free_mat(struct mat *mat)
{
	if (mat==NULL)
		return;
	lrs_clear_mp_matrix(mat->num, mat->m, mat->n);
	lrs_clear_mp_matrix(mat->den, mat->m, mat->n);
	free(mat->linearities);
	free(mat->linmask);
	free(mat->varmask);
	free(mat->rays);
	free(mat->raymask);
	free(mat->vars);
	free(mat);
	return;
}

/* vmat is a V-representation, write the formula
 * defining the polytope (free variables x_%u
 */
void printcombo(struct mat *vmat)
{
	unsigned int i, j, m=vmat->m, n=vmat->n;

	printf("    (exists (");
	for (i=1; i<=m; i++)
		printf("(a_%u Real) ", i);
	printf(")\n      "); /* close variable list */
	printf("(and\n");

	/* x is the appropriate combo of vertices, rays, lines */
	for (j=1; j<n; j++)
	{ /* x_j is the right value */
		printf("        (= x_%u (+ ", j);
		for (i=1; i<=m; i++)
		{
			printf("(* a_%u ", i);
			print_prat(vmat->num[i][j], vmat->den[i][j]);
			printf(") "); /* close * */
		}
		printf("))\n"); /* close +, close = */
	}

	/* these a_i satisfy the necessary properties */
	/* first, a_i for vertices and rays are >= 0 */
	printf("       ");
	for (i=1; i<=m; i++)
		if (!vmat->linmask[i]) /* not lines */
			printf(" (>= a_%u 0)", i);
	printf("\n        ");
	/* next, convex combo of vertices (those a_i sum to 1 */
	if (vmat->nrays+vmat->nlinearity!=m)
	{ /* only if some vertex exists */
		printf("(= 1 (+");
		for (i=1; i<=m; i++)
			if (!vmat->linmask[i] && !vmat->raymask[i])
				printf(" a_%u", i);
		printf("))\n"); /* close +, close = */
	}
	printf("      )\n    )\n"); /* close and, close exists */

	return;
}

/* make an array of new variable names for mat1 so that the variables projecting
 * onto are 1..., and the others mat->n-1-mat->nv...mat->n-1
 */
unsigned int *make_hh_newn(struct mat *mat)
{
	unsigned int i, n=mat->n;
	unsigned int nk=1, ne=n-mat->nv; /* first new names for kept,elim vars*/
	unsigned int *newn = malloc(sizeof(unsigned int)*n+1);

	newn[0]=0;

	for (i=1; i<n; i++)
	{
		if (projvar(mat, i))
			newn[i]=ne++;
		else
			newn[i]=nk++;
	}

	return newn;
}

void print_hv_body(struct mat *hmat, struct mat *vmat)
{
	unsigned int i, n=hmat->n, n2=vmat->n;
	unsigned int *newn = NULL;

	for (i=1; i<n2; i++) /* if nv==0, hmat->n==n2, no projecting out */
		printf("(declare-const x_%u Real)\n", i);

	printf("(assert (not\n  (= \n");

	if (hmat->nv>0) /* projecting out some variables */
	{
		newn = make_hh_newn(hmat);
		printf("   (exists (");
		for (i=1; i<n; i++)
			if (projvar(hmat, i))
				printf("(x_%u Real) ", newn[i]);
		printf(")\n"); /* close variable list */
	}
	printsystem(hmat, 0, newn);
	if (hmat->nv>0)
		printf("   )\n"); /* close exists */
	printcombo(vmat);
	printf("  )))"); /* close =, close not, close assert */
	free(newn);
	return;
}

void print_vv_body(struct mat *mat1, struct mat *mat2)
{
	unsigned int i, n=mat1->n;
	for (i=1; i<n; i++)
		printf("(declare-const x_%u Real)\n", i);

	printf("(assert (=\n");
	printcombo(mat1);
	printf("  (not\n");
	printcombo(mat2);
	printf("  ) ))\n"); /* close not, =, assert */
	return;
}

int run_hvcheck(struct mat *mat1, struct mat *mat2, int argc, char **argv)
{
	struct mat *vmat=NULL, *hmat=NULL;
	int abort=0;

	if (mat1->vrep && !mat2->vrep)
		vmat=mat1,hmat=mat2;
	else if (mat2->vrep && !mat1->vrep)
		vmat=mat2,hmat=mat1;
	else
	{
		fprintf(stderr, "Error: must give one H and one V representation\n");
		abort=1;
	}

	if (abort==0 && vmat->n != hmat->n-hmat->nv)
	{
		fprintf(stderr, "Dimensions don't match\n");
		abort=1;
	}

	if (abort==1)
	{
		free_mat(mat1);
		free_mat(mat2);
		return 0;
	}

	print_header(argc, argv, "H-V representation");
	print_hv_body(hmat, vmat);
	printf("\n(check-sat)\n(get-model)\n");

	free_mat(mat1);
	free_mat(mat2);
	return 0;
}

int run_vvcheck(struct mat *mat1, struct mat *mat2, int argc, char **argv)
{
	if (mat1->n != mat2->n)
	{
		fprintf(stderr, "Dimensions don't match\n");
		free_mat(mat1);
		free_mat(mat2);
		return 0;
	}

	print_header(argc, argv, "V-V representation");
	print_vv_body(mat1, mat2);
	printf("\n(check-sat)\n(get-model)\n");

	free_mat(mat1);
	free_mat(mat2);
	return 0;
}

void print_hh_bit(struct mat *mat1, struct mat *mat2, unsigned int *newn1,
		  unsigned int *newn2)
{
	unsigned int d = mat1->n - mat1->nv;
	unsigned int i;

	printf("  (and\n");
	if (mat1->nv > 0)
	{ 
		printf("    (exists (");
		for (i=d; i<mat1->n; i++)
			printf("(x_%u Real)", i);
		printf(")\n"); /* close variable list */
	}
	printsystem(mat1, 0, newn1);
	printf("%c\n    (not\n", (mat1->nv>0?')':' ')); /* close exists */

	if (mat2->nv > 0)
	{
		printf("      (exists (");
		for (i=d; i<mat2->n; i++)
			printf("(x_%u Real)", i);
		printf(")\n"); /* close variable list */
	}
	printsystem(mat2, 0, newn2);
	printf("%c\n    )\n  )\n", (mat2->nv>0?')':' ')); /* close exists,
							   * not, and
							   */
	return;
}

void print_hh_body(struct mat *mat1, struct mat *mat2)
{
	unsigned int *newn1=make_hh_newn(mat1),
		     *newn2=make_hh_newn(mat2);

	unsigned int i, d = mat1->n - mat1->nv;

	for (i=1; i<d; i++)
		printf("(declare-const x_%u Real)\n", i);

	printf("(assert (or\n");

	print_hh_bit(mat1, mat2, newn1, newn2);
	print_hh_bit(mat2, mat1, newn2, newn1);
	printf("))\n"); /* close or, close assert */

	free(newn1);
	free(newn2);
	return;
}

int run_hhcheck(struct mat *mat1, struct mat *mat2, int argc, char **argv)
{
	int abort=0;

	if (mat1->vrep || mat2->vrep)
	{
		fprintf(stderr, "Error: must give two H representations\n");
		abort=1;
	}
	if (mat1->n - mat1->nv != mat2->n - mat2->nv)
	{
		fprintf(stderr, "Dimensions of projections don't match\n");
		abort=1;
	}

	if (abort==1)
	{
		free_mat(mat1);
		free_mat(mat2);
		return 0;
	}

	print_header(argc, argv, "equivalency of projected H representations");
	print_hh_body(mat1, mat2);
	printf("\n(check-sat)\n(get-model)\n");

	free_mat(mat1);
	free_mat(mat2);
	return 0;
}

void print_intbody(struct mat *mat1, struct mat *mat2, struct mat *mat3)
{
	unsigned int *newn1=NULL, *newn2=NULL, *newn3=NULL;
	unsigned int i, nres = mat1->n - (mat1->vrep? 0 : mat1->nv);

	/* nres the resulting n, mat1,mat2,mat3 already checked for agreement*/

	for (i=1; i<nres; i++)
		printf("(declare-const x_%u Real)\n", i);

	/* when projecting out columns, rename variables so x_1..x_nres-1
	 * are the columns remaining
	 */
	if (!mat1->vrep && mat1->nv>0)
		newn1 = make_hh_newn(mat1);
	if (!mat2->vrep && mat2->nv>0)
		newn2 = make_hh_newn(mat2);
	if (!mat3->vrep && mat3->nv>0)
		newn3 = make_hh_newn(mat3);

	printf("(assert (not\n  (=\n");

	/* left side, mat1 interset mat2 */
	printf("   (and\n");
	if (mat1->vrep)
		printcombo(mat1);
	else
	{
		if (mat1->nv>0) /* TODO refactor this out, merge with
				 * relevant bit of print_hv_body
				 */
		{
			printf("    (exists (");
			for (i=1; i<mat1->n; i++)
				if (projvar(mat1, i))
					printf("(x_%u Real) ",newn1[i]);
			printf(")\n"); /* close variable list */
		} 
		printsystem(mat1, 0, newn1);
		if (mat1->nv>0)
			printf("    )\n"); /* close exists */
	}
	if (mat2->vrep)
		printcombo(mat2);
	else
	{
		if (mat2->nv>0) /* TODO as above */
		{
			printf("    (exists (");
			for (i=1; i<mat2->n; i++)
				if (projvar(mat2, i))
					printf("(x_%u Real) ", newn2[i]);
			printf(")\n"); /* close variable list */
		}
		printsystem(mat2, 0, newn2);
		if (mat2->nv>0)
			printf("    )\n"); /* close exists */
	}
	printf("   )\n"); /* close and */

	/*print right side */
	if (mat3->vrep)
		printcombo(mat3);
	else
	{
		if (mat3->nv>0) /* TODO as above */
		{
			printf("    (exists (");
			for (i=1; i<mat3->n; i++)
				if (projvar(mat3, i))
					printf("(x_%u Real) ", newn3[i]);
			printf(")\n"); /* close variable list */
		}
		printsystem(mat3, 0, newn3);
		if (mat3->nv>0)
			printf("    )\n"); /* close exists */
	}

	printf("   )\n))"); /* close =, not, assert */
	free(newn1);
	free(newn2);
	free(newn3);
	return;
}

int run_intcheck(struct mat *mat1, struct mat *mat2, struct mat *mat3, 
		 int argc, char **argv)
{
	int dim1, dim2, dim3;

	dim1=mat1->n - (mat1->vrep? 0 : mat1->nv);
	dim2=mat2->n - (mat2->vrep? 0 : mat2->nv);
	dim3=mat3->n - (mat3->vrep? 0 : mat3->nv);

	if (dim1!=dim2 || dim1!=dim3)
	{
		fprintf(stderr, "Dimensions don't match\n");
		free_mat(mat1);
		free_mat(mat2);
		free_mat(mat3);
		return 0;
	}

	fprintf(stderr, "*polyv: checking if %s intersect %s is not %s\n",
		filename1, filename2, filename3);
	print_header(argc, argv, "intersection checking");
	print_intbody(mat1, mat2, mat3);
	printf("\n(check-sat)\n(get-model)\n");

	free_mat(mat1);
	free_mat(mat2);
	free_mat(mat3);
	return 0;
}