# Changeset 15813d in git

Ignore:
Timestamp:
Aug 8, 2016, 2:16:11 PM (7 years ago)
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
Children:
7cc3fdae6091a6fab65626857c1ed2c59f027dcd
Parents:
47a774fa4174b578843d43227fcc351165473747
Message:
`format`
Location:
gfanlib
Files:
13 edited

Unmodified
Removed
• ## gfanlib/gfanlib_circuittableint.h

 r47a774 class CircuitTableInt32{ public: class Divisor{ public: int32_t v; int shift; int32_t multiplicativeInverse; Divisor(CircuitTableInt32 const &a)// for exact division { // A good place to read about these tricks seems to be the book "Hacker's Delight" by Warren. v=a.v; shift=0; uint32_t t=v; assert(t); while(!(t&1)){t>>=      1;shift++;} uint32_t inverse=t; while(t*inverse!=1)inverse*=2-t*inverse; multiplicativeInverse=inverse; } }; class Double{ public: int64_t v; Double():v(0){}; Double(int64_t a):v(a){}; Double &operator+=(Double a){v+=a.v;return *this;} Double &operator-=(Double a){v-=a.v;return *this;} CircuitTableInt32 castToSingle()const; bool isZero()const{return v==0;} bool isNegative()const{return v<0;} Double operator-()const{Double ret;ret.v=-v;return ret;} friend Double operator-(Double const &a, Double const &b){return Double(a.v-b.v);} friend Double operator+(Double const &a, Double const &b){return Double(a.v+b.v);} Double &addWithOverflowCheck(Double const &a) { if(a.v<0 || v<0 || a.v>1000000000000000000 || v>1000000000000000000) throw MVMachineIntegerOverflow; v+=a.v; return *this; } std::string toString()const{std::stringstream s;s<>=        1;shift++;} uint32_t inverse=t; while(t*inverse!=1)inverse*=2-t*inverse; multiplicativeInverse=inverse; } }; class Double{ public: int64_t v; Double():v(0){}; Double(int64_t a):v(a){}; Double &operator+=(Double a){v+=a.v;return *this;} Double &operator-=(Double a){v-=a.v;return *this;} CircuitTableInt32 castToSingle()const; bool isZero()const{return v==0;} bool isNegative()const{return v<0;} Double operator-()const{Double ret;ret.v=-v;return ret;} friend Double operator-(Double const &a, Double const &b){return Double(a.v-b.v);} friend Double operator+(Double const &a, Double const &b){return Double(a.v+b.v);} Double &addWithOverflowCheck(Double const &a) { if(a.v<0 || v<0 || a.v>1000000000000000000 || v>1000000000000000000) throw MVMachineIntegerOverflow; v+=a.v; return *this; } std::string toString()const{std::stringstream s;s<-30000 && v<30000;}// Is it better to allow only non-negative numbers? static bool isSuitableSummationNumber(int numberOfSummands){return numberOfSummands<30000;} CircuitTableInt32(){v=0;}; CircuitTableInt32(CircuitTableInt32 const &m):v(m.v){} CircuitTableInt32(int32_t val):v(val){}; class CircuitTableInt32 &operator=(int32_t a){v=a;return *this;} CircuitTableInt32 operator-()const{CircuitTableInt32 ret;ret.v=-v;return ret;} CircuitTableInt32 &operator-=(CircuitTableInt32 a){v-=a.v;return *this;} CircuitTableInt32 &operator+=(CircuitTableInt32 a){v+=a.v;return *this;} CircuitTableInt32 &operator*=(CircuitTableInt32 a){v*=a.v;return *this;} friend bool operator<=(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return a.v<=b.v;} friend bool operator==(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return a.v==b.v;} bool isZero()const{return v==0;} bool isOne()const{return v==1;} bool isNonZero()const{return v!=0;} bool isNegative()const{return v<0;} bool isPositive()const{return v>0;} friend int determinantSign1(CircuitTableInt32 const &a, CircuitTableInt32 const &c, CircuitTableInt32 const &b, CircuitTableInt32 const &d)//NOTICE MIXED ORDER. MUST WORK BY EXTENDING { int64_t r=((int64_t)a.v)*((int64_t)c.v)-((int64_t)b.v)*((int64_t)d.v); if(r>0)return 1; if(r<0)return -1; return 0; } friend int determinantSign2(CircuitTableInt32::Double const &a, CircuitTableInt32 const &c, CircuitTableInt32::Double const &b, CircuitTableInt32 const &d)//NOTICE MIXED ORDER. MUST WORK BY EXTENDING { int64_t r1=(a.v)*((int64_t)c.v); int64_t r2=(b.v)*((int64_t)d.v); if(r1>r2)return 1; if(r1extend();t+=extendedMultiplication(a,b);*this=t.castToSingle();return *this;} CircuitTableInt32 &msubWithOverflowChecking(CircuitTableInt32 const &a, CircuitTableInt32 const&b){Double s=this->extend();s-=extendedMultiplication(a,b);*this=s.castToSingle();return *this;} CircuitTableInt32 &subWithOverflowChecking(CircuitTableInt32 const &a){Double t=this->extend();t-=a.extend();*this=t.castToSingle();return *this;} CircuitTableInt32 &addWithOverflowChecking(CircuitTableInt32 const &a){Double t=this->extend();t+=a.extend();*this=t.castToSingle();return *this;} CircuitTableInt32 negatedWithOverflowChecking()const{return (-extend()).castToSingle();} friend Double extendedMultiplication(CircuitTableInt32 const &a, CircuitTableInt32 const &b){Double ret;ret.v=((int64_t)a.v)*((int64_t)b.v);return ret;} friend Double extendedMultiplication(CircuitTableInt32 const &a, int32_t b){Double ret;ret.v=((int64_t)a.v)*((int64_t)b);return ret;}//to be removed? friend CircuitTableInt32 min(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return (a.v>=b.v)?b:a;} friend CircuitTableInt32 max(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return (a.v<=b.v)?b:a;} friend CircuitTableInt32 negabs(CircuitTableInt32 const &a){return min(a,-a);} friend CircuitTableInt32 dotDiv(CircuitTableInt32 const &s, CircuitTableInt32 const &a, CircuitTableInt32 const &t, CircuitTableInt32 const &b, CircuitTableInt32::Divisor const &q) { return CircuitTableInt32(((((int32_t)(((uint64_t)(((int64_t)s.v)*((int64_t)a.v)+((int64_t)t.v)*((int64_t)b.v)))>>q.shift)))*q.multiplicativeInverse)); } friend void dotDivAssign(CircuitTableInt32 &s, CircuitTableInt32 const &a, CircuitTableInt32 const &t, CircuitTableInt32 const &b, CircuitTableInt32::Divisor const &q)//as above but assigning to first argument. This helps the vectorizer. { s.v=((((int32_t)(((uint64_t)(((int64_t)s.v)*((int64_t)a.v)+((int64_t)t.v)*((int64_t)b.v)))>>q.shift)))*q.multiplicativeInverse); } static int MIN(int32_t a, int32_t b) { return (ab)?a:b; } static CircuitTableInt32 computeNegativeBound(CircuitTableInt32 * __restrict__ Ai, int w) { CircuitTableInt32 M=0; CircuitTableInt32 m=0; for(int j=0;j0){s.v>>=1;t.v>>=1;denominatorDivisor.shift--;denominatorDivisor.v>>=1;} CircuitTableInt32 max=0; CircuitTableInt32 min=0; // This is a bound on the absolute value of the sums of products before dividing by the denominator // The calculation can overflow, but since we are casting to unsigned long, this is not a problem. uint64_t positiveResultBoundTimesD=(negabs(t).v*((int64_t)boundA.v)+negabs(s).v*((int64_t)boundB.v)); /* The first case is the case where we know that the intermediate results fit in 32bit before shifting down. * In other words, the shifted-in bits of the results are equal. * In that case the intermediate result can be computed with 32bits. */ if(positiveResultBoundTimesD<((((int64_t)0x40000000)*denominatorDivisor.v)>>denominatorDivisor.shift))//check this carefully {//FAST VERSION // debug<>denominatorDivisor.shift; min.v=MIN(min.v,aa[i].v); max.v=MAX(max.v,aa[i].v); } if(-max<=min)min=-max; return min; } else { /* The next case would be to check if the result is known to fit in 32bit. * In that case intermediate results would be handled in 64 bit, * but the final result would not have to be checked for 32 bit overflows. */ if(positiveResultBoundTimesD<((((int64_t)0x40000000)*denominatorDivisor.v))) { for(int i=0;i>denominatorDivisor.shift)))*denominatorDivisor.multiplicativeInverse); if(max<=aa[i])max=aa[i]; if(aa[i]<=min)min=aa[i]; } if(-max<=min)min=-max; } /* The last case would be to where we don't know that the results will fit in 32 bit. * Since we need to compute max and min anyway, we can compute these quantities in 64bit * and just check if they fit in 32bit at the end. */ else { bool doesOverflow=(((uint32_t)t.v)==0x80000000); int64_t min64=0; int64_t max64=0; for(int i=0;i-30000 && v<30000;}// Is it better to allow only non-negative numbers? static bool isSuitableSummationNumber(int numberOfSummands){return numberOfSummands<30000;} CircuitTableInt32(){v=0;}; CircuitTableInt32(CircuitTableInt32 const &m):v(m.v){} CircuitTableInt32(int32_t val):v(val){}; class CircuitTableInt32 &operator=(int32_t a){v=a;return *this;} CircuitTableInt32 operator-()const{CircuitTableInt32 ret;ret.v=-v;return ret;} CircuitTableInt32 &operator-=(CircuitTableInt32 a){v-=a.v;return *this;} CircuitTableInt32 &operator+=(CircuitTableInt32 a){v+=a.v;return *this;} CircuitTableInt32 &operator*=(CircuitTableInt32 a){v*=a.v;return *this;} friend bool operator<=(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return a.v<=b.v;} friend bool operator==(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return a.v==b.v;} bool isZero()const{return v==0;} bool isOne()const{return v==1;} bool isNonZero()const{return v!=0;} bool isNegative()const{return v<0;} bool isPositive()const{return v>0;} friend int determinantSign1(CircuitTableInt32 const &a, CircuitTableInt32 const &c, CircuitTableInt32 const &b, CircuitTableInt32 const &d)//NOTICE MIXED ORDER. MUST WORK BY EXTENDING { int64_t r=((int64_t)a.v)*((int64_t)c.v)-((int64_t)b.v)*((int64_t)d.v); if(r>0)return 1; if(r<0)return -1; return 0; } friend int determinantSign2(CircuitTableInt32::Double const &a, CircuitTableInt32 const &c, CircuitTableInt32::Double const &b, CircuitTableInt32 const &d)//NOTICE MIXED ORDER. MUST WORK BY EXTENDING { int64_t r1=(a.v)*((int64_t)c.v); int64_t r2=(b.v)*((int64_t)d.v); if(r1>r2)return 1; if(r1extend();t+=extendedMultiplication(a,b);*this=t.castToSingle();return *this;} CircuitTableInt32 &msubWithOverflowChecking(CircuitTableInt32 const &a, CircuitTableInt32 const&b){Double s=this->extend();s-=extendedMultiplication(a,b);*this=s.castToSingle();return *this;} CircuitTableInt32 &subWithOverflowChecking(CircuitTableInt32 const &a){Double t=this->extend();t-=a.extend();*this=t.castToSingle();return *this;} CircuitTableInt32 &addWithOverflowChecking(CircuitTableInt32 const &a){Double t=this->extend();t+=a.extend();*this=t.castToSingle();return *this;} CircuitTableInt32 negatedWithOverflowChecking()const{return (-extend()).castToSingle();} friend Double extendedMultiplication(CircuitTableInt32 const &a, CircuitTableInt32 const &b){Double ret;ret.v=((int64_t)a.v)*((int64_t)b.v);return ret;} friend Double extendedMultiplication(CircuitTableInt32 const &a, int32_t b){Double ret;ret.v=((int64_t)a.v)*((int64_t)b);return ret;}//to be removed? friend CircuitTableInt32 min(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return (a.v>=b.v)?b:a;} friend CircuitTableInt32 max(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return (a.v<=b.v)?b:a;} friend CircuitTableInt32 negabs(CircuitTableInt32 const &a){return min(a,-a);} friend CircuitTableInt32 dotDiv(CircuitTableInt32 const &s, CircuitTableInt32 const &a, CircuitTableInt32 const &t, CircuitTableInt32 const &b, CircuitTableInt32::Divisor const &q) { return CircuitTableInt32(((((int32_t)(((uint64_t)(((int64_t)s.v)*((int64_t)a.v)+((int64_t)t.v)*((int64_t)b.v)))>>q.shift)))*q.multiplicativeInverse)); } friend void dotDivAssign(CircuitTableInt32 &s, CircuitTableInt32 const &a, CircuitTableInt32 const &t, CircuitTableInt32 const &b, CircuitTableInt32::Divisor const &q)//as above but assigning to first argument. This helps the vectorizer. { s.v=((((int32_t)(((uint64_t)(((int64_t)s.v)*((int64_t)a.v)+((int64_t)t.v)*((int64_t)b.v)))>>q.shift)))*q.multiplicativeInverse); } static int MIN(int32_t a, int32_t b) { return (ab)?a:b; } static CircuitTableInt32 computeNegativeBound(CircuitTableInt32 * __restrict__ Ai, int w) { CircuitTableInt32 M=0; CircuitTableInt32 m=0; for(int j=0;j0){s.v>>=1;t.v>>=1;denominatorDivisor.shift--;denominatorDivisor.v>>=1;} CircuitTableInt32 max=0; CircuitTableInt32 min=0; // This is a bound on the absolute value of the sums of products before dividing by the denominator // The calculation can overflow, but since we are casting to unsigned long, this is not a problem. uint64_t positiveResultBoundTimesD=(negabs(t).v*((int64_t)boundA.v)+negabs(s).v*((int64_t)boundB.v)); /* The first case is the case where we know that the intermediate results fit in 32bit before shifting down. * In other words, the shifted-in bits of the results are equal. * In that case the intermediate result can be computed with 32bits. */ if(positiveResultBoundTimesD<((((int64_t)0x40000000)*denominatorDivisor.v)>>denominatorDivisor.shift))//check this carefully {//FAST VERSION // debug<>denominatorDivisor.shift; min.v=MIN(min.v,aa[i].v); max.v=MAX(max.v,aa[i].v); } if(-max<=min)min=-max; return min; } else { /* The next case would be to check if the result is known to fit in 32bit. * In that case intermediate results would be handled in 64 bit, * but the final result would not have to be checked for 32 bit overflows. */ if(positiveResultBoundTimesD<((((int64_t)0x40000000)*denominatorDivisor.v))) { for(int i=0;i>denominatorDivisor.shift)))*denominatorDivisor.multiplicativeInverse); if(max<=aa[i])max=aa[i]; if(aa[i]<=min)min=aa[i]; } if(-max<=min)min=-max; } /* The last case would be to where we don't know that the results will fit in 32 bit. * Since we need to compute max and min anyway, we can compute these quantities in 64bit * and just check if they fit in 32bit at the end. */ else { bool doesOverflow=(((uint32_t)t.v)==0x80000000); int64_t min64=0; int64_t max64=0; for(int i=0;i=0x7fffffff || -v>=0x7fffffff) throw MVMachineIntegerOverflow; return CircuitTableInt32(v); if(v>=0x7fffffff || -v>=0x7fffffff) throw MVMachineIntegerOverflow; return CircuitTableInt32(v); } }

• ## gfanlib/gfanlib_mixedvolume.cpp

 r47a774 static Integer stringToInteger(char const *s) { Integer ret; while(*s)ret=Integer(10)*ret+((*(s++))-'0'); return ret; Integer ret; while(*s)ret=Integer(10)*ret+((*(s++))-'0'); return ret; } virtual const char* what() const throw() { return "Exception: The number of polytopes passed to the mixed volume computation does not match the ambient dimension."; return "Exception: The number of polytopes passed to the mixed volume computation does not match the ambient dimension."; } }MVAmbientDimensionMismatch; static std::vector > convertTuple(std::vector const &tuple) { vector > ret; ret.reserve(tuple.size()); for(int i=0;i a(tuple[i].getHeight(),tuple[i].getWidth()); for(int j=0;j > ret; ret.reserve(tuple.size()); for(int i=0;i a(tuple[i].getHeight(),tuple[i].getWidth()); for(int j=0;j const &tuple, int nthreads, int steps) { //Zero-dimensional volumes are obtained by taking determinants of 0x0 matrices. Such determinants are 1. Hence the volume polynomial is constant 1. The coefficient of the product of the variables is 1. //Alternative explanation. A generic polynomial system with zero variables has exactly 1 solution. if(tuple.size()==0){return Integer(1);} // For now, let's make in an error to pass the wrong number of polytopes. Of course we could define the // mixed volume to be 0 in that case, but the code calling will almost certainly be wrong. It is safer to throw. // for(const auto &x:tuple) // Range based for is not supported until gcc 4.6/ for(auto x=tuple.begin();x!=tuple.end();x++) if(x->getHeight()!=tuple.size()) {throw MVAmbientDimensionMismatch;} std::vector > tuple2=convertTuple(tuple); //              assert(tuple2.size()); typedef SpecializedRTraverser MySpecializedRTraverser; vector T1; int N=nthreads; if(N==0)N=1;// Even if we do not parallelize, we still need one traverser. T1.reserve(N); vector I; for(int i=0;iaborting; if(aborted) { throw MVExceptionRethrow; } mvtyp::Double total; int totalSteps=0; for(int i=0;imixedVolume.toString()<<"Steps:"<<((MySpecializedRTraverser*)(I[i]))->numberOfExpensiveSteps<<"\n"; //debug<<((SpecializedMTraverser*)(I[i]))->T.counter; total.addWithOverflowCheck(((MySpecializedRTraverser*)(I[i]))->mixedVolume); totalSteps+=((MySpecializedRTraverser*)(I[i]))->numberOfExpensiveSteps; } //              debug<<"Total:"<<(int)total.v<<"\n"; //              debug<<"Totalsteps:"< cyclic(int n) { vector ret; for(int i=1;i noon(int n) { vector ret; for(int i=0;i chandra(int n) { vector ret; for(int i=0;i katsura(int n) { n++; vector ret; for(int i=0;i0 ? (n-1-y-i) : -(n-1-y-i)][y]+=1; } m[i][m.getWidth()-1]=1; ret.push_back(m); } ret.push_back(combineLeftRight(IntMatrix::identity(n),IntMatrix(n,1))); return ret; } vector gaukwa(int n) { vector ret; for(int i=0;i<2*n;i++) ret.push_back(combineLeftRight(combineOnTop(IntMatrix::identity(n),i*IntMatrix::identity(n)),IntMatrix(2*n,1))); return ret; } vector eco(int n) { vector ret; for(int i=0;i const &tuple, int nthreads, int steps) { //Zero-dimensional volumes are obtained by taking determinants of 0x0 matrices. Such determinants are 1. Hence the volume polynomial is constant 1. The coefficient of the product of the variables is 1. //Alternative explanation. A generic polynomial system with zero variables has exactly 1 solution. if(tuple.size()==0){return Integer(1);} // For now, let's make in an error to pass the wrong number of polytopes. Of course we could define the // mixed volume to be 0 in that case, but the code calling will almost certainly be wrong. It is safer to throw. // for(const auto &x:tuple) // Range based for is not supported until gcc 4.6/ for(auto x=tuple.begin();x!=tuple.end();x++) if(x->getHeight()!=tuple.size()) {throw MVAmbientDimensionMismatch;} std::vector > tuple2=convertTuple(tuple); //                assert(tuple2.size()); typedef SpecializedRTraverser MySpecializedRTraverser; vector T1; int N=nthreads; if(N==0)N=1;// Even if we do not parallelize, we still need one traverser. T1.reserve(N); vector I; for(int i=0;iaborting; if(aborted) { throw MVExceptionRethrow; } mvtyp::Double total; int totalSteps=0; for(int i=0;imixedVolume.toString()<<"Steps:"<<((MySpecializedRTraverser*)(I[i]))->numberOfExpensiveSteps<<"\n"; //debug<<((SpecializedMTraverser*)(I[i]))->T.counter; total.addWithOverflowCheck(((MySpecializedRTraverser*)(I[i]))->mixedVolume); totalSteps+=((MySpecializedRTraverser*)(I[i]))->numberOfExpensiveSteps; } //                debug<<"Total:"<<(int)total.v<<"\n"; //                debug<<"Totalsteps:"< cyclic(int n) { vector ret; for(int i=1;i noon(int n) { vector ret; for(int i=0;i chandra(int n) { vector ret; for(int i=0;i katsura(int n) { n++; vector ret; for(int i=0;i0 ? (n-1-y-i) : -(n-1-y-i)][y]+=1; } m[i][m.getWidth()-1]=1; ret.push_back(m); } ret.push_back(combineLeftRight(IntMatrix::identity(n),IntMatrix(n,1))); return ret; } vector gaukwa(int n) { vector ret; for(int i=0;i<2*n;i++) ret.push_back(combineLeftRight(combineOnTop(IntMatrix::identity(n),i*IntMatrix::identity(n)),IntMatrix(2*n,1))); return ret; } vector eco(int n) { vector ret; for(int i=0;i
• ## gfanlib/gfanlib_mixedvolume.h

 r47a774 namespace gfan{ Integer mixedVolume(std::vector const &tuple, int nthreads=0, int steps=500); namespace MixedVolumeExamples{ std::vector cyclic(int n); std::vector noon(int n); std::vector chandra(int n); std::vector katsura(int n); std::vector gaukwa(int n); std::vector eco(int n); } Integer mixedVolume(std::vector const &tuple, int nthreads=0, int steps=500); namespace MixedVolumeExamples{ std::vector cyclic(int n); std::vector noon(int n); std::vector chandra(int n); std::vector katsura(int n); std::vector gaukwa(int n); std::vector eco(int n); } }

• ## gfanlib/gfanlib_paralleltraverser.h

 r47a774 { public: bool aborting;                                  // Added by Anders void abort(){aborting=true;}    // Added by Anders Traverser():aborting(false){}   // Added by Anders bool aborting;                      // Added by Anders void abort(){aborting=true;}        // Added by Anders Traverser():aborting(false){}       // Added by Anders // Virtual destructor virtual ~Traverser( void ) {};
• ## gfanlib/gfanlib_symmetriccomplex.cpp

 r47a774 int SymmetricComplex::getLinDim()const { return linealitySpace.getHeight(); return linealitySpace.getHeight(); }

• ## gfanlib/gfanlib_traversal.h

 r47a774 { public: /** * Go to the cone which is connected to the current facet through the ridge in direction ray. * The "ridge" is a relative interior point of the ridge. */ virtual void changeCone(ZVector const &ridgeVector, ZVector const &rayVector)=0; /** * Go to the cone which is connected to the current facet through the ridge in direction ray. * The "ridge" is a relative interior point of the ridge. */ virtual void changeCone(ZVector const &ridgeVector, ZVector const &rayVector)=0; /** * Compute the link of the fan in the ridge given by the vector ridge IS THIS A FACET NORMAL OR AN INTERIOR POINT? * This gives a list of symmetry invariant points under the actions keeping the link fixed. */ virtual std::list link(ZVector const &ridgeVector)=0; virtual ZCone & refToPolyhedralCone()=0; virtual std::list link(ZVector const &ridgeVector)=0; virtual ZCone & refToPolyhedralCone()=0; /** * If there is no cone state data for the traverser, half of the changeCone() calls can be avoided. * That this is a valid of optimization for the ConeTraverser is indicated returning true in the following function. */ virtual bool hasNoState()const; virtual bool hasNoState()const; }; { public: virtual bool process(FanTraverser &traverser)=0; virtual bool process(FanTraverser &traverser)=0; }; class FanBuilder : public Target { ZFan coneCollection; ZFan coneCollection; public: ZFan const &getFanRef(){return coneCollection;} FanBuilder(int n, SymmetryGroup const &sym); bool process(FanTraverser &Traverser); ZFan const &getFanRef(){return coneCollection;} FanBuilder(int n, SymmetryGroup const &sym); bool process(FanTraverser &Traverser); };
• ## gfanlib/gfanlib_tropicalhomotopy.h

 r47a774 class Flags{ public: static const bool computeDotProductInMatrix=true; static const bool computeDotProductInMatrix=true; }; static Matrix simplex(int n, mvtyp d) { Matrix ret(n,n+1); for(int i=0;i ret(n,n+1); for(int i=0;i class SingleTropicalHomotopyTraverser{ class InequalityComparisonResult{//actual comparison functions were moved to the InequalityTable public: bool empty; int configurationIndex;//used for finding best int columnIndex;//used for finding best }; class InequalityTable //circuit table // This table has been moved inside the IntegersectionTraverser simply because it is used nowhere else and is specific to mixed cells in Cayley configurations. { /* All methods are marked to show if they can overflow without throwing/asserting. * Overflowing methods at the moment are: *  replaceFirstOrSecond:       subroutine calls may overflow (dotDivVector) *  compareInequalities:           only if target entries are too big *  dotVector:                     only if target entries are too big *  setChoicesFromEarlierHomotopy: only if tuple entries are too big */ std::vector > tuple; std::vector offsets; std::vector > choices; Matrix A;//one row for each potential inequality, to store entries with indices in chosen Vector tempA; Vector Abounds;// a negative bound for each row of A, bounding the absolute value of the rows; std::vector svec;//used locally int subconfigurationIndex; mvtyp denominator; int m; int k; bool isLegalIndex(int subconfigurationIndex, int columnIndex)const {// Cannot overflow return choices[subconfigurationIndex].first!=columnIndex && choices[subconfigurationIndex].second!=columnIndex; } mvtyp dotVector(int subconfigurationIndex, int columnIndex, Vector const &target, int onlyK=-1)const {   // May overflow if entries of target are too big. //if onlyK!=-1 then only the onlyKth subconfiguration is considered mvtypDouble ret; if(onlyK!=-1) { if(onlyK==subconfigurationIndex) { int i=subconfigurationIndex; ret+=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].second+offsets[i])); ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i])); ret-=extendedMultiplication(denominator,target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));// the multiplication can be merged with multiplication above except that that could cause and overflow. ret+=extendedMultiplication(denominator,target.UNCHECKEDACCESS(columnIndex+offsets[i])); return ret.castToSingle(); } else { int i=onlyK; if(target.UNCHECKEDACCESS((choices)[i].first+offsets[i]).isNonZero()) { ret+=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].second+offsets[i])); ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i])); } return ret.castToSingle(); } } for(int i=0;i const &target, int onlyK=-1) {// Cannot overflow int J=0; for(int i=0;ib]=a; v[b>a]=b; return 2; } int sort3uniquely(int *v, int a, int b, int c)const//a and b and c different {// Cannot overflow v[(a>b)+int(a>c)]=a; v[(b>a)+int(b>c)]=b; v[(c>a)+int(c>b)]=c; return 3; } int sort4uniquely(int *v, int a, int b, int c, int d)const// a and b different and different from c and d, but c may equal d {// Cannot overflow if(c!=d) { v[(a>b)+int(a>c)+int(a>d)]=a; v[(b>a)+int(b>c)+int(b>d)]=b; v[(c>a)+int(c>b)+int(c>d)]=c; v[(d>a)+int(d>b)+int(d>c)]=d; return 4; } else return sort3uniquely(v,a,b,c); } bool compareReverseLexicographicInverted(int i1, int j1, int i2, int j2, mvtyp s1, mvtyp s2)const//s1 and s2 are always negative {// cannot overflow for(int i=0;i0; } } int indices[4]; int F=choices[i].first; int S=choices[i].second; int toCheck; if(i1==i) if(i2==i) toCheck=sort4uniquely(indices,F,S,j1,j2); else toCheck=sort3uniquely(indices,F,S,j1); else if(i2==i) toCheck=sort3uniquely(indices,F,S,j2); else toCheck=sort2uniquely(indices,F,S); for(int J=0;J0) return true; else if(temp<0) return false; } } return false; } mvtyp getVolume() {// Cannot overflow return denominator; } void replaceFirstOrSecond(bool first, int subconfigurationIndex, int newIndex, Vector const &target)__attribute__ ((always_inline))//updates the inequality table according to replacing first or second by newIndex in the subconfigurationIndex'th configuration {// Cannot overflow int newIndex2=newIndex;for(int i=0;i const &target){replaceFirstOrSecond(true,subconfigurationIndex,newIndex,target);} void replaceSecond(int subconfigurationIndex, int newIndex, Vector const &target){replaceFirstOrSecond(false,subconfigurationIndex,newIndex,target);} InequalityTable(std::vector > const &tuple_, int subconfigurationIndex_): tempA(tuple_.size()+Flags::computeDotProductInMatrix), tuple(tuple_), choices(tuple_.size()), subconfigurationIndex(subconfigurationIndex_), offsets(tuple_.size()) {// Cannot overflow k=tuple.size(); m=0; for(int i=0;i(k+Flags::computeDotProductInMatrix,m); {int offset=0;for(int i=0;i(k+Flags::computeDotProductInMatrix); } void setChoicesInitially() {// Cannot overflow //THIS WILL ONLY WORK FOR THE STARTING CONFIGURATION //sets denominators,A and choices (these must have been initialized to the right sizes) for(int i=0;i (i+0,i+1); for(int i=0;ii+1) for(int ii=i;ii const &target, int onlyK=-1) {// Can only overflow if target entries are too big. Actually it seems that it is not this function that will overflow but dotVector. bool empty=true; int bestConfigurationIndex=-1; int bestColumnIndex=-1; mvtyp targetDotBest=0; for(int i=0;i::const_RowRef Ak=const_cast&>(A)[k]; gfan::Matrix::const_RowRef Ak=const_cast&>(A)[k]; int offsetsi=offsets[i]; int tupleiwidth=tuple[i].getWidth(); if(onlyK!=-1)if(i!=onlyK)continue; for(int j=0;j0); } void setChoicesFromEarlierHomotopy(InequalityTable const &parent, mvtyp degreeScaling, Vector const &target) {       // Overflows may happen, but only if the tuple entries are too big. // Notice that the code below has overflow checks everywhere, except in the loop marked with a comment. //sets denominators,A and choices (these must have been initialized to the right sizes //columns have been added to configuration this->subconfigurationIndex //for that reason we need to introduce new circuits. Old circuits are preserved. //chioices are "relative" so no update is needed. choices=parent.choices; int numberToDrop=(subconfigurationIndex!=0) ? numberToDrop=k+1 : 0; choices[subconfigurationIndex-1].first-=numberToDrop; choices[subconfigurationIndex-1].second-=numberToDrop; denominator=parent.denominator; int offsetOld=0; int offsetNew=0; for(int i=0;i=0) A.UNCHECKEDACCESS(subconfigurationIndex,offsetNew+j).msubWithOverflowChecking(tuple[i].UNCHECKEDACCESS(b,j),denominator); } for(int a=0;a > choices; Vector target; bool useFirstChanged; bool useSecondChanged; std::vector stack; int eliminatedK; int eliminatedKOffset; std::vector > tuple; std::vector offsets; int m; InequalityComparisonResult result; InequalityTable inequalityTable; void constructInequalityTableFromParent(InequalityTable const &parentTable, mvtyp degreeScaling) { inequalityTable.setChoicesFromEarlierHomotopy(parentTable, degreeScaling, target); } void constructInequalityTableInitially(mvtyp degreeScaling) { std::vector > tempTuple;for(int i=0;i(tuple.size(),1)); InequalityTable tempTable(tempTuple,-1); tempTable.setChoicesInitially(); constructInequalityTableFromParent(tempTable,degreeScaling); } SingleTropicalHomotopyTraverser(std::vector > const &tuple_, int m_, std::vector > const &choices_, Vector const &target_, int eliminatedK_): choices(choices_), target(target_), eliminatedK(eliminatedK_), tuple(tuple_), m(m_), inequalityTable(tuple,eliminatedK), offsets(tuple_.size()) { eliminatedKOffset=0;for(int i=0;i > tuple; std::vector offsets; std::vector > choices; Matrix A;//one row for each potential inequality, to store entries with indices in chosen Vector tempA; Vector Abounds;// a negative bound for each row of A, bounding the absolute value of the rows; std::vector svec;//used locally int subconfigurationIndex; mvtyp denominator; int m; int k; bool isLegalIndex(int subconfigurationIndex, int columnIndex)const {// Cannot overflow return choices[subconfigurationIndex].first!=columnIndex && choices[subconfigurationIndex].second!=columnIndex; } mvtyp dotVector(int subconfigurationIndex, int columnIndex, Vector const &target, int onlyK=-1)const {   // May overflow if entries of target are too big. //if onlyK!=-1 then only the onlyKth subconfiguration is considered mvtypDouble ret; if(onlyK!=-1) { if(onlyK==subconfigurationIndex) { int i=subconfigurationIndex; ret+=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].second+offsets[i])); ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i])); ret-=extendedMultiplication(denominator,target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));// the multiplication can be merged with multiplication above except that that could cause and overflow. ret+=extendedMultiplication(denominator,target.UNCHECKEDACCESS(columnIndex+offsets[i])); return ret.castToSingle(); } else { int i=onlyK; if(target.UNCHECKEDACCESS((choices)[i].first+offsets[i]).isNonZero()) { ret+=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].second+offsets[i])); ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i])); } return ret.castToSingle(); } } for(int i=0;i const &target, int onlyK=-1) {// Cannot overflow int J=0; for(int i=0;ib]=a; v[b>a]=b; return 2; } int sort3uniquely(int *v, int a, int b, int c)const//a and b and c different {// Cannot overflow v[(a>b)+int(a>c)]=a; v[(b>a)+int(b>c)]=b; v[(c>a)+int(c>b)]=c; return 3; } int sort4uniquely(int *v, int a, int b, int c, int d)const// a and b different and different from c and d, but c may equal d {// Cannot overflow if(c!=d) { v[(a>b)+int(a>c)+int(a>d)]=a; v[(b>a)+int(b>c)+int(b>d)]=b; v[(c>a)+int(c>b)+int(c>d)]=c; v[(d>a)+int(d>b)+int(d>c)]=d; return 4; } else return sort3uniquely(v,a,b,c); } bool compareReverseLexicographicInverted(int i1, int j1, int i2, int j2, mvtyp s1, mvtyp s2)const//s1 and s2 are always negative {// cannot overflow for(int i=0;i0; } } int indices[4]; int F=choices[i].first; int S=choices[i].second; int toCheck; if(i1==i) if(i2==i) toCheck=sort4uniquely(indices,F,S,j1,j2); else toCheck=sort3uniquely(indices,F,S,j1); else if(i2==i) toCheck=sort3uniquely(indices,F,S,j2); else toCheck=sort2uniquely(indices,F,S); for(int J=0;J0) return true; else if(temp<0) return false; } } return false; } mvtyp getVolume() {// Cannot overflow return denominator; } void replaceFirstOrSecond(bool first, int subconfigurationIndex, int newIndex, Vector const &target)__attribute__ ((always_inline))//updates the inequality table according to replacing first or second by newIndex in the subconfigurationIndex'th configuration {// Cannot overflow int newIndex2=newIndex;for(int i=0;i const &target){replaceFirstOrSecond(true,subconfigurationIndex,newIndex,target);} void replaceSecond(int subconfigurationIndex, int newIndex, Vector const &target){replaceFirstOrSecond(false,subconfigurationIndex,newIndex,target);} InequalityTable(std::vector > const &tuple_, int subconfigurationIndex_): tempA(tuple_.size()+Flags::computeDotProductInMatrix), tuple(tuple_), choices(tuple_.size()), subconfigurationIndex(subconfigurationIndex_), offsets(tuple_.size()) {// Cannot overflow k=tuple.size(); m=0; for(int i=0;i(k+Flags::computeDotProductInMatrix,m); {int offset=0;for(int i=0;i(k+Flags::computeDotProductInMatrix); } void setChoicesInitially() {// Cannot overflow //THIS WILL ONLY WORK FOR THE STARTING CONFIGURATION //sets denominators,A and choices (these must have been initialized to the right sizes) for(int i=0;i (i+0,i+1); for(int i=0;ii+1) for(int ii=i;ii const &target, int onlyK=-1) {// Can only overflow if target entries are too big. Actually it seems that it is not this function that will overflow but dotVector. bool empty=true; int bestConfigurationIndex=-1; int bestColumnIndex=-1; mvtyp targetDotBest=0; for(int i=0;i::const_RowRef Ak=const_cast&>(A)[k]; gfan::Matrix::const_RowRef Ak=const_cast&>(A)[k]; int offsetsi=offsets[i]; int tupleiwidth=tuple[i].getWidth(); if(onlyK!=-1)if(i!=onlyK)continue; for(int j=0;j0); } void setChoicesFromEarlierHomotopy(InequalityTable const &parent, mvtyp degreeScaling, Vector const &target) {        // Overflows may happen, but only if the tuple entries are too big. // Notice that the code below has overflow checks everywhere, except in the loop marked with a comment. //sets denominators,A and choices (these must have been initialized to the right sizes //columns have been added to configuration this->subconfigurationIndex //for that reason we need to introduce new circuits. Old circuits are preserved. //chioices are "relative" so no update is needed. choices=parent.choices; int numberToDrop=(subconfigurationIndex!=0) ? numberToDrop=k+1 : 0; choices[subconfigurationIndex-1].first-=numberToDrop; choices[subconfigurationIndex-1].second-=numberToDrop; denominator=parent.denominator; int offsetOld=0; int offsetNew=0; for(int i=0;i=0) A.UNCHECKEDACCESS(subconfigurationIndex,offsetNew+j).msubWithOverflowChecking(tuple[i].UNCHECKEDACCESS(b,j),denominator); } for(int a=0;a > choices; Vector target; bool useFirstChanged; bool useSecondChanged; std::vector stack; int eliminatedK; int eliminatedKOffset; std::vector > tuple; std::vector offsets; int m; InequalityComparisonResult result; InequalityTable inequalityTable; void constructInequalityTableFromParent(InequalityTable const &parentTable, mvtyp degreeScaling) { inequalityTable.setChoicesFromEarlierHomotopy(parentTable, degreeScaling, target); } void constructInequalityTableInitially(mvtyp degreeScaling) { std::vector > tempTuple;for(int i=0;i(tuple.size(),1)); InequalityTable tempTable(tempTuple,-1); tempTable.setChoicesInitially(); constructInequalityTableFromParent(tempTable,degreeScaling); } SingleTropicalHomotopyTraverser(std::vector > const &tuple_, int m_, std::vector > const &choices_, Vector const &target_, int eliminatedK_): choices(choices_), target(target_), eliminatedK(eliminatedK_), tuple(tuple_), m(m_), inequalityTable(tuple,eliminatedK), offsets(tuple_.size()) { eliminatedKOffset=0;for(int i=0;i class TropicalRegenerationTraverser{ // The following class is an attempt to separate homotopy data from the traversal logic. class Data{ static mvtyp degree(Matrix const &m)//assumes entries of m are non-negative { mvtyp ret=0; for(int i=0;i > targets; std::vector > tuple; std::vector > > tuples; Vector degrees; bool isFiniteIndex(int level, int index) { return index>=tuple[0].getHeight()+1; } std::vector > produceIthTuple(int i) { int n=tuple[0].getHeight(); std::vector > ret; for(int j=0;j(n,degree(tuple[j])),tuple[j])); if(j>i)ret.push_back(simplex(n,1)); } return ret; } Data(std::vector > const &tuple_):tuple(tuple_) { int n=tuple[0].getHeight(); {//adjusting to positive orthant for(int i=0;i targ; for(int j=0;j::allOnes(n+1),Vector(tuple[i].getWidth()))); else targ=concatenation(targ,Vector(tuples[i][j].getWidth())); } targets.push_back(targ); } }; std::vector > firstIntersection() { std::vector > ret; for(int i=0;i(i+0,i+1)); return ret; } void castToNextLevel(std::vector > const &choices, int i, int S, std::vector > &ret) { assert(ret.size()==choices.size()); for(int j=0;j=S); assert(ret[i].second>=S); ret[i].first-=S; ret[i].second-=S; } }; static int cayleyConfigurationWidth(std::vector > const &tuple) { int m2=0; for(int i=0;i MySingleTropicalHomotopyTraverser; std::vector traversers; Data fullData; int level; bool deadEnd; bool isLevelLeaf; bool isSolutionVertex; std::vector isLevelLeafStack; TropicalRegenerationTraverser(std::vector > const &tuple_): fullData(tuple_),counter(0),depth(0) { assert(tuple_.size()); for(int i=0;i class SpecializedRTraverser: public Traverser { public: typedef TropicalRegenerationTraverser MyTropicalRegenerationTraverser; MyTropicalRegenerationTraverser T; mvtypDouble mixedVolume; int numberOfExpensiveSteps; SpecializedRTraverser(std::vector > const &tuple_): mixedVolume(), numberOfExpensiveSteps(0), T(tuple_) //Constructor my throw if input is too big. { numberOfExpensiveSteps++; T.findOutgoingAndProcess(false); } int  getEdgeCountNext( void ) { if(!aborting) { try{ return T.numberOfChildren(); } catch(...){abort();} } return 0; } int  moveToNext( int   index, bool  collect_info ) { if(!aborting) { try{ T.goToNthChild(index); numberOfExpensiveSteps++; T.findOutgoingAndProcess(false); } catch(...){abort();} } return 0; } void  moveToPrev( int  index ) { if(!aborting) { try{ T.goBack(); //index ignored } catch(...){abort();} } } void  collectInfo( void ) { if(!aborting) { try{ if(T.isSolutionVertex) mixedVolume.addWithOverflowCheck(T.traversers[T.level].inequalityTable.getVolume().extend()); } catch(...){abort();} } } void  printState( void ) { T.print(); } }; class TropicalRegenerationTraverser{ // The following class is an attempt to separate homotopy data from the traversal logic. class Data{ static mvtyp degree(Matrix const &m)//assumes entries of m are non-negative { mvtyp ret=0; for(int i=0;i > targets; std::vector > tuple; std::vector > > tuples; Vector degrees; bool isFiniteIndex(int level, int index) { return index>=tuple[0].getHeight()+1; } std::vector > produceIthTuple(int i) { int n=tuple[0].getHeight(); std::vector > ret; for(int j=0;j(n,degree(tuple[j])),tuple[j])); if(j>i)ret.push_back(simplex(n,1)); } return ret; } Data(std::vector > const &tuple_):tuple(tuple_) { int n=tuple[0].getHeight(); {//adjusting to positive orthant for(int i=0;i targ; for(int j=0;j::allOnes(n+1),Vector(tuple[i].getWidth()))); else targ=concatenation(targ,Vector(tuples[i][j].getWidth())); } targets.push_back(targ); } }; std::vector > firstIntersection() { std::vector > ret; for(int i=0;i(i+0,i+1)); return ret; } void castToNextLevel(std::vector > const &choices, int i, int S, std::vector > &ret) { assert(ret.size()==choices.size()); for(int j=0;j=S); assert(ret[i].second>=S); ret[i].first-=S; ret[i].second-=S; } }; static int cayleyConfigurationWidth(std::vector > const &tuple) { int m2=0; for(int i=0;i MySingleTropicalHomotopyTraverser; std::vector traversers; Data fullData; int level; bool deadEnd; bool isLevelLeaf; bool isSolutionVertex; std::vector isLevelLeafStack; TropicalRegenerationTraverser(std::vector > const &tuple_): fullData(tuple_),counter(0),depth(0) { assert(tuple_.size()); for(int i=0;i class SpecializedRTraverser: public Traverser { public: typedef TropicalRegenerationTraverser MyTropicalRegenerationTraverser; MyTropicalRegenerationTraverser T; mvtypDouble mixedVolume; int numberOfExpensiveSteps; SpecializedRTraverser(std::vector > const &tuple_): mixedVolume(), numberOfExpensiveSteps(0), T(tuple_) //Constructor my throw if input is too big. { numberOfExpensiveSteps++; T.findOutgoingAndProcess(false); } int  getEdgeCountNext( void ) { if(!aborting) { try{ return T.numberOfChildren(); } catch(...){abort();} } return 0; } int  moveToNext( int   index, bool  collect_info ) { if(!aborting) { try{ T.goToNthChild(index); numberOfExpensiveSteps++; T.findOutgoingAndProcess(false); } catch(...){abort();} } return 0; } void  moveToPrev( int  index ) { if(!aborting) { try{ T.goBack(); //index ignored } catch(...){abort();} } } void  collectInfo( void ) { if(!aborting) { try{ if(T.isSolutionVertex) mixedVolume.addWithOverflowCheck(T.traversers[T.level].inequalityTable.getVolume().extend()); } catch(...){abort();} } } void  printState( void ) { T.print(); } }; }
• ## gfanlib/gfanlib_vector.h

 r47a774 std::string toString()const { std::stringstream f; f<<*this; return f.str(); std::stringstream f; f<<*this; return f.str(); }
• ## gfanlib/gfanlib_zcone.h

 r47a774 * Returns true if cddlib is needed for the ZCone implementation. */ bool isCddlibRequired(); /** * Only call this function if gfanlib is the only code in your program using cddlib. * Should be paired with a deinitializeCddlibIfRequired() call. * Calling the function repeatedly may cause memory leaks even if deinitializeCddlibIfRequired() is also called. */ void initializeCddlibIfRequired(); /** * This function may do nothing. */ void deinitializeCddlibIfRequired(); bool isCddlibRequired(); /** * Only call this function if gfanlib is the only code in your program using cddlib. * Should be paired with a deinitializeCddlibIfRequired() call. * Calling the function repeatedly may cause memory leaks even if deinitializeCddlibIfRequired() is also called. */ void initializeCddlibIfRequired(); /** * This function may do nothing. */ void deinitializeCddlibIfRequired(); /**
• ## gfanlib/test.cpp

 r47a774 catch (...) { cerr<<"Error - most likely an integer overflow."<
Note: See TracChangeset for help on using the changeset viewer.