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