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 #include <vector>
  11 #include <cmath>
  12 #include <cassert> // it's going to be ok
  13 #include <set>
  14
  15 #define PARANOID_DIGIT_T double // we could theoretically replace this with a template
  16                                                                 // but let's not do that...
  17
  18
  19 //#define PARANOID_CACHE_RESULTS
  20
  21 //#define PARANOID_USE_ARENA
  22 #define PARANOID_SIZE_LIMIT 0
  23
  24
  25 // Define to compare all ops against double ops and check within epsilon
  26 #define PARANOID_COMPARE_EPSILON 1e-6
  27 #define CompareForSanity(...) this->ParanoidNumber::CompareForSanityEx(__func__, __FILE__, __LINE__, __VA_ARGS__)
  28
  29 namespace IPDF
  30 {
  31         typedef enum {ADD, SUBTRACT, MULTIPLY, DIVIDE, NOP} Optype;
  32         inline Optype InverseOp(Optype op)
  33         {
  34                 return ((op == ADD) ? SUBTRACT :
  35                                 (op == SUBTRACT) ? ADD :
  36                                 (op == MULTIPLY) ? DIVIDE :
  37                                 (op == DIVIDE) ? MULTIPLY :
  38                                 (op == NOP) ? NOP : NOP);
  39         }
  40
  41
  42         inline char OpChar(int op)
  43         {
  44                 static char opch[] = {'+','-','*','/'};
  45                 return (op < NOP && op >= 0) ? opch[op] : '?';
  46         }
  47
  48
  49         /** Performs an operation, returning if the result was exact **/
  50         // NOTE: DIFFERENT to ParanoidOp (although that wraps to this...)
  51         template <class T> bool TrustingOp(T & a, const T & b, Optype op);
  52
  53         /** Performs an operation _only_ if the result would be exact **/
  54         template <class T> bool ParanoidOp(T & a, const T & b, Optype op)
  55         {
  56                 T cpy(a);
  57                 if (TrustingOp<T>(cpy, b, op))
  58                 {
  59                         a = cpy;
  60                         return true;
  61                 }
  62                 return false;
  63         }
  64         template <> bool TrustingOp<float>(float & a, const float & b, Optype op);
  65         template <> bool TrustingOp<double>(double & a, const double & b, Optype op);
  66         template <> bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op);
  67
  68         /**
  69          * A ParanoidNumber
  70          * Idea: Perform regular floating point arithmetic but rearrange operations to only ever use exact results
  71          * Memory Usage: O(all of it)
  72          * CPU Usage: O(all of it)
  73          * Accuracy: O(gives better result for 0.3+0.3+0.3, gives same result for everything else, or worse result)
  74          *
  75          * The ParanoidNumber basically stores 4 linked lists which can be split into two "dimensions"
  76          *  1. Terms to ADD and terms to SUBTRACT
  77          *  2. Factors to MULTIPLY and DIVIDE
  78          * Because ADD and SUBTRACT are inverse operations and MULTIPLY and DIVIDE are inverse operations
  79          * See paranoidnumber.cpp and the ParanoidNumber::Operation function
  80          */
  81         class ParanoidNumber
  82         {
  83
  84                 public:
  85                         typedef PARANOID_DIGIT_T digit_t;
  86
  87                         ParanoidNumber(PARANOID_DIGIT_T value=0) : m_value(value), m_next()
  88                         {
  89                                 #ifdef PARANOID_SIZE_LIMIT
  90                                         m_size = 0;
  91                                 #endif
  92                                 #ifdef PARANOID_CACHE_RESULTS
  93                                         m_cached_result = value;
  94                                 #endif
  95                         }
  96
  97                         ParanoidNumber(const ParanoidNumber & cpy) : m_value(cpy.m_value), m_next()
  98                         {
  99
 100                                 #ifdef PARANOID_SIZE_LIMIT
 101                                         m_size = cpy.m_size;
 102                                 #endif
 103                                 #ifdef PARANOID_CACHE_RESULTS
 104                                         m_cached_result = cpy.m_cached_result;
 105                                 #endif
 106                                 for (int i = 0; i < NOP; ++i)
 107                                 {
 108                                         for (auto next : cpy.m_next[i])
 109                                         {
 110                                                 if (next != NULL) // why would this ever be null
 111                                                         m_next[i].push_back(new ParanoidNumber(*next)); // famous last words...
 112                                         }
 113                                 }
 114                                 //assert(SanityCheck());
 115                         }
 116
 117                         //ParanoidNumber(const char * str);
 118                         ParanoidNumber(const std::string & str);// : ParanoidNumber(str.c_str()) {}
 119
 120                         virtual ~ParanoidNumber();
 121
 122
 123                         bool SanityCheck(std::set<ParanoidNumber*> & visited) const;
 124                         bool SanityCheck() const
 125                         {
 126                                 std::set<ParanoidNumber*> s;
 127                                 return SanityCheck(s);
 128                         }
 129
 130                         template <class T> T Convert() const;
 131                         digit_t GetFactors() const;
 132                         digit_t GetTerms() const;
 133
 134                         // This function is declared const purely to trick the compiler.
 135                         // It is not actually const, and therefore, none of the other functions that call it are const either.
 136                         digit_t Digit() const;
 137
 138                         // Like this one. It isn't const.
 139                         double ToDouble() const {return (double)Digit();}
 140
 141                         // This one is probably const.
 142                         bool Floating() const
 143                         {
 144                                 return NoFactors() && NoTerms();
 145                         }
 146                         bool Sunken() const {return !Floating();} // I could not resist...
 147
 148                         bool NoFactors() const {return (m_next[MULTIPLY].size() == 0 && m_next[DIVIDE].size() == 0);}
 149                         bool NoTerms() const {return (m_next[ADD].size() == 0 && m_next[SUBTRACT].size() == 0);}
 150
 151
 152                         ParanoidNumber & operator+=(const ParanoidNumber & a);
 153                         ParanoidNumber & operator-=(const ParanoidNumber & a);
 154                         ParanoidNumber & operator*=(const ParanoidNumber & a);
 155                         ParanoidNumber & operator/=(const ParanoidNumber & a);
 156                         ParanoidNumber & operator=(const ParanoidNumber & a);
 157
 158                         ParanoidNumber & operator+=(const digit_t & a);
 159                         ParanoidNumber & operator-=(const digit_t & a);
 160                         ParanoidNumber & operator*=(const digit_t & a);
 161                         ParanoidNumber & operator/=(const digit_t & a);
 162                         ParanoidNumber & operator=(const digit_t & a);
 163
 164
 165                         ParanoidNumber * OperationTerm(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point = NULL, Optype * mop = NULL);
 166                         ParanoidNumber * OperationFactor(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point = NULL, Optype * mop = NULL);
 167                         ParanoidNumber * TrivialOp(ParanoidNumber * b, Optype op);
 168                         ParanoidNumber * Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point = NULL, Optype * mop = NULL);
 169                         bool Simplify(Optype op);
 170                         bool FullSimplify();
 171
 172
 173                         // None of these are actually const
 174                         bool operator<(const ParanoidNumber & a) const {return ToDouble() < a.ToDouble();}
 175                         bool operator<=(const ParanoidNumber & a) const {return this->operator<(a) || this->operator==(a);}
 176                         bool operator>(const ParanoidNumber & a) const {return !(this->operator<=(a));}
 177                         bool operator>=(const ParanoidNumber & a) const {return !(this->operator<(a));}
 178                         bool operator==(const ParanoidNumber & a) const {return ToDouble() == a.ToDouble();}
 179                         bool operator!=(const ParanoidNumber & a) const {return !(this->operator==(a));}
 180
 181                         ParanoidNumber operator-() const
 182                         {
 183                                 ParanoidNumber neg(0);
 184                                 neg -= *this;
 185                                 return neg;
 186                         }
 187
 188                         ParanoidNumber operator+(const ParanoidNumber & a) const
 189                         {
 190                                 ParanoidNumber result(*this);
 191                                 a.SanityCheck();
 192                                 result += a;
 193                                 return result;
 194                         }
 195                         ParanoidNumber operator-(const ParanoidNumber & a) const
 196                         {
 197                                 ParanoidNumber result(*this);
 198                                 result -= a;
 199                                 return result;
 200                         }
 201                         ParanoidNumber operator*(const ParanoidNumber & a) const
 202                         {
 203                                 ParanoidNumber result(*this);
 204                                 result *= a;
 205                                 if (!result.SanityCheck())
 206                                 {
 207                                         Fatal("Blargh");
 208                                 }
 209                                 return result;
 210                         }
 211                         ParanoidNumber operator/(const ParanoidNumber & a) const
 212                         {
 213                                 ParanoidNumber result(*this);
 214                                 result /= a;
 215                                 return result;
 216                         }
 217
 218                         std::string Str() const;
 219
 220                         inline void CompareForSanityEx(const char * func, const char * file, int line, const digit_t & compare, const digit_t & arg, const digit_t & eps = PARANOID_COMPARE_EPSILON)
 221                         {
 222                                 if (fabs(Digit() - compare) > eps)
 223                                 {
 224                                         Error("Called via %s(%lf) (%s:%d)", func, arg, file, line);
 225                                         Error("Failed: %s", Str().c_str());
 226                                         Fatal("This: %.30lf vs Expected: %.30lf", Digit(), compare);
 227                                 }
 228                         }
 229
 230
 231                         std::string PStr() const;
 232
 233                         #ifdef PARANOID_USE_ARENA
 234                         void * operator new(size_t byes);
 235                         void operator delete(void * p);
 236                         #endif //PARANOID_USE_ARENA
 237
 238                 private:
 239
 240                         void Simplify();
 241                         void SimplifyTerms();
 242                         void SimplifyFactors();
 243
 244                         digit_t m_value;
 245                         #ifdef PARANOID_CACHE_RESULTS
 246                                 digit_t m_cached_result;
 247                         #endif
 248                         std::vector<ParanoidNumber*> m_next[4];
 249                         #ifdef PARANOID_SIZE_LIMIT
 250                                 int64_t m_size;
 251                         #endif //PARANOID_SIZE_LIMIT
 252
 253                         #ifdef PARANOID_USE_ARENA
 254                         class Arena
 255                         {
 256                                 public:
 257                                         Arena(int64_t block_size = 10000000);
 258                                         ~Arena();
 259
 260                                         void * allocate(size_t bytes);
 261                                         void deallocate(void * p);
 262
 263                                 private:
 264                                         struct Block
 265                                         {
 266                                                 void * memory;
 267                                                 int64_t used;
 268                                         };
 269
 270                                         std::vector<Block> m_blocks;
 271                                         int64_t m_block_size;
 272
 273                                         void * m_spare;
 274
 275                         };
 276
 277                         static Arena g_arena;
 278                         #endif //PARANOID_USE_ARENA
 279
 280
 281         };
 282
 283
 284
 285
 286 template <class T>
 287 T ParanoidNumber::Convert() const
 288 {
 289         #ifdef PARANOID_CACHE_RESULTS
 290         if (!isnan((float(m_cached_result))))
 291                 return (T)m_cached_result;
 292         #endif
 293         T value(m_value);
 294         for (auto mul : m_next[MULTIPLY])
 295         {
 296                 value *= mul->Convert<T>();
 297         }
 298         for (auto div : m_next[DIVIDE])
 299         {
 300                 value /= div->Convert<T>();
 301         }
 302         for (auto add : m_next[ADD])
 303                 value += add->Convert<T>();
 304         for (auto sub : m_next[SUBTRACT])
 305                 value -= sub->Convert<T>();
 306         return value;
 307 }
 308
 309
 310
 311 }
 312
 313 #endif //_PARANOIDNUMBER_H
 314