source: git/gfanlib/gfanlib_circuittableint.h @ 5cea7a

spielwiese
Last change on this file since 5cea7a was 5cea7a, checked in by Yue <ren@…>, 8 years ago
chg: new gfanlib version 0.6
  • Property mode set to 100644
File size: 11.0 KB
Line 
1/*
2 * gfanlib_circuittabletypeint.h
3 *
4 *  Created on: Apr 10, 2016
5 *      Author: anders
6 */
7
8#ifndef GFANLIB_CIRCUITTABLEINT_H_
9#define GFANLIB_CIRCUITTABLEINT_H_
10
11#include <cstdint>
12#include <exception>
13#include <sstream>
14#include <assert.h>
15
16namespace gfan{
17
18class MVMachineIntegerOverflow: public std::exception
19{
20  virtual const char* what() const throw()
21  {
22    return "Exception: Overflow (or possible future overflow) in 32/64 bit integers in tropical homotopy.";
23  }
24};
25
26extern MVMachineIntegerOverflow MVMachineIntegerOverflow;
27
28/*
29 * The philosophy here is that if this class overflows, then the computation needs to be restarted. Therefore
30 * all overflows must be caught.
31 */
32class CircuitTableInt32{
33 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        };
72 private:
73public:
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);}
78 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        }
239};
240
241inline CircuitTableInt32 CircuitTableInt32::Double::castToSingle()const//casts and checks precission
242{
243        if(v>=0x7fffffff || -v>=0x7fffffff) throw MVMachineIntegerOverflow;
244        return CircuitTableInt32(v);
245}
246}
247
248
249#endif /* GFANLIB_CIRCUITTABLETYPEINT_H_ */
Note: See TracBrowser for help on using the repository browser.