Changeset 15813d in git


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

Legend:

Unmodified
Added
Removed

  • gfanlib/gfanlib_circuittableint.h

    r47a774 r15813d  
    3232class CircuitTableInt32{
    3333 public:
    34         class Divisor{
    35         public:
    36                 int32_t v;
    37                 int shift;
    38                 int32_t multiplicativeInverse;
    39                 Divisor(CircuitTableInt32 const &a)// for exact division
    40                 { // A good place to read about these tricks seems to be the book "Hacker's Delight" by Warren.
    41                         v=a.v;
    42                         shift=0;
    43                         uint32_t t=v;
    44                         assert(t);
    45                         while(!(t&1)){t>>=      1;shift++;}
    46                         uint32_t inverse=t;
    47                         while(t*inverse!=1)inverse*=2-t*inverse;
    48                         multiplicativeInverse=inverse;
    49                 }
    50         };
    51         class Double{
    52         public:
    53                 int64_t v;
    54                 Double():v(0){};
    55                 Double(int64_t a):v(a){};
    56                 Double &operator+=(Double a){v+=a.v;return *this;}
    57                 Double &operator-=(Double a){v-=a.v;return *this;}
    58                 CircuitTableInt32 castToSingle()const;
    59                 bool isZero()const{return v==0;}
    60                 bool isNegative()const{return v<0;}
    61                 Double operator-()const{Double ret;ret.v=-v;return ret;}
    62                 friend Double operator-(Double const &a, Double const &b){return Double(a.v-b.v);}
    63                 friend Double operator+(Double const &a, Double const &b){return Double(a.v+b.v);}
    64                 Double &addWithOverflowCheck(Double const &a)
    65                 {
    66                         if(a.v<0 || v<0 || a.v>1000000000000000000 || v>1000000000000000000) throw MVMachineIntegerOverflow;
    67                         v+=a.v;
    68                         return *this;
    69                 }
    70                 std::string toString()const{std::stringstream s;s<<v;return s.str();}
    71         };
     34        class Divisor{
     35        public:
     36                int32_t v;
     37                int shift;
     38                int32_t multiplicativeInverse;
     39                Divisor(CircuitTableInt32 const &a)// for exact division
     40                { // A good place to read about these tricks seems to be the book "Hacker's Delight" by Warren.
     41                        v=a.v;
     42                        shift=0;
     43                        uint32_t t=v;
     44                        assert(t);
     45                        while(!(t&1)){t>>=        1;shift++;}
     46                        uint32_t inverse=t;
     47                        while(t*inverse!=1)inverse*=2-t*inverse;
     48                        multiplicativeInverse=inverse;
     49                }
     50        };
     51        class Double{
     52        public:
     53                int64_t v;
     54                Double():v(0){};
     55                Double(int64_t a):v(a){};
     56                Double &operator+=(Double a){v+=a.v;return *this;}
     57                Double &operator-=(Double a){v-=a.v;return *this;}
     58                CircuitTableInt32 castToSingle()const;
     59                bool isZero()const{return v==0;}
     60                bool isNegative()const{return v<0;}
     61                Double operator-()const{Double ret;ret.v=-v;return ret;}
     62                friend Double operator-(Double const &a, Double const &b){return Double(a.v-b.v);}
     63                friend Double operator+(Double const &a, Double const &b){return Double(a.v+b.v);}
     64                Double &addWithOverflowCheck(Double const &a)
     65                {
     66                        if(a.v<0 || v<0 || a.v>1000000000000000000 || v>1000000000000000000) throw MVMachineIntegerOverflow;
     67                        v+=a.v;
     68                        return *this;
     69                }
     70                std::string toString()const{std::stringstream s;s<<v;return s.str();}
     71        };
    7272 private:
    7373public:
    74         int32_t v;
    75         friend CircuitTableInt32 operator+(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return CircuitTableInt32(a.v+b.v);}
    76         friend CircuitTableInt32 operator-(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return CircuitTableInt32(a.v-b.v);}
    77         friend CircuitTableInt32 operator*(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return CircuitTableInt32(a.v*b.v);}
     74        int32_t v;
     75        friend CircuitTableInt32 operator+(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return CircuitTableInt32(a.v+b.v);}
     76        friend CircuitTableInt32 operator-(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return CircuitTableInt32(a.v-b.v);}
     77        friend CircuitTableInt32 operator*(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return CircuitTableInt32(a.v*b.v);}
    7878 public:
    79         // In the code products of CircuitTableInt32 objects are summed. To avoid overflows one of the factors must "fit in half" and the "number of summands" may not exceed a certain number. The bounds are specified in the following functions:
    80         bool fitsInHalf()const{return v>-30000 && v<30000;}// Is it better to allow only non-negative numbers?
    81         static bool isSuitableSummationNumber(int numberOfSummands){return numberOfSummands<30000;}
    82         CircuitTableInt32(){v=0;};
    83         CircuitTableInt32(CircuitTableInt32 const &m):v(m.v){}
    84         CircuitTableInt32(int32_t val):v(val){};
    85         class CircuitTableInt32 &operator=(int32_t a){v=a;return *this;}
    86         CircuitTableInt32 operator-()const{CircuitTableInt32 ret;ret.v=-v;return ret;}
    87         CircuitTableInt32 &operator-=(CircuitTableInt32 a){v-=a.v;return *this;}
    88         CircuitTableInt32 &operator+=(CircuitTableInt32 a){v+=a.v;return *this;}
    89         CircuitTableInt32 &operator*=(CircuitTableInt32 a){v*=a.v;return *this;}
    90         friend bool operator<=(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return a.v<=b.v;}
    91         friend bool operator==(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return a.v==b.v;}
    92         bool isZero()const{return v==0;}
    93         bool isOne()const{return v==1;}
    94         bool isNonZero()const{return v!=0;}
    95         bool isNegative()const{return v<0;}
    96         bool isPositive()const{return v>0;}
    97         friend int determinantSign1(CircuitTableInt32 const &a, CircuitTableInt32 const &c, CircuitTableInt32 const &b, CircuitTableInt32 const &d)//NOTICE MIXED ORDER. MUST WORK BY EXTENDING
    98         {
    99                 int64_t r=((int64_t)a.v)*((int64_t)c.v)-((int64_t)b.v)*((int64_t)d.v);
    100                 if(r>0)return 1;
    101                 if(r<0)return -1;
    102                 return 0;
    103         }
    104         friend int determinantSign2(CircuitTableInt32::Double const &a, CircuitTableInt32 const &c, CircuitTableInt32::Double const &b, CircuitTableInt32 const &d)//NOTICE MIXED ORDER. MUST WORK BY EXTENDING
    105         {
    106                 int64_t r1=(a.v)*((int64_t)c.v);
    107                 int64_t r2=(b.v)*((int64_t)d.v);
    108                 if(r1>r2)return 1;
    109                 if(r1<r2)return -1;
    110                 return 0;
    111         }
    112         std::string toString()const{std::stringstream s;s<<v;return s.str();}
    113         Double extend()const{Double ret;ret.v=v;return ret;}
    114         CircuitTableInt32 &maddWithOverflowChecking(CircuitTableInt32 const &a, CircuitTableInt32 const&b){Double t=this->extend();t+=extendedMultiplication(a,b);*this=t.castToSingle();return *this;}
    115         CircuitTableInt32 &msubWithOverflowChecking(CircuitTableInt32 const &a, CircuitTableInt32 const&b){Double s=this->extend();s-=extendedMultiplication(a,b);*this=s.castToSingle();return *this;}
    116         CircuitTableInt32 &subWithOverflowChecking(CircuitTableInt32 const &a){Double t=this->extend();t-=a.extend();*this=t.castToSingle();return *this;}
    117         CircuitTableInt32 &addWithOverflowChecking(CircuitTableInt32 const &a){Double t=this->extend();t+=a.extend();*this=t.castToSingle();return *this;}
    118         CircuitTableInt32 negatedWithOverflowChecking()const{return (-extend()).castToSingle();}
    119         friend Double extendedMultiplication(CircuitTableInt32 const &a, CircuitTableInt32 const &b){Double ret;ret.v=((int64_t)a.v)*((int64_t)b.v);return ret;}
    120         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?
    121         friend CircuitTableInt32 min(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return (a.v>=b.v)?b:a;}
    122         friend CircuitTableInt32 max(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return (a.v<=b.v)?b:a;}
    123         friend CircuitTableInt32 negabs(CircuitTableInt32 const &a){return min(a,-a);}
    124         friend CircuitTableInt32 dotDiv(CircuitTableInt32 const &s, CircuitTableInt32 const &a, CircuitTableInt32 const &t, CircuitTableInt32 const &b, CircuitTableInt32::Divisor const &q)
    125         {
    126                 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));
    127         }
    128         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.
    129         {
    130                 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);
    131         }
    132         static int MIN(int32_t a, int32_t b)
    133         {
    134                 return (a<b)?a:b;
    135         }
    136         static int MAX(int32_t a, int32_t b)
    137         {
    138                 return (a>b)?a:b;
    139         }
    140         static CircuitTableInt32 computeNegativeBound(CircuitTableInt32 * __restrict__ Ai, int w)
    141         {
    142                 CircuitTableInt32 M=0;
    143                 CircuitTableInt32 m=0;
    144                 for(int j=0;j<w;j++)
    145                 {
    146                         m.v=MIN(m.v,Ai[j].v);
    147                         M.v=MAX(M.v,Ai[j].v);
    148                 }
    149                 return min(m,-M);
    150         }
    151         static CircuitTableInt32 quickNoShiftBounded(CircuitTableInt32 * __restrict__ a, CircuitTableInt32 * __restrict__ b, CircuitTableInt32 s, CircuitTableInt32 t, CircuitTableInt32::Divisor denominatorDivisor,int c)
    152         {
    153                 CircuitTableInt32 *aa=a;
    154                 CircuitTableInt32 *bb=b;
    155                 CircuitTableInt32 max=0;
    156             CircuitTableInt32 min=0;
    157             CircuitTableInt32 ret=0;
    158             for(int i=0;i<c;i++)
    159               {
    160                 aa[i].v=((s.v*denominatorDivisor.multiplicativeInverse)*aa[i].v+
    161                                 (t.v*denominatorDivisor.multiplicativeInverse)*bb[i].v);
    162                 min.v=MIN(min.v,aa[i].v);
    163                 max.v=MAX(max.v,aa[i].v);
    164               }
    165             if(-max<=min)min=-max;
    166             ret=min;
    167             return ret;
    168         }
    169         static CircuitTableInt32 dotDivVector(CircuitTableInt32 * __restrict__ aa, CircuitTableInt32 * __restrict__ bb, CircuitTableInt32 s, CircuitTableInt32 t, CircuitTableInt32::Divisor denominatorDivisor,int c, CircuitTableInt32 boundA, CircuitTableInt32 boundB)__attribute__ ((always_inline))//assigning to first argument. The return negative of the return value is an upper bound on the absolute values of the produced entries
    170         {
    171                 while(((s.v|t.v)&1)==0 && denominatorDivisor.shift>0){s.v>>=1;t.v>>=1;denominatorDivisor.shift--;denominatorDivisor.v>>=1;}
    172 
    173                 CircuitTableInt32 max=0;
    174                 CircuitTableInt32 min=0;
    175 
    176           // This is a bound on the absolute value of the sums of products before dividing by the denominator
    177           // The calculation can overflow, but since we are casting to unsigned long, this is not a problem.
    178           uint64_t positiveResultBoundTimesD=(negabs(t).v*((int64_t)boundA.v)+negabs(s).v*((int64_t)boundB.v));
    179 
    180           /* The first case is the case where we know that the intermediate results fit in 32bit before shifting down.
    181            * In other words, the shifted-in bits of the results are equal.
    182            * In that case the intermediate result can be computed with 32bits.
    183            */
    184           if(positiveResultBoundTimesD<((((int64_t)0x40000000)*denominatorDivisor.v)>>denominatorDivisor.shift))//check this carefully
    185           {//FAST VERSION
    186                   // debug<<s.v<<" "<<t.v<<" "<<denominatorDivisor.v<<" "<<denominatorDivisor.shift<<"\n";
    187                   if(denominatorDivisor.shift==0)return quickNoShiftBounded(aa,bb,s,t,denominatorDivisor,c);
    188                   for(int i=0;i<c;i++)
    189                   {
    190                           aa[i].v=((s.v*denominatorDivisor.multiplicativeInverse)*aa[i].v+
    191                                           (t.v*denominatorDivisor.multiplicativeInverse)*bb[i].v)>>denominatorDivisor.shift;
    192                           min.v=MIN(min.v,aa[i].v);
    193                           max.v=MAX(max.v,aa[i].v);
    194                   }
    195                   if(-max<=min)min=-max;
    196                   return min;
    197           }
    198           else
    199           {
    200                   /* The next case would be to check if the result is known to fit in 32bit.
    201                    * In that case intermediate results would be handled in 64 bit,
    202                    * but the final result would not have to be checked for 32 bit overflows.
    203                    */
    204                   if(positiveResultBoundTimesD<((((int64_t)0x40000000)*denominatorDivisor.v)))
    205                   {
    206                           for(int i=0;i<c;i++)
    207                           {
    208                                   aa[i].v=((((int32_t)(((uint64_t)(((int64_t)s.v)*((int64_t)aa[i].v)+((int64_t)t.v)*((int64_t)bb[i].v)))>>denominatorDivisor.shift)))*denominatorDivisor.multiplicativeInverse);
    209                                   if(max<=aa[i])max=aa[i];
    210                                   if(aa[i]<=min)min=aa[i];
    211                           }
    212                           if(-max<=min)min=-max;
    213                   }
    214                   /* The last case would be to where we don't know that the results will fit in 32 bit.
    215                    * Since we need to compute max and min anyway, we can compute these quantities in 64bit
    216                    * and just check if they fit in 32bit at the end.
    217                    */
    218                   else
    219                   {
    220                           bool doesOverflow=(((uint32_t)t.v)==0x80000000);
    221                           int64_t min64=0;
    222                           int64_t max64=0;
    223                           for(int i=0;i<c;i++)
    224                           {
    225                                   // In the unlikely event that all the numbers to be multiplied are -2^31, the following code will overflow. Therefore the overflow flag was set as above.
    226                                   int64_t temp=(((int64_t)s.v)*((int64_t)aa[i].v)+((int64_t)t.v)*((int64_t)bb[i].v))/denominatorDivisor.v;
    227                                   if(max64<=temp)max64=temp;
    228                                   if(temp<=min64)min64=temp;
    229                                   aa[i].v=temp;
    230                           }
    231                           if(-max64<=min64)min64=-max64;
    232                           if(min64<=-0x7fffffff)doesOverflow=true;
    233                           if(doesOverflow)throw MVMachineIntegerOverflow;
    234                           min=min64;
    235                   }
    236           }
    237           return min;
    238         }
     79        // In the code products of CircuitTableInt32 objects are summed. To avoid overflows one of the factors must "fit in half" and the "number of summands" may not exceed a certain number. The bounds are specified in the following functions:
     80        bool fitsInHalf()const{return v>-30000 && v<30000;}// Is it better to allow only non-negative numbers?
     81        static bool isSuitableSummationNumber(int numberOfSummands){return numberOfSummands<30000;}
     82        CircuitTableInt32(){v=0;};
     83        CircuitTableInt32(CircuitTableInt32 const &m):v(m.v){}
     84        CircuitTableInt32(int32_t val):v(val){};
     85        class CircuitTableInt32 &operator=(int32_t a){v=a;return *this;}
     86        CircuitTableInt32 operator-()const{CircuitTableInt32 ret;ret.v=-v;return ret;}
     87        CircuitTableInt32 &operator-=(CircuitTableInt32 a){v-=a.v;return *this;}
     88        CircuitTableInt32 &operator+=(CircuitTableInt32 a){v+=a.v;return *this;}
     89        CircuitTableInt32 &operator*=(CircuitTableInt32 a){v*=a.v;return *this;}
     90        friend bool operator<=(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return a.v<=b.v;}
     91        friend bool operator==(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return a.v==b.v;}
     92        bool isZero()const{return v==0;}
     93        bool isOne()const{return v==1;}
     94        bool isNonZero()const{return v!=0;}
     95        bool isNegative()const{return v<0;}
     96        bool isPositive()const{return v>0;}
     97        friend int determinantSign1(CircuitTableInt32 const &a, CircuitTableInt32 const &c, CircuitTableInt32 const &b, CircuitTableInt32 const &d)//NOTICE MIXED ORDER. MUST WORK BY EXTENDING
     98        {
     99                int64_t r=((int64_t)a.v)*((int64_t)c.v)-((int64_t)b.v)*((int64_t)d.v);
     100                if(r>0)return 1;
     101                if(r<0)return -1;
     102                return 0;
     103        }
     104        friend int determinantSign2(CircuitTableInt32::Double const &a, CircuitTableInt32 const &c, CircuitTableInt32::Double const &b, CircuitTableInt32 const &d)//NOTICE MIXED ORDER. MUST WORK BY EXTENDING
     105        {
     106                int64_t r1=(a.v)*((int64_t)c.v);
     107                int64_t r2=(b.v)*((int64_t)d.v);
     108                if(r1>r2)return 1;
     109                if(r1<r2)return -1;
     110                return 0;
     111        }
     112        std::string toString()const{std::stringstream s;s<<v;return s.str();}
     113        Double extend()const{Double ret;ret.v=v;return ret;}
     114        CircuitTableInt32 &maddWithOverflowChecking(CircuitTableInt32 const &a, CircuitTableInt32 const&b){Double t=this->extend();t+=extendedMultiplication(a,b);*this=t.castToSingle();return *this;}
     115        CircuitTableInt32 &msubWithOverflowChecking(CircuitTableInt32 const &a, CircuitTableInt32 const&b){Double s=this->extend();s-=extendedMultiplication(a,b);*this=s.castToSingle();return *this;}
     116        CircuitTableInt32 &subWithOverflowChecking(CircuitTableInt32 const &a){Double t=this->extend();t-=a.extend();*this=t.castToSingle();return *this;}
     117        CircuitTableInt32 &addWithOverflowChecking(CircuitTableInt32 const &a){Double t=this->extend();t+=a.extend();*this=t.castToSingle();return *this;}
     118        CircuitTableInt32 negatedWithOverflowChecking()const{return (-extend()).castToSingle();}
     119        friend Double extendedMultiplication(CircuitTableInt32 const &a, CircuitTableInt32 const &b){Double ret;ret.v=((int64_t)a.v)*((int64_t)b.v);return ret;}
     120        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?
     121        friend CircuitTableInt32 min(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return (a.v>=b.v)?b:a;}
     122        friend CircuitTableInt32 max(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return (a.v<=b.v)?b:a;}
     123        friend CircuitTableInt32 negabs(CircuitTableInt32 const &a){return min(a,-a);}
     124        friend CircuitTableInt32 dotDiv(CircuitTableInt32 const &s, CircuitTableInt32 const &a, CircuitTableInt32 const &t, CircuitTableInt32 const &b, CircuitTableInt32::Divisor const &q)
     125        {
     126                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));
     127        }
     128        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.
     129        {
     130                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);
     131        }
     132        static int MIN(int32_t a, int32_t b)
     133        {
     134                return (a<b)?a:b;
     135        }
     136        static int MAX(int32_t a, int32_t b)
     137        {
     138                return (a>b)?a:b;
     139        }
     140        static CircuitTableInt32 computeNegativeBound(CircuitTableInt32 * __restrict__ Ai, int w)
     141        {
     142                CircuitTableInt32 M=0;
     143                CircuitTableInt32 m=0;
     144                for(int j=0;j<w;j++)
     145                {
     146                        m.v=MIN(m.v,Ai[j].v);
     147                        M.v=MAX(M.v,Ai[j].v);
     148                }
     149                return min(m,-M);
     150        }
     151        static CircuitTableInt32 quickNoShiftBounded(CircuitTableInt32 * __restrict__ a, CircuitTableInt32 * __restrict__ b, CircuitTableInt32 s, CircuitTableInt32 t, CircuitTableInt32::Divisor denominatorDivisor,int c)
     152        {
     153                CircuitTableInt32 *aa=a;
     154                CircuitTableInt32 *bb=b;
     155                CircuitTableInt32 max=0;
     156            CircuitTableInt32 min=0;
     157            CircuitTableInt32 ret=0;
     158            for(int i=0;i<c;i++)
     159              {
     160                    aa[i].v=((s.v*denominatorDivisor.multiplicativeInverse)*aa[i].v+
     161                                    (t.v*denominatorDivisor.multiplicativeInverse)*bb[i].v);
     162                    min.v=MIN(min.v,aa[i].v);
     163                    max.v=MAX(max.v,aa[i].v);
     164              }
     165            if(-max<=min)min=-max;
     166            ret=min;
     167            return ret;
     168        }
     169        static CircuitTableInt32 dotDivVector(CircuitTableInt32 * __restrict__ aa, CircuitTableInt32 * __restrict__ bb, CircuitTableInt32 s, CircuitTableInt32 t, CircuitTableInt32::Divisor denominatorDivisor,int c, CircuitTableInt32 boundA, CircuitTableInt32 boundB)__attribute__ ((always_inline))//assigning to first argument. The return negative of the return value is an upper bound on the absolute values of the produced entries
     170        {
     171                while(((s.v|t.v)&1)==0 && denominatorDivisor.shift>0){s.v>>=1;t.v>>=1;denominatorDivisor.shift--;denominatorDivisor.v>>=1;}
     172
     173                CircuitTableInt32 max=0;
     174                CircuitTableInt32 min=0;
     175
     176          // This is a bound on the absolute value of the sums of products before dividing by the denominator
     177          // The calculation can overflow, but since we are casting to unsigned long, this is not a problem.
     178          uint64_t positiveResultBoundTimesD=(negabs(t).v*((int64_t)boundA.v)+negabs(s).v*((int64_t)boundB.v));
     179
     180          /* The first case is the case where we know that the intermediate results fit in 32bit before shifting down.
     181           * In other words, the shifted-in bits of the results are equal.
     182           * In that case the intermediate result can be computed with 32bits.
     183           */
     184          if(positiveResultBoundTimesD<((((int64_t)0x40000000)*denominatorDivisor.v)>>denominatorDivisor.shift))//check this carefully
     185          {//FAST VERSION
     186                  // debug<<s.v<<" "<<t.v<<" "<<denominatorDivisor.v<<" "<<denominatorDivisor.shift<<"\n";
     187                  if(denominatorDivisor.shift==0)return quickNoShiftBounded(aa,bb,s,t,denominatorDivisor,c);
     188                  for(int i=0;i<c;i++)
     189                  {
     190                          aa[i].v=((s.v*denominatorDivisor.multiplicativeInverse)*aa[i].v+
     191                                          (t.v*denominatorDivisor.multiplicativeInverse)*bb[i].v)>>denominatorDivisor.shift;
     192                          min.v=MIN(min.v,aa[i].v);
     193                          max.v=MAX(max.v,aa[i].v);
     194                  }
     195                  if(-max<=min)min=-max;
     196                  return min;
     197          }
     198          else
     199          {
     200                  /* The next case would be to check if the result is known to fit in 32bit.
     201                   * In that case intermediate results would be handled in 64 bit,
     202                   * but the final result would not have to be checked for 32 bit overflows.
     203                   */
     204                  if(positiveResultBoundTimesD<((((int64_t)0x40000000)*denominatorDivisor.v)))
     205                  {
     206                          for(int i=0;i<c;i++)
     207                          {
     208                                  aa[i].v=((((int32_t)(((uint64_t)(((int64_t)s.v)*((int64_t)aa[i].v)+((int64_t)t.v)*((int64_t)bb[i].v)))>>denominatorDivisor.shift)))*denominatorDivisor.multiplicativeInverse);
     209                                  if(max<=aa[i])max=aa[i];
     210                                  if(aa[i]<=min)min=aa[i];
     211                          }
     212                          if(-max<=min)min=-max;
     213                  }
     214                  /* The last case would be to where we don't know that the results will fit in 32 bit.
     215                   * Since we need to compute max and min anyway, we can compute these quantities in 64bit
     216                   * and just check if they fit in 32bit at the end.
     217                   */
     218                  else
     219                  {
     220                          bool doesOverflow=(((uint32_t)t.v)==0x80000000);
     221                          int64_t min64=0;
     222                          int64_t max64=0;
     223                          for(int i=0;i<c;i++)
     224                          {
     225                                  // In the unlikely event that all the numbers to be multiplied are -2^31, the following code will overflow. Therefore the overflow flag was set as above.
     226                                  int64_t temp=(((int64_t)s.v)*((int64_t)aa[i].v)+((int64_t)t.v)*((int64_t)bb[i].v))/denominatorDivisor.v;
     227                                  if(max64<=temp)max64=temp;
     228                                  if(temp<=min64)min64=temp;
     229                                  aa[i].v=temp;
     230                          }
     231                          if(-max64<=min64)min64=-max64;
     232                          if(min64<=-0x7fffffff)doesOverflow=true;
     233                          if(doesOverflow)throw MVMachineIntegerOverflow;
     234                          min=min64;
     235                  }
     236          }
     237          return min;
     238        }
    239239};
    240240
    241241inline CircuitTableInt32 CircuitTableInt32::Double::castToSingle()const//casts and checks precission
    242242{
    243         if(v>=0x7fffffff || -v>=0x7fffffff) throw MVMachineIntegerOverflow;
    244         return CircuitTableInt32(v);
     243        if(v>=0x7fffffff || -v>=0x7fffffff) throw MVMachineIntegerOverflow;
     244        return CircuitTableInt32(v);
    245245}
    246246}

  • gfanlib/gfanlib_matrix.h

    r47a774 r15813d  
    5454      for(int i=0;i<getWidth();i++)
    5555          for(int j=0;j<getHeight();j++)
    56                   ret[i][j]=(*this)[j][i];
     56                  ret[i][j]=(*this)[j][i];
    5757      return ret;
    5858    }
    … …  
    6666    {
    6767      assert(m.getWidth()==width);
    68           data.resize((height+m.height)*width);
    69           int oldHeight=height;
     68          data.resize((height+m.height)*width);
     69          int oldHeight=height;
    7070      height+=m.height;
    7171      for(int i=0;i<m.height;i++)
    7272        {
    7373          for(int j=0;j<m.width;j++)
    74                   (*this)[i+oldHeight][j]=m[i][j];
     74                  (*this)[i+oldHeight][j]=m[i][j];
    7575        }
    7676    }
    7777  void appendRow(Vector<typ> const &v)
    7878    {
    79           assert(v.size()==width);
    80           data.resize((height+1)*width);
    81           height++;
    82           for(int j=0;j<width;j++)
    83                   (*this)[height-1][j]=v[j];
     79          assert(v.size()==width);
     80          data.resize((height+1)*width);
     81          height++;
     82          for(int j=0;j<width;j++)
     83                  (*this)[height-1][j]=v[j];
    8484    }
    8585  void eraseLastRow()
    … …  
    113113      for(int i=0;i<q.height;i++)
    114114          for(int j=0;j<q.width;j++)
    115                   p[i][j]=s*(q[i][j]);
     115                  p[i][j]=s*(q[i][j]);
    116116      return p;
    117117    }
    … …  
    166166  inline typ const &operator[](int j)const __attribute__((always_inline))
    167167    {
    168         assert(j>=0);
    169         assert(j<matrix.width);
    170         return matrix.data[rowNumM+j];
     168        assert(j>=0);
     169        assert(j<matrix.width);
     170        return matrix.data[rowNumM+j];
    171171    }
    172172  inline typ const &UNCHECKEDACCESS(int j)const __attribute__((always_inline))
    173173    {
    174         return matrix.data[rowNumM+j];
     174        return matrix.data[rowNumM+j];
    175175    }
    176176    const Vector<typ> toVector()const
    … …  
    178178      Vector<typ> ret(matrix.width);
    179179      for(int j=0;j<matrix.width;j++)
    180           ret[j]=matrix.data[rowNumM+j];
     180              ret[j]=matrix.data[rowNumM+j];
    181181      return ret;
    182182    }
    183183    operator Vector<typ>()const
    184                 {
    185                         return toVector();
    186                 }
     184                {
     185                        return toVector();
     186                }
    187187    bool operator==(Vector<typ> const &b)const
    188                 {
    189                         return toVector()==b;
    190                 }
     188                {
     189                        return toVector()==b;
     190                }
    191191/*    typ dot(Vector<typ> const &b)const
    192                 {
    193                         return dot(toVector(),b);
    194                 }*/
     192                {
     193                        return dot(toVector(),b);
     194                }*/
    195195    Vector<typ> operator-()const
    196196    {
    197         return -toVector();
     197            return -toVector();
    198198    }
    199199  };
    … …  
    209209    inline typ &operator[](int j) __attribute__((always_inline))
    210210      {
    211         assert(j>=0);
    212         assert(j<matrix.width);
    213         return matrix.data[rowNumM+j];
     211            assert(j>=0);
     212            assert(j<matrix.width);
     213            return matrix.data[rowNumM+j];
    214214      }
    215215    inline typ &UNCHECKEDACCESS(int j)
    216216      {
    217         return matrix.data[rowNumM+j];
     217            return matrix.data[rowNumM+j];
    218218      }
    219219    RowRef &operator=(Vector<typ> const &v)
    … …  
    221221        assert(v.size()==matrix.width);
    222222        for(int j=0;j<matrix.width;j++)
    223                 matrix.data[rowNumM+j]=v[j];
    224 
    225         return *this;
     223                matrix.data[rowNumM+j]=v[j];
     224
     225            return *this;
    226226    }
    227227    RowRef &operator=(RowRef const &v)
    … …  
    229229        assert(v.matrix.width==matrix.width);
    230230        for(int j=0;j<matrix.width;j++)
    231                 matrix.data[rowNumM+j]=v.matrix.data[v.rowNumM+j];
    232 
    233         return *this;
     231                matrix.data[rowNumM+j]=v.matrix.data[v.rowNumM+j];
     232
     233            return *this;
    234234    }
    235235/*    RowRef &operator+=(Vector<typ> const &v)
    … …  
    237237        assert(v.size()==matrix.width);
    238238        for(int j=0;j<matrix.width;j++)
    239                 matrix.data[rowNumM+j]+=v.v[j];
    240 
    241         return *this;
     239                matrix.data[rowNumM+j]+=v.v[j];
     240
     241            return *this;
    242242    }*/
    243243    RowRef &operator+=(RowRef const &v)
    … …  
    245245        assert(v.matrix.width==matrix.width);
    246246        for(int j=0;j<matrix.width;j++)
    247                 matrix.data[rowNumM+j]+=v.matrix.data[v.rowNumM+j];
    248 
    249         return *this;
     247                matrix.data[rowNumM+j]+=v.matrix.data[v.rowNumM+j];
     248
     249            return *this;
    250250    }
    251251    RowRef &operator+=(const_RowRef const &v)
    … …  
    253253        assert(v.matrix.width==matrix.width);
    254254        for(int j=0;j<matrix.width;j++)
    255                 matrix.data[rowNumM+j]+=v.matrix.data[v.rowNumM+j];
    256 
    257         return *this;
     255                matrix.data[rowNumM+j]+=v.matrix.data[v.rowNumM+j];
     256
     257            return *this;
    258258    }
    259259    RowRef &operator=(const_RowRef const &v)
    … …  
    261261        assert(v.matrix.width==matrix.width);
    262262        for(int j=0;j<matrix.width;j++)
    263                 matrix.data[rowNumM+j]=v.matrix.data[v.rowNumM+j];
    264 
    265         return *this;
     263                matrix.data[rowNumM+j]=v.matrix.data[v.rowNumM+j];
     264
     265            return *this;
    266266    }
    267267    const Vector<typ> toVector()const
    … …  
    269269      Vector<typ> ret(matrix.width);
    270270      for(int j=0;j<matrix.width;j++)
    271           ret[j]=matrix.data[rowNumM+j];
     271              ret[j]=matrix.data[rowNumM+j];
    272272      return ret;
    273273    }
    274274    operator Vector<typ>()const
    275                 {
    276                         return toVector();
    277                 }
     275                {
     276                        return toVector();
     277                }
    278278/*    typ dot(Vector<typ> const &b)const
    279                 {
    280                         return dot(toVector(),b);
    281                 }*/
     279                {
     280                        return dot(toVector(),b);
     281                }*/
    282282    bool isZero()const
    283283        {
    … …  
    303303
    304304  const typ& UNCHECKEDACCESS(int i,int j)const __attribute__((always_inline)){
    305 /*          assert(i>=0);
    306             assert(i<height);
    307             assert(j>=0);
    308             assert(j<width);*/
    309           return data[i*width+j];}
     305/*            assert(i>=0);
     306            assert(i<height);
     307            assert(j>=0);
     308            assert(j<width);*/
     309          return data[i*width+j];}
    310310  typ& UNCHECKEDACCESS(int i,int j) __attribute__((always_inline)){
    311 /*          assert(i>=0);
    312             assert(i<height);
    313             assert(j>=0);
    314             assert(j<width);*/
    315             return data[i*width+j];}
     311/*            assert(i>=0);
     312            assert(i<height);
     313            assert(j>=0);
     314            assert(j<width);*/
     315            return data[i*width+j];}
    316316
    317317
    … …  
    343343    for(int k=0;k<width;k++)
    344344      if(!(*this)[i][k].isZero())
    345           (*this)[j][k].madd((*this)[i][k],a);
     345              (*this)[j][k].madd((*this)[i][k],a);
    346346  }
    347347
    … …  
    359359  std::string toString()const
    360360  {
    361           std::stringstream f;
    362           f<<*this;
    363           return f.str();
     361          std::stringstream f;
     362          f<<*this;
     363          return f.str();
    364364  }
    365365
    … …  
    574574    while(nextPivot(pivotI,pivotJ))
    575575      {
    576         if(scalePivotsToOne)
     576            if(scalePivotsToOne)
    577577          (*this)[pivotI]=(*this)[pivotI].toVector()/(*this)[pivotI][pivotJ];
    578578        for(int i=0;i<pivotI;i++)
    … …  
    713713  void sortRows()
    714714  {
    715           std::vector<std::pair<Matrix*,int> > v;
    716           for(int i=0;i<height;i++)v.push_back(std::pair<Matrix*,int>(this,i));
    717           std::sort(v.begin(),v.end(),theRowComparer);
    718           Matrix result(height,width);
    719           for(int i=0;i<height;i++)
    720                   result[i]=(*this)[v[i].second].toVector();
    721           data=result.data;
     715          std::vector<std::pair<Matrix*,int> > v;
     716          for(int i=0;i<height;i++)v.push_back(std::pair<Matrix*,int>(this,i));
     717          std::sort(v.begin(),v.end(),theRowComparer);
     718          Matrix result(height,width);
     719          for(int i=0;i<height;i++)
     720                  result[i]=(*this)[v[i].second].toVector();
     721          data=result.data;
    722722  }
    723723  /**

  • gfanlib/gfanlib_mixedvolume.cpp

    r47a774 r15813d  
    2121static Integer stringToInteger(char const *s)
    2222{
    23         Integer ret;
    24         while(*s)ret=Integer(10)*ret+((*(s++))-'0');
    25         return ret;
     23        Integer ret;
     24        while(*s)ret=Integer(10)*ret+((*(s++))-'0');
     25        return ret;
    2626}
    2727
    … …  
    4040  virtual const char* what() const throw()
    4141  {
    42         return "Exception: The number of polytopes passed to the mixed volume computation does not match the ambient dimension.";
     42        return "Exception: The number of polytopes passed to the mixed volume computation does not match the ambient dimension.";
    4343  }
    4444}MVAmbientDimensionMismatch;
    … …  
    4747static std::vector<Matrix<mvtyp> > convertTuple(std::vector<IntMatrix> const &tuple)
    4848{
    49         vector<Matrix<mvtyp> > ret;
    50         ret.reserve(tuple.size());
    51         for(int i=0;i<tuple.size();i++)
    52         {
    53                 Matrix<mvtyp> a(tuple[i].getHeight(),tuple[i].getWidth());
    54                 for(int j=0;j<a.getHeight();j++)
    55                         for(int k=0;k<a.getWidth();k++)
    56                                 a[j][k]=tuple[i][j][k];
    57                 ret.push_back(a);
    58         }
    59         return ret;
     49        vector<Matrix<mvtyp> > ret;
     50        ret.reserve(tuple.size());
     51        for(int i=0;i<tuple.size();i++)
     52        {
     53                Matrix<mvtyp> a(tuple[i].getHeight(),tuple[i].getWidth());
     54                for(int j=0;j<a.getHeight();j++)
     55                        for(int k=0;k<a.getWidth();k++)
     56                                a[j][k]=tuple[i][j][k];
     57                ret.push_back(a);
     58        }
     59        return ret;
    6060}
    6161
    62 //      mvtyp::Double
    63         Integer mixedVolume(vector<IntMatrix> const &tuple, int nthreads, int steps)
    64         {
    65                 //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.
    66                 //Alternative explanation. A generic polynomial system with zero variables has exactly 1 solution.
    67                 if(tuple.size()==0){return Integer(1);}
    68 
    69                 // For now, let's make in an error to pass the wrong number of polytopes. Of course we could define the
    70                 // mixed volume to be 0 in that case, but the code calling will almost certainly be wrong. It is safer to throw.
    71                 // for(const auto &x:tuple) // Range based for is not supported until gcc 4.6/
    72                 for(auto x=tuple.begin();x!=tuple.end();x++)
    73                         if(x->getHeight()!=tuple.size())
    74                                 {throw MVAmbientDimensionMismatch;}
    75 
    76                 std::vector<Matrix<mvtyp> > tuple2=convertTuple(tuple);
    77                 //              assert(tuple2.size());
    78                 typedef SpecializedRTraverser<mvtyp,mvtyp::Double,mvtyp::Divisor> MySpecializedRTraverser;
    79                 vector<MySpecializedRTraverser> T1;
    80                 int N=nthreads;
    81                 if(N==0)N=1;// Even if we do not parallelize, we still need one traverser.
    82                 T1.reserve(N);
    83                 vector<Traverser*> I;
    84                 for(int i=0;i<N;i++)T1.emplace_back(tuple2);
    85                 for(int i=0;i<N;i++)I.push_back(&(T1[i]));
    86 
    87 //              debug<<"Initialized.\n";
    88 
    89                 if(nthreads)
    90                         traverse_threaded(&(I[0]),N,steps);
    91                 else
    92                         traverse_simple(I[0]);
    93 
    94                 bool aborted=false;
    95                 for(int i=0;i<N;i++)aborted|=I[i]->aborting;
    96 
    97                 if(aborted)
    98                 {
    99                         throw MVExceptionRethrow;
    100                 }
    101 
    102                 mvtyp::Double total;
    103                 int totalSteps=0;
    104                 for(int i=0;i<N;i++)
    105                 {
    106 //                      debug<<"#"<<((MySpecializedRTraverser*)(I[i]))->mixedVolume.toString()<<"Steps:"<<((MySpecializedRTraverser*)(I[i]))->numberOfExpensiveSteps<<"\n";
    107                         //debug<<((SpecializedMTraverser*)(I[i]))->T.counter;
    108                         total.addWithOverflowCheck(((MySpecializedRTraverser*)(I[i]))->mixedVolume);
    109                         totalSteps+=((MySpecializedRTraverser*)(I[i]))->numberOfExpensiveSteps;
    110                 }
    111 //              debug<<"Total:"<<(int)total.v<<"\n";
    112 //              debug<<"Totalsteps:"<<totalSteps<<"\n";
    113 
    114 
    115                 return stringToInteger(total.toString().c_str());
    116         }
    117 
    118 
    119         namespace MixedVolumeExamples{
    120                 vector<IntMatrix> cyclic(int n)
    121                 {
    122                         vector<IntMatrix> ret;
    123                         for(int i=1;i<n;i++)
    124                         {
    125                                 IntMatrix m(n,n);
    126                                 for(int y=0;y<n;y++)
    127                                         for(int x=0;x<n;x++)
    128                                                 m[y][x]=((x-y+n)%n)<i;
    129                                 ret.push_back(m);
    130                         }
    131 
    132                         IntMatrix last(n,2);
    133                         for(int y=0;y<n;y++)last[y][0]=1;
    134                         ret.push_back(last);
    135 
    136                         return ret;
    137                 }
    138                 vector<IntMatrix> noon(int n)
    139                 {
    140                         vector<IntMatrix> ret;
    141                         for(int i=0;i<n;i++)
    142                         {
    143                                 IntMatrix m(n,n+1);
    144                                 for(int y=0;y<n-1;y++)
    145                                         m[y+(i<=y)][y]=2;
    146                                 for(int x=0;x<n;x++)
    147                                         m[i][x]=1;
    148                                 ret.push_back(m);
    149                         }
    150                         return ret;
    151                 }
    152                 vector<IntMatrix> chandra(int n)
    153                 {
    154                         vector<IntMatrix> ret;
    155                         for(int i=0;i<n;i++)
    156                         {
    157                                 IntMatrix m(n,n+1);
    158                                 for(int y=0;y<n-1;y++)
    159                                         m[y][y+1]=1;
    160                                 for(int x=0;x<n;x++)
    161                                         m[i][x]+=1;
    162                                 ret.push_back(m);
    163                         }
    164                         return ret;
    165                 }
    166                 vector<IntMatrix> katsura(int n)
    167                 {
    168                         n++;
    169                         vector<IntMatrix> ret;
    170                         for(int i=0;i<n-1;i++)
    171                         {
    172                                 IntMatrix m(n,n+1-((i+1)/2));
    173                                 for(int y=0;y<n-((i+1)/2);y++)
    174                                 {
    175                                         m[n-1-y][y]=1;
    176                                         m[(n-1-y-i)>0 ? (n-1-y-i) : -(n-1-y-i)][y]+=1;
    177                                 }
    178                                 m[i][m.getWidth()-1]=1;
    179                                 ret.push_back(m);
    180                         }
    181                         ret.push_back(combineLeftRight(IntMatrix::identity(n),IntMatrix(n,1)));
    182                         return ret;
    183                 }
    184                 vector<IntMatrix> gaukwa(int n)
    185                 {
    186                         vector<IntMatrix> ret;
    187                         for(int i=0;i<2*n;i++)
    188                                 ret.push_back(combineLeftRight(combineOnTop(IntMatrix::identity(n),i*IntMatrix::identity(n)),IntMatrix(2*n,1)));
    189                         return ret;
    190                 }
    191                 vector<IntMatrix> eco(int n)
    192                 {
    193                         vector<IntMatrix> ret;
    194                         for(int i=0;i<n-1;i++)
    195                         {
    196                                 IntMatrix m(n,n-i);
    197                                 for(int y=0;y<n-(i+1);y++)
    198                                 {
    199                                         m[y+i][y]=1;
    200                                         m[n-1][y]=1;
    201                                         if(y)m[y-1][y]=1;
    202                                 }
    203                                 ret.push_back(m);
    204                         }
    205                         ret.push_back(combineLeftRight(combineOnTop(IntMatrix::identity(n-1),IntMatrix(1,n-1)),IntMatrix(n,1)));
    206                         return ret;
    207                 }
    208         }
     62//        mvtyp::Double
     63        Integer mixedVolume(vector<IntMatrix> const &tuple, int nthreads, int steps)
     64        {
     65                //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.
     66                //Alternative explanation. A generic polynomial system with zero variables has exactly 1 solution.
     67                if(tuple.size()==0){return Integer(1);}
     68
     69                // For now, let's make in an error to pass the wrong number of polytopes. Of course we could define the
     70                // mixed volume to be 0 in that case, but the code calling will almost certainly be wrong. It is safer to throw.
     71                // for(const auto &x:tuple) // Range based for is not supported until gcc 4.6/
     72                for(auto x=tuple.begin();x!=tuple.end();x++)
     73                        if(x->getHeight()!=tuple.size())
     74                                {throw MVAmbientDimensionMismatch;}
     75
     76                std::vector<Matrix<mvtyp> > tuple2=convertTuple(tuple);
     77                //                assert(tuple2.size());
     78                typedef SpecializedRTraverser<mvtyp,mvtyp::Double,mvtyp::Divisor> MySpecializedRTraverser;
     79                vector<MySpecializedRTraverser> T1;
     80                int N=nthreads;
     81                if(N==0)N=1;// Even if we do not parallelize, we still need one traverser.
     82                T1.reserve(N);
     83                vector<Traverser*> I;
     84                for(int i=0;i<N;i++)T1.emplace_back(tuple2);
     85                for(int i=0;i<N;i++)I.push_back(&(T1[i]));
     86
     87//                debug<<"Initialized.\n";
     88
     89                if(nthreads)
     90                        traverse_threaded(&(I[0]),N,steps);
     91                else
     92                        traverse_simple(I[0]);
     93
     94                bool aborted=false;
     95                for(int i=0;i<N;i++)aborted|=I[i]->aborting;
     96
     97                if(aborted)
     98                {
     99                        throw MVExceptionRethrow;
     100                }
     101
     102                mvtyp::Double total;
     103                int totalSteps=0;
     104                for(int i=0;i<N;i++)
     105                {
     106//                        debug<<"#"<<((MySpecializedRTraverser*)(I[i]))->mixedVolume.toString()<<"Steps:"<<((MySpecializedRTraverser*)(I[i]))->numberOfExpensiveSteps<<"\n";
     107                        //debug<<((SpecializedMTraverser*)(I[i]))->T.counter;
     108                        total.addWithOverflowCheck(((MySpecializedRTraverser*)(I[i]))->mixedVolume);
     109                        totalSteps+=((MySpecializedRTraverser*)(I[i]))->numberOfExpensiveSteps;
     110                }
     111//                debug<<"Total:"<<(int)total.v<<"\n";
     112//                debug<<"Totalsteps:"<<totalSteps<<"\n";
     113
     114
     115                return stringToInteger(total.toString().c_str());
     116        }
     117
     118
     119        namespace MixedVolumeExamples{
     120                vector<IntMatrix> cyclic(int n)
     121                {
     122                        vector<IntMatrix> ret;
     123                        for(int i=1;i<n;i++)
     124                        {
     125                                IntMatrix m(n,n);
     126                                for(int y=0;y<n;y++)
     127                                        for(int x=0;x<n;x++)
     128                                                m[y][x]=((x-y+n)%n)<i;
     129                                ret.push_back(m);
     130                        }
     131
     132                        IntMatrix last(n,2);
     133                        for(int y=0;y<n;y++)last[y][0]=1;
     134                        ret.push_back(last);
     135
     136                        return ret;
     137                }
     138                vector<IntMatrix> noon(int n)
     139                {
     140                        vector<IntMatrix> ret;
     141                        for(int i=0;i<n;i++)
     142                        {
     143                                IntMatrix m(n,n+1);
     144                                for(int y=0;y<n-1;y++)
     145                                        m[y+(i<=y)][y]=2;
     146                                for(int x=0;x<n;x++)
     147                                        m[i][x]=1;
     148                                ret.push_back(m);
     149                        }
     150                        return ret;
     151                }
     152                vector<IntMatrix> chandra(int n)
     153                {
     154                        vector<IntMatrix> ret;
     155                        for(int i=0;i<n;i++)
     156                        {
     157                                IntMatrix m(n,n+1);
     158                                for(int y=0;y<n-1;y++)
     159                                        m[y][y+1]=1;
     160                                for(int x=0;x<n;x++)
     161                                        m[i][x]+=1;
     162                                ret.push_back(m);
     163                        }
     164                        return ret;
     165                }
     166                vector<IntMatrix> katsura(int n)
     167                {
     168                        n++;
     169                        vector<IntMatrix> ret;
     170                        for(int i=0;i<n-1;i++)
     171                        {
     172                                IntMatrix m(n,n+1-((i+1)/2));
     173                                for(int y=0;y<n-((i+1)/2);y++)
     174                                {
     175                                        m[n-1-y][y]=1;
     176                                        m[(n-1-y-i)>0 ? (n-1-y-i) : -(n-1-y-i)][y]+=1;
     177                                }
     178                                m[i][m.getWidth()-1]=1;
     179                                ret.push_back(m);
     180                        }
     181                        ret.push_back(combineLeftRight(IntMatrix::identity(n),IntMatrix(n,1)));
     182                        return ret;
     183                }
     184                vector<IntMatrix> gaukwa(int n)
     185                {
     186                        vector<IntMatrix> ret;
     187                        for(int i=0;i<2*n;i++)
     188                                ret.push_back(combineLeftRight(combineOnTop(IntMatrix::identity(n),i*IntMatrix::identity(n)),IntMatrix(2*n,1)));
     189                        return ret;
     190                }
     191                vector<IntMatrix> eco(int n)
     192                {
     193                        vector<IntMatrix> ret;
     194                        for(int i=0;i<n-1;i++)
     195                        {
     196                                IntMatrix m(n,n-i);
     197                                for(int y=0;y<n-(i+1);y++)
     198                                {
     199                                        m[y+i][y]=1;
     200                                        m[n-1][y]=1;
     201                                        if(y)m[y-1][y]=1;
     202                                }
     203                                ret.push_back(m);
     204                        }
     205                        ret.push_back(combineLeftRight(combineOnTop(IntMatrix::identity(n-1),IntMatrix(1,n-1)),IntMatrix(n,1)));
     206                        return ret;
     207                }
     208        }
    209209
    210210}

  • gfanlib/gfanlib_mixedvolume.h

    r47a774 r15813d  
    1515
    1616namespace gfan{
    17         Integer mixedVolume(std::vector<IntMatrix> const &tuple, int nthreads=0, int steps=500);
    18         namespace MixedVolumeExamples{
    19                 std::vector<IntMatrix> cyclic(int n);
    20                 std::vector<IntMatrix> noon(int n);
    21                 std::vector<IntMatrix> chandra(int n);
    22                 std::vector<IntMatrix> katsura(int n);
    23                 std::vector<IntMatrix> gaukwa(int n);
    24                 std::vector<IntMatrix> eco(int n);
    25         }
     17        Integer mixedVolume(std::vector<IntMatrix> const &tuple, int nthreads=0, int steps=500);
     18        namespace MixedVolumeExamples{
     19                std::vector<IntMatrix> cyclic(int n);
     20                std::vector<IntMatrix> noon(int n);
     21                std::vector<IntMatrix> chandra(int n);
     22                std::vector<IntMatrix> katsura(int n);
     23                std::vector<IntMatrix> gaukwa(int n);
     24                std::vector<IntMatrix> eco(int n);
     25        }
    2626
    2727}

  • gfanlib/gfanlib_paralleltraverser.cpp

    r47a774 r15813d  
    159159
    160160public:
    161         bool aborting;                                                                                                          // Added by Anders
     161        bool aborting;                                                                                                                // Added by Anders
    162162  // Create a new Job. if first_split is not set (or it is -2), the
    163163  // first split will be found.
    164164  Job( vector<TraverseState>*  stack = new vector<TraverseState>(),
    165165       int                     first_split = -2 )
    166   :aborting(false)                                                                                                              // Added by Anders
     166  :aborting(false)                                                                                                                // Added by Anders
    167167  {
    168168    if (first_split == -2) {
    … …  
    255255
    256256        int  prev_index  =  traverser->moveToNext(state.next_index, true);
    257         aborting=traverser->aborting;                                           /* Added by Anders */
     257        aborting=traverser->aborting;                                                /* Added by Anders */
    258258        steps++;
    259259
    … …  
    409409
    410410public:
    411         bool aborting;                                                                  // Added by Anders
     411        bool aborting;                                                                        // Added by Anders
    412412  // step_count is the number of algorithm steps taken between
    413413  // possible job transfers. This value should be high (e.g. 100) if
    … …  
    417417              int          count,
    418418              int          step_count )
    419   :aborting(false)                                                                      // Added by Anders
     419  :aborting(false)                                                                        // Added by Anders
    420420  {
    421421    context_count  =  count;
    … …  
    463463  bool  hasTransfer( void )
    464464  {
    465         return !transfers->empty();
     465        return !transfers->empty();
    466466  }
    467467
    … …  
    517517
    518518      for(deque<JobTransfer*>::const_iterator tr=transfers->begin();tr!=transfers->end();tr++){
    519         //      for (JobTransfer*  tr : *transfers) {
     519        //      for (JobTransfer*  tr : *transfers) {
    520520        (*tr)->setJob(NULL);
    521521      }
    … …  
    542542    int  step_count  =  central->hasTransfer() ? 1 : context->step_count;
    543543
    544     if(central->aborting)job->aborting=true;                    // Added by Anders
     544    if(central->aborting)job->aborting=true;                        // Added by Anders
    545545
    546546    while (job->step(step_count)) {
    547       if(job->aborting)central->aborting=true;                  // Added by Anders
     547      if(job->aborting)central->aborting=true;                        // Added by Anders
    548548      JobTransfer*  transfer  =  central->getTransfer();
    549549

  • gfanlib/gfanlib_paralleltraverser.h

    r47a774 r15813d  
    2727{
    2828public:
    29         bool aborting;                                  // Added by Anders
    30         void abort(){aborting=true;}    // Added by Anders
    31         Traverser():aborting(false){}   // Added by Anders
     29        bool aborting;                      // Added by Anders
     30        void abort(){aborting=true;}        // Added by Anders
     31        Traverser():aborting(false){}       // Added by Anders
    3232  // Virtual destructor
    3333  virtual ~Traverser( void ) {};

  • gfanlib/gfanlib_symmetriccomplex.cpp

    r47a774 r15813d  
    185185int SymmetricComplex::getLinDim()const
    186186{
    187         return linealitySpace.getHeight();
     187        return linealitySpace.getHeight();
    188188}
    189189

  • gfanlib/gfanlib_traversal.cpp

    r47a774 r15813d  
    1313
    1414static list<ZVector> rowsToList(ZMatrix const &m)
    15                 {
    16                         list<ZVector> ret;
    17                         for(int i=0;i<m.getHeight();i++)
    18                                 ret.push_back(m[i]);
    19                         return ret;
    20                 }
     15                {
     16                        list<ZVector> ret;
     17                        for(int i=0;i<m.getHeight();i++)
     18                                ret.push_back(m[i]);
     19                        return ret;
     20                }
    2121
    2222bool FanTraverser::hasNoState()const
    … …  
    3131
    3232FanBuilder::FanBuilder(int n, SymmetryGroup const &sym):
    33         coneCollection(n)
     33        coneCollection(n)
    3434{
    3535}
    … …  
    3838bool FanBuilder::process(FanTraverser &traverser)
    3939{
    40         ZCone cone2=traverser.refToPolyhedralCone();
    41         cone2.canonicalize();
    42         coneCollection.insert(cone2);
    43         return true;
     40        ZCone cone2=traverser.refToPolyhedralCone();
     41        cone2.canonicalize();
     42        coneCollection.insert(cone2);
     43        return true;
    4444}
    4545
    … …  
    128128  void removeDuplicates(ZVector const &ridge, list<ZVector> &rays, list<ZVector> *normals=0)const
    129129  {
    130           list<ZVector> ret;
    131           list<ZVector> normalsRet;
    132           set<ZVector> representatives;
    133           list<ZVector>::const_iterator I;
     130          list<ZVector> ret;
     131          list<ZVector> normalsRet;
     132          set<ZVector> representatives;
     133          list<ZVector>::const_iterator I;
    134134    if(normals)I=normals->begin();
    135135    for(list<ZVector>::const_iterator i=rays.begin();i!=rays.end();i++)
    136136      {
    137         ZVector rep=sym.orbitRepresentativeFixing(*i,ridge);
    138         if(representatives.count(rep)==0)
    139           {
    140             representatives.insert(rep);
    141             ret.push_back(*i);
    142             if(normals)normalsRet.push_back(*I);
    143           }
    144         if(normals)I++;
     137        ZVector rep=sym.orbitRepresentativeFixing(*i,ridge);
     138        if(representatives.count(rep)==0)
     139          {
     140            representatives.insert(rep);
     141            ret.push_back(*i);
     142            if(normals)normalsRet.push_back(*I);
     143          }
     144        if(normals)I++;
    145145      }
    146146    rays=ret;
    … …  
    152152    for(map<EFirst,ESecond>::const_iterator i=theSet.begin();i!=theSet.end();i++)
    153153      {
    154         cerr << i->first.first << i->first.second;
    155         cerr << endl;
     154            cerr << i->first.first << i->first.second;
     155            cerr << endl;
    156156      }
    157157    cerr<<endl<<endl;
    … …  
    182182struct pathStepFacet
    183183{
    184         list<ZVector> ridges;
    185         list<ZVector> ridgesRayUniqueVector;//stores the ray of the link that we came from
     184        list<ZVector> ridges;
     185        list<ZVector> ridgesRayUniqueVector;//stores the ray of the link that we came from
    186186};
    187187
    … …  
    263263  //Alternative the links should be computed and stored the first time a facet is seen.
    264264  //Or the conetraverse should be given more info about the ridge to make computations quicker.
    265         int lastNumberOfEdges=0;
    266         float averageEdge=0;
    267         int n=traverser.refToPolyhedralCone().ambientDimension();//symmetryGroup->sizeOfBaseSet();
    268         SymmetryGroup localSymmetryGroup(n);
    269         if(!symmetryGroup)symmetryGroup=&localSymmetryGroup;
    270 
    271         ZMatrix linealitySpaceGenerators=traverser.refToPolyhedralCone().generatorsOfLinealitySpace();
    272 
    273         int d=traverser.refToPolyhedralCone().dimension();
    274 
    275           Boundary boundary(*symmetryGroup);
    276           list<pathStepFacet> facetStack;
    277           list<pathStepRidge> ridgeStack;
    278 
    279           int numberOfCompletedFacets=0;
    280           int numberOfCompletedRidges=0;
    281           int stackSize=0;
    282 
    283           ZVector facetUniqueVector;
    284           goto entry;
    285           while(1)
    286             {
    287             L1:
    288 //          printStack(facetStack,ridgeStack);
    289             //if we have more ProcessRidge calls to do
    290               if(!facetStack.front().ridges.empty())
    291                 {
    292                   //ProcessRidge "called"
    293                   pathStepRidge top;
     265        int lastNumberOfEdges=0;
     266        float averageEdge=0;
     267        int n=traverser.refToPolyhedralCone().ambientDimension();//symmetryGroup->sizeOfBaseSet();
     268        SymmetryGroup localSymmetryGroup(n);
     269        if(!symmetryGroup)symmetryGroup=&localSymmetryGroup;
     270
     271        ZMatrix linealitySpaceGenerators=traverser.refToPolyhedralCone().generatorsOfLinealitySpace();
     272
     273        int d=traverser.refToPolyhedralCone().dimension();
     274
     275          Boundary boundary(*symmetryGroup);
     276          list<pathStepFacet> facetStack;
     277          list<pathStepRidge> ridgeStack;
     278
     279          int numberOfCompletedFacets=0;
     280          int numberOfCompletedRidges=0;
     281          int stackSize=0;
     282
     283          ZVector facetUniqueVector;
     284          goto entry;
     285          while(1)
     286            {
     287            L1:
     288//            printStack(facetStack,ridgeStack);
     289            //if we have more ProcessRidge calls to do
     290              if(!facetStack.front().ridges.empty())
     291                {
     292                      //ProcessRidge "called"
     293                  pathStepRidge top;
    294294
    295295
    … …  
    303303                    }
    304304
    305                   top.parentRidge=facetStack.front().ridges.front();
    306 //                AsciiPrinter(Stderr)<<top.parentRay<<"--------------------------------++++\n";
    307                   list<ZVector> rays;
     305                  top.parentRidge=facetStack.front().ridges.front();
     306//                  AsciiPrinter(Stderr)<<top.parentRay<<"--------------------------------++++\n";
     307                  list<ZVector> rays;
    308308                  if(traverser.hasNoState())
    309309                    {
    … …  
    314314                    rays=traverser.link(facetStack.front().ridges.front());
    315315
    316                   assert(!rays.empty());
    317                   boundary.removeDuplicates(top.parentRidge,rays);
    318                   ridgeStack.push_front(top);stackSize++;
    319                   ZVector ridgeRidgeRidge=facetStack.front().ridges.front();
    320                   for(list<ZVector>::const_iterator i=rays.begin();i!=rays.end();i++)
    321                     {
    322                       ridgeStack.front().rays.push_front(*i);
    323                       if(boundary.containsFlip(ridgeRidgeRidge,*i,&ridgeStack.front().rays,ridgeStack.front().rays.begin(),0,ridgeStack.front().rays.begin()))
    324                         ridgeStack.front().rays.pop_front();
    325                     }
    326                   // "state saved" ready to do calls to ProcessFacet
    327                   numberOfCompletedRidges++;
    328                 }
    329               else
    330                 {
    331                   // No more calls to do - we now return from ProcessFacet
    332                   //THIS IS THE PLACE TO CHANGE CONE BACK
    333                   facetStack.pop_front();stackSize--;
    334                   if(facetStack.empty())break;
    335 //      log1 cerr<<"BACK"<<endl;
    336           if(!traverser.hasNoState())traverser.changeCone(ridgeStack.front().parentRidge,ridgeStack.front().parentRay);
    337                 }
    338             L2:
    339 //          printStack(facetStack,ridgeStack);
    340             //check if ProcessRidge needs to make more ProcessFacet calls
    341               if(!ridgeStack.front().rays.empty())
    342                 {
    343 //                log1 cerr<<"FORWARD"<<endl;
    344                   traverser.changeCone(ridgeStack.front().parentRidge,ridgeStack.front().rays.front());
    345                 entry:
    346                 //ProcessFacet()
    347                 averageEdge=0.99*averageEdge+0.01*(boundary.size()-lastNumberOfEdges);
    348 //              log1 fprintf(Stderr,"\n-------------------------------------\n");
    349 //                log1 fprintf(Stderr,"Boundary edges in bipartite graph: %i, Completed ridges: %i, Completed facets: %i, Recursion depth:%i Average new edge/facet:%0.2f\n",boundary.size(),numberOfCompletedRidges,numberOfCompletedFacets,stackSize,averageEdge);
    350 //                log1 fprintf(Stderr,"-------------------------------------\n");
    351                   lastNumberOfEdges=boundary.size();
    352 
    353 //                target.process(traverser);//Postponed until extrem rays have been computed
    354                   ZMatrix extremeRays=traverser.refToPolyhedralCone().extremeRays(&linealitySpaceGenerators);
     316                  assert(!rays.empty());
     317                  boundary.removeDuplicates(top.parentRidge,rays);
     318                  ridgeStack.push_front(top);stackSize++;
     319                  ZVector ridgeRidgeRidge=facetStack.front().ridges.front();
     320                  for(list<ZVector>::const_iterator i=rays.begin();i!=rays.end();i++)
     321                    {
     322                      ridgeStack.front().rays.push_front(*i);
     323                      if(boundary.containsFlip(ridgeRidgeRidge,*i,&ridgeStack.front().rays,ridgeStack.front().rays.begin(),0,ridgeStack.front().rays.begin()))
     324                        ridgeStack.front().rays.pop_front();
     325                    }
     326                  // "state saved" ready to do calls to ProcessFacet
     327                  numberOfCompletedRidges++;
     328                }
     329              else
     330                {
     331                      // No more calls to do - we now return from ProcessFacet
     332                      //THIS IS THE PLACE TO CHANGE CONE BACK
     333                  facetStack.pop_front();stackSize--;
     334                  if(facetStack.empty())break;
     335//            log1 cerr<<"BACK"<<endl;
     336          if(!traverser.hasNoState())traverser.changeCone(ridgeStack.front().parentRidge,ridgeStack.front().parentRay);
     337                }
     338            L2:
     339//            printStack(facetStack,ridgeStack);
     340            //check if ProcessRidge needs to make more ProcessFacet calls
     341              if(!ridgeStack.front().rays.empty())
     342                {
     343//                      log1 cerr<<"FORWARD"<<endl;
     344                      traverser.changeCone(ridgeStack.front().parentRidge,ridgeStack.front().rays.front());
     345                entry:
     346                //ProcessFacet()
     347                averageEdge=0.99*averageEdge+0.01*(boundary.size()-lastNumberOfEdges);
     348//                log1 fprintf(Stderr,"\n-------------------------------------\n");
     349//                  log1 fprintf(Stderr,"Boundary edges in bipartite graph: %i, Completed ridges: %i, Completed facets: %i, Recursion depth:%i Average new edge/facet:%0.2f\n",boundary.size(),numberOfCompletedRidges,numberOfCompletedFacets,stackSize,averageEdge);
     350//                  log1 fprintf(Stderr,"-------------------------------------\n");
     351                  lastNumberOfEdges=boundary.size();
     352
     353//                  target.process(traverser);//Postponed until extrem rays have been computed
     354                  ZMatrix extremeRays=traverser.refToPolyhedralCone().extremeRays(&linealitySpaceGenerators);
    355355                  target.process(traverser);
    356356
    357 //                IntegerVectorList inequalities=traverser.refToPolyhedralCone().getHalfSpaces();
     357//                  IntegerVectorList inequalities=traverser.refToPolyhedralCone().getHalfSpaces();
    358358                  ZMatrix equations=traverser.refToPolyhedralCone().getEquations();
    359 //                facetUniqueVector=traverser.refToPolyhedralCone().getUniquePoint();
    360                   facetUniqueVector=traverser.refToPolyhedralCone().getUniquePointFromExtremeRays(extremeRays);
    361                   list<ZVector> facetNormals=rowsToList(traverser.refToPolyhedralCone().getFacets());
    362 
    363                   pathStepFacet stepFacet;
    364                   list<ZVector> ridges;
    365 
    366                   for(list<ZVector>::iterator i=facetNormals.begin();i!=facetNormals.end();i++)
    367                     {
    368                           ZVector v(n);
    369 //                        for(IntegerVectorList::const_iterator j=extremeRays.begin();j!=extremeRays.end();j++)
    370                           for(int j=0;j<extremeRays.getHeight();j++)
    371                                   if(dot(*i,extremeRays[j]).isZero())v+=extremeRays[j];
    372                           ridges.push_back(v);
    373                     }
    374 
    375                   ZVector temp(n);
    376 //                boundary.removeDuplicates(temp,ridges);//use facetUniqueVector instead
    377                   boundary.removeDuplicates(facetUniqueVector,ridges,&facetNormals);//use facetUniqueVector instead
    378 
    379                   facetStack.push_front(stepFacet);stackSize++;
    380                   list<ZVector>::const_iterator I=facetNormals.begin();
    381                   for(list<ZVector>::const_iterator i=ridges.begin();i!=ridges.end();i++,I++)
    382                     {
    383                           ZVector rayUniqueVector;
    384 
    385                           if(d==n)
    386                           {
    387                                 rayUniqueVector =I->normalized();
    388 //                              if(dotLong(rayUniqueVector,*I)
    389                           }
    390                           else
    391                           {
    392                                   ZCone rayCone=traverser.refToPolyhedralCone().link(*i);
    393                                   rayCone.canonicalize();
    394                                   rayUniqueVector=rayCone.getUniquePoint();
    395 //                                debug<<traverser.refToPolyhedralCone();
    396 //                                debug<<rayCone;
    397                           }
    398                       facetStack.front().ridges.push_front(*i);
    399                       if(traverser.hasNoState())facetStack.front().ridgesRayUniqueVector.push_front(rayUniqueVector);
    400 
    401 
    402                       if(!traverser.hasNoState())
    403                         {
    404                         if(boundary.containsFlip(*i,rayUniqueVector,&facetStack.front().ridges,facetStack.front().ridges.begin(),0,facetStack.front().ridges.begin()))
     359//                  facetUniqueVector=traverser.refToPolyhedralCone().getUniquePoint();
     360                  facetUniqueVector=traverser.refToPolyhedralCone().getUniquePointFromExtremeRays(extremeRays);
     361                  list<ZVector> facetNormals=rowsToList(traverser.refToPolyhedralCone().getFacets());
     362
     363                  pathStepFacet stepFacet;
     364                  list<ZVector> ridges;
     365
     366                  for(list<ZVector>::iterator i=facetNormals.begin();i!=facetNormals.end();i++)
     367                    {
     368                          ZVector v(n);
     369//                          for(IntegerVectorList::const_iterator j=extremeRays.begin();j!=extremeRays.end();j++)
     370                          for(int j=0;j<extremeRays.getHeight();j++)
     371                                  if(dot(*i,extremeRays[j]).isZero())v+=extremeRays[j];
     372                          ridges.push_back(v);
     373                    }
     374
     375                  ZVector temp(n);
     376//                  boundary.removeDuplicates(temp,ridges);//use facetUniqueVector instead
     377                  boundary.removeDuplicates(facetUniqueVector,ridges,&facetNormals);//use facetUniqueVector instead
     378
     379                  facetStack.push_front(stepFacet);stackSize++;
     380                  list<ZVector>::const_iterator I=facetNormals.begin();
     381                  for(list<ZVector>::const_iterator i=ridges.begin();i!=ridges.end();i++,I++)
     382                    {
     383                          ZVector rayUniqueVector;
     384
     385                          if(d==n)
     386                          {
     387                                rayUniqueVector =I->normalized();
     388//                                if(dotLong(rayUniqueVector,*I)
     389                          }
     390                          else
     391                          {
     392                                  ZCone rayCone=traverser.refToPolyhedralCone().link(*i);
     393                                  rayCone.canonicalize();
     394                                  rayUniqueVector=rayCone.getUniquePoint();
     395//                                  debug<<traverser.refToPolyhedralCone();
     396//                                  debug<<rayCone;
     397                          }
     398                      facetStack.front().ridges.push_front(*i);
     399                      if(traverser.hasNoState())facetStack.front().ridgesRayUniqueVector.push_front(rayUniqueVector);
     400
     401
     402                      if(!traverser.hasNoState())
     403                        {
     404                        if(boundary.containsFlip(*i,rayUniqueVector,&facetStack.front().ridges,facetStack.front().ridges.begin(),0,facetStack.front().ridges.begin()))
    405405                        {
    406406                          facetStack.front().ridges.pop_front();
    407407                        }
    408                         }
    409                       else
    410                         {
    411                       if(boundary.containsFlip(*i,rayUniqueVector,&facetStack.front().ridges,facetStack.front().ridges.begin(),&facetStack.front().ridgesRayUniqueVector,facetStack.front().ridgesRayUniqueVector.begin()))
    412                         {
    413                           facetStack.front().ridges.pop_front();
    414                           facetStack.front().ridgesRayUniqueVector.pop_front();
    415                         }
    416                         }
    417                     }
    418                   //"State pushed" ready to call ProcessRidge
    419 
    420                   numberOfCompletedFacets++;
    421                 }
    422               else
    423                 {
    424                   //ProcessRidge is done making its calls to ProcessFacet so we can return from ProcessRidge
    425 //                cerr<<"BACK"<<endl;
    426 //                traverser.changeCone(ridgeStack.front().parentRidge,ridgeStack.front().parentRay);
    427                   ridgeStack.pop_front();stackSize--;
    428                 }
    429             }//goto L1
    430           //      log1 fprintf(Stderr,"\n-------------------------------------\n");
    431 //        log1 fprintf(Stderr,"Boundary edges in bipartite graph: %i, Completed ridges: %i, Completed facets: %i, Recursion depth:%i\n",boundary.size(),numberOfCompletedRidges,numberOfCompletedFacets,stackSize);
    432 //        log1 fprintf(Stderr,"-------------------------------------\n");
     408                        }
     409                      else
     410                        {
     411                      if(boundary.containsFlip(*i,rayUniqueVector,&facetStack.front().ridges,facetStack.front().ridges.begin(),&facetStack.front().ridgesRayUniqueVector,facetStack.front().ridgesRayUniqueVector.begin()))
     412                        {
     413                          facetStack.front().ridges.pop_front();
     414                          facetStack.front().ridgesRayUniqueVector.pop_front();
     415                        }
     416                        }
     417                    }
     418                  //"State pushed" ready to call ProcessRidge
     419
     420                  numberOfCompletedFacets++;
     421                }
     422              else
     423                {
     424                      //ProcessRidge is done making its calls to ProcessFacet so we can return from ProcessRidge
     425//                      cerr<<"BACK"<<endl;
     426//                  traverser.changeCone(ridgeStack.front().parentRidge,ridgeStack.front().parentRay);
     427                  ridgeStack.pop_front();stackSize--;
     428                }
     429            }//goto L1
     430          //          log1 fprintf(Stderr,"\n-------------------------------------\n");
     431//          log1 fprintf(Stderr,"Boundary edges in bipartite graph: %i, Completed ridges: %i, Completed facets: %i, Recursion depth:%i\n",boundary.size(),numberOfCompletedRidges,numberOfCompletedFacets,stackSize);
     432//          log1 fprintf(Stderr,"-------------------------------------\n");
    433433}

  • gfanlib/gfanlib_traversal.h

    r47a774 r15813d  
    1111{
    1212public:
    13         /**
    14         * Go to the cone which is connected to the current facet through the ridge in direction ray.
    15         * The "ridge" is a relative interior point of the ridge.
    16         */
    17         virtual void changeCone(ZVector const &ridgeVector, ZVector const &rayVector)=0;
     13        /**
     14        * Go to the cone which is connected to the current facet through the ridge in direction ray.
     15        * The "ridge" is a relative interior point of the ridge.
     16        */
     17        virtual void changeCone(ZVector const &ridgeVector, ZVector const &rayVector)=0;
    1818/**
    1919 * Compute the link of the fan in the ridge given by the vector ridge IS THIS A FACET NORMAL OR AN INTERIOR POINT?
    2020 * This gives a list of symmetry invariant points under the actions keeping the link fixed.
    2121 */
    22         virtual std::list<ZVector> link(ZVector const &ridgeVector)=0;
    23         virtual ZCone & refToPolyhedralCone()=0;
     22        virtual std::list<ZVector> link(ZVector const &ridgeVector)=0;
     23        virtual ZCone & refToPolyhedralCone()=0;
    2424/**
    2525 * If there is no cone state data for the traverser, half of the changeCone() calls can be avoided.
    2626 * That this is a valid of optimization for the ConeTraverser is indicated returning true in the following function.
    2727 */
    28         virtual bool hasNoState()const;
     28        virtual bool hasNoState()const;
    2929};
    3030
    … …  
    3232{
    3333public:
    34         virtual bool process(FanTraverser &traverser)=0;
     34        virtual bool process(FanTraverser &traverser)=0;
    3535};
    3636
    3737class FanBuilder : public Target
    3838{
    39         ZFan coneCollection;
     39        ZFan coneCollection;
    4040public:
    41         ZFan const &getFanRef(){return coneCollection;}
    42         FanBuilder(int n, SymmetryGroup const &sym);
    43         bool process(FanTraverser &Traverser);
     41        ZFan const &getFanRef(){return coneCollection;}
     42        FanBuilder(int n, SymmetryGroup const &sym);
     43        bool process(FanTraverser &Traverser);
    4444};
    4545

  • gfanlib/gfanlib_tropicalhomotopy.h

    r47a774 r15813d  
    1616class Flags{
    1717public:
    18         static const bool computeDotProductInMatrix=true;
     18        static const bool computeDotProductInMatrix=true;
    1919 };
    2020
    … …  
    5858static Matrix<mvtyp> simplex(int n, mvtyp d)
    5959{
    60           Matrix<mvtyp> ret(n,n+1);
    61           for(int i=0;i<n;i++)ret[i][i+1]=d;
    62           return ret;
     60          Matrix<mvtyp> ret(n,n+1);
     61          for(int i=0;i<n;i++)ret[i][i+1]=d;
     62          return ret;
    6363}
    6464
    … …  
    6666
    6767template<class mvtyp, class mvtypDouble, class mvtypDivisor>
    68         class SingleTropicalHomotopyTraverser{
    69                 class InequalityComparisonResult{//actual comparison functions were moved to the InequalityTable
    70                 public:
    71                           bool empty;
    72                           int configurationIndex;//used for finding best
    73                           int columnIndex;//used for finding best
    74                 };
    75         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.
    76         {
    77                 /* All methods are marked to show if they can overflow without throwing/asserting.
    78                 * Overflowing methods at the moment are:
    79                 *  replaceFirstOrSecond:       subroutine calls may overflow (dotDivVector)
    80                 *  compareInequalities:           only if target entries are too big
    81                 *  dotVector:                     only if target entries are too big
    82                 *  setChoicesFromEarlierHomotopy: only if tuple entries are too big
    83                 */
    84                 std::vector<Matrix<mvtyp> > tuple;
    85                 std::vector<int> offsets;
    86                 std::vector<std::pair<int,int> > choices;
    87                 Matrix<mvtyp> A;//one row for each potential inequality, to store entries with indices in chosen
    88                 Vector<mvtyp> tempA;
    89                 Vector<mvtyp> Abounds;// a negative bound for each row of A, bounding the absolute value of the rows;
    90                 std::vector<mvtyp> svec;//used locally
    91                 int subconfigurationIndex;
    92                 mvtyp denominator;
    93                 int m;
    94                 int k;
    95                 bool isLegalIndex(int subconfigurationIndex, int columnIndex)const
    96                 {// Cannot overflow
    97                         return choices[subconfigurationIndex].first!=columnIndex && choices[subconfigurationIndex].second!=columnIndex;
    98                 }
    99                   mvtyp dotVector(int subconfigurationIndex, int columnIndex, Vector<mvtyp> const &target, int onlyK=-1)const
    100                   {   // May overflow if entries of target are too big.
    101                           //if onlyK!=-1 then only the onlyKth subconfiguration is considered
    102                           mvtypDouble ret;
    103                           if(onlyK!=-1)
    104                           {
    105                                   if(onlyK==subconfigurationIndex)
    106                                   {
    107                                           int i=subconfigurationIndex;
    108                                           ret+=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].second+offsets[i]));
    109                                           ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));
    110                                           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.
    111                                           ret+=extendedMultiplication(denominator,target.UNCHECKEDACCESS(columnIndex+offsets[i]));
    112                                           return ret.castToSingle();
    113                                   }
    114                                   else
    115                                   {
    116                                           int i=onlyK;
    117                                           if(target.UNCHECKEDACCESS((choices)[i].first+offsets[i]).isNonZero())
    118                                           {
    119                                                   ret+=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].second+offsets[i]));
    120                                                   ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));
    121                                           }
    122                                           return ret.castToSingle();
    123                                   }
    124                           }
    125                           for(int i=0;i<tuple.size();i++)
    126                           {
    127                                   ret+=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].second+offsets[i]));
    128                                   if(i==subconfigurationIndex)
    129                                   {
    130                                           ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));
    131                                           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.
    132                                           ret+=extendedMultiplication(denominator,target.UNCHECKEDACCESS(columnIndex+offsets[i]));
    133                                   }
    134                                   else
    135                                   {
    136                                           ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));
    137                                   }
    138                           }
    139                           return ret.castToSingle();
    140                   }
    141                   void assignDotProducts(Vector<mvtyp> const &target, int onlyK=-1)
    142                   {// Cannot overflow
    143                           int J=0;
    144                           for(int i=0;i<k;i++)
    145                                   for(int j=0;j<tuple[i].getWidth();j++,J++)
    146                                           A[k][J]=dotVector(i,j,target,onlyK);
    147                   }
    148                   bool isReverseLexInvertedLessThanZero(int subconfigurationIndex, int columnIndex)const __attribute__ ((always_inline))//As in ReverseLexicographicInvertedTermOrder. Compare against zero
    149                   {// Cannot overflow
    150                           int i;
    151                           int index=columnIndex+offsets[subconfigurationIndex];
    152                           for(i=0;i<subconfigurationIndex;i++)
    153                                   if(A.UNCHECKEDACCESS(i,index).isNonZero())
    154                                   {
    155                                           if(choices[i].first<choices[i].second)
    156                                                   return A.UNCHECKEDACCESS(i,index).isNegative();
    157                                           else
    158                                                   return A.UNCHECKEDACCESS(i,index).isPositive();
    159                                   }
    160 
    161                                   mvtyp a=A.UNCHECKEDACCESS(i,index);
    162                                   {
    163                                           int firstIndex=choices[i].first;
    164                                           int secondIndex=choices[i].second;
    165                                           int thirdIndex=columnIndex;
    166                                           mvtyp firstValue=-a-denominator;
    167                                           mvtyp secondValue=a;
    168                                           mvtyp thirdValue=denominator;
    169 
    170                                           // Bubble sort
    171                                           if(secondIndex<firstIndex)
    172                                           {
    173                                                   std::swap(secondIndex,firstIndex);
    174                                                   std::swap(secondValue,firstValue);
    175                                           }
    176                                           if(thirdIndex<secondIndex)
    177                                           {
    178                                                   std::swap(secondIndex,thirdIndex);
    179                                                   std::swap(secondValue,thirdValue);
    180                                           }
    181                                           if(secondIndex<firstIndex)
    182                                           {
    183                                                   std::swap(secondIndex,firstIndex);
    184                                                   std::swap(secondValue,firstValue);
    185                                           }
    186 
    187                                           if(firstValue.isNonZero())
    188                                                   return firstValue.isPositive();
    189                                           if(secondValue.isNonZero())
    190                                                   return secondValue.isPositive();
    191                                           if(thirdValue.isNonZero())
    192                                                   return thirdValue.isPositive();
    193                                   }
    194                                   i++;
    195                                   for(;i<k;i++)
    196                                           if(A.UNCHECKEDACCESS(i,index).isNonZero())
    197                                           {
    198                                                   if(choices[i].first<choices[i].second)
    199                                                           return A.UNCHECKEDACCESS(i,index).isNegative();
    200                                                   else
    201                                                           return A.UNCHECKEDACCESS(i,index).isPositive();
    202                                           }
    203 
    204 
    205                                   return false;
    206                   }
    207         public:
    208                         void computeABounds()
    209                         {//Cannot overflow
    210                                 for(int i=0;i<A.getHeight();i++)
    211                                         Abounds[i]=mvtyp::computeNegativeBound(&(A[i][0]),A.getWidth());
    212                         }
    213                         void checkABounds()const//remove?
    214                         {//Cannot overflow
    215                                 for(int i=0;i<A.getHeight();i++)
    216                                 {
    217                                         mvtyp M=0;
    218                                         mvtyp m=0;
    219                                         for(int j=0;j<A.getWidth();j++)
    220                                         {
    221                                                 if(M<=A[i][j])M=A[i][j];
    222                                                 if(A[i][j]<=m)m=A[i][j];
    223                                         }
    224                                         assert(Abounds[i]<=m);
    225                                         assert(Abounds[i]<=-M);
    226                                 }
    227                         }
    228                   mvtypDouble getCoordinateOfInequality(int subconfigurationIndex, int columnIndex, int i, int j)const
    229                   {// Cannot overflows
    230                           //get (i,j)th coordinate of (subconfigurationIndex,columnIndex)th inequality
    231                           if(i==subconfigurationIndex)
    232                           {
    233                                   if(choices[i].first==j)return -A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+columnIndex).extend()-denominator.extend();
    234                                   else if(choices[i].second==j)return A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+columnIndex).extend();
    235                                   else if(j==columnIndex)return denominator.extend();
    236                                   else return mvtypDouble(0);
    237                           }
    238                           else
    239                                   if(choices[i].first==j)return -A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+columnIndex).extend();
    240                                   else if(choices[i].second==j)return A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+columnIndex).extend();
    241                                   else return mvtypDouble(0);
    242                   }
    243                   int sort2uniquely(int *v, int a, int b)const//a and b different
    244                   {// Cannot overflow
    245                           v[a>b]=a;
    246                           v[b>a]=b;
    247                           return 2;
    248                   }
    249                   int sort3uniquely(int *v, int a, int b, int c)const//a and b and c different
    250                   {// Cannot overflow
    251                           v[(a>b)+int(a>c)]=a;
    252                           v[(b>a)+int(b>c)]=b;
    253                           v[(c>a)+int(c>b)]=c;
    254                           return 3;
    255                   }
    256                   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
    257                   {// Cannot overflow
    258                           if(c!=d)
    259                           {
    260                           v[(a>b)+int(a>c)+int(a>d)]=a;
    261                           v[(b>a)+int(b>c)+int(b>d)]=b;
    262                           v[(c>a)+int(c>b)+int(c>d)]=c;
    263                           v[(d>a)+int(d>b)+int(d>c)]=d;
    264                           return 4;
    265                           }
    266                           else return sort3uniquely(v,a,b,c);
    267                   }
    268                   bool compareReverseLexicographicInverted(int i1, int j1, int i2, int j2, mvtyp s1, mvtyp s2)const//s1 and s2 are always negative
    269                   {// cannot overflow
    270                           for(int i=0;i<k;i++)
    271                           {
    272                                           if(i1!=i && i2!=i)
    273                                           {
    274                                                   int temp=determinantSign1(A.UNCHECKEDACCESS(i,offsets[i1]+j1),s1,A.UNCHECKEDACCESS(i,offsets[i2]+j2),s2);
    275                                                   if(temp)
    276                                                   {
    277                                                           if(choices[i].first<choices[i].second)
    278                                                                   return temp<0;
    279                                                           else
    280                                                                   return temp>0;
    281                                                   }
    282                                           }
    283                                   int indices[4];
    284                                   int F=choices[i].first;
    285                                   int S=choices[i].second;
    286                                   int toCheck;
    287                                   if(i1==i)
    288                                           if(i2==i)
    289                                                   toCheck=sort4uniquely(indices,F,S,j1,j2);
    290                                           else
    291                                                   toCheck=sort3uniquely(indices,F,S,j1);
    292                                   else
    293                                           if(i2==i)
    294                                                   toCheck=sort3uniquely(indices,F,S,j2);
    295                                           else
    296                                                   toCheck=sort2uniquely(indices,F,S);
    297 
    298                                   for(int J=0;J<toCheck;J++)
    299                                   {
    300                                           int j=indices[J];
    301                                           int temp=determinantSign2(getCoordinateOfInequality(i1,j1,i,j),s1,getCoordinateOfInequality(i2,j2,i,j),s2);
    302                                           if(temp>0)
    303                                                   return true;
    304                                           else if(temp<0)
    305                                                   return false;
    306                                   }
    307                           }
    308                           return false;
    309                   }
    310                   mvtyp getVolume()
    311                   {// Cannot overflow
    312                           return denominator;
    313                   }
    314                 void replaceFirstOrSecond(bool first, int subconfigurationIndex, int newIndex, Vector<mvtyp> const &target)__attribute__ ((always_inline))//updates the inequality table according to replacing first or second by newIndex in the subconfigurationIndex'th configuration
    315                 {// Cannot overflow
    316                         int newIndex2=newIndex;for(int i=0;i<subconfigurationIndex;i++)newIndex2+=tuple[i].getWidth();
    317                         int oldIndex=first?choices[subconfigurationIndex].first:choices[subconfigurationIndex].second;
    318                         (first?choices[subconfigurationIndex].first:choices[subconfigurationIndex].second)=newIndex;
    319 
    320                         mvtyp nextDenominator=(first?A[subconfigurationIndex][newIndex2].extend()+denominator.extend():-A[subconfigurationIndex][newIndex2].extend()).castToSingle();
    321                         mvtyp t=nextDenominator;
    322                         mvtypDivisor denominatorDivisor(denominator);
    323                         for(int c=0;c<k+Flags::computeDotProductInMatrix;c++)tempA.UNCHECKEDACCESS(c)=A.UNCHECKEDACCESS(c,newIndex2);
    324 
    325                         (-Abounds[subconfigurationIndex].extend()).castToSingle();//This conversion will fail if the following fails
    326                         for(int u=0;u<m;u++)
    327                                 svec[u]=first?-A.UNCHECKEDACCESS(subconfigurationIndex,u):A.UNCHECKEDACCESS(subconfigurationIndex,u);
    328                         mvtyp svecBound=Abounds[subconfigurationIndex];
    329 
    330                         if(first)
    331                                 for(int j=0;j<tuple[subconfigurationIndex].getWidth();j++)
    332                                         svec[offsets[subconfigurationIndex]+j].subWithOverflowChecking(denominator);
    333                         svecBound.subWithOverflowChecking(denominator);
    334 
    335 
    336                         for(int a=0;a<k+Flags::computeDotProductInMatrix;a++)
    337                                 if(first||a!=subconfigurationIndex)
    338                                 Abounds.UNCHECKEDACCESS(a)=mvtyp::dotDivVector(&A.UNCHECKEDACCESS(a,0),&(svec[0]),t,tempA.UNCHECKEDACCESS(a),denominatorDivisor,m,Abounds[a],svecBound);
    339         //                      for(int u=0;u<m;u++)
    340         //                              *A.UNCHECKEDACCESSRESTRICT(a,u)=dotDiv(svec[u],tempA.UNCHECKEDACCESS(a),t,A.UNCHECKEDACCESS(a,u),denominatorDivisor);
    341 
    342 
    343                         {
    344                                 for(int a=0;a<k+Flags::computeDotProductInMatrix;a++)
    345                                         {
    346                                                 A[a][offsets[subconfigurationIndex]+oldIndex]=tempA[a].negatedWithOverflowChecking();
    347                                                 Abounds[a]=min(Abounds[a],negabs(tempA[a]));
    348                                         }
    349 
    350                                 if(!first)
    351                                 {
    352                                         A[subconfigurationIndex][offsets[subconfigurationIndex]+oldIndex]=denominator.negatedWithOverflowChecking();
    353                                         Abounds[subconfigurationIndex].subWithOverflowChecking(denominator);
    354                                 }
    355                         }
    356                         denominator=nextDenominator;
    357 
    358                         // We clear these unused entries of A to ensure that these columns are not chosen when comparing inequalities
    359                         for(int i=0;i<k+Flags::computeDotProductInMatrix;i++)
    360                         {
    361                                 A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+choices[subconfigurationIndex].first)=0;
    362                                 A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+choices[subconfigurationIndex].second)=0;
    363                         }
    364                 }
    365                 void replaceFirst(int subconfigurationIndex, int newIndex, Vector<mvtyp> const &target){replaceFirstOrSecond(true,subconfigurationIndex,newIndex,target);}
    366                 void replaceSecond(int subconfigurationIndex, int newIndex, Vector<mvtyp> const &target){replaceFirstOrSecond(false,subconfigurationIndex,newIndex,target);}
    367 
    368                 InequalityTable(std::vector<Matrix<mvtyp> > const &tuple_, int subconfigurationIndex_):
    369                         tempA(tuple_.size()+Flags::computeDotProductInMatrix),
    370                         tuple(tuple_),
    371                         choices(tuple_.size()),
    372                         subconfigurationIndex(subconfigurationIndex_),
    373                         offsets(tuple_.size())
    374                 {// Cannot overflow
    375                         k=tuple.size();
    376                         m=0;
    377                         for(int i=0;i<tuple.size();i++)m+=tuple[i].getWidth();
    378                         svec.resize(m);
    379                         A=Matrix<mvtyp>(k+Flags::computeDotProductInMatrix,m);
    380                         {int offset=0;for(int i=0;i<tuple.size();i++){offsets[i]=offset;offset+=tuple[i].getWidth();}}
    381                         Abounds=Vector<mvtyp>(k+Flags::computeDotProductInMatrix);
    382                 }
    383                 void setChoicesInitially()
    384                 {// Cannot overflow
    385                         //THIS WILL ONLY WORK FOR THE STARTING CONFIGURATION
    386                         //sets denominators,A and choices (these must have been initialized to the right sizes)
    387                         for(int i=0;i<k;i++)
    388                                 choices[i]=std::pair<int, int> (i+0,i+1);
    389                         for(int i=0;i<m;i++)
    390                                 for(int j=0;j<k;j++)
    391                                         A[j][i]=0;
    392                         //we follow the lemma in the article. Why does the proof of the lemma split into 3 cases and not just 2?
    393                         int a=0;
    394                         for(int i=0;i<k;i++)
    395                                 for(int gamma=0;gamma<tuple[i].getWidth();gamma++,a++)
    396                                         if(gamma>i+1)
    397                                                 for(int ii=i;ii<gamma;ii++)
    398                                                         A[ii][a]=-1;
    399                                         else
    400                                                 for(int ii=gamma;ii<i;ii++)
    401                                                         A[ii][a]=1;
    402                         denominator=1;
    403                         for(int i=0;i<k;i++)Abounds[i]=-1;
    404         //              assignDotProducts(target);
    405                 }
    406                 void compareInequalities(InequalityComparisonResult &result, Vector<mvtyp> const &target, int onlyK=-1)
    407                 {// Can only overflow if target entries are too big. Actually it seems that it is not this function that will overflow but dotVector.
    408                         bool empty=true;
    409                         int bestConfigurationIndex=-1;
    410                         int bestColumnIndex=-1;
    411                         mvtyp targetDotBest=0;
    412 
    413                         for(int i=0;i<k;i++)
    414                         {
    415                                 // TO DO: exchange the line with the template parameter version:
    416 //                  gfan::Matrix<mvtyp>::const_RowRef Ak=const_cast<const gfan::Matrix<mvtyp>&>(A)[k];
    417                     gfan::Matrix<CircuitTableInt32>::const_RowRef Ak=const_cast<const gfan::Matrix<mvtyp>&>(A)[k];
    418                     int offsetsi=offsets[i];
    419                     int tupleiwidth=tuple[i].getWidth();
    420 
    421                                 if(onlyK!=-1)if(i!=onlyK)continue;
    422                                 for(int j=0;j<tupleiwidth;j++)
    423                                         if(Flags::computeDotProductInMatrix || isLegalIndex(i,j))//unused inequalities will have value 0. Therefore isLegalIndex(i,j) is not required if values are stored.
    424                                         {
    425                                                 mvtyp ineqDotTarget=Flags::computeDotProductInMatrix?Ak.UNCHECKEDACCESS(offsetsi+j):dotVector(i,j,target,onlyK);
    426                                                 if(ineqDotTarget.isNegative())
    427                                                 {
    428                                                         if(!isReverseLexInvertedLessThanZero(i,j))
    429                                                         {
    430                                                                 if(empty||compareReverseLexicographicInverted(bestConfigurationIndex,bestColumnIndex,i,j,ineqDotTarget,targetDotBest))
    431                                                                 {
    432                                                                         targetDotBest=ineqDotTarget;
    433                                                                         empty=false;
    434                                                                         bestConfigurationIndex=i;
    435                                                                         bestColumnIndex=j;
    436                                                                 }
    437                                                         }
    438                 //                                      assert(!ineq.isZero());
    439                                                 }
    440                                         }
    441                         }
    442                         result.empty=empty;
    443                         result.configurationIndex=bestConfigurationIndex;
    444                         result.columnIndex=bestColumnIndex;
    445         //              assert(denominator>0);
    446                 }
    447                 void setChoicesFromEarlierHomotopy(InequalityTable const &parent, mvtyp degreeScaling, Vector<mvtyp> const &target)
    448                 {       // Overflows may happen, but only if the tuple entries are too big.
    449                         // Notice that the code below has overflow checks everywhere, except in the loop marked with a comment.
    450 
    451                         //sets denominators,A and choices (these must have been initialized to the right sizes
    452                         //columns have been added to configuration this->subconfigurationIndex
    453                         //for that reason we need to introduce new circuits. Old circuits are preserved.
    454                         //chioices are "relative" so no update is needed.
    455 
    456                         choices=parent.choices;
    457                         int numberToDrop=(subconfigurationIndex!=0) ? numberToDrop=k+1 : 0;
    458 
    459                         choices[subconfigurationIndex-1].first-=numberToDrop;
    460                         choices[subconfigurationIndex-1].second-=numberToDrop;
    461 
    462                         denominator=parent.denominator;
    463                         int offsetOld=0;
    464                         int offsetNew=0;
    465                         for(int i=0;i<k;i++)
    466                         {
    467                                 int localNumberToDrop=0;
    468                                 if(i==subconfigurationIndex-1)
    469                                         localNumberToDrop=numberToDrop;
    470                                 for(int a=0;a<A.getHeight()-Flags::computeDotProductInMatrix;a++)
    471                                         if(a==subconfigurationIndex)
    472                                                 for(int j=0;j<parent.tuple[i].getWidth()-localNumberToDrop;j++)
    473                                                         A.UNCHECKEDACCESS(a,offsetNew+j)=parent.A.UNCHECKEDACCESS(a,offsetOld+j+localNumberToDrop);
    474                                         else
    475                                         {
    476                                                 // We can get an estimated 3 percent speed up by using known bounds to check if any of the multiplications below can overflow.
    477                                                 // Moreover, these multiplications can be moved outside the outer loop, to give better vectorisation.
    478                                                 for(int j=0;j<parent.tuple[i].getWidth()-localNumberToDrop;j++)
    479                                                         A.UNCHECKEDACCESS(a,offsetNew+j)=extendedMultiplication(degreeScaling,parent.A.UNCHECKEDACCESS(a,offsetOld+j+localNumberToDrop)).castToSingle();
    480                                         }
    481                                 if(i==subconfigurationIndex)
    482                                 {
    483                                         for(int j=parent.tuple[i].getWidth();j<tuple[i].getWidth();j++)
    484                                         {
    485                                                 for(int a=0;a<A.getHeight()-Flags::computeDotProductInMatrix;a++)
    486                                                         A.UNCHECKEDACCESS(a,offsetNew+j)=0;
    487                                                 {
    488                                                         int b=choices[subconfigurationIndex].second-1;
    489                                                         if(b>=0)
    490                                                                 A.UNCHECKEDACCESS(subconfigurationIndex,offsetNew+j).msubWithOverflowChecking(tuple[i].UNCHECKEDACCESS(b,j),denominator);
    491                                                 }
    492                                                 for(int a=0;a<k+Flags::computeDotProductInMatrix;a++)
    493                                                 {   // If tuple entries are not bounded, overflow can happen in the following loop
    494                                                         mvtypDouble tempDouble=A.UNCHECKEDACCESS(a,offsetNew+j).extend();
    495                                                         for(int b=0;b<k;b++)
    496                                                                 if(choices[subconfigurationIndex].second!=b+1 &&choices[subconfigurationIndex].first!=b+1)
    497                                                                         tempDouble+=extendedMultiplication(tuple[i].UNCHECKEDACCESS(b,j),A.UNCHECKEDACCESS(a,offsetNew+b+1));
    498                                                         A.UNCHECKEDACCESS(a,offsetNew+j)=tempDouble.castToSingle();
    499                                                 }
    500                                                 mvtypDouble left=degreeScaling.extend();
    501                                                 for(int b=0;b<k;b++)
    502                                                         left-=mvtyp(tuple[i].UNCHECKEDACCESS(b,j)).extend();
    503 
    504                                                 mvtyp leftS=left.castToSingle();
    505 
    506                                                 if(choices[subconfigurationIndex].second==0)
    507                                                         A.UNCHECKEDACCESS(subconfigurationIndex,offsetNew+j).msubWithOverflowChecking(leftS,denominator);
    508                                                 else if(choices[subconfigurationIndex].first!=0)
    509                                                         for(int a=0;a<k+Flags::computeDotProductInMatrix;a++)
    510                                                                 A.UNCHECKEDACCESS(a,offsetNew+j).maddWithOverflowChecking(leftS,mvtyp(A.UNCHECKEDACCESS(a,offsetNew)));
    511                                         }
    512                                         for(int j=0;j<parent.tuple[i].getWidth();j++)// <-- this loop has side effects on addExpressionOfCeb() above. Therefore it is placed after the loop above.
    513                                                 for(int a=0;a<A.getHeight()-Flags::computeDotProductInMatrix;a++)
    514                                                         A.UNCHECKEDACCESS(a,offsetNew+j)=extendedMultiplication(A.UNCHECKEDACCESS(a,offsetNew+j),degreeScaling).castToSingle();
    515                                 }
    516                                 offsetOld+=parent.tuple[i].getWidth();
    517                                 offsetNew+=tuple[i].getWidth();
    518                         }
    519                         denominator=extendedMultiplication(denominator,degreeScaling).castToSingle();
    520                         if(Flags::computeDotProductInMatrix)assignDotProducts(target,subconfigurationIndex);
    521                         computeABounds();
    522                         for(int a=0;a<k;a++)
    523                         for(int i=0;i<k+Flags::computeDotProductInMatrix;i++)
    524                         {
    525                                 A[i][offsets[a]+choices[a].first]=0;
    526                                 A[i][offsets[a]+choices[a].second]=0;
    527                         }
    528                 }
    529         };
    530                 struct StackItem{
    531                         int columnIndex;
    532                         int configurationIndex;
    533                         bool b;
    534                         int choice;
    535                         bool useFirstChanged,useSecondChanged;
    536                         StackItem(int columnIndex_, int configurationIndex_, bool b_, int choice_, bool useFirstChanged_, bool useSecondChanged_):
    537                                 columnIndex(columnIndex_),
    538                                 configurationIndex(configurationIndex_),
    539                                 b(b_),
    540                                 choice(choice_),
    541                                 useFirstChanged(useFirstChanged_),
    542                                 useSecondChanged(useSecondChanged_)
    543                         {
    544                         }
    545                 };
    546         public:
    547                 std::vector<std::pair<int,int> > choices;
    548                 Vector<mvtyp> target;
    549                 bool useFirstChanged;
    550                 bool useSecondChanged;
    551                 std::vector<StackItem> stack;
    552                 int eliminatedK;
    553                 int eliminatedKOffset;
    554                 std::vector<Matrix<mvtyp> > tuple;
    555                 std::vector<int> offsets;
    556                 int m;
    557                 InequalityComparisonResult result;
    558                 InequalityTable inequalityTable;
    559                 void constructInequalityTableFromParent(InequalityTable const &parentTable, mvtyp degreeScaling)
    560                 {
    561                         inequalityTable.setChoicesFromEarlierHomotopy(parentTable, degreeScaling, target);
    562                 }
    563                 void constructInequalityTableInitially(mvtyp degreeScaling)
    564                 {
    565                         std::vector<Matrix<mvtyp> > tempTuple;for(int i=0;i<tuple.size();i++)tempTuple.push_back(simplex<mvtyp>(tuple.size(),1));
    566                         InequalityTable tempTable(tempTuple,-1);
    567                         tempTable.setChoicesInitially();
    568                         constructInequalityTableFromParent(tempTable,degreeScaling);
    569                 }
    570                 SingleTropicalHomotopyTraverser(std::vector<Matrix<mvtyp> > const &tuple_, int m_, std::vector<std::pair<int,int> > const &choices_, Vector<mvtyp> const &target_, int eliminatedK_):
    571                         choices(choices_),
    572                         target(target_),
    573                         eliminatedK(eliminatedK_),
    574                         tuple(tuple_),
    575                         m(m_),
    576                         inequalityTable(tuple,eliminatedK),
    577                         offsets(tuple_.size())
    578                 {
    579                         eliminatedKOffset=0;for(int i=0;i<eliminatedK;i++)eliminatedKOffset+=tuple_[i].getWidth();
    580                         {int offset=0;for(int i=0;i<tuple.size();i++){offsets[i]=offset;offset+=tuple[i].getWidth();}}
    581 
    582                         /* We also need to check that the input is within the required bounds
    583                         * This is important to not have overflows in:
    584                         *   dotVector
    585                         *   setChoicesFromEarlierHomotopy
    586                         * The number of summands is bounded by 3k+4*/
    587 
    588                         bool isOK=mvtyp::isSuitableSummationNumber(3*tuple.size()+4);
    589                         for(int i=0;i<target.size();i++)isOK&=mvtyp(target[i]).fitsInHalf();
    590                         for(int i=0;i<tuple.size();i++)
    591                                 for(int j=0;j<tuple[i].getHeight();j++)
    592                                         for(int k=0;k<tuple[i].getWidth();k++)
    593                                                 isOK&=mvtyp(tuple[i][j][k]).fitsInHalf();
    594                         if(!isOK)throw MVMachineIntegerOverflow;
    595                 }
    596                 virtual void process()
    597                 {
    598                 }
    599                 bool findOutgoingAndProcess(bool doProcess)//sets up useFirstChanged and useSecondChanged and processes if process is true
    600                 {//returns true if we are at a leaf
    601 //                      inequalityTable.checkABounds();
    602                         useFirstChanged=false;
    603                         useSecondChanged=false;
    604                         int onlyK=-1;
     68        class SingleTropicalHomotopyTraverser{
     69                class InequalityComparisonResult{//actual comparison functions were moved to the InequalityTable
     70                public:
     71                          bool empty;
     72                          int configurationIndex;//used for finding best
     73                          int columnIndex;//used for finding best
     74                };
     75        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.
     76        {
     77                /* All methods are marked to show if they can overflow without throwing/asserting.
     78                * Overflowing methods at the moment are:
     79                *  replaceFirstOrSecond:       subroutine calls may overflow (dotDivVector)
     80                *  compareInequalities:           only if target entries are too big
     81                *  dotVector:                     only if target entries are too big
     82                *  setChoicesFromEarlierHomotopy: only if tuple entries are too big
     83                */
     84                std::vector<Matrix<mvtyp> > tuple;
     85                std::vector<int> offsets;
     86                std::vector<std::pair<int,int> > choices;
     87                Matrix<mvtyp> A;//one row for each potential inequality, to store entries with indices in chosen
     88                Vector<mvtyp> tempA;
     89                Vector<mvtyp> Abounds;// a negative bound for each row of A, bounding the absolute value of the rows;
     90                std::vector<mvtyp> svec;//used locally
     91                int subconfigurationIndex;
     92                mvtyp denominator;
     93                int m;
     94                int k;
     95                bool isLegalIndex(int subconfigurationIndex, int columnIndex)const
     96                {// Cannot overflow
     97                        return choices[subconfigurationIndex].first!=columnIndex && choices[subconfigurationIndex].second!=columnIndex;
     98                }
     99                  mvtyp dotVector(int subconfigurationIndex, int columnIndex, Vector<mvtyp> const &target, int onlyK=-1)const
     100                  {   // May overflow if entries of target are too big.
     101                          //if onlyK!=-1 then only the onlyKth subconfiguration is considered
     102                          mvtypDouble ret;
     103                          if(onlyK!=-1)
     104                          {
     105                                  if(onlyK==subconfigurationIndex)
     106                                  {
     107                                          int i=subconfigurationIndex;
     108                                          ret+=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].second+offsets[i]));
     109                                          ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));
     110                                          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.
     111                                          ret+=extendedMultiplication(denominator,target.UNCHECKEDACCESS(columnIndex+offsets[i]));
     112                                          return ret.castToSingle();
     113                                  }
     114                                  else
     115                                  {
     116                                          int i=onlyK;
     117                                          if(target.UNCHECKEDACCESS((choices)[i].first+offsets[i]).isNonZero())
     118                                          {
     119                                                  ret+=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].second+offsets[i]));
     120                                                  ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));
     121                                          }
     122                                          return ret.castToSingle();
     123                                  }
     124                          }
     125                          for(int i=0;i<tuple.size();i++)
     126                          {
     127                                  ret+=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].second+offsets[i]));
     128                                  if(i==subconfigurationIndex)
     129                                  {
     130                                          ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));
     131                                          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.
     132                                          ret+=extendedMultiplication(denominator,target.UNCHECKEDACCESS(columnIndex+offsets[i]));
     133                                  }
     134                                  else
     135                                  {
     136                                          ret-=extendedMultiplication(A.UNCHECKEDACCESS(i,columnIndex+offsets[subconfigurationIndex]),target.UNCHECKEDACCESS((choices)[i].first+offsets[i]));
     137                                  }
     138                          }
     139                          return ret.castToSingle();
     140                  }
     141                  void assignDotProducts(Vector<mvtyp> const &target, int onlyK=-1)
     142                  {// Cannot overflow
     143                          int J=0;
     144                          for(int i=0;i<k;i++)
     145                                  for(int j=0;j<tuple[i].getWidth();j++,J++)
     146                                          A[k][J]=dotVector(i,j,target,onlyK);
     147                  }
     148                  bool isReverseLexInvertedLessThanZero(int subconfigurationIndex, int columnIndex)const __attribute__ ((always_inline))//As in ReverseLexicographicInvertedTermOrder. Compare against zero
     149                  {// Cannot overflow
     150                          int i;
     151                          int index=columnIndex+offsets[subconfigurationIndex];
     152                          for(i=0;i<subconfigurationIndex;i++)
     153                                  if(A.UNCHECKEDACCESS(i,index).isNonZero())
     154                                  {
     155                                          if(choices[i].first<choices[i].second)
     156                                                  return A.UNCHECKEDACCESS(i,index).isNegative();
     157                                          else
     158                                                  return A.UNCHECKEDACCESS(i,index).isPositive();
     159                                  }
     160
     161                                  mvtyp a=A.UNCHECKEDACCESS(i,index);
     162                                  {
     163                                          int firstIndex=choices[i].first;
     164                                          int secondIndex=choices[i].second;
     165                                          int thirdIndex=columnIndex;
     166                                          mvtyp firstValue=-a-denominator;
     167                                          mvtyp secondValue=a;
     168                                          mvtyp thirdValue=denominator;
     169
     170                                          // Bubble sort
     171                                          if(secondIndex<firstIndex)
     172                                          {
     173                                                  std::swap(secondIndex,firstIndex);
     174                                                  std::swap(secondValue,firstValue);
     175                                          }
     176                                          if(thirdIndex<secondIndex)
     177                                          {
     178                                                  std::swap(secondIndex,thirdIndex);
     179                                                  std::swap(secondValue,thirdValue);
     180                                          }
     181                                          if(secondIndex<firstIndex)
     182                                          {
     183                                                  std::swap(secondIndex,firstIndex);
     184                                                  std::swap(secondValue,firstValue);
     185                                          }
     186
     187                                          if(firstValue.isNonZero())
     188                                                  return firstValue.isPositive();
     189                                          if(secondValue.isNonZero())
     190                                                  return secondValue.isPositive();
     191                                          if(thirdValue.isNonZero())
     192                                                  return thirdValue.isPositive();
     193                                  }
     194                                  i++;
     195                                  for(;i<k;i++)
     196                                          if(A.UNCHECKEDACCESS(i,index).isNonZero())
     197                                          {
     198                                                  if(choices[i].first<choices[i].second)
     199                                                          return A.UNCHECKEDACCESS(i,index).isNegative();
     200                                                  else
     201                                                          return A.UNCHECKEDACCESS(i,index).isPositive();
     202                                          }
     203
     204
     205                                  return false;
     206                  }
     207        public:
     208                        void computeABounds()
     209                        {//Cannot overflow
     210                                for(int i=0;i<A.getHeight();i++)
     211                                        Abounds[i]=mvtyp::computeNegativeBound(&(A[i][0]),A.getWidth());
     212                        }
     213                        void checkABounds()const//remove?
     214                        {//Cannot overflow
     215                                for(int i=0;i<A.getHeight();i++)
     216                                {
     217                                        mvtyp M=0;
     218                                        mvtyp m=0;
     219                                        for(int j=0;j<A.getWidth();j++)
     220                                        {
     221                                                if(M<=A[i][j])M=A[i][j];
     222                                                if(A[i][j]<=m)m=A[i][j];
     223                                        }
     224                                        assert(Abounds[i]<=m);
     225                                        assert(Abounds[i]<=-M);
     226                                }
     227                        }
     228                  mvtypDouble getCoordinateOfInequality(int subconfigurationIndex, int columnIndex, int i, int j)const
     229                  {// Cannot overflows
     230                          //get (i,j)th coordinate of (subconfigurationIndex,columnIndex)th inequality
     231                          if(i==subconfigurationIndex)
     232                          {
     233                                  if(choices[i].first==j)return -A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+columnIndex).extend()-denominator.extend();
     234                                  else if(choices[i].second==j)return A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+columnIndex).extend();
     235                                  else if(j==columnIndex)return denominator.extend();
     236                                  else return mvtypDouble(0);
     237                          }
     238                          else
     239                                  if(choices[i].first==j)return -A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+columnIndex).extend();
     240                                  else if(choices[i].second==j)return A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+columnIndex).extend();
     241                                  else return mvtypDouble(0);
     242                  }
     243                  int sort2uniquely(int *v, int a, int b)const//a and b different
     244                  {// Cannot overflow
     245                          v[a>b]=a;
     246                          v[b>a]=b;
     247                          return 2;
     248                  }
     249                  int sort3uniquely(int *v, int a, int b, int c)const//a and b and c different
     250                  {// Cannot overflow
     251                          v[(a>b)+int(a>c)]=a;
     252                          v[(b>a)+int(b>c)]=b;
     253                          v[(c>a)+int(c>b)]=c;
     254                          return 3;
     255                  }
     256                  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
     257                  {// Cannot overflow
     258                          if(c!=d)
     259                          {
     260                          v[(a>b)+int(a>c)+int(a>d)]=a;
     261                          v[(b>a)+int(b>c)+int(b>d)]=b;
     262                          v[(c>a)+int(c>b)+int(c>d)]=c;
     263                          v[(d>a)+int(d>b)+int(d>c)]=d;
     264                          return 4;
     265                          }
     266                          else return sort3uniquely(v,a,b,c);
     267                  }
     268                  bool compareReverseLexicographicInverted(int i1, int j1, int i2, int j2, mvtyp s1, mvtyp s2)const//s1 and s2 are always negative
     269                  {// cannot overflow
     270                          for(int i=0;i<k;i++)
     271                          {
     272                                          if(i1!=i && i2!=i)
     273                                          {
     274                                                  int temp=determinantSign1(A.UNCHECKEDACCESS(i,offsets[i1]+j1),s1,A.UNCHECKEDACCESS(i,offsets[i2]+j2),s2);
     275                                                  if(temp)
     276                                                  {
     277                                                          if(choices[i].first<choices[i].second)
     278                                                                  return temp<0;
     279                                                          else
     280                                                                  return temp>0;
     281                                                  }
     282                                          }
     283                                  int indices[4];
     284                                  int F=choices[i].first;
     285                                  int S=choices[i].second;
     286                                  int toCheck;
     287                                  if(i1==i)
     288                                          if(i2==i)
     289                                                  toCheck=sort4uniquely(indices,F,S,j1,j2);
     290                                          else
     291                                                  toCheck=sort3uniquely(indices,F,S,j1);
     292                                  else
     293                                          if(i2==i)
     294                                                  toCheck=sort3uniquely(indices,F,S,j2);
     295                                          else
     296                                                  toCheck=sort2uniquely(indices,F,S);
     297
     298                                  for(int J=0;J<toCheck;J++)
     299                                  {
     300                                          int j=indices[J];
     301                                          int temp=determinantSign2(getCoordinateOfInequality(i1,j1,i,j),s1,getCoordinateOfInequality(i2,j2,i,j),s2);
     302                                          if(temp>0)
     303                                                  return true;
     304                                          else if(temp<0)
     305                                                  return false;
     306                                  }
     307                          }
     308                          return false;
     309                  }
     310                  mvtyp getVolume()
     311                  {// Cannot overflow
     312                          return denominator;
     313                  }
     314                void replaceFirstOrSecond(bool first, int subconfigurationIndex, int newIndex, Vector<mvtyp> const &target)__attribute__ ((always_inline))//updates the inequality table according to replacing first or second by newIndex in the subconfigurationIndex'th configuration
     315                {// Cannot overflow
     316                        int newIndex2=newIndex;for(int i=0;i<subconfigurationIndex;i++)newIndex2+=tuple[i].getWidth();
     317                        int oldIndex=first?choices[subconfigurationIndex].first:choices[subconfigurationIndex].second;
     318                        (first?choices[subconfigurationIndex].first:choices[subconfigurationIndex].second)=newIndex;
     319
     320                        mvtyp nextDenominator=(first?A[subconfigurationIndex][newIndex2].extend()+denominator.extend():-A[subconfigurationIndex][newIndex2].extend()).castToSingle();
     321                        mvtyp t=nextDenominator;
     322                        mvtypDivisor denominatorDivisor(denominator);
     323                        for(int c=0;c<k+Flags::computeDotProductInMatrix;c++)tempA.UNCHECKEDACCESS(c)=A.UNCHECKEDACCESS(c,newIndex2);
     324
     325                        (-Abounds[subconfigurationIndex].extend()).castToSingle();//This conversion will fail if the following fails
     326                        for(int u=0;u<m;u++)
     327                                svec[u]=first?-A.UNCHECKEDACCESS(subconfigurationIndex,u):A.UNCHECKEDACCESS(subconfigurationIndex,u);
     328                        mvtyp svecBound=Abounds[subconfigurationIndex];
     329
     330                        if(first)
     331                                for(int j=0;j<tuple[subconfigurationIndex].getWidth();j++)
     332                                        svec[offsets[subconfigurationIndex]+j].subWithOverflowChecking(denominator);
     333                        svecBound.subWithOverflowChecking(denominator);
     334
     335
     336                        for(int a=0;a<k+Flags::computeDotProductInMatrix;a++)
     337                                if(first||a!=subconfigurationIndex)
     338                                Abounds.UNCHECKEDACCESS(a)=mvtyp::dotDivVector(&A.UNCHECKEDACCESS(a,0),&(svec[0]),t,tempA.UNCHECKEDACCESS(a),denominatorDivisor,m,Abounds[a],svecBound);
     339        //                        for(int u=0;u<m;u++)
     340        //                                *A.UNCHECKEDACCESSRESTRICT(a,u)=dotDiv(svec[u],tempA.UNCHECKEDACCESS(a),t,A.UNCHECKEDACCESS(a,u),denominatorDivisor);
     341
     342
     343                        {
     344                                for(int a=0;a<k+Flags::computeDotProductInMatrix;a++)
     345                                        {
     346                                                A[a][offsets[subconfigurationIndex]+oldIndex]=tempA[a].negatedWithOverflowChecking();
     347                                                Abounds[a]=min(Abounds[a],negabs(tempA[a]));
     348                                        }
     349
     350                                if(!first)
     351                                {
     352                                        A[subconfigurationIndex][offsets[subconfigurationIndex]+oldIndex]=denominator.negatedWithOverflowChecking();
     353                                        Abounds[subconfigurationIndex].subWithOverflowChecking(denominator);
     354                                }
     355                        }
     356                        denominator=nextDenominator;
     357
     358                        // We clear these unused entries of A to ensure that these columns are not chosen when comparing inequalities
     359                        for(int i=0;i<k+Flags::computeDotProductInMatrix;i++)
     360                        {
     361                                A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+choices[subconfigurationIndex].first)=0;
     362                                A.UNCHECKEDACCESS(i,offsets[subconfigurationIndex]+choices[subconfigurationIndex].second)=0;
     363                        }
     364                }
     365                void replaceFirst(int subconfigurationIndex, int newIndex, Vector<mvtyp> const &target){replaceFirstOrSecond(true,subconfigurationIndex,newIndex,target);}
     366                void replaceSecond(int subconfigurationIndex, int newIndex, Vector<mvtyp> const &target){replaceFirstOrSecond(false,subconfigurationIndex,newIndex,target);}
     367
     368                InequalityTable(std::vector<Matrix<mvtyp> > const &tuple_, int subconfigurationIndex_):
     369                        tempA(tuple_.size()+Flags::computeDotProductInMatrix),
     370                        tuple(tuple_),
     371                        choices(tuple_.size()),
     372                        subconfigurationIndex(subconfigurationIndex_),
     373                        offsets(tuple_.size())
     374                {// Cannot overflow
     375                        k=tuple.size();
     376                        m=0;
     377                        for(int i=0;i<tuple.size();i++)m+=tuple[i].getWidth();
     378                        svec.resize(m);
     379                        A=Matrix<mvtyp>(k+Flags::computeDotProductInMatrix,m);
     380                        {int offset=0;for(int i=0;i<tuple.size();i++){offsets[i]=offset;offset+=tuple[i].getWidth();}}
     381                        Abounds=Vector<mvtyp>(k+Flags::computeDotProductInMatrix);
     382                }
     383                void setChoicesInitially()
     384                {// Cannot overflow
     385                        //THIS WILL ONLY WORK FOR THE STARTING CONFIGURATION
     386                        //sets denominators,A and choices (these must have been initialized to the right sizes)
     387                        for(int i=0;i<k;i++)
     388                                choices[i]=std::pair<int, int> (i+0,i+1);
     389                        for(int i=0;i<m;i++)
     390                                for(int j=0;j<k;j++)
     391                                        A[j][i]=0;
     392                        //we follow the lemma in the article. Why does the proof of the lemma split into 3 cases and not just 2?
     393                        int a=0;
     394                        for(int i=0;i<k;i++)
     395                                for(int gamma=0;gamma<tuple[i].getWidth();gamma++,a++)
     396                                        if(gamma>i+1)
     397                                                for(int ii=i;ii<gamma;ii++)
     398                                                        A[ii][a]=-1;
     399                                        else
     400                                                for(int ii=gamma;ii<i;ii++)
     401                                                        A[ii][a]=1;
     402                        denominator=1;
     403                        for(int i=0;i<k;i++)Abounds[i]=-1;
     404        //                assignDotProducts(target);
     405                }
     406                void compareInequalities(InequalityComparisonResult &result, Vector<mvtyp> const &target, int onlyK=-1)
     407                {// Can only overflow if target entries are too big. Actually it seems that it is not this function that will overflow but dotVector.
     408                        bool empty=true;
     409                        int bestConfigurationIndex=-1;
     410                        int bestColumnIndex=-1;
     411                        mvtyp targetDotBest=0;
     412
     413                        for(int i=0;i<k;i++)
     414                        {
     415                                // TO DO: exchange the line with the template parameter version:
     416//                    gfan::Matrix<mvtyp>::const_RowRef Ak=const_cast<const gfan::Matrix<mvtyp>&>(A)[k];
     417                    gfan::Matrix<CircuitTableInt32>::const_RowRef Ak=const_cast<const gfan::Matrix<mvtyp>&>(A)[k];
     418                    int offsetsi=offsets[i];
     419                    int tupleiwidth=tuple[i].getWidth();
     420
     421                                if(onlyK!=-1)if(i!=onlyK)continue;
     422                                for(int j=0;j<tupleiwidth;j++)
     423                                        if(Flags::computeDotProductInMatrix || isLegalIndex(i,j))//unused inequalities will have value 0. Therefore isLegalIndex(i,j) is not required if values are stored.
     424                                        {
     425                                                mvtyp ineqDotTarget=Flags::computeDotProductInMatrix?Ak.UNCHECKEDACCESS(offsetsi+j):dotVector(i,j,target,onlyK);
     426                                                if(ineqDotTarget.isNegative())
     427                                                {
     428                                                        if(!isReverseLexInvertedLessThanZero(i,j))
     429                                                        {
     430                                                                if(empty||compareReverseLexicographicInverted(bestConfigurationIndex,bestColumnIndex,i,j,ineqDotTarget,targetDotBest))
     431                                                                {
     432                                                                        targetDotBest=ineqDotTarget;
     433                                                                        empty=false;
     434                                                                        bestConfigurationIndex=i;
     435                                                                        bestColumnIndex=j;
     436                                                                }
     437                                                        }
     438                //                                        assert(!ineq.isZero());
     439                                                }
     440                                        }
     441                        }
     442                        result.empty=empty;
     443                        result.configurationIndex=bestConfigurationIndex;
     444                        result.columnIndex=bestColumnIndex;
     445        //                assert(denominator>0);
     446                }
     447                void setChoicesFromEarlierHomotopy(InequalityTable const &parent, mvtyp degreeScaling, Vector<mvtyp> const &target)
     448                {        // Overflows may happen, but only if the tuple entries are too big.
     449                        // Notice that the code below has overflow checks everywhere, except in the loop marked with a comment.
     450
     451                        //sets denominators,A and choices (these must have been initialized to the right sizes
     452                        //columns have been added to configuration this->subconfigurationIndex
     453                        //for that reason we need to introduce new circuits. Old circuits are preserved.
     454                        //chioices are "relative" so no update is needed.
     455
     456                        choices=parent.choices;
     457                        int numberToDrop=(subconfigurationIndex!=0) ? numberToDrop=k+1 : 0;
     458
     459                        choices[subconfigurationIndex-1].first-=numberToDrop;
     460                        choices[subconfigurationIndex-1].second-=numberToDrop;
     461
     462                        denominator=parent.denominator;
     463                        int offsetOld=0;
     464                        int offsetNew=0;
     465                        for(int i=0;i<k;i++)
     466                        {
     467                                int localNumberToDrop=0;
     468                                if(i==subconfigurationIndex-1)
     469                                        localNumberToDrop=numberToDrop;
     470                                for(int a=0;a<A.getHeight()-Flags::computeDotProductInMatrix;a++)
     471                                        if(a==subconfigurationIndex)
     472                                                for(int j=0;j<parent.tuple[i].getWidth()-localNumberToDrop;j++)
     473                                                        A.UNCHECKEDACCESS(a,offsetNew+j)=parent.A.UNCHECKEDACCESS(a,offsetOld+j+localNumberToDrop);
     474                                        else
     475                                        {
     476                                                // We can get an estimated 3 percent speed up by using known bounds to check if any of the multiplications below can overflow.
     477                                                // Moreover, these multiplications can be moved outside the outer loop, to give better vectorisation.
     478                                                for(int j=0;j<parent.tuple[i].getWidth()-localNumberToDrop;j++)
     479                                                        A.UNCHECKEDACCESS(a,offsetNew+j)=extendedMultiplication(degreeScaling,parent.A.UNCHECKEDACCESS(a,offsetOld+j+localNumberToDrop)).castToSingle();
     480                                        }
     481                                if(i==subconfigurationIndex)
     482                                {
     483                                        for(int j=parent.tuple[i].getWidth();j<tuple[i].getWidth();j++)
     484                                        {
     485                                                for(int a=0;a<A.getHeight()-Flags::computeDotProductInMatrix;a++)
     486                                                        A.UNCHECKEDACCESS(a,offsetNew+j)=0;
     487                                                {
     488                                                        int b=choices[subconfigurationIndex].second-1;
     489                                                        if(b>=0)
     490                                                                A.UNCHECKEDACCESS(subconfigurationIndex,offsetNew+j).msubWithOverflowChecking(tuple[i].UNCHECKEDACCESS(b,j),denominator);
     491                                                }
     492                                                for(int a=0;a<k+Flags::computeDotProductInMatrix;a++)
     493                                                {   // If tuple entries are not bounded, overflow can happen in the following loop
     494                                                        mvtypDouble tempDouble=A.UNCHECKEDACCESS(a,offsetNew+j).extend();
     495                                                        for(int b=0;b<k;b++)
     496                                                                if(choices[subconfigurationIndex].second!=b+1 &&choices[subconfigurationIndex].first!=b+1)
     497                                                                        tempDouble+=extendedMultiplication(tuple[i].UNCHECKEDACCESS(b,j),A.UNCHECKEDACCESS(a,offsetNew+b+1));
     498                                                        A.UNCHECKEDACCESS(a,offsetNew+j)=tempDouble.castToSingle();
     499                                                }
     500                                                mvtypDouble left=degreeScaling.extend();
     501                                                for(int b=0;b<k;b++)
     502                                                        left-=mvtyp(tuple[i].UNCHECKEDACCESS(b,j)).extend();
     503
     504                                                mvtyp leftS=left.castToSingle();
     505
     506                                                if(choices[subconfigurationIndex].second==0)
     507                                                        A.UNCHECKEDACCESS(subconfigurationIndex,offsetNew+j).msubWithOverflowChecking(leftS,denominator);
     508                                                else if(choices[subconfigurationIndex].first!=0)
     509                                                        for(int a=0;a<k+Flags::computeDotProductInMatrix;a++)
     510                                                                A.UNCHECKEDACCESS(a,offsetNew+j).maddWithOverflowChecking(leftS,mvtyp(A.UNCHECKEDACCESS(a,offsetNew)));
     511                                        }
     512                                        for(int j=0;j<parent.tuple[i].getWidth();j++)// <-- this loop has side effects on addExpressionOfCeb() above. Therefore it is placed after the loop above.
     513                                                for(int a=0;a<A.getHeight()-Flags::computeDotProductInMatrix;a++)
     514                                                        A.UNCHECKEDACCESS(a,offsetNew+j)=extendedMultiplication(A.UNCHECKEDACCESS(a,offsetNew+j),degreeScaling).castToSingle();
     515                                }
     516                                offsetOld+=parent.tuple[i].getWidth();
     517                                offsetNew+=tuple[i].getWidth();
     518                        }
     519                        denominator=extendedMultiplication(denominator,degreeScaling).castToSingle();
     520                        if(Flags::computeDotProductInMatrix)assignDotProducts(target,subconfigurationIndex);
     521                        computeABounds();
     522                        for(int a=0;a<k;a++)
     523                        for(int i=0;i<k+Flags::computeDotProductInMatrix;i++)
     524                        {
     525                                A[i][offsets[a]+choices[a].first]=0;
     526                                A[i][offsets[a]+choices[a].second]=0;
     527                        }
     528                }
     529        };
     530                struct StackItem{
     531                        int columnIndex;
     532                        int configurationIndex;
     533                        bool b;
     534                        int choice;
     535                        bool useFirstChanged,useSecondChanged;
     536                        StackItem(int columnIndex_, int configurationIndex_, bool b_, int choice_, bool useFirstChanged_, bool useSecondChanged_):
     537                                columnIndex(columnIndex_),
     538                                configurationIndex(configurationIndex_),
     539                                b(b_),
     540                                choice(choice_),
     541                                useFirstChanged(useFirstChanged_),
     542                                useSecondChanged(useSecondChanged_)
     543                        {
     544                        }
     545                };
     546        public:
     547                std::vector<std::pair<int,int> > choices;
     548                Vector<mvtyp> target;
     549                bool useFirstChanged;
     550                bool useSecondChanged;
     551                std::vector<StackItem> stack;
     552                int eliminatedK;
     553                int eliminatedKOffset;
     554                std::vector<Matrix<mvtyp> > tuple;
     555                std::vector<int> offsets;
     556                int m;
     557                InequalityComparisonResult result;
     558                InequalityTable inequalityTable;
     559                void constructInequalityTableFromParent(InequalityTable const &parentTable, mvtyp degreeScaling)
     560                {
     561                        inequalityTable.setChoicesFromEarlierHomotopy(parentTable, degreeScaling, target);
     562                }
     563                void constructInequalityTableInitially(mvtyp degreeScaling)
     564                {
     565                        std::vector<Matrix<mvtyp> > tempTuple;for(int i=0;i<tuple.size();i++)tempTuple.push_back(simplex<mvtyp>(tuple.size(),1));
     566                        InequalityTable tempTable(tempTuple,-1);
     567                        tempTable.setChoicesInitially();
     568                        constructInequalityTableFromParent(tempTable,degreeScaling);
     569                }
     570                SingleTropicalHomotopyTraverser(std::vector<Matrix<mvtyp> > const &tuple_, int m_, std::vector<std::pair<int,int> > const &choices_, Vector<mvtyp> const &target_, int eliminatedK_):
     571                        choices(choices_),
     572                        target(target_),
     573                        eliminatedK(eliminatedK_),
     574                        tuple(tuple_),
     575                        m(m_),
     576                        inequalityTable(tuple,eliminatedK),
     577                        offsets(tuple_.size())
     578                {
     579                        eliminatedKOffset=0;for(int i=0;i<eliminatedK;i++)eliminatedKOffset+=tuple_[i].getWidth();
     580                        {int offset=0;for(int i=0;i<tuple.size();i++){offsets[i]=offset;offset+=tuple[i].getWidth();}}
     581
     582                        /* We also need to check that the input is within the required bounds
     583                        * This is important to not have overflows in:
     584                        *   dotVector
     585                        *   setChoicesFromEarlierHomotopy
     586                        * The number of summands is bounded by 3k+4*/
     587
     588                        bool isOK=mvtyp::isSuitableSummationNumber(3*tuple.size()+4);
     589                        for(int i=0;i<target.size();i++)isOK&=mvtyp(target[i]).fitsInHalf();
     590                        for(int i=0;i<tuple.size();i++)
     591                                for(int j=0;j<tuple[i].getHeight();j++)
     592                                        for(int k=0;k<tuple[i].getWidth();k++)
     593                                                isOK&=mvtyp(tuple[i][j][k]).fitsInHalf();
     594                        if(!isOK)throw MVMachineIntegerOverflow;
     595                }
     596                virtual void process()
     597                {
     598                }
     599                bool findOutgoingAndProcess(bool doProcess)//sets up useFirstChanged and useSecondChanged and processes if process is true
     600                {//returns true if we are at a leaf
     601//                        inequalityTable.checkABounds();
     602                        useFirstChanged=false;
     603                        useSecondChanged=false;
     604                        int onlyK=-1;
    605605#if 1
    606                         if(eliminatedK!=-1)
    607                                 if(target[choices[eliminatedK].first+eliminatedKOffset]==target[choices[eliminatedK].second+eliminatedKOffset])
    608                                         onlyK=eliminatedK;
     606                        if(eliminatedK!=-1)
     607                                if(target[choices[eliminatedK].first+eliminatedKOffset]==target[choices[eliminatedK].second+eliminatedKOffset])
     608                                        onlyK=eliminatedK;
    609609#endif
    610610
    611                         inequalityTable.compareInequalities(result,target,onlyK);
    612 
    613                         if(result.empty)
    614                         {
    615                                 if(doProcess)process();
    616                                 return true;
    617                         }
    618 
    619                         // reverse search rule:
    620 
    621                         mvtypDouble bestAtFirst=inequalityTable.getCoordinateOfInequality(result.configurationIndex,result.columnIndex,result.configurationIndex,choices[result.configurationIndex].first);
    622                         mvtypDouble bestAtSecond=inequalityTable.getCoordinateOfInequality(result.configurationIndex,result.columnIndex,result.configurationIndex,choices[result.configurationIndex].second);
    623                         if(bestAtFirst.isNegative())
    624                         {
    625                                 if(bestAtSecond.isNegative())
    626                                 {
    627                                         useFirstChanged=true;
    628                                         useSecondChanged=true;
    629                                 }
    630                                 else
    631                                 {
    632                                         if(bestAtSecond.isZero())
    633                                         useFirstChanged=true;
    634                                         else
    635                                         if(choices[result.configurationIndex].second<result.columnIndex)
    636                                         useFirstChanged=true;
    637                                 }
    638                         }
    639                         else
    640                         {
    641                                 if(bestAtSecond.isNegative())
    642                                 {
    643                                         if(bestAtFirst.isZero())
    644                                         useSecondChanged=true;
    645                                         else
    646                                         if(choices[result.configurationIndex].first<result.columnIndex)
    647                                         useSecondChanged=true;
    648                                 }
    649                                 else
    650                                 {
    651                                         assert(0);
    652                                 }
    653                         }
    654                         return false;
    655                 }
    656                 void goToFirstChild()
    657                 {
    658 //                      debug<<"First:  configuration index:"<<data.configurationIndex<<" column index:"<<data.columnIndex<<"\n";
    659                         assert(useFirstChanged);
    660                         {
    661                                 stack.push_back(StackItem(
    662                                                 result.columnIndex,
    663                                                 result.configurationIndex,
    664                                                 false,
    665                                                 choices[result.configurationIndex].first,
    666                                                 useFirstChanged,
    667                                                 useSecondChanged));
    668                                 choices[result.configurationIndex].first=result.columnIndex;
    669                                 inequalityTable.replaceFirst(result.configurationIndex,result.columnIndex,target);
    670                         }
    671                 }
    672                 void goToSecondChild()
    673                 {
    674 //                      debug<<"Second: configuration index:"<<data.configurationIndex<<" column index:"<<data.columnIndex<<"\n";
    675                         assert(useSecondChanged);
    676                         {
    677                                 stack.emplace_back(StackItem(
    678                                                 result.columnIndex,
    679                                                 result.configurationIndex,
    680                                                 true,
    681                                                 choices[result.configurationIndex].second,
    682                                                 useFirstChanged,
    683                                                 useSecondChanged));
    684                                 choices[result.configurationIndex].second=result.columnIndex;
    685                                 inequalityTable.replaceSecond(result.configurationIndex,result.columnIndex,target);
    686                         }
    687                 }
    688                 int numberOfChildren()
    689                 {
    690                         return int(useFirstChanged)+int(useSecondChanged);
    691                 }
    692                 void goToNthChild(int n)
    693                 {
    694                         if(n==0)
    695                                 if(useFirstChanged)
    696                                         goToFirstChild();
    697                                 else
    698                                         goToSecondChild();
    699                         else
    700                                 goToSecondChild();
    701                 }
    702                 void goBack()
    703                 {
    704                         StackItem &B=stack.back();
    705                         result.columnIndex=B.columnIndex;
    706                         result.configurationIndex=B.configurationIndex;
    707                         if(B.b)
    708                         {
    709                                 choices[result.configurationIndex].second=B.choice;
    710                                 inequalityTable.replaceSecond(result.configurationIndex,B.choice,target);
    711                         }
    712                         else
    713                         {
    714                                 choices[result.configurationIndex].first=B.choice;
    715                                 inequalityTable.replaceFirst(result.configurationIndex,B.choice,target);
    716                         }
    717                         useFirstChanged=B.useFirstChanged;
    718                         useSecondChanged=B.useSecondChanged;
    719                         stack.pop_back();
    720                 }
    721                 bool atRoot()
    722                 {
    723                         return stack.empty();
    724                 }
    725         };
     611                        inequalityTable.compareInequalities(result,target,onlyK);
     612
     613                        if(result.empty)
     614                        {
     615                                if(doProcess)process();
     616                                return true;
     617                        }
     618
     619                        // reverse search rule:
     620
     621                        mvtypDouble bestAtFirst=inequalityTable.getCoordinateOfInequality(result.configurationIndex,result.columnIndex,result.configurationIndex,choices[result.configurationIndex].first);
     622                        mvtypDouble bestAtSecond=inequalityTable.getCoordinateOfInequality(result.configurationIndex,result.columnIndex,result.configurationIndex,choices[result.configurationIndex].second);
     623                        if(bestAtFirst.isNegative())
     624                        {
     625                                if(bestAtSecond.isNegative())
     626                                {
     627                                        useFirstChanged=true;
     628                                        useSecondChanged=true;
     629                                }
     630                                else
     631                                {
     632                                        if(bestAtSecond.isZero())
     633                                        useFirstChanged=true;
     634                                        else
     635                                        if(choices[result.configurationIndex].second<result.columnIndex)
     636                                        useFirstChanged=true;
     637                                }
     638                        }
     639                        else
     640                        {
     641                                if(bestAtSecond.isNegative())
     642                                {
     643                                        if(bestAtFirst.isZero())
     644                                        useSecondChanged=true;
     645                                        else
     646                                        if(choices[result.configurationIndex].first<result.columnIndex)
     647                                        useSecondChanged=true;
     648                                }
     649                                else
     650                                {
     651                                        assert(0);
     652                                }
     653                        }
     654                        return false;
     655                }
     656                void goToFirstChild()
     657                {
     658//                        debug<<"First:  configuration index:"<<data.configurationIndex<<" column index:"<<data.columnIndex<<"\n";
     659                        assert(useFirstChanged);
     660                        {
     661                                stack.push_back(StackItem(
     662                                                result.columnIndex,
     663                                                result.configurationIndex,
     664                                                false,
     665                                                choices[result.configurationIndex].first,
     666                                                useFirstChanged,
     667                                                useSecondChanged));
     668                                choices[result.configurationIndex].first=result.columnIndex;
     669                                inequalityTable.replaceFirst(result.configurationIndex,result.columnIndex,target);
     670                        }
     671                }
     672                void goToSecondChild()
     673                {
     674//                        debug<<"Second: configuration index:"<<data.configurationIndex<<" column index:"<<data.columnIndex<<"\n";
     675                        assert(useSecondChanged);
     676                        {
     677                                stack.emplace_back(StackItem(
     678                                                result.columnIndex,
     679                                                result.configurationIndex,
     680                                                true,
     681                                                choices[result.configurationIndex].second,
     682                                                useFirstChanged,
     683                                                useSecondChanged));
     684                                choices[result.configurationIndex].second=result.columnIndex;
     685                                inequalityTable.replaceSecond(result.configurationIndex,result.columnIndex,target);
     686                        }
     687                }
     688                int numberOfChildren()
     689                {
     690                        return int(useFirstChanged)+int(useSecondChanged);
     691                }
     692                void goToNthChild(int n)
     693                {
     694                        if(n==0)
     695                                if(useFirstChanged)
     696                                        goToFirstChild();
     697                                else
     698                                        goToSecondChild();
     699                        else
     700                                goToSecondChild();
     701                }
     702                void goBack()
     703                {
     704                        StackItem &B=stack.back();
     705                        result.columnIndex=B.columnIndex;
     706                        result.configurationIndex=B.configurationIndex;
     707                        if(B.b)
     708                        {
     709                                choices[result.configurationIndex].second=B.choice;
     710                                inequalityTable.replaceSecond(result.configurationIndex,B.choice,target);
     711                        }
     712                        else
     713                        {
     714                                choices[result.configurationIndex].first=B.choice;
     715                                inequalityTable.replaceFirst(result.configurationIndex,B.choice,target);
     716                        }
     717                        useFirstChanged=B.useFirstChanged;
     718                        useSecondChanged=B.useSecondChanged;
     719                        stack.pop_back();
     720                }
     721                bool atRoot()
     722                {
     723                        return stack.empty();
     724                }
     725        };
    726726
    727727
    … …  
    730730 */
    731731template<class mvtyp, class mvtypDouble, class mvtypDivisor>
    732         class TropicalRegenerationTraverser{
    733                 // The following class is an attempt to separate homotopy data from the traversal logic.
    734                 class Data{
    735                         static mvtyp degree(Matrix<mvtyp> const &m)//assumes entries of m are non-negative
    736                         {
    737                                   mvtyp ret=0;
    738                                   for(int i=0;i<m.getWidth();i++)
    739                                   {
    740                                           mvtyp s(0);
    741                                           for(int j=0;j<m.getHeight();j++)
    742                                                   s.addWithOverflowChecking(m[j][i]);
    743                                           ret=max(s,ret);
    744                                   }
    745                                   return ret;
    746                         }
    747                   public:
    748                           std::vector<Vector<mvtyp> > targets;
    749                           std::vector<Matrix<mvtyp> > tuple;
    750                           std::vector<std::vector<Matrix<mvtyp> > > tuples;
    751                           Vector<mvtyp> degrees;
    752                           bool isFiniteIndex(int level, int index)
    753                           {
    754                                   return index>=tuple[0].getHeight()+1;
    755                           }
    756 
    757                           std::vector<Matrix<mvtyp> > produceIthTuple(int i)
    758                                 {
    759                                   int n=tuple[0].getHeight();
    760                                   std::vector<Matrix<mvtyp> > ret;
    761                                   for(int j=0;j<tuple.size();j++)
    762                                   {
    763                                           if(j<i)ret.push_back(tuple[j]);
    764                                           if(j==i)ret.push_back(combineLeftRight(simplex<mvtyp>(n,degree(tuple[j])),tuple[j]));
    765                                           if(j>i)ret.push_back(simplex<mvtyp>(n,1));
    766                                   }
    767                                   return ret;
    768                                 }
    769                           Data(std::vector<Matrix<mvtyp> > const &tuple_):tuple(tuple_)
    770                           {
    771                                   int n=tuple[0].getHeight();
    772 
    773                                   {//adjusting to positive orthant
    774                                           for(int i=0;i<tuple.size();i++)
    775                                                   for(int j=0;j<tuple[i].getHeight();j++)
    776                                                   {
    777                                                           mvtyp m;
    778                                                           if(tuple[i].getWidth()==0)
    779                                                                   m=0;
    780                                                           else
    781                                                                   m=tuple[i][j][0];
    782                                                           for(int k=0;k<tuple[i].getWidth();k++)m=min(m,tuple[i][j][k]);
    783                                                           for(int k=0;k<tuple[i].getWidth();k++)tuple[i][j][k].subWithOverflowChecking(m);
    784                                                   }
    785                                   }
    786 
    787                                   for(int i=0;i<tuple.size();i++)
    788                                           degrees.push_back(degree(tuple[i]));
    789 
    790                                   for(int i=0;i<tuple.size();i++)
    791                                           tuples.push_back(produceIthTuple(i));
    792 
    793                                   for(int i=0;i<tuple.size();i++)
    794                                   {
    795                                           Vector<mvtyp> targ;
    796                                           for(int j=0;j<tuple.size();j++)
    797                                           {
    798                                                   if(j==i)
    799                                                           targ=concatenation(targ,concatenation(Vector<mvtyp>::allOnes(n+1),Vector<mvtyp>(tuple[i].getWidth())));
    800                                                   else
    801                                                           targ=concatenation(targ,Vector<mvtyp>(tuples[i][j].getWidth()));
    802                                           }
    803                                           targets.push_back(targ);
    804                                   }
    805                           };
    806 
    807                           std::vector<std::pair<int,int> > firstIntersection()
    808                           {
    809                                   std::vector<std::pair<int,int> > ret;
    810                                   for(int i=0;i<tuple.size();i++)
    811                                           ret.push_back(std::pair<int,int>(i+0,i+1));
    812                                   return ret;
    813                           }
    814 
    815                           void castToNextLevel(std::vector<std::pair<int,int> > const &choices, int i, int S, std::vector<std::pair<int,int> > &ret)
    816                           {
    817                                   assert(ret.size()==choices.size());
    818                                   for(int j=0;j<choices.size();j++)
    819                                           ret[j]=choices[j];
    820 
    821                                   assert(ret[i].first>=S);
    822                                   assert(ret[i].second>=S);
    823                                   ret[i].first-=S;
    824                                   ret[i].second-=S;
    825                           }
    826                   };
    827                 static int cayleyConfigurationWidth(std::vector<Matrix<mvtyp> > const &tuple)
    828                 {
    829                           int m2=0;
    830                           for(int i=0;i<tuple.size();i++)
    831                                   m2+=tuple[i].getWidth();
    832                           return m2;
    833                 }
    834         public:
    835                 int depth;
    836                 int counter;
    837                 typedef SingleTropicalHomotopyTraverser<mvtyp,mvtypDouble,mvtypDivisor> MySingleTropicalHomotopyTraverser;
    838                 std::vector<MySingleTropicalHomotopyTraverser> traversers;
    839                 Data fullData;
    840                 int level;
    841                 bool deadEnd;
    842                 bool isLevelLeaf;
    843                 bool isSolutionVertex;
    844                 std::vector<bool> isLevelLeafStack;
    845                 TropicalRegenerationTraverser(std::vector<Matrix<mvtyp> > const &tuple_):
    846                         fullData(tuple_),counter(0),depth(0)
    847                 {
    848                         assert(tuple_.size());
    849                         for(int i=0;i<tuple_.size();i++)
    850                                 traversers.push_back(MySingleTropicalHomotopyTraverser(fullData.tuples[i],cayleyConfigurationWidth(fullData.tuples[i]),fullData.firstIntersection(),fullData.targets[i],i));
    851                         traversers[0].constructInequalityTableInitially(fullData.degrees[0]);
    852                         level=0;
    853                 }
    854                 virtual void process()
    855                 {
    856                 }
    857                 bool findOutgoingAndProcess(bool doProcess)
    858                 {
    859                         isSolutionVertex=false;
    860                         deadEnd=false;
    861                         isLevelLeaf=traversers[level].findOutgoingAndProcess(false);
    862                         if(isLevelLeaf)
    863                         {//leaf
    864                                 bool isFinite=fullData.isFiniteIndex(level,traversers[level].choices[level].first)&&fullData.isFiniteIndex(level,traversers[level].choices[level].second);
    865                                 deadEnd=!isFinite;
    866                                 if(isFinite && (level==fullData.tuple.size()-1))
    867                                 {
    868                                         isSolutionVertex=true;
    869                                         if(doProcess){process();}
    870                                         return true;
    871                                 }
    872                         }
    873                         return false;
    874                 }
    875                 int numberOfChildren()
    876                 {
    877                         if(isLevelLeaf&&(level==fullData.tuple.size()-1))return 0;
    878                         if(!isLevelLeaf)
    879                                 return traversers[level].numberOfChildren();
    880                         else
    881                                 return 1-deadEnd;
    882                 }
    883                 void goToNthChild(int n)
    884                 {
    885                         depth++;
    886                         isLevelLeafStack.push_back(isLevelLeaf);
    887                         if(!isLevelLeaf)
    888                                 traversers[level].goToNthChild(n);
    889                         else
    890                         {
    891                                 fullData.castToNextLevel(traversers[level].choices,level,fullData.tuples[level][level].getWidth()-fullData.tuples[level+1][level].getWidth(),traversers[level+1].choices);
    892                                 traversers[level+1].constructInequalityTableFromParent(traversers[level].inequalityTable,fullData.degrees[level+1]);
    893                                 level++;
    894                         }
    895                 }
    896                 void print()
    897                 {
    898                 }
    899                 void goBack()
    900                 {
    901                         depth--;
    902                         counter++;
    903                         deadEnd=false;
    904                         if(traversers[level].atRoot())
    905                                 level--;
    906                         else
    907                                 traversers[level].goBack();
    908                         isLevelLeaf=isLevelLeafStack.back();
    909                         isLevelLeafStack.pop_back();
    910                 }
    911         };
    912 
    913         /*
    914         * This class glues Bjarne Knudsen's parallel traverser to our MultiLevelIntersectionTraverser.
    915         * This class should eventually be written so that it works on any homotopy.
    916         *
    917         * Since we are not sure about how exceptions would be handled by the parallel traverser,
    918         * we write ugly aborting code, which may actually work.
    919         */
    920         template<class mvtyp, class mvtypDouble, class mvtypDivisor>
    921         class SpecializedRTraverser: public Traverser
    922         {
    923         public:
    924                 typedef TropicalRegenerationTraverser<mvtyp,mvtypDouble,mvtypDivisor> MyTropicalRegenerationTraverser;
    925                 MyTropicalRegenerationTraverser T;
    926                 mvtypDouble mixedVolume;
    927                 int numberOfExpensiveSteps;
    928                 SpecializedRTraverser(std::vector<Matrix<mvtyp> > const &tuple_):
    929                         mixedVolume(),
    930                         numberOfExpensiveSteps(0),
    931                         T(tuple_) //Constructor my throw if input is too big.
    932                 {
    933                         numberOfExpensiveSteps++;
    934                         T.findOutgoingAndProcess(false);
    935                 }
    936                 int  getEdgeCountNext( void )
    937                 {
    938                         if(!aborting)
    939                         {
    940                                 try{
    941                                         return T.numberOfChildren();
    942                                 }
    943                                 catch(...){abort();}
    944                         }
    945                         return 0;
    946                 }
    947 
    948                 int  moveToNext( int   index,
    949                                    bool  collect_info )
    950                 {
    951                         if(!aborting)
    952                         {
    953                                 try{
    954                                         T.goToNthChild(index);
    955                                         numberOfExpensiveSteps++;
    956                                         T.findOutgoingAndProcess(false);
    957                                 }
    958                                 catch(...){abort();}
    959                         }
    960                         return 0;
    961                 }
    962 
    963           void  moveToPrev( int  index )
    964           {
    965                   if(!aborting)
    966                   {
    967                           try{
    968                                   T.goBack(); //index ignored
    969                           }
    970                           catch(...){abort();}
    971                   }
    972           }
    973 
    974           void  collectInfo( void )
    975           {
    976                   if(!aborting)
    977                   {
    978                           try{
    979                                   if(T.isSolutionVertex)
    980                                           mixedVolume.addWithOverflowCheck(T.traversers[T.level].inequalityTable.getVolume().extend());
    981                           }
    982                           catch(...){abort();}
    983                   }
    984           }
    985 
    986           void  printState( void )
    987           {
    988                   T.print();
    989           }
    990         };
     732        class TropicalRegenerationTraverser{
     733                // The following class is an attempt to separate homotopy data from the traversal logic.
     734                class Data{
     735                        static mvtyp degree(Matrix<mvtyp> const &m)//assumes entries of m are non-negative
     736                        {
     737                                  mvtyp ret=0;
     738                                  for(int i=0;i<m.getWidth();i++)
     739                                  {
     740                                          mvtyp s(0);
     741                                          for(int j=0;j<m.getHeight();j++)
     742                                                  s.addWithOverflowChecking(m[j][i]);
     743                                          ret=max(s,ret);
     744                                  }
     745                                  return ret;
     746                        }
     747                  public:
     748                          std::vector<Vector<mvtyp> > targets;
     749                          std::vector<Matrix<mvtyp> > tuple;
     750                          std::vector<std::vector<Matrix<mvtyp> > > tuples;
     751                          Vector<mvtyp> degrees;
     752                          bool isFiniteIndex(int level, int index)
     753                          {
     754                                  return index>=tuple[0].getHeight()+1;
     755                          }
     756
     757                          std::vector<Matrix<mvtyp> > produceIthTuple(int i)
     758                                {
     759                                  int n=tuple[0].getHeight();
     760                                  std::vector<Matrix<mvtyp> > ret;
     761                                  for(int j=0;j<tuple.size();j++)
     762                                  {
     763                                          if(j<i)ret.push_back(tuple[j]);
     764                                          if(j==i)ret.push_back(combineLeftRight(simplex<mvtyp>(n,degree(tuple[j])),tuple[j]));
     765                                          if(j>i)ret.push_back(simplex<mvtyp>(n,1));
     766                                  }
     767                                  return ret;
     768                                }
     769                          Data(std::vector<Matrix<mvtyp> > const &tuple_):tuple(tuple_)
     770                          {
     771                                  int n=tuple[0].getHeight();
     772
     773                                  {//adjusting to positive orthant
     774                                          for(int i=0;i<tuple.size();i++)
     775                                                  for(int j=0;j<tuple[i].getHeight();j++)
     776                                                  {
     777                                                          mvtyp m;
     778                                                          if(tuple[i].getWidth()==0)
     779                                                                  m=0;
     780                                                          else
     781                                                                  m=tuple[i][j][0];
     782                                                          for(int k=0;k<tuple[i].getWidth();k++)m=min(m,tuple[i][j][k]);
     783                                                          for(int k=0;k<tuple[i].getWidth();k++)tuple[i][j][k].subWithOverflowChecking(m);
     784                                                  }
     785                                  }
     786
     787                                  for(int i=0;i<tuple.size();i++)
     788                                          degrees.push_back(degree(tuple[i]));
     789
     790                                  for(int i=0;i<tuple.size();i++)
     791                                          tuples.push_back(produceIthTuple(i));
     792
     793                                  for(int i=0;i<tuple.size();i++)
     794                                  {
     795                                          Vector<mvtyp> targ;
     796                                          for(int j=0;j<tuple.size();j++)
     797                                          {
     798                                                  if(j==i)
     799                                                          targ=concatenation(targ,concatenation(Vector<mvtyp>::allOnes(n+1),Vector<mvtyp>(tuple[i].getWidth())));
     800                                                  else
     801                                                          targ=concatenation(targ,Vector<mvtyp>(tuples[i][j].getWidth()));
     802                                          }
     803                                          targets.push_back(targ);
     804                                  }
     805                          };
     806
     807                          std::vector<std::pair<int,int> > firstIntersection()
     808                          {
     809                                  std::vector<std::pair<int,int> > ret;
     810                                  for(int i=0;i<tuple.size();i++)
     811                                          ret.push_back(std::pair<int,int>(i+0,i+1));
     812                                  return ret;
     813                          }
     814
     815                          void castToNextLevel(std::vector<std::pair<int,int> > const &choices, int i, int S, std::vector<std::pair<int,int> > &ret)
     816                          {
     817                                  assert(ret.size()==choices.size());
     818                                  for(int j=0;j<choices.size();j++)
     819                                          ret[j]=choices[j];
     820
     821                                  assert(ret[i].first>=S);
     822                                  assert(ret[i].second>=S);
     823                                  ret[i].first-=S;
     824                                  ret[i].second-=S;
     825                          }
     826                  };
     827                static int cayleyConfigurationWidth(std::vector<Matrix<mvtyp> > const &tuple)
     828                {
     829                          int m2=0;
     830                          for(int i=0;i<tuple.size();i++)
     831                                  m2+=tuple[i].getWidth();
     832                          return m2;
     833                }
     834        public:
     835                int depth;
     836                int counter;
     837                typedef SingleTropicalHomotopyTraverser<mvtyp,mvtypDouble,mvtypDivisor> MySingleTropicalHomotopyTraverser;
     838                std::vector<MySingleTropicalHomotopyTraverser> traversers;
     839                Data fullData;
     840                int level;
     841                bool deadEnd;
     842                bool isLevelLeaf;
     843                bool isSolutionVertex;
     844                std::vector<bool> isLevelLeafStack;
     845                TropicalRegenerationTraverser(std::vector<Matrix<mvtyp> > const &tuple_):
     846                        fullData(tuple_),counter(0),depth(0)
     847                {
     848                        assert(tuple_.size());
     849                        for(int i=0;i<tuple_.size();i++)
     850                                traversers.push_back(MySingleTropicalHomotopyTraverser(fullData.tuples[i],cayleyConfigurationWidth(fullData.tuples[i]),fullData.firstIntersection(),fullData.targets[i],i));
     851                        traversers[0].constructInequalityTableInitially(fullData.degrees[0]);
     852                        level=0;
     853                }
     854                virtual void process()
     855                {
     856                }
     857                bool findOutgoingAndProcess(bool doProcess)
     858                {
     859                        isSolutionVertex=false;
     860                        deadEnd=false;
     861                        isLevelLeaf=traversers[level].findOutgoingAndProcess(false);
     862                        if(isLevelLeaf)
     863                        {//leaf
     864                                bool isFinite=fullData.isFiniteIndex(level,traversers[level].choices[level].first)&&fullData.isFiniteIndex(level,traversers[level].choices[level].second);
     865                                deadEnd=!isFinite;
     866                                if(isFinite && (level==fullData.tuple.size()-1))
     867                                {
     868                                        isSolutionVertex=true;
     869                                        if(doProcess){process();}
     870                                        return true;
     871                                }
     872                        }
     873                        return false;
     874                }
     875                int numberOfChildren()
     876                {
     877                        if(isLevelLeaf&&(level==fullData.tuple.size()-1))return 0;
     878                        if(!isLevelLeaf)
     879                                return traversers[level].numberOfChildren();
     880                        else
     881                                return 1-deadEnd;
     882                }
     883                void goToNthChild(int n)
     884                {
     885                        depth++;
     886                        isLevelLeafStack.push_back(isLevelLeaf);
     887                        if(!isLevelLeaf)
     888                                traversers[level].goToNthChild(n);
     889                        else
     890                        {
     891                                fullData.castToNextLevel(traversers[level].choices,level,fullData.tuples[level][level].getWidth()-fullData.tuples[level+1][level].getWidth(),traversers[level+1].choices);
     892                                traversers[level+1].constructInequalityTableFromParent(traversers[level].inequalityTable,fullData.degrees[level+1]);
     893                                level++;
     894                        }
     895                }
     896                void print()
     897                {
     898                }
     899                void goBack()
     900                {
     901                        depth--;
     902                        counter++;
     903                        deadEnd=false;
     904                        if(traversers[level].atRoot())
     905                                level--;
     906                        else
     907                                traversers[level].goBack();
     908                        isLevelLeaf=isLevelLeafStack.back();
     909                        isLevelLeafStack.pop_back();
     910                }
     911        };
     912
     913        /*
     914        * This class glues Bjarne Knudsen's parallel traverser to our MultiLevelIntersectionTraverser.
     915        * This class should eventually be written so that it works on any homotopy.
     916        *
     917        * Since we are not sure about how exceptions would be handled by the parallel traverser,
     918        * we write ugly aborting code, which may actually work.
     919        */
     920        template<class mvtyp, class mvtypDouble, class mvtypDivisor>
     921        class SpecializedRTraverser: public Traverser
     922        {
     923        public:
     924                typedef TropicalRegenerationTraverser<mvtyp,mvtypDouble,mvtypDivisor> MyTropicalRegenerationTraverser;
     925                MyTropicalRegenerationTraverser T;
     926                mvtypDouble mixedVolume;
     927                int numberOfExpensiveSteps;
     928                SpecializedRTraverser(std::vector<Matrix<mvtyp> > const &tuple_):
     929                        mixedVolume(),
     930                        numberOfExpensiveSteps(0),
     931                        T(tuple_) //Constructor my throw if input is too big.
     932                {
     933                        numberOfExpensiveSteps++;
     934                        T.findOutgoingAndProcess(false);
     935                }
     936                int  getEdgeCountNext( void )
     937                {
     938                        if(!aborting)
     939                        {
     940                                try{
     941                                        return T.numberOfChildren();
     942                                }
     943                                catch(...){abort();}
     944                        }
     945                        return 0;
     946                }
     947
     948                int  moveToNext( int   index,
     949                                   bool  collect_info )
     950                {
     951                        if(!aborting)
     952                        {
     953                                try{
     954                                        T.goToNthChild(index);
     955                                        numberOfExpensiveSteps++;
     956                                        T.findOutgoingAndProcess(false);
     957                                }
     958                                catch(...){abort();}
     959                        }
     960                        return 0;
     961                }
     962
     963          void  moveToPrev( int  index )
     964          {
     965                  if(!aborting)
     966                  {
     967                          try{
     968                                  T.goBack(); //index ignored
     969                          }
     970                          catch(...){abort();}
     971                  }
     972          }
     973
     974          void  collectInfo( void )
     975          {
     976                  if(!aborting)
     977                  {
     978                          try{
     979                                  if(T.isSolutionVertex)
     980                                          mixedVolume.addWithOverflowCheck(T.traversers[T.level].inequalityTable.getVolume().extend());
     981                          }
     982                          catch(...){abort();}
     983                  }
     984          }
     985
     986          void  printState( void )
     987          {
     988                  T.print();
     989          }
     990        };
    991991}
    992992

  • gfanlib/gfanlib_vector.h

    r47a774 r15813d  
    294294  std::string toString()const
    295295  {
    296           std::stringstream f;
    297           f<<*this;
    298           return f.str();
     296          std::stringstream f;
     297          f<<*this;
     298          return f.str();
    299299  }
    300300

  • gfanlib/gfanlib_zcone.h

    r47a774 r15813d  
    1515 * Returns true if cddlib is needed for the ZCone implementation.
    1616 */
    17         bool isCddlibRequired();
    18         /**
    19         * Only call this function if gfanlib is the only code in your program using cddlib.
    20         * Should be paired with a deinitializeCddlibIfRequired() call.
    21         * Calling the function repeatedly may cause memory leaks even if deinitializeCddlibIfRequired() is also called.
    22         */
    23         void initializeCddlibIfRequired();
    24         /**
    25         * This function may do nothing.
    26         */
    27         void deinitializeCddlibIfRequired();
     17        bool isCddlibRequired();
     18        /**
     19        * Only call this function if gfanlib is the only code in your program using cddlib.
     20        * Should be paired with a deinitializeCddlibIfRequired() call.
     21        * Calling the function repeatedly may cause memory leaks even if deinitializeCddlibIfRequired() is also called.
     22        */
     23        void initializeCddlibIfRequired();
     24        /**
     25        * This function may do nothing.
     26        */
     27        void deinitializeCddlibIfRequired();
    2828
    2929/**

  • gfanlib/test.cpp

    r47a774 r15813d  
    1212    catch (...)
    1313    {
    14         cerr<<"Error - most likely an integer overflow."<<endl;
     14        cerr<<"Error - most likely an integer overflow."<<endl;
    1515    }
    1616    return 0;
Note: See TracChangeset for help on using the changeset viewer.