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