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 = 0;
  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 >= PARANOID_SIZE_LIMIT)
 434                 {
 435                         if (op == ADD)
 436                         {
 437                                 m_value += b->Digit() / GetFactors();
 438                         }
 439                         else
 440                         {
 441                                 m_value -= b->Digit() / GetFactors();
 442                         }
 443                         return b;
 444                 }
 445                 //Debug("At size limit %d", m_size);
 446         #endif
 447
 448
 449         if (Floating() && m_value == 0) // 0 + b = b
 450         {
 451                 m_value = b->m_value;
 452                 if (op == SUBTRACT)
 453                 {
 454                         m_value = -m_value;
 455                         swap(b->m_next[ADD], b->m_next[SUBTRACT]);
 456                 }
 457
 458                 for (int i = 0; i < NOP; ++i)
 459                 {
 460                         m_next[i] = b->m_next[i];
 461                         b->m_next[i].clear();
 462                 }
 463
 464                 //assert(SanityCheck());
 465                 return b;
 466         }
 467         if (b->Floating() && b->m_value == 0) // a + 0 = a
 468                 return b;
 469
 470
 471
 472         if ((NoFactors() && b->NoFactors())
 473                 || (GetFactors() == b->GetFactors()))
 474         {
 475                 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
 476                 {
 477                         Optype addop = (op == ADD) ? ADD : SUBTRACT;
 478                         for (auto add : b->m_next[ADD])
 479                         {
 480                                 delete (OperationTerm(add, addop));
 481                         }
 482                         Optype subop = (op == ADD) ? SUBTRACT : ADD;
 483                         for (auto sub : b->m_next[SUBTRACT])
 484                                 delete (OperationTerm(sub, subop));
 485
 486                         b->m_next[ADD].clear();
 487                         b->m_next[SUBTRACT].clear();
 488                         //assert(SanityCheck());
 489                         return b;
 490                 }
 491         }
 492
 493
 494
 495         bool parent = (merge_point == NULL);
 496         ParanoidNumber * merge = this;
 497         Optype mop = op;
 498         //assert(mop != NOP); // silence compiler warning
 499         if (parent)
 500         {
 501                 merge_point = &merge;
 502                 merge_op = &mop;
 503         }
 504         else
 505         {
 506                 merge = *merge_point;
 507                 mop = *merge_op;
 508         }
 509
 510         Optype invop = InverseOp(op); // inverse of p
 511         Optype fwd = op;
 512         Optype rev = invop;
 513         if (op == SUBTRACT)
 514         {
 515                 fwd = ADD;
 516                 rev = SUBTRACT;
 517         }
 518
 519         for (auto prev : m_next[invop])
 520         {
 521                 if (prev->OperationTerm(b, rev, merge_point, merge_op) == b)
 522                 {
 523                         //assert(SanityCheck());
 524                         return b;
 525                 }
 526
 527         }
 528         for (auto next : m_next[op])
 529         {
 530                 if (next->OperationTerm(b, fwd, merge_point, merge_op) == b)
 531                 {
 532                         //assert(SanityCheck());
 533                         return b;
 534                 }
 535         }
 536
 537
 538
 539
 540         if (parent)
 541         {
 542                 //merge->m_next[*merge_op].push_back(b);
 543                 m_next[op].push_back(b);
 544                 #ifdef PARANOID_SIZE_LIMIT
 545                         m_size += 1+b->m_size;
 546                 #endif
 547         }
 548         else
 549         {
 550                 if (m_next[op].size() == 0)
 551                 {
 552                         *merge_point = this;
 553                         *merge_op = op;
 554                 }
 555         }
 556
 557         //assert(SanityCheck());
 558
 559         return NULL;
 560 }
 561
 562 ParanoidNumber * ParanoidNumber::OperationFactor(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
 563 {
 564         #ifdef PARANOID_CACHE_RESULTS
 565         m_cached_result = NAN;
 566         #endif
 567         #ifdef PARANOID_SIZE_LIMIT
 568                 if (m_size >= PARANOID_SIZE_LIMIT)
 569                 {
 570                         if (op == MULTIPLY)
 571                                 m_value *= b->Digit();
 572                         else
 573                                 m_value /= b->Digit();
 574
 575                         for (auto n : m_next[ADD])
 576                                 delete n->OperationFactor(new ParanoidNumber(*b), op);
 577                         for (auto n : m_next[SUBTRACT])
 578                                 delete n->OperationFactor(new ParanoidNumber(*b), op);
 579                         return b;
 580
 581                 }
 582         #endif
 583
 584         if (Floating() && m_value == 0)
 585         {
 586                 return b;
 587         }
 588
 589         if (Floating() && m_value == 1 && op == MULTIPLY)
 590         {
 591                 m_value = b->m_value;
 592                 for (int i = 0; i < NOP; ++i)
 593                 {
 594                         for (auto n : m_next[i])
 595                                 delete (n);
 596                         m_next[i].clear();
 597                         swap(m_next[i], b->m_next[i]);
 598                 }
 599                 //assert(SanityCheck());
 600                 return b;
 601         }
 602         if (b->Floating() && b->m_value == 1)
 603                 return b;
 604         if (b->Floating() && b->m_value == 0 && op == MULTIPLY)
 605         {
 606                 operator=(*b);
 607                 return b;
 608         }
 609
 610
 611         if (NoTerms() && b->NoTerms())
 612         {
 613                 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
 614                 {
 615                         Optype mulop = (op == MULTIPLY) ? MULTIPLY : DIVIDE;
 616                         for (auto mul : b->m_next[MULTIPLY])
 617                         {
 618                                 delete(OperationFactor(mul, mulop));
 619                         }
 620                         Optype divop = (op == MULTIPLY) ? DIVIDE : MULTIPLY;
 621                         for (auto div : b->m_next[DIVIDE])
 622                                 delete(OperationFactor(div, divop));
 623
 624                         b->m_next[DIVIDE].clear();
 625                         b->m_next[MULTIPLY].clear();
 626                         return b;
 627                 }
 628         }
 629
 630
 631         bool parent = (merge_point == NULL);
 632         ParanoidNumber * merge = this;
 633         Optype mop = op;
 634         if (parent)
 635         {
 636                 merge_point = &merge;
 637                 merge_op = &mop;
 638         }
 639         else
 640         {
 641                 merge = *merge_point;
 642                 mop = *merge_op;
 643         }
 644
 645         Optype invop = InverseOp(op); // inverse of p
 646         Optype fwd = op;
 647         Optype rev = invop;
 648         if (op == DIVIDE)
 649         {
 650                 fwd = MULTIPLY;
 651                 rev = DIVIDE;
 652         }
 653
 654         ParanoidNumber * cpy_b = new ParanoidNumber(*b);
 655         for (auto prev : m_next[invop])
 656         {
 657                 if (prev->OperationFactor(b, rev, merge_point, merge_op) == b)
 658                 {
 659                         for (auto add : m_next[ADD])
 660                                 delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
 661                         for (auto sub : m_next[SUBTRACT])
 662                                 delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
 663
 664                         delete(cpy_b);
 665                         return b;
 666                 }
 667         }
 668         for (auto next : m_next[op])
 669         {
 670                 if (next->OperationFactor(b, fwd, merge_point, merge_op) == b)
 671                 {
 672                         for (auto add : m_next[ADD])
 673                         {
 674                                 delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
 675                         }
 676                         for (auto sub : m_next[SUBTRACT])
 677                         {
 678                                 delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
 679                         }
 680                         delete(cpy_b);
 681                         return b;
 682                 }
 683         }
 684
 685         if (parent)
 686         {
 687                 //assert(b != NULL);
 688                 m_next[op].push_back(b);
 689                 for (auto add : m_next[ADD])
 690                         delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
 691                 for (auto sub : m_next[SUBTRACT])
 692                         delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
 693
 694                 #ifdef PARANOID_SIZE_LIMIT
 695                         m_size += 1+b->m_size;
 696                 #endif
 697         }
 698         //assert(SanityCheck());
 699
 700
 701
 702         return NULL;
 703 }
 704
 705
 706
 707 /**
 708  * Performs the operation on a with argument b (a += b, a -= b, a *= b, a /= b)
 709  * @returns b if b can safely be deleted
 710  * @returns NULL if b has been merged with a
 711  * append indicates that b should be merged
 712  */
 713 ParanoidNumber * ParanoidNumber::Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
 714 {
 715
 716         if (b == NULL)
 717                 return NULL;
 718
 719
 720         if (op == SUBTRACT || op == ADD)
 721                 return OperationTerm(b, op, merge_point, merge_op);
 722         if (op == MULTIPLY || op == DIVIDE)
 723                 return OperationFactor(b, op, merge_point, merge_op);
 724         return b;
 725 }
 726
 727
 728
 729 string ParanoidNumber::PStr() const
 730 {
 731         stringstream s;
 732         for (int i = 0; i < NOP; ++i)
 733         {
 734                 Optype f = Optype(i);
 735                 s << this;
 736                 for (auto n : m_next[f])
 737                 {
 738                         s << OpChar(f) << n->PStr();
 739                 }
 740         }
 741         return s.str();
 742 }
 743
 744 bool ParanoidNumber::Simplify(Optype op)
 745 {
 746
 747         if (Floating())
 748                 return false;
 749
 750         //assert(SanityCheck());
 751         vector<ParanoidNumber*> next;
 752         next.clear();
 753         swap(m_next[op], next);
 754         m_next[op].clear();
 755         //assert(m_next[op].size() == 0);
 756         //assert(SanityCheck());
 757         Optype fwd = op;
 758         if (op == DIVIDE)
 759                 fwd = MULTIPLY;
 760         else if (op == SUBTRACT)
 761                 fwd = ADD;
 762
 763
 764         vector<ParanoidNumber*> hold[2];
 765         if (op == MULTIPLY || op == DIVIDE)
 766         {
 767                 swap(m_next[ADD], hold[0]);
 768                 swap(m_next[SUBTRACT], hold[1]);
 769         }
 770
 771         for (vector<ParanoidNumber*>::iterator n = next.begin(); n != next.end(); ++n)
 772         {
 773                 if (*n == NULL)
 774                         continue;
 775                 for (vector<ParanoidNumber*>::iterator m = n; m != next.end(); ++m)
 776                 {
 777                         if ((*m) == (*n))
 778                                 continue;
 779                         if (*m == NULL)
 780                                 continue;
 781
 782                         ParanoidNumber * parent = this;
 783                         Optype mop = op;
 784                         ParanoidNumber * result = (*n)->Operation(*m, fwd, &parent, &mop);
 785                         if (result != NULL)
 786                         {
 787                                 #ifdef PARANOID_SIZE_LIMIT
 788                                         m_size -= (1+result->m_size);
 789                                 #endif
 790                                 *m = NULL;
 791                                 delete(result);
 792                         }
 793                 }
 794         }
 795
 796
 797
 798         for (auto n : next)
 799         {
 800                 if (n != NULL)
 801                 {
 802                         #ifdef PARANOID_SIZE_LIMIT
 803                                 if (Operation(n, op) == n)
 804                                 {
 805                                         m_size -= (1+n->m_size);
 806                                         delete n;
 807                                 }
 808                         #else
 809                                 delete(Operation(n, op));
 810                         #endif
 811                 }
 812         }
 813
 814         if (op == MULTIPLY || op == DIVIDE)
 815         {
 816                 swap(m_next[ADD], hold[0]);
 817                 swap(m_next[SUBTRACT], hold[1]);
 818         }
 819
 820         set<ParanoidNumber*> s;
 821         //if (!SanityCheck(s))
 822         //{
 823         //      Error("Simplify broke Sanity");
 824         //}
 825         return (next.size() > m_next[op].size());
 826 }
 827
 828 bool ParanoidNumber::FullSimplify()
 829 {
 830         bool result = false;
 831         result |= Simplify(MULTIPLY);
 832         result |= Simplify(DIVIDE);
 833         result |= Simplify(ADD);
 834         result |= Simplify(SUBTRACT);
 835         return result;
 836 }
 837
 838 ParanoidNumber::digit_t ParanoidNumber::Digit() const
 839 {
 840
 841         // Get around the absurd requirement that const correctness be observed.
 842         #ifdef PARANOID_CACHE_RESULTS
 843         digit_t & result = ((ParanoidNumber*)(this))->m_cached_result;
 844
 845         if (!isnan(float(result))) // le sigh ambiguous function compiler warnings
 846                 return result;
 847         #else
 848                 digit_t result;
 849         #endif
 850         result = m_value;
 851         for (auto mul : m_next[MULTIPLY])
 852         {
 853                 result *= mul->Digit();
 854         }
 855         for (auto div : m_next[DIVIDE])
 856         {
 857                 result /= div->Digit();
 858         }
 859         for (auto add : m_next[ADD])
 860                 result += add->Digit();
 861         for (auto sub : m_next[SUBTRACT])
 862                 result -= sub->Digit();
 863         return result;
 864
 865 }
 866
 867 ParanoidNumber::digit_t ParanoidNumber::GetFactors() const
 868 {
 869         digit_t value = 1;
 870         for (auto mul : m_next[MULTIPLY])
 871                 value *= mul->Digit();
 872         for (auto div : m_next[DIVIDE])
 873                 value /= div->Digit();
 874         return value;
 875 }
 876
 877
 878 ParanoidNumber::digit_t ParanoidNumber::GetTerms() const
 879 {
 880         digit_t value = 0;
 881         for (auto add : m_next[ADD])
 882                 value += add->Digit();
 883         for (auto sub : m_next[SUBTRACT])
 884                 value -= sub->Digit();
 885         return value;
 886 }
 887
 888 bool ParanoidNumber::SanityCheck(set<ParanoidNumber*> & visited) const
 889 {
 890         if (this == NULL)
 891         {
 892                 Error("NULL pointer in tree");
 893                 return false;
 894         }
 895
 896         if (visited.find((ParanoidNumber*)this) != visited.end())
 897         {
 898                 Error("I think I've seen this tree before...");
 899                 return false;
 900         }
 901
 902         visited.insert((ParanoidNumber*)this);
 903
 904         for (auto add : m_next[ADD])
 905         {
 906                 if (!add->SanityCheck(visited))
 907                         return false;
 908         }
 909         for (auto sub : m_next[SUBTRACT])
 910         {
 911                 if (!sub->SanityCheck(visited))
 912                         return false;
 913         }
 914         for (auto mul : m_next[MULTIPLY])
 915         {
 916                 if (!mul->SanityCheck(visited))
 917                         return false;
 918         }
 919
 920         for (auto div : m_next[DIVIDE])
 921         {
 922                 if (!div->SanityCheck(visited))
 923                         return false;
 924                 if (div->Digit() == 0)
 925                 {
 926                         Error("Divide by zero");
 927                         return false;
 928                 }
 929         }
 930         return true;
 931 }
 932
 933 void ParanoidNumber::Negate()
 934 {
 935         swap(m_next[ADD], m_next[SUBTRACT]);
 936         m_value = -m_value;
 937 }
 938
 939 #ifdef PARANOID_USE_ARENA
 940
 941 void * ParanoidNumber::operator new(size_t s)
 942 {
 943         return g_arena.allocate(s);
 944 }
 945
 946 void ParanoidNumber::operator delete(void * p)
 947 {
 948         g_arena.deallocate(p);
 949 }
 950
 951 ParanoidNumber::Arena::Arena(int64_t block_size) : m_block_size(block_size), m_spare(NULL)
 952 {
 953         m_blocks.push_back({malloc(block_size*sizeof(ParanoidNumber)),0});
 954 }
 955
 956 ParanoidNumber::Arena::~Arena()
 957 {
 958         for (auto block : m_blocks)
 959         {
 960                 free(block.memory);
 961         }
 962         m_blocks.clear();
 963 }
 964
 965 void * ParanoidNumber::Arena::allocate(size_t s)
 966 {
 967         if (m_spare != NULL)
 968         {
 969                 void * result = m_spare;
 970                 m_spare = NULL;
 971                 return result;
 972         }
 973
 974         Block & b = m_blocks.back();
 975         void * result = (ParanoidNumber*)(b.memory) + (b.used++);
 976         if (b.used >= m_block_size)
 977         {
 978                 m_block_size *= 2;
 979                 Debug("Add block of size %d", m_block_size);
 980                 m_blocks.push_back({malloc(m_block_size*sizeof(ParanoidNumber)), 0});
 981         }
 982         return result;
 983 }
 984
 985 void ParanoidNumber::Arena::deallocate(void * p)
 986 {
 987         m_spare = p;
 988 }
 989 #endif //PARANOID_USE_ARENA
 990
 991 }