src/paranoidnumber.h

   1 #ifndef _PARANOIDNUMBER_H
   2 #define _PARANOIDNUMBER_H
   3
   4 #include <list>
   5 #include <cfloat>
   6 #include <map>
   7 #include <string>
   8 #include "log.h"
   9 #include <fenv.h>
  10
  11 #define PARANOID_DIGIT_T float // we could theoretically replace this with a template
  12                                                                 // but let's not do that...
  13
  14 namespace IPDF
  15 {
  16         typedef enum {ADD, SUBTRACT, MULTIPLY, DIVIDE, NOP} Optype;
  17         inline Optype InverseOp(Optype op)
  18         {
  19                 return ((op == ADD) ? SUBTRACT :
  20                                 (op == SUBTRACT) ? ADD :
  21                                 (op == MULTIPLY) ? DIVIDE :
  22                                 (op == DIVIDE) ? MULTIPLY :
  23                                 (op == NOP) ? NOP : NOP);
  24         }
  25         inline Optype AdjacentOp(Optype op)
  26         {
  27                 return ((op == ADD) ? MULTIPLY :
  28                                 (op == SUBTRACT) ? DIVIDE :
  29                                 (op == MULTIPLY) ? ADD :
  30                                 (op == DIVIDE) ? SUBTRACT :
  31                                 (op == NOP) ? NOP : NOP);
  32         }
  33
  34         inline char OpChar(int op)
  35         {
  36                 static char opch[] = {'+','-','*','/'};
  37                 return (op < NOP && op >= 0) ? opch[op] : '?';
  38         }
  39
  40
  41         /** Performs an operation, returning if the result was exact **/
  42         // NOTE: DIFFERENT to ParanoidOp (although that wraps to this...)
  43         template <class T> bool TrustingOp(T & a, const T & b, Optype op);
  44
  45         /** Performs an operation _only_ if the result would be exact **/
  46         template <class T> bool ParanoidOp(T & a, const T & b, Optype op)
  47         {
  48                 T cpy(a);
  49                 if (TrustingOp<T>(cpy, b, op))
  50                 {
  51                         a = cpy;
  52                         return true;
  53                 }
  54                 return false;
  55         }
  56         template <> bool TrustingOp<float>(float & a, const float & b, Optype op);
  57         template <> bool TrustingOp<double>(double & a, const double & b, Optype op);
  58         template <> bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op);
  59
  60         /**
  61          * A ParanoidNumber
  62          * Idea: Perform regular floating point arithmetic but rearrange operations to only ever use exact results
  63          * Memory Usage: O(all of it)
  64          * CPU Usage: O(all of it)
  65          * Accuracy: O(gives better result for 0.3+0.3+0.3, gives same result for everything else, or worse result)
  66          *
  67          * The ParanoidNumber basically stores 4 linked lists which can be split into two "dimensions"
  68          *  1. Terms to ADD and terms to SUBTRACT
  69          *  2. Factors to MULTIPLY and DIVIDE
  70          * Because ADD and SUBTRACT are inverse operations and MULTIPLY and DIVIDE are inverse operations
  71          * See paranoidnumber.cpp and the ParanoidNumber::Operation function
  72          */
  73         class ParanoidNumber
  74         {
  75
  76                 public:
  77                         typedef PARANOID_DIGIT_T digit_t;
  78
  79                         ParanoidNumber(digit_t value=0) : m_value(value)
  80                         {
  81                                 Construct();
  82                         }
  83
  84                         ParanoidNumber(const ParanoidNumber & cpy) : m_value(cpy.m_value)
  85                         {
  86                                 Construct();
  87                                 for (int i = 0; i < NOP; ++i)
  88                                 {
  89                                         if (cpy.m_next[i] != NULL)
  90                                                 m_next[i] = new ParanoidNumber(*(cpy.m_next[i]));
  91                                 }
  92                         }
  93
  94                         ParanoidNumber(const char * str);
  95                         ParanoidNumber(const std::string & str) : ParanoidNumber(str.c_str()) {Construct();}
  96
  97                         virtual ~ParanoidNumber();
  98
  99                         inline void Construct()
 100                         {
 101                                 for (int i = 0; i < NOP; ++i)
 102                                         m_next[i] = NULL;
 103                                 g_count++;
 104                         }
 105
 106
 107                         template <class T> T Convert() const;
 108                         template <class T> T AddTerms(T value = T(0)) const;
 109                         template <class T> T MultiplyFactors(T value = T(1)) const;
 110                         template <class T> T Head() const {return (m_op == SUBTRACT) ? T(-m_value) : T(m_value);}
 111
 112
 113
 114
 115                         double ToDouble() const {return Convert<double>();}
 116                         float ToFloat() const {return Convert<float>();}
 117                         digit_t Digit() const {return Convert<digit_t>();}
 118
 119                         bool Floating() const
 120                         {
 121                                 for (int i = 0; i < NOP; ++i)
 122                                 {
 123                                         if (m_next[i] != NULL)
 124                                                 return false;
 125                                 }
 126                                 return true;
 127                         }
 128                         bool Sunken() const {return !Floating();} // I could not resist...
 129
 130                         bool Pure(Optype op) const
 131                         {
 132                                 if (op == ADD || op == SUBTRACT)
 133                                         return (m_next[MULTIPLY] == NULL && m_next[DIVIDE] == NULL);
 134                                 return (m_next[ADD] == NULL && m_next[SUBTRACT] == NULL);
 135                         }
 136
 137                         ParanoidNumber & operator+=(const ParanoidNumber & a);
 138                         ParanoidNumber & operator-=(const ParanoidNumber & a);
 139                         ParanoidNumber & operator*=(const ParanoidNumber & a);
 140                         ParanoidNumber & operator/=(const ParanoidNumber & a);
 141                         ParanoidNumber & operator=(const ParanoidNumber & a);
 142
 143
 144                         ParanoidNumber * Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** parent = NULL);
 145                         bool Simplify(Optype op);
 146
 147
 148                         bool operator<(const ParanoidNumber & a) const {return ToDouble() < a.ToDouble();}
 149                         bool operator<=(const ParanoidNumber & a) const {return this->operator<(a) || this->operator==(a);}
 150                         bool operator>(const ParanoidNumber & a) const {return !(this->operator<=(a));}
 151                         bool operator>=(const ParanoidNumber & a) const {return !(this->operator<(a));}
 152                         bool operator==(const ParanoidNumber & a) const {return ToDouble() == a.ToDouble();}
 153                         bool operator!=(const ParanoidNumber & a) const {return !(this->operator==(a));}
 154
 155                         ParanoidNumber operator+(const ParanoidNumber & a) const
 156                         {
 157                                 ParanoidNumber result(*this);
 158                                 result += a;
 159                                 return result;
 160                         }
 161                         ParanoidNumber operator-(const ParanoidNumber & a) const
 162                         {
 163                                 ParanoidNumber result(*this);
 164                                 result -= a;
 165                                 return result;
 166                         }
 167                         ParanoidNumber operator*(const ParanoidNumber & a) const
 168                         {
 169                                 ParanoidNumber result(*this);
 170                                 result *= a;
 171                                 return result;
 172                         }
 173                         ParanoidNumber operator/(const ParanoidNumber & a) const
 174                         {
 175                                 ParanoidNumber result(*this);
 176                                 result /= a;
 177                                 return result;
 178                         }
 179
 180                         std::string Str() const;
 181
 182
 183                         static int64_t Paranoia() {return g_count;}
 184
 185                         std::string PStr() const;
 186
 187                 private:
 188                         static int64_t g_count;
 189                         void Simplify();
 190                         void SimplifyTerms();
 191                         void SimplifyFactors();
 192
 193
 194                         digit_t m_value;
 195                         Optype m_op;
 196                         ParanoidNumber * m_next[4]; // Next by Operation
 197         };
 198
 199 template <class T>
 200 T ParanoidNumber::AddTerms(T value) const
 201 {
 202         ParanoidNumber * add = m_next[ADD];
 203         ParanoidNumber * sub = m_next[SUBTRACT];
 204         while (add != NULL && sub != NULL)
 205         {
 206                 value += add->m_value * add->MultiplyFactors<T>();
 207                 value -= sub->m_value * sub->MultiplyFactors<T>();
 208                 add = add->m_next[ADD];
 209                 sub = sub->m_next[SUBTRACT];
 210         }
 211         while (add != NULL)
 212         {
 213                 value += add->m_value * add->MultiplyFactors<T>();
 214                 add = add->m_next[ADD];
 215         }
 216         while (sub != NULL)
 217         {
 218                 value -= sub->m_value * sub->MultiplyFactors<T>();
 219                 sub = sub->m_next[SUBTRACT];;
 220         }
 221         return value;
 222 }
 223
 224 template <class T>
 225 T ParanoidNumber::MultiplyFactors(T value) const
 226 {
 227         ParanoidNumber * mul = m_next[MULTIPLY];
 228         ParanoidNumber * div = m_next[DIVIDE];
 229         while (mul != NULL && div != NULL)
 230         {
 231                 value *= (mul->m_value + mul->AddTerms<T>());
 232                 value /= (div->m_value + div->AddTerms<T>());
 233                 mul = mul->m_next[MULTIPLY];
 234                 div = div->m_next[DIVIDE];
 235         }
 236         while (mul != NULL)
 237         {
 238                 value *= (mul->m_value + mul->AddTerms<T>());
 239                 mul = mul->m_next[MULTIPLY];
 240         }
 241         while (div != NULL)
 242         {
 243                 value /= (div->m_value + div->AddTerms<T>());
 244                 div = div->m_next[DIVIDE];
 245         }
 246         return value;
 247 }
 248
 249
 250
 251 template <class T>
 252 T ParanoidNumber::Convert() const
 253 {
 254         return MultiplyFactors<T>(m_value) + AddTerms<T>(0);
 255 }
 256
 257
 258
 259 }
 260
 261 #endif //_PARANOIDNUMBER_H
 262