source: git/gfanlib/gfanlib_circuittableint.h

spielwiese
Last change on this file was aa4d31, checked in by Reimer Behrends <behrends@…>, 5 years ago
parallel preprocessor-related fixes.
  • Property mode set to 100644
File size: 13.7 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 "../factory/globaldefs.h"
12
13#include <cstdint>
14#include <exception>
15#include <sstream>
16#include <assert.h>
17
18namespace gfan{
19
20class MVMachineIntegerOverflow: public std::exception
21{
22  virtual const char* what() const throw()
23  {
24    return "Exception: Overflow (or possible future overflow) in 32/64 bit integers in tropical homotopy.";
25  }
26};
27
28typedef MVMachineIntegerOverflow MVMachineIntegerOverflowType;
29
30EXTERN_INST_VAR MVMachineIntegerOverflowType MVMachineIntegerOverflow;
31
32/*
33 * The philosophy here is that if this class overflows, then the computation needs to be restarted. Therefore
34 * all overflows must be caught.
35 */
36class CircuitTableInt32{
37 public:
38        class Divisor{
39        public:
40                int32_t v;
41                int shift;
42                int32_t multiplicativeInverse;
43                Divisor(CircuitTableInt32 const &a)// for exact division
44                { // A good place to read about these tricks seems to be the book "Hacker's Delight" by Warren.
45                        v=a.v;
46                        shift=0;
47                        int32_t t=v;
48//                        uint32_t t=v;
49                        assert(t);
50                        while(!(t&1)){t>>=        1;shift++;}
51                        uint32_t inverse=t;
52                        while(t*inverse!=1)inverse*=2-t*inverse;
53                        multiplicativeInverse=inverse;
54                }
55        };
56        class Double{
57        public:
58                int64_t v;
59                Double():v(0){};
60                Double(int64_t a):v(a){};
61                Double &operator+=(Double a){v+=a.v;return *this;}
62                Double &operator-=(Double a){v-=a.v;return *this;}
63                CircuitTableInt32 castToSingle()const;
64                bool isZero()const{return v==0;}
65                bool isNegative()const{return v<0;}
66                Double operator-()const{Double ret;ret.v=-v;return ret;}
67                friend Double operator-(Double const &a, Double const &b){return Double(a.v-b.v);}
68                friend Double operator+(Double const &a, Double const &b){return Double(a.v+b.v);}
69                Double &addWithOverflowCheck(Double const &a)
70                {
71                        if(a.v<0 || v<0 || a.v>1000000000000000000 || v>1000000000000000000) throw MVMachineIntegerOverflow;
72                        v+=a.v;
73                        return *this;
74                }
75                std::string toString()const{std::stringstream s;s<<v;return s.str();}
76        };
77 private:
78public:
79        int32_t v;
80        friend CircuitTableInt32 operator+(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return CircuitTableInt32(a.v+b.v);}
81        friend CircuitTableInt32 operator-(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return CircuitTableInt32(a.v-b.v);}
82        friend CircuitTableInt32 operator*(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return CircuitTableInt32(a.v*b.v);}
83 public:
84        // 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:
85        bool fitsInHalf()const{return v>-30000 && v<30000;}// Is it better to allow only non-negative numbers?
86        static bool isSuitableSummationNumber(int numberOfSummands){return numberOfSummands<30000;}
87        CircuitTableInt32(){v=0;};
88        CircuitTableInt32(CircuitTableInt32 const &m):v(m.v){}
89        CircuitTableInt32(int32_t val):v(val){};
90        class CircuitTableInt32 &operator=(int32_t a){v=a;return *this;}
91        CircuitTableInt32 operator-()const{CircuitTableInt32 ret;ret.v=-v;return ret;}
92        CircuitTableInt32 &operator-=(CircuitTableInt32 a){v-=a.v;return *this;}
93        CircuitTableInt32 &operator+=(CircuitTableInt32 a){v+=a.v;return *this;}
94        CircuitTableInt32 &operator*=(CircuitTableInt32 a){v*=a.v;return *this;}
95        friend bool operator<=(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return a.v<=b.v;}
96        friend bool operator==(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return a.v==b.v;}
97        bool isZero()const{return v==0;}
98        bool isOne()const{return v==1;}
99        bool isNonZero()const{return v!=0;}
100        bool isNegative()const{return v<0;}
101        bool isPositive()const{return v>0;}
102        friend int determinantSign1(CircuitTableInt32 const &a, CircuitTableInt32 const &c, CircuitTableInt32 const &b, CircuitTableInt32 const &d)//NOTICE MIXED ORDER. MUST WORK BY EXTENDING
103        {
104                int64_t r=((int64_t)a.v)*((int64_t)c.v)-((int64_t)b.v)*((int64_t)d.v);
105                if(r>0)return 1;
106                if(r<0)return -1;
107                return 0;
108        }
109        friend int determinantSign2(CircuitTableInt32::Double const &a, CircuitTableInt32 const &c, CircuitTableInt32::Double const &b, CircuitTableInt32 const &d)//NOTICE MIXED ORDER. MUST WORK BY EXTENDING
110        {
111                int64_t r1=(a.v)*((int64_t)c.v);
112                int64_t r2=(b.v)*((int64_t)d.v);
113                if(r1>r2)return 1;
114                if(r1<r2)return -1;
115                return 0;
116        }
117        std::string toString()const{std::stringstream s;s<<v;return s.str();}
118        friend std::ostream &operator<<(std::ostream &f, CircuitTableInt32 const &a){f<<(int)a.v;return f;}
119        Double extend()const{Double ret;ret.v=v;return ret;}
120        CircuitTableInt32 &maddWithOverflowChecking(CircuitTableInt32 const &a, CircuitTableInt32 const&b){Double t=this->extend();t+=extendedMultiplication(a,b);*this=t.castToSingle();return *this;}
121        CircuitTableInt32 &msubWithOverflowChecking(CircuitTableInt32 const &a, CircuitTableInt32 const&b){Double s=this->extend();s-=extendedMultiplication(a,b);*this=s.castToSingle();return *this;}
122        CircuitTableInt32 &subWithOverflowChecking(CircuitTableInt32 const &a){Double t=this->extend();t-=a.extend();*this=t.castToSingle();return *this;}
123        CircuitTableInt32 &addWithOverflowChecking(CircuitTableInt32 const &a){Double t=this->extend();t+=a.extend();*this=t.castToSingle();return *this;}
124        CircuitTableInt32 negatedWithOverflowChecking()const{return (-extend()).castToSingle();}
125        friend Double extendedMultiplication(CircuitTableInt32 const &a, CircuitTableInt32 const &b){Double ret;ret.v=((int64_t)a.v)*((int64_t)b.v);return ret;}
126        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?
127        friend CircuitTableInt32 min(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return (a.v>=b.v)?b:a;}
128        friend CircuitTableInt32 max(CircuitTableInt32 const &a, CircuitTableInt32 const &b){return (a.v<=b.v)?b:a;}
129        friend CircuitTableInt32 negabs(CircuitTableInt32 const &a){return min(a,-a);}
130        friend CircuitTableInt32 dotDiv(CircuitTableInt32 const &s, CircuitTableInt32 const &a, CircuitTableInt32 const &t, CircuitTableInt32 const &b, CircuitTableInt32::Divisor const &q)
131        {
132                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));
133        }
134        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.
135        {
136                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);
137        }
138        static int MIN(int32_t a, int32_t b)
139        {
140                return (a<b)?a:b;
141        }
142        static int MAX(int32_t a, int32_t b)
143        {
144                return (a>b)?a:b;
145        }
146        static CircuitTableInt32 computeNegativeBound(CircuitTableInt32 * __restrict__ Ai, int w)
147        {
148                CircuitTableInt32 M=0;
149                CircuitTableInt32 m=0;
150                for(int j=0;j<w;j++)
151                {
152                        m.v=MIN(m.v,Ai[j].v);
153                        M.v=MAX(M.v,Ai[j].v);
154                }
155                return min(m,-M);
156        }
157        static CircuitTableInt32 quickNoShiftBounded(CircuitTableInt32 * __restrict__ a, CircuitTableInt32 * __restrict__ b, CircuitTableInt32 s, CircuitTableInt32 t, CircuitTableInt32::Divisor denominatorDivisor,int c)
158        {
159                CircuitTableInt32 *aa=a;
160                CircuitTableInt32 *bb=b;
161                CircuitTableInt32 max=0;
162            CircuitTableInt32 min=0;
163            CircuitTableInt32 ret=0;
164            for(int i=0;i<c;i++)
165              {
166                    aa[i].v=((s.v*denominatorDivisor.multiplicativeInverse)*aa[i].v+
167                                    (t.v*denominatorDivisor.multiplicativeInverse)*bb[i].v);
168                    min.v=MIN(min.v,aa[i].v);
169                    max.v=MAX(max.v,aa[i].v);
170              }
171            if(-max<=min)min=-max;
172            ret=min;
173            return ret;
174        }
175        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
176        {
177                while(((s.v|t.v)&1)==0 && denominatorDivisor.shift>0){s.v>>=1;t.v>>=1;denominatorDivisor.shift--;denominatorDivisor.v>>=1;}
178
179                CircuitTableInt32 max=0;
180                CircuitTableInt32 min=0;
181
182          // This is a bound on the absolute value of the sums of products before dividing by the denominator
183          // The calculation can overflow, but since we are casting to unsigned long, this is not a problem.
184          uint64_t positiveResultBoundTimesD=(negabs(t).v*((int64_t)boundA.v)+negabs(s).v*((int64_t)boundB.v));
185
186          /* The first case is the case where we know that the intermediate results fit in 32bit before shifting down.
187           * In other words, the shifted-in bits of the results are equal.
188           * In that case the intermediate result can be computed with 32bits.
189           */
190          if(positiveResultBoundTimesD<((((int64_t)0x40000000)*denominatorDivisor.v)>>denominatorDivisor.shift))//check this carefully
191          {//FAST VERSION
192                  // debug<<s.v<<" "<<t.v<<" "<<denominatorDivisor.v<<" "<<denominatorDivisor.shift<<"\n";
193                  if(denominatorDivisor.shift==0)return quickNoShiftBounded(aa,bb,s,t,denominatorDivisor,c);
194                  for(int i=0;i<c;i++)
195                  {
196                          aa[i].v=((s.v*denominatorDivisor.multiplicativeInverse)*aa[i].v+
197                                          (t.v*denominatorDivisor.multiplicativeInverse)*bb[i].v)>>denominatorDivisor.shift;
198                          min.v=MIN(min.v,aa[i].v);
199                          max.v=MAX(max.v,aa[i].v);
200                  }
201                  if(-max<=min)min=-max;
202                  return min;
203          }
204          else
205          {
206                  /* The next case would be to check if the result is known to fit in 32bit.
207                   * In that case intermediate results would be handled in 64 bit,
208                   * but the final result would not have to be checked for 32 bit overflows.
209                   */
210                  if(positiveResultBoundTimesD<((((int64_t)0x40000000)*denominatorDivisor.v)))
211                  {
212                          for(int i=0;i<c;i++)
213                          {
214                                  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);
215                                  if(max<=aa[i])max=aa[i];
216                                  if(aa[i]<=min)min=aa[i];
217                          }
218                          if(-max<=min)min=-max;
219                  }
220                  /* The last case would be to where we don't know that the results will fit in 32 bit.
221                   * Since we need to compute max and min anyway, we can compute these quantities in 64bit
222                   * and just check if they fit in 32bit at the end.
223                   */
224                  else
225                  {
226                          bool doesOverflow=(((uint32_t)t.v)==0x80000000);
227                          int64_t min64=0;
228                          int64_t max64=0;
229                          for(int i=0;i<c;i++)
230                          {
231                                  // 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.
232                                  int64_t temp=(((int64_t)s.v)*((int64_t)aa[i].v)+((int64_t)t.v)*((int64_t)bb[i].v))/denominatorDivisor.v;
233                                  if(max64<=temp)max64=temp;
234                                  if(temp<=min64)min64=temp;
235                                  aa[i].v=temp;
236                          }
237                          if(-max64<=min64)min64=-max64;
238                          if(min64<=-0x7fffffff)doesOverflow=true;
239                          if(doesOverflow)throw MVMachineIntegerOverflow;
240                          min=min64;
241                  }
242          }
243          return min;
244        }
245};
246
247inline CircuitTableInt32 CircuitTableInt32::Double::castToSingle()const//casts and checks precission
248{
249        if(v>=0x7fffffff || -v>=0x7fffffff) throw MVMachineIntegerOverflow;
250        return CircuitTableInt32(v);
251}
252}
253
254
255#endif /* GFANLIB_CIRCUITTABLETYPEINT_H_ */
Note: See TracBrowser for help on using the repository browser.