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 string & str) : m_value(0), m_cached_result(0), m_cache_valid(false), m_next()
  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         assert(this != NULL);
  58         m_value = a.m_value;
  59         m_cached_result = a.m_cached_result;
  60         for (int i = 0; i < NOP; ++i)
  61         {
  62                 for (auto n : m_next[i])
  63                         delete n;
  64                 m_next[i].clear();
  65                 for (auto n : a.m_next[i])
  66                         m_next[i].push_back(new ParanoidNumber(*n));
  67         }
  68                 /*
  69                 for (unsigned j = 0; j < m_next[i].size() && j < a.m_next[i].size(); ++j)
  70                 {
  71                         if (a.m_next[i][j] != NULL)
  72                                 m_next[i][j]->operator=(*(a.m_next[i][j]));
  73                 }
  74
  75                 for (unsigned j = a.m_next[i].size(); j < m_next[i].size(); ++j)
  76                 {
  77                         delete m_next[i][j];
  78                 }
  79                 m_next[i].resize(a.m_next[i].size());
  80                 */
  81         //}
  82         return *this;
  83 }
  84
  85
  86 string ParanoidNumber::Str() const
  87 {
  88
  89         assert(this != NULL);
  90         string result("");
  91         stringstream s;
  92         s << (double)m_value;
  93         result += s.str();
  94         for (auto mul : m_next[MULTIPLY])
  95         {
  96                 result += "*";
  97                 if (!mul->Floating())
  98                         result += "(" + mul->Str() + ")";
  99                 else
 100                         result += mul->Str();
 101         }
 102         for (auto div : m_next[DIVIDE])
 103         {
 104                 result += "/";
 105                 if (!div->Floating())
 106                         result += "(" + div->Str() + ")";
 107                 else
 108                         result += div->Str();
 109         }
 110
 111         for (auto add : m_next[ADD])
 112         {
 113                 result += "+";
 114                 if (!add->Floating())
 115                         result += "(" + add->Str() + ")";
 116                 else
 117                         result += add->Str();
 118         }
 119         for (auto sub : m_next[SUBTRACT])
 120         {
 121                 result += "-";
 122                 if (!sub->Floating())
 123                         result += "(" + sub->Str() + ")";
 124                 else
 125                         result += sub->Str();
 126         }
 127
 128
 129         return result;
 130 }
 131
 132 template <>
 133 bool TrustingOp<float>(float & a, const float & b, Optype op)
 134 {
 135
 136         feclearexcept(FE_ALL_EXCEPT);
 137         switch (op)
 138         {
 139                 case ADD:
 140                         a += b;
 141                         break;
 142                 case SUBTRACT:
 143                         a -= b;
 144                         break;
 145                 case MULTIPLY:
 146                         a *= b;
 147                         break;
 148                 case DIVIDE:
 149                         if (b == 0)
 150                         {
 151                                 a = (a >= 0) ? INFINITY : -INFINITY;
 152                                 return false;
 153                         }
 154
 155                         a /= b;
 156                         break;
 157                 case NOP:
 158                         break;
 159         }
 160         return !fetestexcept(FE_ALL_EXCEPT);
 161 }
 162
 163 template <>
 164 bool TrustingOp<double>(double & a, const double & b, Optype op)
 165 {
 166         feclearexcept(FE_ALL_EXCEPT);
 167         switch (op)
 168         {
 169                 case ADD:
 170                         a += b;
 171                         break;
 172                 case SUBTRACT:
 173                         a -= b;
 174                         break;
 175                 case MULTIPLY:
 176                         a *= b;
 177                         break;
 178                 case DIVIDE:
 179                         if (b == 0)
 180                         {
 181                                 a = (a >= 0) ? INFINITY : -INFINITY;
 182                                 return false;
 183                         }
 184                         a /= b;
 185                         break;
 186                 case NOP:
 187                         break;
 188         }
 189         return !fetestexcept(FE_ALL_EXCEPT);
 190 }
 191
 192 template <>
 193 bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op)
 194 {
 195         int16_t sa(a);
 196         bool exact = true;
 197         switch (op)
 198         {
 199                 case ADD:
 200                         sa += b;
 201                         exact = (abs(sa) <= 127);
 202                         break;
 203                 case SUBTRACT:
 204                         sa -= b;
 205                         exact = (abs(sa) <= 127);
 206                         break;
 207                 case MULTIPLY:
 208                         sa *= b;
 209                         exact = (abs(sa) <= 127);
 210                         break;
 211                 case DIVIDE:
 212                         exact = (b != 0 && sa > b && sa % b == 0);
 213                         sa /= b;
 214                         break;
 215                 case NOP:
 216                         break;
 217         }
 218         a = (int8_t)(sa);
 219         return exact;
 220 }
 221
 222
 223 ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
 224 {
 225
 226         assert(this != NULL);
 227         delete Operation(SafeConstruct(a), ADD);
 228         Simplify(ADD);
 229         Simplify(SUBTRACT);
 230         return *this;
 231 }
 232
 233
 234 ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
 235 {
 236         delete Operation(SafeConstruct(a), SUBTRACT);
 237         Simplify(SUBTRACT);
 238         Simplify(ADD);
 239         return *this;
 240 }
 241
 242 ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
 243 {
 244         delete Operation(SafeConstruct(a), MULTIPLY);
 245         Simplify(MULTIPLY);
 246         Simplify(DIVIDE);
 247         return *this;
 248 }
 249
 250
 251 ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
 252 {
 253         delete Operation(SafeConstruct(a), DIVIDE);
 254         Simplify(MULTIPLY);
 255         Simplify(DIVIDE);
 256         return *this;
 257 }
 258
 259 // a + b
 260 ParanoidNumber * ParanoidNumber::OperationTerm(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
 261 {
 262         if (!SanityCheck())
 263         {
 264                 Fatal("What...");
 265         }
 266         assert(b->SanityCheck());
 267
 268         m_cached_result = nan("");
 269         if (Floating() && m_value == 0) // 0 + b = b
 270         {
 271                 m_value = b->m_value;
 272                 if (op == SUBTRACT)
 273                 {
 274                         m_value = -m_value;
 275                         swap(b->m_next[ADD], b->m_next[SUBTRACT]);
 276                 }
 277
 278                 for (int i = 0; i < NOP; ++i)
 279                 {
 280                         m_next[i] = b->m_next[i];
 281                         b->m_next[i].clear();
 282                 }
 283
 284                 assert(SanityCheck());
 285                 return b;
 286         }
 287         if (b->Floating() && b->m_value == 0) // a + 0 = a
 288                 return b;
 289
 290
 291
 292         if ((NoFactors() && b->NoFactors())
 293                 || (GetFactors() == b->GetFactors()))
 294         {
 295                 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
 296                 {
 297                         Optype addop = (op == ADD) ? ADD : SUBTRACT;
 298                         for (auto add : b->m_next[ADD])
 299                         {
 300                                 delete OperationTerm(add, addop);
 301                         }
 302                         Optype subop = (op == ADD) ? SUBTRACT : ADD;
 303                         for (auto sub : b->m_next[SUBTRACT])
 304                                 delete OperationTerm(sub, subop);
 305
 306                         b->m_next[ADD].clear();
 307                         b->m_next[SUBTRACT].clear();
 308                         assert(SanityCheck());
 309                         return b;
 310                 }
 311         }
 312
 313
 314
 315         bool parent = (merge_point == NULL);
 316         ParanoidNumber * merge = this;
 317         Optype mop = op;
 318         assert(mop != NOP); // silence compiler warning
 319         if (parent)
 320         {
 321                 merge_point = &merge;
 322                 merge_op = &mop;
 323         }
 324         else
 325         {
 326                 merge = *merge_point;
 327                 mop = *merge_op;
 328         }
 329
 330         Optype invop = InverseOp(op); // inverse of p
 331         Optype fwd = op;
 332         Optype rev = invop;
 333         if (op == SUBTRACT)
 334         {
 335                 fwd = ADD;
 336                 rev = SUBTRACT;
 337         }
 338
 339         for (auto prev : m_next[invop])
 340         {
 341                 if (prev->OperationTerm(b, rev, merge_point, merge_op) == b)
 342                 {
 343                         assert(SanityCheck());
 344                         return b;
 345                 }
 346
 347         }
 348         for (auto next : m_next[op])
 349         {
 350                 if (next->OperationTerm(b, fwd, merge_point, merge_op) == b)
 351                 {
 352                         assert(SanityCheck());
 353                         return b;
 354                 }
 355         }
 356
 357
 358
 359
 360         if (parent)
 361         {
 362                 //merge->m_next[*merge_op].push_back(b);
 363                 m_next[op].push_back(b);
 364         }
 365         else
 366         {
 367                 if (m_next[op].size() == 0)
 368                 {
 369                         *merge_point = this;
 370                         *merge_op = op;
 371                 }
 372         }
 373
 374         assert(SanityCheck());
 375         return NULL;
 376 }
 377
 378 ParanoidNumber * ParanoidNumber::OperationFactor(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
 379 {
 380         assert(SanityCheck());
 381         assert(b->SanityCheck());
 382         m_cached_result = nan("");
 383         if (Floating() && m_value == 0)
 384         {
 385                 return b;
 386         }
 387
 388         if (Floating() && m_value == 1 && op == MULTIPLY)
 389         {
 390                 m_value = b->m_value;
 391                 for (int i = 0; i < NOP; ++i)
 392                 {
 393                         for (auto n : m_next[i])
 394                                 delete n;
 395                         m_next[i].clear();
 396                         swap(m_next[i], b->m_next[i]);
 397                 }
 398                 assert(SanityCheck());
 399                 return b;
 400         }
 401         if (b->Floating() && b->m_value == 1)
 402                 return b;
 403
 404
 405
 406         if (NoTerms() && b->NoTerms())
 407         {
 408                 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
 409                 {
 410                         Optype mulop = (op == MULTIPLY) ? MULTIPLY : DIVIDE;
 411                         for (auto mul : b->m_next[MULTIPLY])
 412                         {
 413                                 delete OperationFactor(mul, mulop);
 414                         }
 415                         Optype divop = (op == MULTIPLY) ? DIVIDE : MULTIPLY;
 416                         for (auto div : b->m_next[DIVIDE])
 417                                 delete OperationFactor(div, divop);
 418
 419                         b->m_next[DIVIDE].clear();
 420                         b->m_next[MULTIPLY].clear();
 421
 422
 423                         assert(SanityCheck());
 424                         return b;
 425                 }
 426         }
 427
 428
 429         bool parent = (merge_point == NULL);
 430         ParanoidNumber * merge = this;
 431         Optype mop = op;
 432         if (parent)
 433         {
 434                 merge_point = &merge;
 435                 merge_op = &mop;
 436         }
 437         else
 438         {
 439                 merge = *merge_point;
 440                 mop = *merge_op;
 441         }
 442
 443         Optype invop = InverseOp(op); // inverse of p
 444         Optype fwd = op;
 445         Optype rev = invop;
 446         if (op == DIVIDE)
 447         {
 448                 fwd = MULTIPLY;
 449                 rev = DIVIDE;
 450         }
 451
 452         ParanoidNumber * cpy_b = NULL;
 453
 454         if (m_next[ADD].size() > 0 || m_next[SUBTRACT].size() > 0)
 455         {
 456                 cpy_b = SafeConstruct(*b);
 457         }
 458
 459         for (auto prev : m_next[invop])
 460         {
 461                 if (prev->OperationFactor(b, rev, merge_point, merge_op) == b)
 462                 {
 463                         for (auto add : m_next[ADD])
 464                                 delete add->OperationFactor(SafeConstruct(*cpy_b), op);
 465                         for (auto sub : m_next[SUBTRACT])
 466                                 delete sub->OperationFactor(SafeConstruct(*cpy_b), op);
 467
 468                         delete cpy_b;
 469                         assert(SanityCheck());
 470                         return b;
 471                 }
 472         }
 473         for (auto next : m_next[op])
 474         {
 475                 if (next->OperationFactor(b, fwd, merge_point, merge_op) == b)
 476                 {
 477                         for (auto add : m_next[ADD])
 478                         {
 479                                 delete add->OperationFactor(SafeConstruct(*cpy_b), op);
 480                         }
 481                         for (auto sub : m_next[SUBTRACT])
 482                         {
 483                                 delete sub->OperationFactor(SafeConstruct(*cpy_b), op);
 484                         }
 485                         delete cpy_b;
 486                         assert(SanityCheck());
 487                         return b;
 488                 }
 489         }
 490
 491         if (parent)
 492         {
 493                 assert(b != NULL);
 494                 m_next[op].push_back(b);
 495                 for (auto add : m_next[ADD])
 496                         delete add->OperationFactor(SafeConstruct(*cpy_b), op);
 497                 for (auto sub : m_next[SUBTRACT])
 498                         delete sub->OperationFactor(SafeConstruct(*cpy_b), op);
 499         }
 500         assert(SanityCheck());
 501         return NULL;
 502 }
 503
 504
 505
 506 /**
 507  * Performs the operation on a with argument b (a += b, a -= b, a *= b, a /= b)
 508  * @returns b if b can safely be deleted
 509  * @returns NULL if b has been merged with a
 510  * append indicates that b should be merged
 511  */
 512 ParanoidNumber * ParanoidNumber::Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
 513 {
 514
 515         if (b == NULL)
 516                 return NULL;
 517
 518
 519         if (op == SUBTRACT || op == ADD)
 520                 return OperationTerm(b, op, merge_point, merge_op);
 521         if (op == MULTIPLY || op == DIVIDE)
 522                 return OperationFactor(b, op, merge_point, merge_op);
 523         return b;
 524 }
 525
 526
 527
 528 string ParanoidNumber::PStr() const
 529 {
 530         stringstream s;
 531         for (int i = 0; i < NOP; ++i)
 532         {
 533                 Optype f = Optype(i);
 534                 s << this;
 535                 for (auto n : m_next[f])
 536                 {
 537                         s << OpChar(f) << n->PStr();
 538                 }
 539         }
 540         return s.str();
 541 }
 542
 543 bool ParanoidNumber::Simplify(Optype op)
 544 {
 545         if (Floating())
 546                 return false;
 547
 548         assert(SanityCheck());
 549         vector<ParanoidNumber*> next;
 550         next.clear();
 551         swap(m_next[op], next);
 552         m_next[op].clear();
 553         assert(m_next[op].size() == 0);
 554         assert(SanityCheck());
 555         Optype fwd = op;
 556         if (op == DIVIDE)
 557                 fwd = MULTIPLY;
 558         else if (op == SUBTRACT)
 559                 fwd = ADD;
 560
 561
 562         for (vector<ParanoidNumber*>::iterator n = next.begin(); n != next.end(); ++n)
 563         {
 564                 if (*n == NULL)
 565                         continue;
 566                 for (vector<ParanoidNumber*>::iterator m = n; m != next.end(); ++m)
 567                 {
 568                         if ((*m) == (*n))
 569                                 continue;
 570                         if (*m == NULL)
 571                                 continue;
 572
 573                         ParanoidNumber * parent = this;
 574                         Optype mop = op;
 575                         ParanoidNumber * result = (*n)->Operation(*m, fwd, &parent, &mop);
 576                         if (result != NULL)
 577                         {
 578                                 *m = NULL;
 579                                 delete result;
 580                         }
 581                 }
 582                 if (*n != NULL)
 583                         delete Operation(*n, op);
 584         }
 585         set<ParanoidNumber*> s;
 586         if (!SanityCheck(s))
 587         {
 588                 Error("Simplify broke Sanity");
 589         }
 590         return (next.size() > m_next[op].size());
 591 }
 592
 593 bool ParanoidNumber::FullSimplify()
 594 {
 595         bool result = false;
 596         result |= Simplify(MULTIPLY);
 597         result |= Simplify(DIVIDE);
 598         result |= Simplify(ADD);
 599         result |= Simplify(SUBTRACT);
 600         return result;
 601 }
 602
 603 ParanoidNumber::digit_t ParanoidNumber::Digit() const
 604 {
 605         if (!SanityCheck())
 606         {
 607                 Fatal("Blargh");
 608         }
 609         //if (!isnan(m_cached_result))
 610         //      return m_cached_result;
 611
 612         // Get around the absurd requirement that const correctness be observed.
 613         digit_t result;// = ((ParanoidNumber*)(this))->m_cached_result;
 614         result = m_value;
 615         for (auto mul : m_next[MULTIPLY])
 616         {
 617                 result *= mul->Digit();
 618         }
 619         for (auto div : m_next[DIVIDE])
 620         {
 621                 result /= div->Digit();
 622         }
 623         for (auto add : m_next[ADD])
 624                 result += add->Digit();
 625         for (auto sub : m_next[SUBTRACT])
 626                 result -= sub->Digit();
 627         return result;
 628
 629 }
 630
 631 ParanoidNumber::digit_t ParanoidNumber::GetFactors() const
 632 {
 633         digit_t value = 1;
 634         for (auto mul : m_next[MULTIPLY])
 635                 value *= mul->Digit();
 636         for (auto div : m_next[DIVIDE])
 637                 value /= div->Digit();
 638         return value;
 639 }
 640
 641
 642 ParanoidNumber::digit_t ParanoidNumber::GetTerms() const
 643 {
 644         digit_t value = 0;
 645         for (auto add : m_next[ADD])
 646                 value += add->Digit();
 647         for (auto sub : m_next[SUBTRACT])
 648                 value -= sub->Digit();
 649         return value;
 650 }
 651
 652 bool ParanoidNumber::SanityCheck(set<ParanoidNumber*> & visited) const
 653 {
 654         if (this == NULL)
 655         {
 656                 Error("NULL pointer in tree");
 657                 return false;
 658         }
 659
 660         if (visited.find((ParanoidNumber*)this) != visited.end())
 661         {
 662                 Error("I think I've seen this tree before...");
 663                 return false;
 664         }
 665
 666         visited.insert((ParanoidNumber*)this);
 667
 668         for (auto add : m_next[ADD])
 669         {
 670                 if (!add->SanityCheck(visited))
 671                         return false;
 672         }
 673         for (auto sub : m_next[SUBTRACT])
 674         {
 675                 if (!sub->SanityCheck(visited))
 676                         return false;
 677         }
 678         for (auto mul : m_next[MULTIPLY])
 679         {
 680                 if (!mul->SanityCheck(visited))
 681                         return false;
 682         }
 683
 684         for (auto div : m_next[DIVIDE])
 685         {
 686                 if (!div->SanityCheck(visited))
 687                         return false;
 688         }
 689         return true;
 690 }
 691
 692 }