Added cddlib-094h from http://www.inf.ethz.ch/personal/fukudak/cdd_home/
Change-Id: I64519509269e434b1b9ea87c3fe0805e711c0ac9
diff --git a/third_party/cddlib/lib-src/cddlib.c b/third_party/cddlib/lib-src/cddlib.c
new file mode 100644
index 0000000..ca9469f
--- /dev/null
+++ b/third_party/cddlib/lib-src/cddlib.c
@@ -0,0 +1,392 @@
+/* cddlib.c: cdd library (library version of cdd)
+ written by Komei Fukuda, fukuda@math.ethz.ch
+ Version 0.94h, April 30, 2015
+ Standard ftp site: ftp.ifor.math.ethz.ch, Directory: pub/fukuda/cdd
+*/
+
+/* cdd : C-Implementation of the double description method for
+ computing all vertices and extreme rays of the polyhedron
+ P= {x : b - A x >= 0}.
+ Please read COPYING (GNU General Public Licence) and
+ the manual cddlibman.tex for detail.
+*/
+
+/* 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., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+/* The first version C0.21 was created on November 10,1993
+ with Dave Gillespie's p2c translator
+ from the Pascal program pdd.p written by Komei Fukuda.
+*/
+
+#include "setoper.h"
+ /* set operation library header (Ver. June 1, 2000 or later) */
+#include "cdd.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <math.h>
+#include <string.h>
+
+/* Global Variables */
+dd_boolean dd_debug =dd_FALSE;
+dd_boolean dd_log =dd_FALSE;
+/* GLOBAL CONSTANTS and STATICS VARIABLES (to be set by dd_set_global_constants() */
+mytype dd_zero;
+mytype dd_one;
+mytype dd_purezero;
+mytype dd_minuszero;
+mytype dd_minusone;
+
+time_t dd_statStartTime; /* cddlib starting time */
+long dd_statBApivots; /* basis finding pivots */
+long dd_statCCpivots; /* criss-cross pivots */
+long dd_statDS1pivots; /* phase 1 pivots */
+long dd_statDS2pivots; /* phase 2 pivots */
+long dd_statACpivots; /* anticycling (cc) pivots */
+#ifdef GMPRATIONAL
+long dd_statBSpivots; /* basis status checking pivots */
+#endif
+dd_LPSolverType dd_choiceLPSolverDefault; /* Default LP solver Algorithm */
+dd_LPSolverType dd_choiceRedcheckAlgorithm; /* Redundancy Checking Algorithm */
+dd_boolean dd_choiceLexicoPivotQ; /* whether to use the lexicographic pivot */
+
+/* #include <profile.h> THINK C PROFILER */
+/* #include <console.h> THINK C PROFILER */
+
+void dd_DDInit(dd_ConePtr cone)
+{
+ cone->Error=dd_NoError;
+ cone->CompStatus=dd_InProgress;
+ cone->RayCount = 0;
+ cone->TotalRayCount = 0;
+ cone->FeasibleRayCount = 0;
+ cone->WeaklyFeasibleRayCount = 0;
+ cone->EdgeCount=0; /* active edge count */
+ cone->TotalEdgeCount=0; /* active edge count */
+ dd_SetInequalitySets(cone);
+ dd_ComputeRowOrderVector(cone);
+ cone->RecomputeRowOrder=dd_FALSE;
+}
+
+void dd_DDMain(dd_ConePtr cone)
+{
+ dd_rowrange hh, itemp, otemp;
+ dd_boolean locallog=dd_log; /* if dd_log=dd_FALSE, no log will be written. */
+
+ if (cone->d<=0){
+ cone->Iteration=cone->m;
+ cone->FeasibleRayCount=0;
+ cone->CompStatus=dd_AllFound;
+ goto _L99;
+ }
+ if (locallog) {
+ fprintf(stderr,"(Initially added rows ) = ");
+ set_fwrite(stderr,cone->InitialHalfspaces);
+ }
+ while (cone->Iteration <= cone->m) {
+ dd_SelectNextHalfspace(cone, cone->WeaklyAddedHalfspaces, &hh);
+ if (set_member(hh,cone->NonequalitySet)){ /* Skip the row hh */
+ if (dd_debug) {
+ fprintf(stderr,"*The row # %3ld should be inactive and thus skipped.\n", hh);
+ }
+ set_addelem(cone->WeaklyAddedHalfspaces, hh);
+ } else {
+ if (cone->PreOrderedRun)
+ dd_AddNewHalfspace2(cone, hh);
+ else{
+ dd_AddNewHalfspace1(cone, hh);
+ }
+ set_addelem(cone->AddedHalfspaces, hh);
+ set_addelem(cone->WeaklyAddedHalfspaces, hh);
+ }
+ if (!cone->PreOrderedRun){
+ for (itemp=1; cone->OrderVector[itemp]!=hh; itemp++);
+ otemp=cone->OrderVector[cone->Iteration];
+ cone->OrderVector[cone->Iteration]=hh;
+ /* store the dynamic ordering in ordervec */
+ cone->OrderVector[itemp]=otemp;
+ /* store the dynamic ordering in ordervec */
+ }
+ if (locallog){
+ fprintf(stderr,"(Iter, Row, #Total, #Curr, #Feas)= %5ld %5ld %9ld %6ld %6ld\n",
+ cone->Iteration, hh, cone->TotalRayCount, cone->RayCount,
+ cone->FeasibleRayCount);
+ }
+ if (cone->CompStatus==dd_AllFound||cone->CompStatus==dd_RegionEmpty) {
+ set_addelem(cone->AddedHalfspaces, hh);
+ goto _L99;
+ }
+ (cone->Iteration)++;
+ }
+ _L99:;
+ if (cone->d<=0 || cone->newcol[1]==0){ /* fixing the number of output */
+ cone->parent->n=cone->LinearityDim + cone->FeasibleRayCount -1;
+ cone->parent->ldim=cone->LinearityDim - 1;
+ } else {
+ cone->parent->n=cone->LinearityDim + cone->FeasibleRayCount;
+ cone->parent->ldim=cone->LinearityDim;
+ }
+}
+
+
+void dd_InitialDataSetup(dd_ConePtr cone)
+{
+ long j, r;
+ dd_rowset ZSet;
+ static dd_Arow Vector1,Vector2;
+ static dd_colrange last_d=0;
+
+ if (last_d < cone->d){
+ if (last_d>0) {
+ for (j=0; j<last_d; j++){
+ dd_clear(Vector1[j]);
+ dd_clear(Vector2[j]);
+ }
+ free(Vector1); free(Vector2);
+ }
+ Vector1=(mytype*)calloc(cone->d,sizeof(mytype));
+ Vector2=(mytype*)calloc(cone->d,sizeof(mytype));
+ for (j=0; j<cone->d; j++){
+ dd_init(Vector1[j]);
+ dd_init(Vector2[j]);
+ }
+ last_d=cone->d;
+ }
+
+ cone->RecomputeRowOrder=dd_FALSE;
+ cone->ArtificialRay = NULL;
+ cone->FirstRay = NULL;
+ cone->LastRay = NULL;
+ set_initialize(&ZSet,cone->m);
+ dd_AddArtificialRay(cone);
+ set_copy(cone->AddedHalfspaces, cone->InitialHalfspaces);
+ set_copy(cone->WeaklyAddedHalfspaces, cone->InitialHalfspaces);
+ dd_UpdateRowOrderVector(cone, cone->InitialHalfspaces);
+ for (r = 1; r <= cone->d; r++) {
+ for (j = 0; j < cone->d; j++){
+ dd_set(Vector1[j], cone->B[j][r-1]);
+ dd_neg(Vector2[j], cone->B[j][r-1]);
+ }
+ dd_Normalize(cone->d, Vector1);
+ dd_Normalize(cone->d, Vector2);
+ dd_ZeroIndexSet(cone->m, cone->d, cone->A, Vector1, ZSet);
+ if (set_subset(cone->EqualitySet, ZSet)){
+ if (dd_debug) {
+ fprintf(stderr,"add an initial ray with zero set:");
+ set_fwrite(stderr,ZSet);
+ }
+ dd_AddRay(cone, Vector1);
+ if (cone->InitialRayIndex[r]==0) {
+ dd_AddRay(cone, Vector2);
+ if (dd_debug) {
+ fprintf(stderr,"and add its negative also.\n");
+ }
+ }
+ }
+ }
+ dd_CreateInitialEdges(cone);
+ cone->Iteration = cone->d + 1;
+ if (cone->Iteration > cone->m) cone->CompStatus=dd_AllFound; /* 0.94b */
+ set_free(ZSet);
+}
+
+dd_boolean dd_CheckEmptiness(dd_PolyhedraPtr poly, dd_ErrorType *err)
+{
+ dd_rowset R, S;
+ dd_MatrixPtr M=NULL;
+ dd_boolean answer=dd_FALSE;
+
+ *err=dd_NoError;
+
+ if (poly->representation==dd_Inequality){
+ M=dd_CopyInequalities(poly);
+ set_initialize(&R, M->rowsize);
+ set_initialize(&S, M->rowsize);
+ if (!dd_ExistsRestrictedFace(M, R, S, err)){
+ poly->child->CompStatus=dd_AllFound;
+ poly->IsEmpty=dd_TRUE;
+ poly->n=0;
+ answer=dd_TRUE;
+ }
+ set_free(R);
+ set_free(S);
+ dd_FreeMatrix(M);
+ } else if (poly->representation==dd_Generator && poly->m<=0){
+ *err=dd_EmptyVrepresentation;
+ poly->IsEmpty=dd_TRUE;
+ poly->child->CompStatus=dd_AllFound;
+ answer=dd_TRUE;
+ poly->child->Error=*err;
+ }
+
+ return answer;
+}
+
+
+dd_boolean dd_DoubleDescription(dd_PolyhedraPtr poly, dd_ErrorType *err)
+{
+ dd_ConePtr cone=NULL;
+ dd_boolean found=dd_FALSE;
+
+ *err=dd_NoError;
+
+ if (poly!=NULL && (poly->child==NULL || poly->child->CompStatus!=dd_AllFound)){
+ cone=dd_ConeDataLoad(poly);
+ /* create a cone associated with poly by homogenization */
+ time(&cone->starttime);
+ dd_DDInit(cone);
+ if (poly->representation==dd_Generator && poly->m<=0){
+ *err=dd_EmptyVrepresentation;
+ cone->Error=*err;
+ goto _L99;
+ }
+ /* Check emptiness of the polyhedron */
+ dd_CheckEmptiness(poly,err);
+
+ if (cone->CompStatus!=dd_AllFound){
+ dd_FindInitialRays(cone, &found);
+ if (found) {
+ dd_InitialDataSetup(cone);
+ if (cone->CompStatus==dd_AllFound) goto _L99;
+ dd_DDMain(cone);
+ if (cone->FeasibleRayCount!=cone->RayCount) *err=dd_NumericallyInconsistent; /* cddlib-093d */
+ }
+ }
+ time(&cone->endtime);
+ }
+
+_L99: ;
+
+ return found;
+}
+
+dd_boolean dd_DoubleDescription2(dd_PolyhedraPtr poly, dd_RowOrderType horder, dd_ErrorType *err)
+{
+ dd_ConePtr cone=NULL;
+ dd_boolean found=dd_FALSE;
+
+ *err=dd_NoError;
+
+ if (poly!=NULL && (poly->child==NULL || poly->child->CompStatus!=dd_AllFound)){
+ cone=dd_ConeDataLoad(poly);
+ /* create a cone associated with poly by homogenization */
+ cone->HalfspaceOrder=horder; /* set the row order */
+ time(&cone->starttime);
+ dd_DDInit(cone);
+ if (poly->representation==dd_Generator && poly->m<=0){
+ *err=dd_EmptyVrepresentation;
+ cone->Error=*err;
+ goto _L99;
+ }
+ /* Check emptiness of the polyhedron */
+ dd_CheckEmptiness(poly,err);
+
+ if (cone->CompStatus!=dd_AllFound){
+ dd_FindInitialRays(cone, &found);
+ if (found) {
+ dd_InitialDataSetup(cone);
+ if (cone->CompStatus==dd_AllFound) goto _L99;
+ dd_DDMain(cone);
+ if (cone->FeasibleRayCount!=cone->RayCount) *err=dd_NumericallyInconsistent; /* cddlib-093d */
+ }
+ }
+ time(&cone->endtime);
+ }
+
+_L99: ;
+
+ return found;
+}
+
+dd_boolean dd_DDInputAppend(dd_PolyhedraPtr *poly, dd_MatrixPtr M,
+ dd_ErrorType *err)
+{
+ /* This is imcomplete. It simply solves the problem from scratch. */
+ dd_boolean found;
+ dd_ErrorType error;
+
+ if ((*poly)->child!=NULL) dd_FreeDDMemory(*poly);
+ dd_AppendMatrix2Poly(poly, M);
+ (*poly)->representation=dd_Inequality;
+ found=dd_DoubleDescription(*poly, &error);
+ *err=error;
+ return found;
+}
+
+dd_boolean dd_DDFile2File(char *ifile, char *ofile, dd_ErrorType *err)
+{
+ /* The representation conversion from an input file to an outfile. */
+ /* modified by D. Avis to allow stdin/stdout */
+ dd_boolean found=dd_TRUE;
+ FILE *reading=NULL,*writing=NULL;
+ dd_PolyhedraPtr poly;
+ dd_MatrixPtr M, A, G;
+
+ if (strcmp(ifile,"**stdin") == 0 )
+ reading = stdin;
+ else if ( ( reading = fopen(ifile, "r") )!= NULL) {
+ fprintf(stderr,"input file %s is open\n", ifile);
+ }
+ else{
+ fprintf(stderr,"The input file %s not found\n",ifile);
+ found=dd_FALSE;
+ *err=dd_IFileNotFound;
+ goto _L99;
+ }
+
+ if (found){
+ if (strcmp(ofile,"**stdout") == 0 )
+ writing = stdout;
+ else if ( (writing = fopen(ofile, "w") ) != NULL){
+ fprintf(stderr,"output file %s is open\n",ofile);
+ found=dd_TRUE;
+ } else {
+ fprintf(stderr,"The output file %s cannot be opened\n",ofile);
+ found=dd_FALSE;
+ *err=dd_OFileNotOpen;
+ goto _L99;
+ }
+ }
+
+ M=dd_PolyFile2Matrix(reading, err);
+ if (*err!=dd_NoError){
+ goto _L99;
+ } poly=dd_DDMatrix2Poly(M, err); /* compute the second representation */ dd_FreeMatrix(M);
+
+ if (*err==dd_NoError) {
+ dd_WriteRunningMode(writing, poly);
+ A=dd_CopyInequalities(poly);
+ G=dd_CopyGenerators(poly);
+
+ if (poly->representation==dd_Inequality) {
+ dd_WriteMatrix(writing,G);
+ } else {
+ dd_WriteMatrix(writing,A);
+ }
+
+ dd_FreePolyhedra(poly);
+ dd_FreeMatrix(A);
+ dd_FreeMatrix(G);
+ }
+
+_L99: ;
+ if (*err!=dd_NoError) dd_WriteErrorMessages(stderr,*err);
+ if (reading!=NULL) fclose(reading);
+ if (writing!=NULL) fclose(writing);
+ return found;
+}
+
+/* end of cddlib.c */