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