src/paranoidnumber.cpp

   1 #include "paranoidnumber.h"
   2
   3 #include <sstream>
   4 #include <fenv.h>
   5 #include "log.h"
   6 #include <cassert>
   7 #include <iostream>
   8
   9 using namespace std;
  10 namespace IPDF
  11 {
  12 int64_t ParanoidNumber::g_count = 0;
  13
  14
  15 ParanoidNumber::~ParanoidNumber()
  16 {
  17         g_count--;
  18         for (int i = 0; i < NOP; ++i)
  19                 delete m_next[i];
  20 }
  21
  22 ParanoidNumber::ParanoidNumber(const char * str) : m_value(0)
  23 {
  24         Construct();
  25         int dp = 0;
  26         int end = 0;
  27         while (str[dp] != '\0' && str[dp] != '.')
  28         {
  29                 ++dp;
  30                 ++end;
  31         }
  32         while (str[end] != '\0')
  33                 ++end;
  34         ParanoidNumber m(1);
  35         for (int i = dp-1; i >= 0; --i)
  36         {
  37                 ParanoidNumber b(str[i]-'0');
  38                 b*=m;
  39                 this->operator+=(b);
  40                 m*=10;
  41         }
  42         ParanoidNumber n(1);
  43         for (int i = dp+1; i < end; ++i)
  44         {
  45                 Debug("{%s} /= 10", n.Str().c_str());
  46                 n/=10;
  47                 Debug("{%s}", n.Str().c_str());
  48                 ParanoidNumber b(str[i]-'0');
  49                 b*=n;
  50                 Debug("{%s} += {%s}", Str().c_str(), b.Str().c_str());
  51                 this->operator+=(b);
  52         }
  53 }
  54
  55 ParanoidNumber & ParanoidNumber::operator=(const ParanoidNumber & a)
  56 {
  57         m_value = a.m_value;
  58         for (int i = 0; i < NOP; ++i)
  59         {
  60                 if (a.m_next[i] == NULL)
  61                 {
  62                         if (m_next[i] != NULL)
  63                                 delete m_next[i];
  64                         m_next[i] = NULL;
  65                         continue;
  66                 }
  67
  68                 if (m_next[i] != NULL)
  69                 {
  70                         m_next[i]->operator=(*(a.m_next[i]));
  71                 }
  72                 else
  73                 {
  74                         m_next[i] = new ParanoidNumber(*(a.m_next[i]));
  75                 }
  76         }
  77         return *this;
  78 }
  79
  80
  81 string ParanoidNumber::Str() const
  82 {
  83         string result("");
  84         stringstream s;
  85         s << (double)m_value;
  86         result += s.str();
  87         if (m_next[MULTIPLY] != NULL)
  88         {
  89                 result += "*";
  90                 if (m_next[MULTIPLY]->m_next[ADD] != NULL || m_next[MULTIPLY]->m_next[SUBTRACT] != NULL)
  91                         result += "(" + m_next[MULTIPLY]->Str() + ")";
  92                 else
  93                         result += m_next[MULTIPLY]->Str();
  94         }
  95         if (m_next[DIVIDE] != NULL)
  96         {
  97                 result += "/";
  98                 if (m_next[DIVIDE]->m_next[ADD] != NULL || m_next[DIVIDE]->m_next[SUBTRACT] != NULL)
  99                         result += "(" + m_next[DIVIDE]->Str() + ")";
 100                 else
 101                         result += m_next[DIVIDE]->Str();
 102         }
 103
 104         if (m_next[ADD] != NULL)
 105         {
 106                 result += "+";
 107                 if (m_next[ADD]->m_next[MULTIPLY] != NULL || m_next[ADD]->m_next[DIVIDE] != NULL)
 108                         result += "(" + m_next[ADD]->Str() + ")";
 109                 else
 110                         result += m_next[ADD]->Str();
 111         }
 112         if (m_next[SUBTRACT] != NULL)
 113         {
 114                 result += "-";
 115                 if (m_next[SUBTRACT]->m_next[MULTIPLY] != NULL || m_next[SUBTRACT]->m_next[DIVIDE] != NULL)
 116                         result += "(" + m_next[SUBTRACT]->Str() + ")";
 117                 else
 118                         result += m_next[SUBTRACT]->Str();
 119         }
 120
 121
 122         return result;
 123 }
 124
 125 template <>
 126 bool TrustingOp<float>(float & a, const float & b, Optype op)
 127 {
 128         feclearexcept(FE_ALL_EXCEPT);
 129         switch (op)
 130         {
 131                 case ADD:
 132                         a += b;
 133                         break;
 134                 case SUBTRACT:
 135                         a -= b;
 136                         break;
 137                 case MULTIPLY:
 138                         a *= b;
 139                         break;
 140                 case DIVIDE:
 141                         a /= b;
 142                         break;
 143                 case NOP:
 144                         break;
 145         }
 146         return !fetestexcept(FE_ALL_EXCEPT);
 147 }
 148
 149 template <>
 150 bool TrustingOp<double>(double & a, const double & b, Optype op)
 151 {
 152         feclearexcept(FE_ALL_EXCEPT);
 153         switch (op)
 154         {
 155                 case ADD:
 156                         a += b;
 157                         break;
 158                 case SUBTRACT:
 159                         a -= b;
 160                         break;
 161                 case MULTIPLY:
 162                         a *= b;
 163                         break;
 164                 case DIVIDE:
 165                         a /= b;
 166                         break;
 167                 case NOP:
 168                         break;
 169         }
 170         return !fetestexcept(FE_ALL_EXCEPT);
 171 }
 172
 173 template <>
 174 bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op)
 175 {
 176         int16_t sa(a);
 177         bool exact = true;
 178         switch (op)
 179         {
 180                 case ADD:
 181                         sa += b;
 182                         exact = (abs(sa) <= 127);
 183                         break;
 184                 case SUBTRACT:
 185                         sa -= b;
 186                         exact = (abs(sa) <= 127);
 187                         break;
 188                 case MULTIPLY:
 189                         sa *= b;
 190                         exact = (abs(sa) <= 127);
 191                         break;
 192                 case DIVIDE:
 193                         exact = (b != 0 && sa > b && sa % b == 0);
 194                         sa /= b;
 195                         break;
 196                 case NOP:
 197                         break;
 198         }
 199         a = (int8_t)(sa);
 200         return exact;
 201 }
 202
 203
 204 ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
 205 {
 206         delete Operation(new ParanoidNumber(a), ADD);
 207         return *this;
 208 }
 209
 210
 211 ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
 212 {
 213         delete Operation(new ParanoidNumber(a), SUBTRACT);
 214         return *this;
 215 }
 216
 217 ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
 218 {
 219         delete Operation(new ParanoidNumber(a), MULTIPLY);
 220         return *this;
 221 }
 222
 223
 224 ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
 225 {
 226         delete Operation(new ParanoidNumber(a), DIVIDE);
 227         return *this;
 228 }
 229
 230 /**
 231  * Performs the operation on a with argument b (a += b, a -= b, a *= b, a /= b)
 232  * @returns b if b can safely be deleted
 233  * @returns NULL if b has been merged with a
 234  * append indicates that b should be merged
 235  */
 236 ParanoidNumber * ParanoidNumber::Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** parent)
 237 {
 238         if (b == NULL)
 239                 return NULL;
 240
 241         Optype invop = InverseOp(op); // inverse of p
 242         ParanoidNumber * append_at = this;
 243
 244         if (Floating())
 245         {
 246                 if ((op == ADD || op == SUBTRACT) && (m_value == 0))
 247                 {
 248                         m_value = b->m_value;
 249                         for (int i = 0; i < NOP; ++i)
 250                         {
 251                                 m_next[i] = b->m_next[i];
 252                                 b->m_next[i] = NULL;
 253                         }
 254                         return b;
 255                 }
 256                 if ((op == MULTIPLY) && (m_value == 1))
 257                 {
 258                         m_value = b->m_value;
 259                         for (int i = 0; i < NOP; ++i)
 260                         {
 261                                 m_next[i] = b->m_next[i];
 262                                 b->m_next[i] = NULL;
 263                         }
 264                         return b;
 265                         return b;
 266                 }
 267
 268         }
 269
 270         if (b->Floating())
 271         {
 272                 if ((op == ADD || op == SUBTRACT) && (b->m_value == 0))
 273                         return b;
 274                 if ((op == MULTIPLY || op == DIVIDE) && (b->m_value == 1))
 275                         return b;
 276         }
 277
 278         // Operation can be applied directly to the m_value of this and b
 279         // ie: op is + or - and this and b have no * or / children
 280         // or: op is * or / and this and b have no + or - children
 281         if (Pure(op) && (b->Pure(op)))
 282         {
 283                 if (ParanoidOp<digit_t>(m_value, b->m_value, op)) // op applied successfully...
 284                 {
 285                         Simplify(op);
 286                         Simplify(invop);
 287                         for (int i = 0; i < NOP; ++i) // Try applying b's children to this
 288                         {
 289                                 delete Operation(b->m_next[i], Optype(i));
 290                                 b->m_next[i] = NULL;
 291                         }
 292                         return b; // can delete b
 293                 }
 294         }
 295
 296         // Try to simplify the cases:
 297         // a + b*c == (a/c + b)*c
 298         // a + b/c == (a*c + b)/c
 299         else if ((op == ADD || op == SUBTRACT) &&
 300                         (Pure(op) || b->Pure(op)))
 301         {
 302
 303                 Debug("Simplify: {%s} %c {%s}", Str().c_str(), OpChar(op), b->Str().c_str());
 304                 Optype adj[] = {MULTIPLY, DIVIDE};
 305                 for (int i = 0; i < 2; ++i)
 306                 {
 307
 308                         Optype f = adj[i];
 309                         Optype invf = InverseOp(f);
 310
 311                         Debug("Try %c", OpChar(f));
 312
 313                         if (m_next[f] == NULL && b->m_next[f] == NULL)
 314                                 continue;
 315
 316                         ParanoidNumber * tmp_a = new ParanoidNumber(*this);
 317                         ParanoidNumber * tmp_b = new ParanoidNumber(*b);
 318
 319
 320                         ParanoidNumber * af = (tmp_a->m_next[f] != NULL) ? new ParanoidNumber(*(tmp_a->m_next[f])) : NULL;
 321                         ParanoidNumber * bf = (tmp_b->m_next[f] != NULL) ? new ParanoidNumber(*(tmp_b->m_next[f])) : NULL;
 322
 323                         Debug("{%s} %c {%s}", tmp_a->Str().c_str(), OpChar(op), tmp_b->Str().c_str());
 324                         Debug("{%s} %c {%s}", tmp_a->Str().c_str(), OpChar(op), tmp_b->Str().c_str());
 325                         if (tmp_a->Operation(af, invf) != af || tmp_b->Operation(bf, invf) != bf)
 326                         {
 327                                 delete af;
 328                                 delete bf;
 329                                 delete tmp_a;
 330                                 delete tmp_b;
 331                                 continue;
 332                         }
 333                         Debug("{%s} %c {%s}", tmp_a->Str().c_str(), OpChar(op), tmp_b->Str().c_str());
 334
 335                         if (tmp_a->Operation(bf, invf) == bf && tmp_b->Operation(af, invf) == af) // a / c simplifies
 336                         {
 337                                 if (tmp_a->Operation(tmp_b, op) != NULL) // (a/c) + b simplifies
 338                                 {
 339                                         this->operator=(*tmp_a);
 340                                         if (bf != NULL)
 341                                                 delete Operation(bf, f);
 342                                         if (af != NULL)
 343                                                 delete Operation(af, f);
 344                                         delete tmp_a;
 345                                         delete tmp_b;
 346                                         return b; // It simplified after all!
 347                                 }
 348                                 else
 349                                 {
 350                                         tmp_b = NULL;
 351                                         delete af;
 352                                         delete bf;
 353                                 }
 354                         }
 355                         //Debug("tmp_a : %s", tmp_a->PStr().c_str());
 356                         //Debug("tmp_b : %s", tmp_b->PStr().c_str());
 357                         delete tmp_a;
 358                         delete tmp_b;
 359                 }
 360         }
 361
 362                 // See if operation can be applied to children of this in the same dimension
 363         {
 364                 // (a / b) / c = a / (b*c)
 365                 // (a * b) / c = a * (b/c)
 366                 // (a / b) * c = a / (b/c)
 367                 // (a * b) * c = a * (b*c)
 368                 // (a + b) + c = a + (b+c)
 369                 // (a - b) + c = a - (b-c)
 370                 // (a + b) - c = a + (b-c)
 371                 // (a - b) - c = a - (b+c)
 372                 Optype fwd(op);
 373                 Optype rev(invop);
 374                 if (op == DIVIDE || op == SUBTRACT)
 375                 {
 376                         fwd = invop;
 377                         rev = op;
 378                 }
 379                 // opposite direction first (because ideally things will cancel each other out...)
 380                 if (m_next[invop] != NULL && m_next[invop]->Operation(b, rev, &append_at) != NULL)
 381                         return b;
 382                 // forward direction
 383                 if (m_next[op] != NULL && m_next[op]->Operation(b, fwd, &append_at) != NULL)
 384                         return b;
 385         }
 386
 387         // At this point, we have no choice but to merge 'b' with this ParanoidNumber
 388
 389         // we are a child; the merge operation needs to be applied by the root, so leave
 390         if (parent != NULL)
 391         {
 392                 if (m_next[op] == NULL)
 393                         *parent = this; // last element in list
 394                 return NULL;
 395         }
 396
 397         append_at->m_next[op] = b; // Merge with b
 398
 399         // MULTIPLY and DIVIDE operations need to be performed on each term in the ADD/SUBTRACT dimension
 400         if (op == DIVIDE || op == MULTIPLY)
 401         {
 402                 // apply the operation to each term
 403                 if (m_next[ADD] != NULL) delete m_next[ADD]->Operation(new ParanoidNumber(*b), op);
 404                 if (m_next[SUBTRACT] != NULL) delete m_next[SUBTRACT]->Operation(new ParanoidNumber(*b), op);
 405
 406                 // try and simplify this by adding the terms (you never know...)
 407                 Simplify(ADD);
 408                 Simplify(SUBTRACT);
 409         }
 410         // failed to simplify
 411         return NULL;
 412 }
 413
 414 bool ParanoidNumber::Simplify(Optype op)
 415 {
 416         ParanoidNumber * n = m_next[op];
 417         m_next[op] = NULL;
 418         if (Operation(n, Optype(op)))
 419         {
 420                 delete n;
 421                 return true;
 422         }
 423         else
 424         {
 425                 m_next[op] = n;
 426                 return false;
 427         }
 428 }
 429
 430 string ParanoidNumber::PStr() const
 431 {
 432         stringstream s;
 433         for (int i = 0; i < NOP; ++i)
 434         {
 435                 Optype f = Optype(i);
 436                 s << this << OpChar(f) << m_next[f] << "\n";
 437         }
 438         return s.str();
 439 }
 440
 441
 442
 443
 444
 445 }