[81384b] | 1 | /* |
---|
| 2 | gfan.h Interface to gfan.cc |
---|
| 3 | |
---|
| 4 | Author: monerjan |
---|
| 5 | */ |
---|
| 6 | #ifndef GFAN_H |
---|
| 7 | #define GFAN_H |
---|
| 8 | |
---|
[b5dd9b] | 9 | #include "kernel/mod2.h" |
---|
[9996416] | 10 | |
---|
| 11 | #if HAVE_GFANLIB |
---|
| 12 | |
---|
[b5dd9b] | 13 | #include "misc/int64vec.h" |
---|
[81384b] | 14 | |
---|
[b5dd9b] | 15 | #include "gfanlib/config.h" |
---|
[888b457] | 16 | #ifdef HAVE_CDD_SETOPER_H |
---|
| 17 | #include <cdd/setoper.h> |
---|
| 18 | #include <cdd/cdd.h> |
---|
| 19 | #include <cdd/cddmp.h> |
---|
[24aeb0f] | 20 | #elif HAVE_CDDLIB_SETOPER_H |
---|
| 21 | #include <cddlib/setoper.h> |
---|
| 22 | #include <cddlib/cdd.h> |
---|
| 23 | #include <cddlib/cddmp.h> |
---|
[888b457] | 24 | #else |
---|
[81384b] | 25 | #include <setoper.h> |
---|
| 26 | #include <cdd.h> |
---|
| 27 | #include <cddmp.h> |
---|
[888b457] | 28 | #endif |
---|
[b5dd9b] | 29 | #include "bbfan.h" |
---|
| 30 | #include "bbcone.h" |
---|
[a3f0fea] | 31 | EXTERN_VAR int gfanHeuristic; |
---|
[81384b] | 32 | |
---|
| 33 | #ifndef USE_ZFAN |
---|
| 34 | #define USE_ZFAN |
---|
| 35 | #endif |
---|
| 36 | #ifndef USE_ZFAN |
---|
| 37 | lists grfan(ideal inputIdeal, int heuristic, bool singleCone); |
---|
| 38 | #else |
---|
[b5dd9b] | 39 | #include "gfanlib/gfanlib.h" |
---|
[81384b] | 40 | gfan::ZFan *grfan(ideal inputIdeal, int h, bool singleCone); |
---|
| 41 | #endif |
---|
| 42 | // lists grcone_by_intvec(ideal inputIdeal); |
---|
| 43 | |
---|
| 44 | class facet |
---|
| 45 | { |
---|
[2bf04b] | 46 | private: |
---|
| 47 | /** \brief Inner normal of the facet, describing it uniquely up to isomorphism */ |
---|
| 48 | int64vec *fNormal; |
---|
[81384b] | 49 | |
---|
[2bf04b] | 50 | /** \brief An interior point of the facet*/ |
---|
| 51 | int64vec *interiorPoint; |
---|
[81384b] | 52 | |
---|
[2bf04b] | 53 | /** \brief Universal Cone Number |
---|
| 54 | * The number of the cone the facet belongs to, Set in getConeNormals() |
---|
| 55 | */ |
---|
| 56 | int UCN; |
---|
[81384b] | 57 | |
---|
[2bf04b] | 58 | /** \brief The codim of the facet |
---|
| 59 | */ |
---|
| 60 | short codim; |
---|
[81384b] | 61 | |
---|
[2bf04b] | 62 | /** \brief The Groebner basis on the other side of a shared facet |
---|
| 63 | * |
---|
| 64 | * In order not to have to compute the flipped GB twice we store the basis we already get |
---|
| 65 | * when identifying search facets. Thus in the next step of the reverse search we can |
---|
| 66 | * just copy the old cone and update the facet and the gcBasis. |
---|
| 67 | * facet::flibGB is set via facet::setFlipGB() and printed via facet::printFlipGB |
---|
| 68 | */ |
---|
| 69 | ideal flipGB; //The Groebner Basis on the other side, computed via gcone::flip |
---|
[81384b] | 70 | |
---|
[2bf04b] | 71 | public: |
---|
| 72 | /** \brief Boolean value to indicate whether a facet is flippable or not |
---|
| 73 | * This is also used to mark facets that nominally are flippable but which do |
---|
| 74 | * not intersect with the positive orthant. This check is done in gcone::getCodim2Normals |
---|
| 75 | */ |
---|
| 76 | bool isFlippable; //**flippable facet? */ |
---|
| 77 | //bool isIncoming; //Is the facet incoming or outgoing in the reverse search? No longer in use |
---|
| 78 | facet *next; //Pointer to next facet |
---|
| 79 | facet *prev; //Pointer to predecessor. Needed for the SearchList in noRevS |
---|
| 80 | facet *codim2Ptr; //Pointer to (codim-2)-facet. Bit of recursion here ;-) |
---|
| 81 | int numCodim2Facets; //#of (codim-2)-facets of this facet. Set in getCodim2Normals() |
---|
| 82 | unsigned numRays; //Number of spanning rays of the facet |
---|
| 83 | ring flipRing; //the ring on the other side of the facet |
---|
| 84 | // int64vec **fRays; |
---|
[81384b] | 85 | |
---|
[2bf04b] | 86 | /** The default constructor. */ |
---|
| 87 | facet(); |
---|
| 88 | /** Constructor for lower dimensional faces*/ |
---|
| 89 | facet(const int &n); |
---|
| 90 | /** The copy constructor */ |
---|
| 91 | facet(const facet& f); |
---|
| 92 | /** A shallow copy of facets*/ |
---|
| 93 | facet* shallowCopy(const facet& f); |
---|
| 94 | void shallowDelete(); |
---|
| 95 | /** The default destructor */ |
---|
| 96 | ~facet(); |
---|
| 97 | /** Comparison operator*/ |
---|
| 98 | // inline bool operator==(const facet *f,const facet *g); |
---|
| 99 | /** \brief Comparison of facets*/ |
---|
| 100 | // inline bool areEqual(facet *f, facet *g);//Now static |
---|
| 101 | /** Stores the facet normal \param int64vec*/ |
---|
| 102 | inline void setFacetNormal(int64vec *iv); |
---|
| 103 | /** Returns the facet normal */ |
---|
| 104 | inline int64vec *getFacetNormal() const; |
---|
| 105 | /** Return a reference to the facet normal*/ |
---|
| 106 | inline const int64vec *getRef2FacetNormal() const; |
---|
| 107 | /** Method to print the facet normal*/ |
---|
| 108 | inline void printNormal() const; |
---|
| 109 | /** Store the flipped GB*/ |
---|
| 110 | inline void setFlipGB(ideal I); |
---|
| 111 | /** Return the flipped GB*/ |
---|
| 112 | inline ideal getFlipGB(); |
---|
| 113 | /** Print the flipped GB*/ |
---|
| 114 | inline void printFlipGB(); |
---|
| 115 | /** Set the UCN */ |
---|
| 116 | inline void setUCN(int n); |
---|
| 117 | /** \brief Get the UCN |
---|
| 118 | * Returns the UCN iff this != NULL, else -1 |
---|
| 119 | */ |
---|
| 120 | inline int getUCN(); |
---|
| 121 | /** Store an interior point of the facet */ |
---|
| 122 | inline void setInteriorPoint(int64vec *iv); |
---|
| 123 | inline int64vec *getInteriorPoint(); |
---|
| 124 | inline const int64vec *getRef2InteriorPoint(); |
---|
| 125 | /** \brief Debugging function |
---|
| 126 | * prints the facet normal an all (codim-2)-facets that belong to it |
---|
| 127 | */ |
---|
| 128 | volatile void fDebugPrint(); |
---|
| 129 | friend class gcone; |
---|
[81384b] | 130 | }; |
---|
| 131 | |
---|
| 132 | |
---|
| 133 | /** |
---|
| 134 | *\brief Implements the cone structure |
---|
| 135 | * |
---|
| 136 | * A cone is represented by a linked list of facet normals |
---|
| 137 | * @see facet |
---|
| 138 | */ |
---|
| 139 | |
---|
| 140 | class gcone |
---|
| 141 | { |
---|
[2bf04b] | 142 | private: |
---|
| 143 | ideal inputIdeal; //the original |
---|
| 144 | ring baseRing; //the basering of the cone |
---|
| 145 | int64vec *ivIntPt; //an interior point of the cone |
---|
| 146 | int UCN; //unique number of the cone |
---|
| 147 | int pred; //UCN of the cone this one is derived from |
---|
[a3f0fea] | 148 | STATIC_VAR int counter; |
---|
[81384b] | 149 | |
---|
[2bf04b] | 150 | public: |
---|
| 151 | /** \brief Pointer to the first facet */ |
---|
| 152 | facet *facetPtr; //Will hold the adress of the first facet; set by gcone::getConeNormals |
---|
[81384b] | 153 | #ifdef gfanp |
---|
[a3f0fea] | 154 | STATIC_VAR float time_getConeNormals; |
---|
| 155 | STATIC_VAR float time_getCodim2Normals; |
---|
| 156 | STATIC_VAR float t_getExtremalRays; |
---|
| 157 | STATIC_VAR float t_ddPolyh; |
---|
| 158 | STATIC_VAR float time_flip; |
---|
| 159 | STATIC_VAR float time_flip2; |
---|
| 160 | STATIC_VAR float t_areEqual; |
---|
| 161 | STATIC_VAR float t_ffG; |
---|
| 162 | STATIC_VAR float t_markings; |
---|
| 163 | STATIC_VAR float t_dd; |
---|
| 164 | STATIC_VAR float t_kStd; |
---|
| 165 | STATIC_VAR float time_enqueue; |
---|
| 166 | STATIC_VAR float time_computeInv; |
---|
| 167 | STATIC_VAR float t_ddMC; |
---|
| 168 | STATIC_VAR float t_mI; |
---|
| 169 | STATIC_VAR float t_iP; |
---|
| 170 | STATIC_VAR float t_isParallel; |
---|
| 171 | STATIC_VAR unsigned parallelButNotEqual; |
---|
| 172 | STATIC_VAR unsigned numberOfFacetChecks; |
---|
[81384b] | 173 | #endif |
---|
[2bf04b] | 174 | /** Matrix to contain the homogeneity/lineality space */ |
---|
[a3f0fea] | 175 | STATIC_VAR dd_MatrixPtr dd_LinealitySpace; |
---|
| 176 | STATIC_VAR int lengthOfSearchList; |
---|
[2bf04b] | 177 | /** Maximum size of the searchlist*/ |
---|
[a3f0fea] | 178 | STATIC_VAR int maxSize; |
---|
[2bf04b] | 179 | /** is the ideal homogeneous? */ |
---|
[a3f0fea] | 180 | STATIC_VAR bool hasHomInput; |
---|
[2bf04b] | 181 | /** # of variables in the ring */ |
---|
[a3f0fea] | 182 | STATIC_VAR int numVars; //#of variables in the ring |
---|
[2bf04b] | 183 | /** The hilbert function - for the homogeneous case*/ |
---|
[a3f0fea] | 184 | STATIC_VAR int64vec *hilbertFunction; |
---|
[2bf04b] | 185 | /** The zero vector. Needed in case of fNormal mismatch*/ |
---|
[a3f0fea] | 186 | STATIC_VAR int64vec *ivZeroVector; |
---|
[81384b] | 187 | |
---|
[2bf04b] | 188 | /** # of facets of the cone |
---|
| 189 | * This value is set by gcone::getConeNormals |
---|
| 190 | */ |
---|
| 191 | int numFacets; //#of facets of the cone |
---|
[81384b] | 192 | |
---|
[2bf04b] | 193 | /** |
---|
| 194 | * At least as a workaround we store the irredundant facets of a matrix here. |
---|
| 195 | * This is needed to compute an interior points of a cone. Note that there |
---|
| 196 | * will be non-flippable facets in it! |
---|
| 197 | */ |
---|
| 198 | dd_MatrixPtr ddFacets; //Matrix to store irredundant facets of the cone |
---|
[81384b] | 199 | |
---|
[2bf04b] | 200 | /** Array of intvecs representing the rays of the cone*/ |
---|
| 201 | int64vec **gcRays; |
---|
| 202 | unsigned numRays; //#rays of the cone |
---|
| 203 | /** Contains the Groebner basis of the cone. Is set by gcone::getGB(ideal I)*/ |
---|
| 204 | ideal gcBasis; //GB of the cone, set by gcone::getGB(); |
---|
| 205 | gcone *next; //Pointer to next cone |
---|
| 206 | gcone *prev; |
---|
[81384b] | 207 | |
---|
[2bf04b] | 208 | gcone(); |
---|
| 209 | gcone(ring r, ideal I); |
---|
| 210 | gcone(const gcone& gc, const facet &f); |
---|
| 211 | ~gcone(); |
---|
| 212 | inline int getCounter(); |
---|
| 213 | inline ring getBaseRing(); |
---|
| 214 | inline ring getRef2BaseRing(); |
---|
| 215 | inline void setBaseRing(ring r); |
---|
| 216 | inline void setIntPoint(int64vec *iv); |
---|
| 217 | inline int64vec *getIntPoint(bool shallow=FALSE); |
---|
| 218 | inline void showIntPoint(); |
---|
| 219 | inline void setNumFacets(); |
---|
| 220 | inline int getNumFacets(); |
---|
| 221 | inline int getUCN(); |
---|
| 222 | inline int getPredUCN(); |
---|
| 223 | volatile void showFacets(short codim=1); |
---|
| 224 | // volatile void showSLA(facet &f); |
---|
| 225 | // void idDebugPrint(const ideal &I); |
---|
| 226 | // void invPrint(const ideal &I); |
---|
| 227 | // bool isMonomial(const ideal &I); |
---|
| 228 | // int64vec *ivNeg(const int64vec *iv); |
---|
| 229 | // inline int dotProduct(int64vec &iva, int64vec &ivb); |
---|
| 230 | // inline int dotProduct(const int64vec &iva, const int64vec &ivb); |
---|
| 231 | // inline bool isParallel(const int64vec &a, const int64vec &b); |
---|
| 232 | void noRevS(gcone &gcRoot, bool usingIntPoint=FALSE); |
---|
| 233 | // inline int intgcd(const int &a, const int &b); |
---|
| 234 | void writeConeToFile(const gcone &gc, bool usingIntPoints=FALSE); |
---|
| 235 | void readConeFromFile(int gcNum, gcone *gc); |
---|
| 236 | int64vec f2M(gcone *gc, facet *f, int n=1); |
---|
| 237 | // inline void sortRays(gcone *gc); |
---|
| 238 | //The real stuff |
---|
| 239 | void getConeNormals(const ideal &I, bool compIntPoint=FALSE); |
---|
| 240 | void getCodim2Normals(const gcone &gc); |
---|
| 241 | void getExtremalRays(const gcone &gc); |
---|
| 242 | void orderRays(); |
---|
| 243 | void flip(ideal gb, facet *f); |
---|
| 244 | void flip2(const ideal &gb, facet *f); |
---|
| 245 | void computeInv(const ideal &gb, ideal &inv, const int64vec &f); |
---|
| 246 | //poly restOfDiv(poly const &f, ideal const &I); removed with r12286 |
---|
| 247 | inline ideal ffG(const ideal &H, const ideal &G); |
---|
| 248 | inline void getGB(ideal const &inputIdeal); |
---|
| 249 | void interiorPoint( dd_MatrixPtr &M, int64vec &iv);//used from flip and optionally from getConeNormals |
---|
| 250 | // void interiorPoint2(); //removed Feb 8th, 2010, new method Feb 19th, 2010, again removed Mar 16th, 2010 |
---|
| 251 | void preprocessInequalities(dd_MatrixPtr &M); |
---|
| 252 | ring rCopyAndAddWeight(const ring &r, int64vec *ivw); |
---|
| 253 | ring rCopyAndAddWeight2(const ring &, const int64vec *, const int64vec *); |
---|
| 254 | // ring rCopyAndChangeWeight(const ring &r, int64vec *ivw); //NOTE remove |
---|
| 255 | // void reverseSearch(gcone *gcAct); //NOTE both removed from r12286 |
---|
| 256 | // bool isSearchFacet(gcone &gcTmp, facet *testfacet); //NOTE remove |
---|
| 257 | void makeInt(const dd_MatrixPtr &M, const int line, int64vec &n); |
---|
| 258 | // void normalize();//NOTE REMOVE |
---|
| 259 | facet * enqueueNewFacets(facet *f); |
---|
| 260 | facet * enqueue2(facet *f); |
---|
| 261 | // dd_MatrixPtr facets2Matrix(const gcone &gc);//NOTE remove |
---|
| 262 | /** Compute the lineality space Ax=0 and return it as dd_MatrixPtr dd_LinealitySpace*/ |
---|
| 263 | dd_MatrixPtr computeLinealitySpace(); |
---|
| 264 | inline bool iv64isStrictlyPositive(const int64vec *); |
---|
| 265 | /** Exchange 2 ordertype_a by just 1 */ |
---|
| 266 | void replaceDouble_ringorder_a_ByASingleOne(); |
---|
| 267 | // static void gcone::idPrint(ideal &I); |
---|
| 268 | // friend class facet; |
---|
[81384b] | 269 | }; |
---|
| 270 | lists lprepareResult(gcone *gc, const int n); |
---|
[69b2c1] | 271 | /* static int64 int64gcd(const int64 &a, const int64 &b); */ |
---|
| 272 | /* static int intgcd(const int &a, const int &b); */ |
---|
| 273 | /* static int dotProduct(const int64vec &iva, const int64vec &ivb); */ |
---|
| 274 | /* static bool isParallel(const int64vec &a, const int64vec &b); */ |
---|
| 275 | /* static int64vec *ivNeg(/\*const*\/ int64vec *iv); */ |
---|
| 276 | /* static void idDebugPrint(const ideal &I); */ |
---|
| 277 | /* static volatile void showSLA(facet &f); */ |
---|
| 278 | /* static bool isMonomial(const ideal &I); */ |
---|
| 279 | /* static bool ivAreEqual(const int64vec &a, const int64vec &b); */ |
---|
| 280 | /* static bool areEqual2(facet *f, facet *g); */ |
---|
| 281 | /* static bool areEqual( facet *f, facet *g); */ |
---|
[81384b] | 282 | // bool iv64isStrictlyPositive(int64vec *); |
---|
| 283 | #endif |
---|
[9996416] | 284 | #endif |
---|