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         {
  20                 for (auto n : m_next[i])
  21                         delete n;
  22         }
  23 }
  24
  25 ParanoidNumber::ParanoidNumber(const char * str) : m_value(0), m_cached_result(0)
  26 {
  27         Construct();
  28         int dp = 0;
  29         int end = 0;
  30         while (str[dp] != '\0' && str[dp] != '.')
  31         {
  32                 ++dp;
  33                 ++end;
  34         }
  35         while (str[end] != '\0')
  36                 ++end;
  37         ParanoidNumber m(1);
  38         for (int i = dp-1; i >= 0; --i)
  39         {
  40                 ParanoidNumber b(str[i]-'0');
  41                 b*=m;
  42                 this->operator+=(b);
  43                 m*=10;
  44         }
  45         ParanoidNumber n(1);
  46         for (int i = dp+1; i < end; ++i)
  47         {
  48                 n/=10;
  49                 ParanoidNumber b(str[i]-'0');
  50                 b*=n;
  51                 this->operator+=(b);
  52         }
  53 }
  54
  55 ParanoidNumber & ParanoidNumber::operator=(const ParanoidNumber & a)
  56 {
  57         m_value = a.m_value;
  58         m_cached_result = a.m_cached_result;
  59         for (int i = 0; i < NOP; ++i)
  60         {
  61                 for (unsigned j = 0; j < m_next[i].size() && j < a.m_next[i].size(); ++j)
  62                 {
  63                         m_next[i][j]->operator=(*(a.m_next[i][j]));
  64                 }
  65
  66                 for (unsigned j = a.m_next[i].size(); j < m_next[i].size(); ++j)
  67                 {
  68                         delete m_next[i][j];
  69                 }
  70                 m_next[i].resize(a.m_next[i].size());
  71         }
  72         return *this;
  73 }
  74
  75
  76 string ParanoidNumber::Str() const
  77 {
  78         string result("");
  79         stringstream s;
  80         s << (double)m_value;
  81         result += s.str();
  82         for (auto mul : m_next[MULTIPLY])
  83         {
  84                 result += "*";
  85                 if (!mul->Floating())
  86                         result += "(" + mul->Str() + ")";
  87                 else
  88                         result += mul->Str();
  89         }
  90         for (auto div : m_next[DIVIDE])
  91         {
  92                 result += "/";
  93                 if (!div->Floating())
  94                         result += "(" + div->Str() + ")";
  95                 else
  96                         result += div->Str();
  97         }
  98
  99         for (auto add : m_next[ADD])
 100         {
 101                 result += "+";
 102                 if (!add->Floating())
 103                         result += "(" + add->Str() + ")";
 104                 else
 105                         result += add->Str();
 106         }
 107         for (auto sub : m_next[SUBTRACT])
 108         {
 109                 result += "-";
 110                 if (!sub->Floating())
 111                         result += "(" + sub->Str() + ")";
 112                 else
 113                         result += sub->Str();
 114         }
 115
 116
 117         return result;
 118 }
 119
 120 template <>
 121 bool TrustingOp<float>(float & a, const float & b, Optype op)
 122 {
 123         feclearexcept(FE_ALL_EXCEPT);
 124         switch (op)
 125         {
 126                 case ADD:
 127                         a += b;
 128                         break;
 129                 case SUBTRACT:
 130                         a -= b;
 131                         break;
 132                 case MULTIPLY:
 133                         a *= b;
 134                         break;
 135                 case DIVIDE:
 136                         a /= b;
 137                         break;
 138                 case NOP:
 139                         break;
 140         }
 141         return !fetestexcept(FE_ALL_EXCEPT);
 142 }
 143
 144 template <>
 145 bool TrustingOp<double>(double & a, const double & b, Optype op)
 146 {
 147         feclearexcept(FE_ALL_EXCEPT);
 148         switch (op)
 149         {
 150                 case ADD:
 151                         a += b;
 152                         break;
 153                 case SUBTRACT:
 154                         a -= b;
 155                         break;
 156                 case MULTIPLY:
 157                         a *= b;
 158                         break;
 159                 case DIVIDE:
 160                         a /= b;
 161                         break;
 162                 case NOP:
 163                         break;
 164         }
 165         return !fetestexcept(FE_ALL_EXCEPT);
 166 }
 167
 168 template <>
 169 bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op)
 170 {
 171         int16_t sa(a);
 172         bool exact = true;
 173         switch (op)
 174         {
 175                 case ADD:
 176                         sa += b;
 177                         exact = (abs(sa) <= 127);
 178                         break;
 179                 case SUBTRACT:
 180                         sa -= b;
 181                         exact = (abs(sa) <= 127);
 182                         break;
 183                 case MULTIPLY:
 184                         sa *= b;
 185                         exact = (abs(sa) <= 127);
 186                         break;
 187                 case DIVIDE:
 188                         exact = (b != 0 && sa > b && sa % b == 0);
 189                         sa /= b;
 190                         break;
 191                 case NOP:
 192                         break;
 193         }
 194         a = (int8_t)(sa);
 195         return exact;
 196 }
 197
 198
 199 ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
 200 {
 201         delete Operation(new ParanoidNumber(a), ADD);
 202         Simplify(ADD);
 203         Simplify(SUBTRACT);
 204         return *this;
 205 }
 206
 207
 208 ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
 209 {
 210         delete Operation(new ParanoidNumber(a), SUBTRACT);
 211         //Simplify(SUBTRACT);
 212         //Simplify(ADD);
 213         return *this;
 214 }
 215
 216 ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
 217 {
 218         delete Operation(new ParanoidNumber(a), MULTIPLY);
 219         return *this;
 220 }
 221
 222
 223 ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
 224 {
 225         delete Operation(new ParanoidNumber(a), DIVIDE);
 226         return *this;
 227 }
 228
 229 // a + b
 230 ParanoidNumber * ParanoidNumber::OperationTerm(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
 231 {
 232         m_cached_result = nan("");
 233         if (Floating() && m_value == 0) // 0 + b = b
 234         {
 235                 m_value = b->m_value;
 236                 if (op == SUBTRACT)
 237                 {
 238                         m_value = -m_value;
 239                         swap(b->m_next[ADD], b->m_next[SUBTRACT]);
 240                 }
 241
 242                 for (int i = 0; i < NOP; ++i)
 243                 {
 244                         m_next[i] = b->m_next[i];
 245                         b->m_next[i].clear();
 246                 }
 247                 return b;
 248         }
 249         if (b->Floating() && b->m_value == 0) // a + 0 = a
 250                 return b;
 251
 252
 253
 254         if ((NoFactors() && b->NoFactors())
 255                 || (GetFactors() == b->GetFactors()))
 256         {
 257                 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
 258                 {
 259                         Optype addop = (op == ADD) ? ADD : SUBTRACT;
 260                         for (auto add : b->m_next[ADD])
 261                         {
 262                                 delete OperationTerm(add, addop);
 263                         }
 264                         Optype subop = (op == ADD) ? SUBTRACT : ADD;
 265                         for (auto sub : b->m_next[SUBTRACT])
 266                                 delete OperationTerm(sub, subop);
 267
 268                         b->m_next[ADD].clear();
 269                         b->m_next[SUBTRACT].clear();
 270                         return b;
 271                 }
 272         }
 273
 274
 275
 276         bool parent = (merge_point == NULL);
 277         ParanoidNumber * merge = this;
 278         Optype mop = op;
 279         assert(mop != NOP); // silence compiler warning
 280         if (parent)
 281         {
 282                 merge_point = &merge;
 283                 merge_op = &mop;
 284         }
 285         else
 286         {
 287                 merge = *merge_point;
 288                 mop = *merge_op;
 289         }
 290
 291         Optype invop = InverseOp(op); // inverse of p
 292         Optype fwd = op;
 293         Optype rev = invop;
 294         if (op == SUBTRACT)
 295         {
 296                 fwd = ADD;
 297                 rev = SUBTRACT;
 298         }
 299
 300         for (auto prev : m_next[invop])
 301         {
 302                 if (prev->OperationTerm(b, rev, merge_point, merge_op) == b)
 303                         return b;
 304
 305         }
 306         for (auto next : m_next[op])
 307         {
 308                 if (next->OperationTerm(b, fwd, merge_point, merge_op) == b)
 309                         return b;
 310         }
 311
 312
 313
 314
 315         if (parent)
 316         {
 317                 //merge->m_next[*merge_op].push_back(b);
 318                 m_next[op].push_back(b);
 319         }
 320         else
 321         {
 322                 if (m_next[op].size() == 0)
 323                 {
 324                         *merge_point = this;
 325                         *merge_op = op;
 326                 }
 327         }
 328         return NULL;
 329 }
 330
 331 ParanoidNumber * ParanoidNumber::OperationFactor(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
 332 {
 333         m_cached_result = nan("");
 334         if (Floating() && m_value == 0)
 335         {
 336                 return b;
 337         }
 338
 339         if (Floating() && m_value == 1 && op == MULTIPLY)
 340         {
 341                 m_value = b->m_value;
 342                 for (int i = 0; i < NOP; ++i)
 343                 {
 344                         for (auto n : m_next[i])
 345                                 delete n;
 346                         m_next[i].clear();
 347                         swap(m_next[i], b->m_next[i]);
 348                 }
 349                 return b;
 350         }
 351         if (b->Floating() && b->m_value == 1)
 352                 return b;
 353
 354
 355
 356         if (NoTerms() && b->NoTerms())
 357         {
 358                 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
 359                 {
 360                         Optype mulop = (op == MULTIPLY) ? MULTIPLY : DIVIDE;
 361                         for (auto mul : b->m_next[MULTIPLY])
 362                         {
 363                                 delete OperationFactor(mul, mulop);
 364                         }
 365                         Optype divop = (op == MULTIPLY) ? DIVIDE : MULTIPLY;
 366                         for (auto div : b->m_next[DIVIDE])
 367                                 delete OperationFactor(div, divop);
 368
 369                         b->m_next[DIVIDE].clear();
 370                         b->m_next[MULTIPLY].clear();
 371
 372
 373
 374                         return b;
 375                 }
 376         }
 377
 378
 379         bool parent = (merge_point == NULL);
 380         ParanoidNumber * merge = this;
 381         Optype mop = op;
 382         if (parent)
 383         {
 384                 merge_point = &merge;
 385                 merge_op = &mop;
 386         }
 387         else
 388         {
 389                 merge = *merge_point;
 390                 mop = *merge_op;
 391         }
 392
 393         Optype invop = InverseOp(op); // inverse of p
 394         Optype fwd = op;
 395         Optype rev = invop;
 396         if (op == DIVIDE)
 397         {
 398                 fwd = MULTIPLY;
 399                 rev = DIVIDE;
 400         }
 401
 402         ParanoidNumber * cpy_b = NULL;
 403
 404         if (m_next[ADD].size() > 0 || m_next[SUBTRACT].size() > 0)
 405         {
 406                 cpy_b = new ParanoidNumber(*b);
 407         }
 408
 409         for (auto prev : m_next[invop])
 410         {
 411                 if (prev->OperationFactor(b, rev, merge_point, merge_op) == b)
 412                 {
 413                         for (auto add : m_next[ADD])
 414                                 delete add->OperationFactor(new ParanoidNumber(*cpy_b), op);
 415                         for (auto sub : m_next[SUBTRACT])
 416                                 delete sub->OperationFactor(new ParanoidNumber(*cpy_b), op);
 417
 418                         delete cpy_b;
 419                         return b;
 420                 }
 421         }
 422         for (auto next : m_next[op])
 423         {
 424                 if (next->OperationFactor(b, fwd, merge_point, merge_op) == b)
 425                 {
 426                         for (auto add : m_next[ADD])
 427                                 delete add->OperationFactor(new ParanoidNumber(*cpy_b), op);
 428                         for (auto sub : m_next[SUBTRACT])
 429                                 delete sub->OperationFactor(new ParanoidNumber(*cpy_b), op);
 430                         delete cpy_b;
 431                         return b;
 432                 }
 433         }
 434
 435         if (parent)
 436         {
 437                 m_next[op].push_back(b);
 438                 for (auto add : m_next[ADD])
 439                         delete add->OperationFactor(new ParanoidNumber(*cpy_b), op);
 440                 for (auto sub : m_next[SUBTRACT])
 441                         delete sub->OperationFactor(new ParanoidNumber(*cpy_b), op);
 442         }
 443         return NULL;
 444 }
 445
 446
 447
 448 /**
 449  * Performs the operation on a with argument b (a += b, a -= b, a *= b, a /= b)
 450  * @returns b if b can safely be deleted
 451  * @returns NULL if b has been merged with a
 452  * append indicates that b should be merged
 453  */
 454 ParanoidNumber * ParanoidNumber::Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
 455 {
 456
 457         if (b == NULL)
 458                 return NULL;
 459
 460
 461         if (op == SUBTRACT || op == ADD)
 462                 return OperationTerm(b, op, merge_point, merge_op);
 463         if (op == MULTIPLY || op == DIVIDE)
 464                 return OperationFactor(b, op, merge_point, merge_op);
 465         return b;
 466 }
 467
 468
 469
 470 string ParanoidNumber::PStr() const
 471 {
 472         stringstream s;
 473         for (int i = 0; i < NOP; ++i)
 474         {
 475                 Optype f = Optype(i);
 476                 s << this;
 477                 for (auto n : m_next[f])
 478                 {
 479                         s << OpChar(f) << n->PStr();
 480                 }
 481         }
 482         return s.str();
 483 }
 484
 485 bool ParanoidNumber::Simplify(Optype op)
 486 {
 487         vector<ParanoidNumber*> next(0);
 488         swap(m_next[op], next);
 489         for (auto n : next)
 490         {
 491                 delete Operation(n, op);
 492         }
 493         return (next.size() > m_next[op].size());
 494 }
 495
 496 bool ParanoidNumber::FullSimplify()
 497 {
 498         bool result = false;
 499         result |= Simplify(MULTIPLY);
 500         result |= Simplify(DIVIDE);
 501         result |= Simplify(ADD);
 502         result |= Simplify(SUBTRACT);
 503         return result;
 504 }
 505
 506
 507 ParanoidNumber::digit_t ParanoidNumber::Digit()
 508 {
 509         if (!isnan(m_cached_result))
 510                 return m_cached_result;
 511         m_cached_result = m_value;
 512         for (auto mul : m_next[MULTIPLY])
 513         {
 514                 m_cached_result *= mul->Digit();
 515         }
 516         for (auto div : m_next[DIVIDE])
 517         {
 518                 m_cached_result /= div->Digit();
 519         }
 520         for (auto add : m_next[ADD])
 521                 m_cached_result += add->Digit();
 522         for (auto sub : m_next[SUBTRACT])
 523                 m_cached_result -= sub->Digit();
 524         return m_cached_result;
 525
 526 }
 527
 528 ParanoidNumber::digit_t ParanoidNumber::GetFactors()
 529 {
 530         digit_t value = 1;
 531         for (auto mul : m_next[MULTIPLY])
 532                 value *= mul->Digit();
 533         for (auto div : m_next[DIVIDE])
 534                 value /= div->Digit();
 535         return value;
 536 }
 537
 538
 539 ParanoidNumber::digit_t ParanoidNumber::GetTerms()
 540 {
 541         digit_t value = 0;
 542         for (auto add : m_next[ADD])
 543                 value += add->Digit();
 544         for (auto sub : m_next[SUBTRACT])
 545                 value -= sub->Digit();
 546         return value;
 547 }
 548
 549
 550 }