Ignore:
Timestamp:
Aug 8, 2016, 2:16:11 PM (8 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b4f17ed1d25f93d46dbe29e4b499baecc2fd51bb')
Children:
7cc3fdae6091a6fab65626857c1ed2c59f027dcd
Parents:
47a774fa4174b578843d43227fcc351165473747
Message:
format
File:
1 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}
Note: See TracChangeset for help on using the changeset viewer.