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 // here be many copy paste bugs
  10
  11 using namespace std;
  12 namespace IPDF
  13 {
  14
  15
  16 #ifdef PARANOID_USE_ARENA
  17 ParanoidNumber::Arena ParanoidNumber::g_arena;
  18 #endif //PARANOID_USE_ARENA
  19
  20 ParanoidNumber::~ParanoidNumber()
  21 {
  22         for (int i = 0; i < NOP; ++i)
  23         {
  24                 for (auto n : m_next[i])
  25                         delete n;
  26         }
  27 }
  28
  29 ParanoidNumber::ParanoidNumber(const string & str) : m_value(0), m_next()
  30 {
  31         #ifdef PARANOID_SIZE_LIMIT
  32                 m_size = 1;
  33         #endif
  34         #ifdef PARANOID_CACHE_RESULTS
  35         m_cached_result = NAN;
  36         #endif
  37
  38         int dp = 0;
  39         int end = 0;
  40         bool negate = str[0] == '-';
  41         if (negate)
  42         {
  43                 dp++;
  44                 end++;
  45         }
  46         while (str[dp] != '\0' && str[dp] != '.')
  47         {
  48                 ++dp;
  49                 ++end;
  50         }
  51         while (str[end] != '\0')
  52                 ++end;
  53         ParanoidNumber m(1);
  54         for (int i = dp-1; i >= negate; --i)
  55         {
  56                 ParanoidNumber b(str[i]-'0');
  57                 b*=m;
  58                 this->operator+=(b);
  59                 m*=10;
  60         }
  61         ParanoidNumber n(1);
  62         for (int i = dp+1; i < end; ++i)
  63         {
  64                 n/=10;
  65                 ParanoidNumber b(str[i]-'0');
  66                 b*=n;
  67                 this->operator+=(b);
  68         }
  69
  70         if (negate)
  71                 Negate();
  72
  73         #ifdef PARANOID_COMPARE_EPSILON
  74                 double d = strtod(str.c_str(), NULL);
  75                 CompareForSanity(d, d);
  76         #endif
  77 }
  78
  79 ParanoidNumber & ParanoidNumber::operator=(const ParanoidNumber & a)
  80 {
  81         //assert(this != NULL);
  82
  83         #ifdef PARANOID_SIZE_LIMIT
  84                 m_size = a.m_size;
  85         #endif
  86
  87         m_value = a.m_value;
  88         #ifdef PARANOID_CACHE_RESULT
  89         m_cached_result = a.m_cached_result;
  90         #endif
  91         for (int i = 0; i < NOP; ++i)
  92         {
  93                 for (auto n : m_next[i])
  94                         delete n;
  95                 m_next[i].clear();
  96                 for (auto n : a.m_next[i])
  97                         m_next[i].push_back(new ParanoidNumber(*n));
  98         }
  99         #ifdef PARANOID_COMPARE_EPSILON
 100                 CompareForSanity(a.Digit(),a.Digit());
 101         #endif
 102         return *this;
 103 }
 104
 105
 106 string ParanoidNumber::Str() const
 107 {
 108
 109         //assert(this != NULL);
 110         string result("");
 111         stringstream s;
 112         s << (double)m_value;
 113         result += s.str();
 114         for (auto mul : m_next[MULTIPLY])
 115         {
 116                 result += "*";
 117                 if (!mul->Floating())
 118                         result += "(" + mul->Str() + ")";
 119                 else
 120                         result += mul->Str();
 121         }
 122         for (auto div : m_next[DIVIDE])
 123         {
 124                 result += "/";
 125                 if (!div->Floating())
 126                         result += "(" + div->Str() + ")";
 127                 else
 128                         result += div->Str();
 129         }
 130
 131         for (auto add : m_next[ADD])
 132         {
 133                 result += "+";
 134                 if (!add->Floating())
 135                         result += "(" + add->Str() + ")";
 136                 else
 137                         result += add->Str();
 138         }
 139         for (auto sub : m_next[SUBTRACT])
 140         {
 141                 result += "-";
 142                 if (!sub->Floating())
 143                         result += "(" + sub->Str() + ")";
 144                 else
 145                         result += sub->Str();
 146         }
 147
 148
 149         return result;
 150 }
 151
 152 template <>
 153 bool TrustingOp<float>(float & a, const float & b, Optype op)
 154 {
 155
 156
 157         feclearexcept(FE_ALL_EXCEPT);
 158         switch (op)
 159         {
 160                 case ADD:
 161                         a += b;
 162                         break;
 163                 case SUBTRACT:
 164                         a -= b;
 165                         break;
 166                 case MULTIPLY:
 167                         a *= b;
 168                         break;
 169                 case DIVIDE:
 170                         if (b == 0)
 171                         {
 172                                 a = (a >= 0) ? INFINITY : -INFINITY;
 173                                 return false;
 174                         }
 175
 176                         a /= b;
 177                         break;
 178                 case NOP:
 179                         break;
 180         }
 181         return !fetestexcept(FE_ALL_EXCEPT);
 182 }
 183
 184 template <>
 185 bool TrustingOp<double>(double & a, const double & b, Optype op)
 186 {
 187         feclearexcept(FE_ALL_EXCEPT);
 188         switch (op)
 189         {
 190                 case ADD:
 191                         a += b;
 192                         break;
 193                 case SUBTRACT:
 194                         a -= b;
 195                         break;
 196                 case MULTIPLY:
 197                         a *= b;
 198                         break;
 199                 case DIVIDE:
 200                         if (b == 0)
 201                         {
 202                                 a = (a >= 0) ? INFINITY : -INFINITY;
 203                                 return false;
 204                         }
 205                         a /= b;
 206                         break;
 207                 case NOP:
 208                         break;
 209         }
 210         return !fetestexcept(FE_ALL_EXCEPT);
 211 }
 212
 213
 214
 215 template <>
 216 bool TrustingOp<long double>(long double & a, const long double & b, Optype op)
 217 {
 218
 219
 220         feclearexcept(FE_ALL_EXCEPT);
 221         switch (op)
 222         {
 223                 case ADD:
 224                         a += b;
 225                         break;
 226                 case SUBTRACT:
 227                         a -= b;
 228                         break;
 229                 case MULTIPLY:
 230                         a *= b;
 231                         break;
 232                 case DIVIDE:
 233                         if (b == 0)
 234                         {
 235                                 a = (a >= 0) ? INFINITY : -INFINITY;
 236                                 return false;
 237                         }
 238
 239                         a /= b;
 240                         break;
 241                 case NOP:
 242                         break;
 243         }
 244         return !fetestexcept(FE_ALL_EXCEPT);
 245 }
 246
 247
 248 template <>
 249 bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op)
 250 {
 251         int16_t sa(a);
 252         bool exact = true;
 253         switch (op)
 254         {
 255                 case ADD:
 256                         sa += b;
 257                         exact = (abs(sa) <= 127);
 258                         break;
 259                 case SUBTRACT:
 260                         sa -= b;
 261                         exact = (abs(sa) <= 127);
 262                         break;
 263                 case MULTIPLY:
 264                         sa *= b;
 265                         exact = (abs(sa) <= 127);
 266                         break;
 267                 case DIVIDE:
 268                         exact = (b != 0 && sa > b && sa % b == 0);
 269                         sa /= b;
 270                         break;
 271                 case NOP:
 272                         break;
 273         }
 274         a = (int8_t)(sa);
 275         return exact;
 276 }
 277
 278
 279 ParanoidNumber & ParanoidNumber::operator+=(const digit_t & a)
 280 {
 281         #ifdef PARANOID_COMPARE_EPSILON
 282                 digit_t compare = Digit();
 283                 compare += a;
 284         #endif
 285         delete Operation(new ParanoidNumber(a), ADD);
 286         Simplify(SUBTRACT);
 287         Simplify(ADD);
 288         #ifdef PARANOID_COMPARE_EPSILON
 289                 CompareForSanity(compare, a);
 290         #endif
 291         return *this;
 292 }
 293
 294
 295 ParanoidNumber & ParanoidNumber::operator-=(const digit_t & a)
 296 {
 297         #ifdef PARANOID_COMPARE_EPSILON
 298                 digit_t compare = Digit();
 299                 compare -= a;
 300         #endif
 301         delete Operation(new ParanoidNumber(a), SUBTRACT);
 302         Simplify(ADD);
 303         Simplify(SUBTRACT);
 304         #ifdef PARANOID_COMPARE_EPSILON
 305                 CompareForSanity(compare, a);
 306         #endif
 307         return *this;
 308 }
 309
 310 ParanoidNumber & ParanoidNumber::operator*=(const digit_t & a)
 311 {
 312         #ifdef PARANOID_COMPARE_EPSILON
 313                 digit_t compare = Digit();
 314                 compare *= a;
 315         #endif
 316         delete Operation(new ParanoidNumber(a), MULTIPLY);
 317         Simplify(DIVIDE);
 318         Simplify(MULTIPLY);
 319         #ifdef PARANOID_COMPARE_EPSILON
 320                 CompareForSanity(compare, a);
 321         #endif
 322         return *this;
 323 }
 324
 325
 326 ParanoidNumber & ParanoidNumber::operator/=(const digit_t & a)
 327 {
 328         #ifdef PARANOID_COMPARE_EPSILON
 329                 digit_t compare = Digit();
 330                 compare /= a;
 331         #endif
 332         delete Operation(new ParanoidNumber(a), DIVIDE);
 333         Simplify(MULTIPLY);
 334         Simplify(DIVIDE);
 335         #ifdef PARANOID_COMPARE_EPSILON
 336                 CompareForSanity(compare, a);
 337         #endif
 338         return *this;
 339 }
 340
 341
 342 ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
 343 {
 344         #ifdef PARANOID_COMPARE_EPSILON
 345                 digit_t compare = Digit();
 346                 compare += a.Digit();
 347         #endif
 348         delete Operation(new ParanoidNumber(a), ADD);
 349         Simplify(SUBTRACT);
 350         Simplify(ADD);
 351         #ifdef PARANOID_COMPARE_EPSILON
 352                 CompareForSanity(compare, a.Digit());
 353         #endif
 354         return *this;
 355 }
 356
 357
 358 ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
 359 {
 360         #ifdef PARANOID_COMPARE_EPSILON
 361                 digit_t compare = Digit();
 362                 compare -= a.Digit();
 363         #endif
 364         delete Operation(new ParanoidNumber(a), SUBTRACT);
 365         Simplify(ADD);
 366         Simplify(SUBTRACT);
 367         #ifdef PARANOID_COMPARE_EPSILON
 368                 CompareForSanity(compare, a.Digit());
 369         #endif
 370         return *this;
 371 }
 372
 373 ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
 374 {
 375         #ifdef PARANOID_COMPARE_EPSILON
 376                 digit_t compare = Digit();
 377                 compare *= a.Digit();
 378         #endif
 379         delete Operation(new ParanoidNumber(a), MULTIPLY);
 380         Simplify(DIVIDE);
 381         Simplify(MULTIPLY);
 382         #ifdef PARANOID_COMPARE_EPSILON
 383                 CompareForSanity(compare, a.Digit());
 384         #endif
 385         return *this;
 386 }
 387
 388
 389 ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
 390 {
 391         #ifdef PARANOID_COMPARE_EPSILON
 392                 digit_t compare = Digit();
 393                 compare /= a.Digit();
 394         #endif
 395         delete Operation(new ParanoidNumber(a), DIVIDE);
 396         Simplify(MULTIPLY);
 397         Simplify(DIVIDE);
 398         #ifdef PARANOID_COMPARE_EPSILON
 399                 CompareForSanity(compare, a.Digit());
 400         #endif
 401         return *this;
 402 }
 403
 404 ParanoidNumber & ParanoidNumber::operator=(const digit_t & a)
 405 {
 406
 407         for (int i = 0; i < NOP; ++i)
 408         {
 409                 for (auto n : m_next[i])
 410                         delete n;
 411                 m_next[i].clear();
 412         }
 413         m_value = a;
 414         #ifdef PARANOID_CACHE_RESULT
 415         m_cached_result = a;
 416         #endif
 417
 418         #ifdef PARANOID_COMPARE_EPSILON
 419                 CompareForSanity(a,a);
 420         #endif
 421
 422         return *this;
 423 }
 424
 425 // a + b
 426 ParanoidNumber * ParanoidNumber::OperationTerm(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
 427 {
 428         ////assert(b->SanityCheck());
 429         #ifdef PARANOID_CACHE_RESULTS
 430         m_cached_result = NAN;
 431         #endif
 432         #ifdef PARANOID_SIZE_LIMIT
 433                 if (m_size + b->m_size >= PARANOID_SIZE_LIMIT)
 434                 {
 435                         this->operator=(this->Digit());
 436                         if (op == ADD)
 437                                 m_value += b->Digit();
 438                         else
 439                                 m_value -= b->Digit();
 440                         m_size = 1;
 441                         //Debug("Cut off %p", this);
 442                         return b;
 443                 }
 444                 //Debug("At size limit %d", m_size);
 445         #endif
 446
 447
 448         if (Floating() && m_value == 0) // 0 + b = b
 449         {
 450                 m_value = b->m_value;
 451                 if (op == SUBTRACT)
 452                 {
 453                         m_value = -m_value;
 454                         swap(b->m_next[ADD], b->m_next[SUBTRACT]);
 455                 }
 456
 457                 for (int i = 0; i < NOP; ++i)
 458                 {
 459                         m_next[i] = b->m_next[i];
 460                         b->m_next[i].clear();
 461                 }
 462
 463                 //assert(SanityCheck());
 464                 return b;
 465         }
 466         if (b->Floating() && b->m_value == 0) // a + 0 = a
 467                 return b;
 468
 469
 470
 471         if ((NoFactors() && b->NoFactors())
 472                 || (GetFactors() == b->GetFactors()))
 473         {
 474                 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
 475                 {
 476                         Optype addop = (op == ADD) ? ADD : SUBTRACT;
 477                         for (auto add : b->m_next[ADD])
 478                         {
 479                                 delete (OperationTerm(add, addop));
 480                         }
 481                         Optype subop = (op == ADD) ? SUBTRACT : ADD;
 482                         for (auto sub : b->m_next[SUBTRACT])
 483                                 delete (OperationTerm(sub, subop));
 484
 485                         b->m_next[ADD].clear();
 486                         b->m_next[SUBTRACT].clear();
 487                         //assert(SanityCheck());
 488                         return b;
 489                 }
 490         }
 491
 492
 493
 494         bool parent = (merge_point == NULL);
 495         ParanoidNumber * merge = this;
 496         Optype mop = op;
 497         //assert(mop != NOP); // silence compiler warning
 498         if (parent)
 499         {
 500                 merge_point = &merge;
 501                 merge_op = &mop;
 502         }
 503         else
 504         {
 505                 merge = *merge_point;
 506                 mop = *merge_op;
 507         }
 508
 509         Optype invop = InverseOp(op); // inverse of p
 510         Optype fwd = op;
 511         Optype rev = invop;
 512         if (op == SUBTRACT)
 513         {
 514                 fwd = ADD;
 515                 rev = SUBTRACT;
 516         }
 517
 518         for (auto prev : m_next[invop])
 519         {
 520                 if (prev->OperationTerm(b, rev, merge_point, merge_op) == b)
 521                 {
 522                         //assert(SanityCheck());
 523                         return b;
 524                 }
 525
 526         }
 527         for (auto next : m_next[op])
 528         {
 529                 if (next->OperationTerm(b, fwd, merge_point, merge_op) == b)
 530                 {
 531                         //assert(SanityCheck());
 532                         return b;
 533                 }
 534         }
 535
 536
 537
 538
 539         if (parent)
 540         {
 541                 //merge->m_next[*merge_op].push_back(b);
 542                 m_next[op].push_back(b);
 543                 #ifdef PARANOID_SIZE_LIMIT
 544                         m_size += b->m_size;
 545                 #endif
 546         }
 547         else
 548         {
 549                 if (m_next[op].size() == 0)
 550                 {
 551                         *merge_point = this;
 552                         *merge_op = op;
 553                 }
 554         }
 555
 556         //assert(SanityCheck());
 557
 558         return NULL;
 559 }
 560
 561 ParanoidNumber * ParanoidNumber::OperationFactor(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
 562 {
 563         #ifdef PARANOID_CACHE_RESULTS
 564         m_cached_result = NAN;
 565         #endif
 566         #ifdef PARANOID_SIZE_LIMIT
 567                 if (m_size + b->m_size >= PARANOID_SIZE_LIMIT)
 568                 {
 569                         this->operator=(this->Digit());
 570                         if (op == MULTIPLY)
 571                                 m_value *= b->Digit();
 572                         else
 573                                 m_value /= b->Digit();
 574                         m_size = 1;
 575
 576                         //Debug("Cut off %p", this);
 577                         return b;
 578
 579                 }
 580         #endif
 581
 582         if (Floating() && m_value == 0)
 583         {
 584                 return b;
 585         }
 586
 587         if (Floating() && m_value == 1 && op == MULTIPLY)
 588         {
 589                 m_value = b->m_value;
 590                 for (int i = 0; i < NOP; ++i)
 591                 {
 592                         for (auto n : m_next[i])
 593                                 delete (n);
 594                         m_next[i].clear();
 595                         swap(m_next[i], b->m_next[i]);
 596                 }
 597                 //assert(SanityCheck());
 598                 return b;
 599         }
 600         if (b->Floating() && b->m_value == 1)
 601                 return b;
 602         if (b->Floating() && b->m_value == 0 && op == MULTIPLY)
 603         {
 604                 operator=(*b);
 605                 return b;
 606         }
 607
 608
 609         if (NoTerms() && b->NoTerms())
 610         {
 611                 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
 612                 {
 613                         Optype mulop = (op == MULTIPLY) ? MULTIPLY : DIVIDE;
 614                         for (auto mul : b->m_next[MULTIPLY])
 615                         {
 616                                 delete(OperationFactor(mul, mulop));
 617                         }
 618                         Optype divop = (op == MULTIPLY) ? DIVIDE : MULTIPLY;
 619                         for (auto div : b->m_next[DIVIDE])
 620                                 delete(OperationFactor(div, divop));
 621
 622                         b->m_next[DIVIDE].clear();
 623                         b->m_next[MULTIPLY].clear();
 624                         return b;
 625                 }
 626         }
 627
 628
 629         bool parent = (merge_point == NULL);
 630         ParanoidNumber * merge = this;
 631         Optype mop = op;
 632         if (parent)
 633         {
 634                 merge_point = &merge;
 635                 merge_op = &mop;
 636         }
 637         else
 638         {
 639                 merge = *merge_point;
 640                 mop = *merge_op;
 641         }
 642
 643         Optype invop = InverseOp(op); // inverse of p
 644         Optype fwd = op;
 645         Optype rev = invop;
 646         if (op == DIVIDE)
 647         {
 648                 fwd = MULTIPLY;
 649                 rev = DIVIDE;
 650         }
 651
 652         ParanoidNumber * cpy_b = new ParanoidNumber(*b);
 653         for (auto prev : m_next[invop])
 654         {
 655                 if (prev->OperationFactor(b, rev, merge_point, merge_op) == b)
 656                 {
 657                         for (auto add : m_next[ADD])
 658                                 delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
 659                         for (auto sub : m_next[SUBTRACT])
 660                                 delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
 661
 662                         delete(cpy_b);
 663                         return b;
 664                 }
 665         }
 666         for (auto next : m_next[op])
 667         {
 668                 if (next->OperationFactor(b, fwd, merge_point, merge_op) == b)
 669                 {
 670                         for (auto add : m_next[ADD])
 671                         {
 672                                 delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
 673                         }
 674                         for (auto sub : m_next[SUBTRACT])
 675                         {
 676                                 delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
 677                         }
 678                         delete(cpy_b);
 679                         return b;
 680                 }
 681         }
 682
 683         if (parent)
 684         {
 685                 //assert(b != NULL);
 686                 m_next[op].push_back(b);
 687                 for (auto add : m_next[ADD])
 688                         delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
 689                 for (auto sub : m_next[SUBTRACT])
 690                         delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
 691
 692                 #ifdef PARANOID_SIZE_LIMIT
 693                         m_size += b->m_size;
 694                 #endif
 695         }
 696         //assert(SanityCheck());
 697
 698
 699
 700         return NULL;
 701 }
 702
 703
 704
 705 /**
 706  * Performs the operation on a with argument b (a += b, a -= b, a *= b, a /= b)
 707  * @returns b if b can safely be deleted
 708  * @returns NULL if b has been merged with a
 709  * append indicates that b should be merged
 710  */
 711 ParanoidNumber * ParanoidNumber::Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
 712 {
 713
 714         if (b == NULL)
 715                 return NULL;
 716
 717
 718         if (op == SUBTRACT || op == ADD)
 719                 return OperationTerm(b, op, merge_point, merge_op);
 720         if (op == MULTIPLY || op == DIVIDE)
 721                 return OperationFactor(b, op, merge_point, merge_op);
 722         return b;
 723 }
 724
 725
 726
 727 string ParanoidNumber::PStr() const
 728 {
 729         stringstream s;
 730         for (int i = 0; i < NOP; ++i)
 731         {
 732                 Optype f = Optype(i);
 733                 s << this;
 734                 for (auto n : m_next[f])
 735                 {
 736                         s << OpChar(f) << n->PStr();
 737                 }
 738         }
 739         return s.str();
 740 }
 741
 742 bool ParanoidNumber::Simplify(Optype op)
 743 {
 744
 745         if (Floating())
 746                 return false;
 747
 748         //assert(SanityCheck());
 749         vector<ParanoidNumber*> next;
 750         next.clear();
 751         swap(m_next[op], next);
 752         m_next[op].clear();
 753         //assert(m_next[op].size() == 0);
 754         //assert(SanityCheck());
 755         Optype fwd = op;
 756         if (op == DIVIDE)
 757                 fwd = MULTIPLY;
 758         else if (op == SUBTRACT)
 759                 fwd = ADD;
 760
 761
 762         vector<ParanoidNumber*> hold[2];
 763         if (op == MULTIPLY || op == DIVIDE)
 764         {
 765                 swap(m_next[ADD], hold[0]);
 766                 swap(m_next[SUBTRACT], hold[1]);
 767         }
 768
 769         for (vector<ParanoidNumber*>::iterator n = next.begin(); n != next.end(); ++n)
 770         {
 771                 if (*n == NULL)
 772                         continue;
 773                 for (vector<ParanoidNumber*>::iterator m = n; m != next.end(); ++m)
 774                 {
 775                         if ((*m) == (*n))
 776                                 continue;
 777                         if (*m == NULL)
 778                                 continue;
 779
 780                         ParanoidNumber * parent = this;
 781                         Optype mop = op;
 782                         ParanoidNumber * result = (*n)->Operation(*m, fwd, &parent, &mop);
 783                         if (result != NULL)
 784                         {
 785                                 #ifdef PARANOID_SIZE_LIMIT
 786                                         m_size -= result->m_size;
 787                                 #endif
 788                                 *m = NULL;
 789                                 delete(result);
 790                         }
 791                 }
 792         }
 793
 794
 795
 796         for (auto n : next)
 797         {
 798                 if (n != NULL)
 799                 {
 800                         #ifdef PARANOID_SIZE_LIMIT
 801                                 if (Operation(n, op) == n)
 802                                 {
 803                                         m_size -= n->m_size;
 804                                         delete n;
 805                                 }
 806                         #else
 807                                 delete(Operation(n, op));
 808                         #endif
 809                 }
 810         }
 811
 812         if (op == MULTIPLY || op == DIVIDE)
 813         {
 814                 swap(m_next[ADD], hold[0]);
 815                 swap(m_next[SUBTRACT], hold[1]);
 816         }
 817
 818         set<ParanoidNumber*> s;
 819         //if (!SanityCheck(s))
 820         //{
 821         //      Error("Simplify broke Sanity");
 822         //}
 823         return (next.size() > m_next[op].size());
 824 }
 825
 826 bool ParanoidNumber::FullSimplify()
 827 {
 828         bool result = false;
 829         result |= Simplify(MULTIPLY);
 830         result |= Simplify(DIVIDE);
 831         result |= Simplify(ADD);
 832         result |= Simplify(SUBTRACT);
 833         return result;
 834 }
 835
 836 ParanoidNumber::digit_t ParanoidNumber::Digit() const
 837 {
 838
 839         // Get around the absurd requirement that const correctness be observed.
 840         #ifdef PARANOID_CACHE_RESULTS
 841         digit_t & result = ((ParanoidNumber*)(this))->m_cached_result;
 842
 843         if (!isnan(float(result))) // le sigh ambiguous function compiler warnings
 844                 return result;
 845         #else
 846                 digit_t result;
 847         #endif
 848         result = m_value;
 849         for (auto mul : m_next[MULTIPLY])
 850         {
 851                 result *= mul->Digit();
 852         }
 853         for (auto div : m_next[DIVIDE])
 854         {
 855                 result /= div->Digit();
 856         }
 857         for (auto add : m_next[ADD])
 858                 result += add->Digit();
 859         for (auto sub : m_next[SUBTRACT])
 860                 result -= sub->Digit();
 861         return result;
 862
 863 }
 864
 865 ParanoidNumber::digit_t ParanoidNumber::GetFactors() const
 866 {
 867         digit_t value = 1;
 868         for (auto mul : m_next[MULTIPLY])
 869                 value *= mul->Digit();
 870         for (auto div : m_next[DIVIDE])
 871                 value /= div->Digit();
 872         return value;
 873 }
 874
 875
 876 ParanoidNumber::digit_t ParanoidNumber::GetTerms() const
 877 {
 878         digit_t value = 0;
 879         for (auto add : m_next[ADD])
 880                 value += add->Digit();
 881         for (auto sub : m_next[SUBTRACT])
 882                 value -= sub->Digit();
 883         return value;
 884 }
 885
 886 bool ParanoidNumber::SanityCheck(set<ParanoidNumber*> & visited) const
 887 {
 888         if (this == NULL)
 889         {
 890                 Error("NULL pointer in tree");
 891                 return false;
 892         }
 893
 894         if (visited.find((ParanoidNumber*)this) != visited.end())
 895         {
 896                 Error("I think I've seen this tree before...");
 897                 return false;
 898         }
 899
 900         visited.insert((ParanoidNumber*)this);
 901
 902         for (auto add : m_next[ADD])
 903         {
 904                 if (!add->SanityCheck(visited))
 905                         return false;
 906         }
 907         for (auto sub : m_next[SUBTRACT])
 908         {
 909                 if (!sub->SanityCheck(visited))
 910                         return false;
 911         }
 912         for (auto mul : m_next[MULTIPLY])
 913         {
 914                 if (!mul->SanityCheck(visited))
 915                         return false;
 916         }
 917
 918         for (auto div : m_next[DIVIDE])
 919         {
 920                 if (!div->SanityCheck(visited))
 921                         return false;
 922                 if (div->Digit() == 0)
 923                 {
 924                         Error("Divide by zero");
 925                         return false;
 926                 }
 927         }
 928         return true;
 929 }
 930
 931 void ParanoidNumber::Negate()
 932 {
 933         swap(m_next[ADD], m_next[SUBTRACT]);
 934         m_value = -m_value;
 935 }
 936
 937 #ifdef PARANOID_USE_ARENA
 938
 939 void * ParanoidNumber::operator new(size_t s)
 940 {
 941         return g_arena.allocate(s);
 942 }
 943
 944 void ParanoidNumber::operator delete(void * p)
 945 {
 946         g_arena.deallocate(p);
 947 }
 948
 949 ParanoidNumber::Arena::Arena(int64_t block_size) : m_block_size(block_size), m_spare(NULL)
 950 {
 951         m_blocks.push_back({malloc(block_size*sizeof(ParanoidNumber)),0});
 952 }
 953
 954 ParanoidNumber::Arena::~Arena()
 955 {
 956         for (auto block : m_blocks)
 957         {
 958                 free(block.memory);
 959         }
 960         m_blocks.clear();
 961 }
 962
 963 void * ParanoidNumber::Arena::allocate(size_t s)
 964 {
 965         if (m_spare != NULL)
 966         {
 967                 void * result = m_spare;
 968                 m_spare = NULL;
 969                 return result;
 970         }
 971
 972         Block & b = m_blocks.back();
 973         void * result = (ParanoidNumber*)(b.memory) + (b.used++);
 974         if (b.used >= m_block_size)
 975         {
 976                 m_block_size *= 2;
 977                 Debug("Add block of size %d", m_block_size);
 978                 m_blocks.push_back({malloc(m_block_size*sizeof(ParanoidNumber)), 0});
 979         }
 980         return result;
 981 }
 982
 983 void ParanoidNumber::Arena::deallocate(void * p)
 984 {
 985         m_spare = p;
 986 }
 987 #endif //PARANOID_USE_ARENA
 988
 989 }