1 #include "paranoidnumber.h"
9 // here be many copy paste bugs
16 #ifdef PARANOID_USE_ARENA
17 ParanoidNumber::Arena ParanoidNumber::g_arena;
18 #endif //PARANOID_USE_ARENA
20 ParanoidNumber::~ParanoidNumber()
22 for (int i = 0; i < NOP; ++i)
24 for (auto n : m_next[i])
29 ParanoidNumber::ParanoidNumber(const string & str) : m_value(0), m_next()
31 #ifdef PARANOID_SIZE_LIMIT
34 #ifdef PARANOID_CACHE_RESULTS
35 m_cached_result = NAN;
40 bool negate = str[0] == '-';
46 while (str[dp] != '\0' && str[dp] != '.')
51 while (str[end] != '\0')
54 for (int i = dp-1; i >= negate; --i)
56 ParanoidNumber b(str[i]-'0');
62 for (int i = dp+1; i < end; ++i)
65 ParanoidNumber b(str[i]-'0');
73 #ifdef PARANOID_COMPARE_EPSILON
74 double d = strtod(str.c_str(), NULL);
75 CompareForSanity(d, d);
79 ParanoidNumber & ParanoidNumber::operator=(const ParanoidNumber & a)
81 //assert(this != NULL);
83 #ifdef PARANOID_SIZE_LIMIT
88 #ifdef PARANOID_CACHE_RESULT
89 m_cached_result = a.m_cached_result;
91 for (int i = 0; i < NOP; ++i)
93 for (auto n : m_next[i])
96 for (auto n : a.m_next[i])
97 m_next[i].push_back(new ParanoidNumber(*n));
99 #ifdef PARANOID_COMPARE_EPSILON
100 CompareForSanity(a.Digit(),a.Digit());
106 string ParanoidNumber::Str() const
109 //assert(this != NULL);
112 s << (double)m_value;
114 for (auto mul : m_next[MULTIPLY])
117 if (!mul->Floating())
118 result += "(" + mul->Str() + ")";
120 result += mul->Str();
122 for (auto div : m_next[DIVIDE])
125 if (!div->Floating())
126 result += "(" + div->Str() + ")";
128 result += div->Str();
131 for (auto add : m_next[ADD])
134 if (!add->Floating())
135 result += "(" + add->Str() + ")";
137 result += add->Str();
139 for (auto sub : m_next[SUBTRACT])
142 if (!sub->Floating())
143 result += "(" + sub->Str() + ")";
145 result += sub->Str();
153 bool TrustingOp<float>(float & a, const float & b, Optype op)
157 feclearexcept(FE_ALL_EXCEPT);
172 a = (a >= 0) ? INFINITY : -INFINITY;
181 return !fetestexcept(FE_ALL_EXCEPT);
185 bool TrustingOp<double>(double & a, const double & b, Optype op)
187 feclearexcept(FE_ALL_EXCEPT);
202 a = (a >= 0) ? INFINITY : -INFINITY;
210 return !fetestexcept(FE_ALL_EXCEPT);
216 bool TrustingOp<long double>(long double & a, const long double & b, Optype op)
220 feclearexcept(FE_ALL_EXCEPT);
235 a = (a >= 0) ? INFINITY : -INFINITY;
244 return !fetestexcept(FE_ALL_EXCEPT);
249 bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op)
257 exact = (abs(sa) <= 127);
261 exact = (abs(sa) <= 127);
265 exact = (abs(sa) <= 127);
268 exact = (b != 0 && sa > b && sa % b == 0);
279 ParanoidNumber & ParanoidNumber::operator+=(const digit_t & a)
281 #ifdef PARANOID_COMPARE_EPSILON
282 digit_t compare = Digit();
285 delete Operation(new ParanoidNumber(a), ADD);
288 #ifdef PARANOID_COMPARE_EPSILON
289 CompareForSanity(compare, a);
295 ParanoidNumber & ParanoidNumber::operator-=(const digit_t & a)
297 #ifdef PARANOID_COMPARE_EPSILON
298 digit_t compare = Digit();
301 delete Operation(new ParanoidNumber(a), SUBTRACT);
304 #ifdef PARANOID_COMPARE_EPSILON
305 CompareForSanity(compare, a);
310 ParanoidNumber & ParanoidNumber::operator*=(const digit_t & a)
312 #ifdef PARANOID_COMPARE_EPSILON
313 digit_t compare = Digit();
316 delete Operation(new ParanoidNumber(a), MULTIPLY);
319 #ifdef PARANOID_COMPARE_EPSILON
320 CompareForSanity(compare, a);
326 ParanoidNumber & ParanoidNumber::operator/=(const digit_t & a)
328 #ifdef PARANOID_COMPARE_EPSILON
329 digit_t compare = Digit();
332 delete Operation(new ParanoidNumber(a), DIVIDE);
335 #ifdef PARANOID_COMPARE_EPSILON
336 CompareForSanity(compare, a);
342 ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
344 #ifdef PARANOID_COMPARE_EPSILON
345 digit_t compare = Digit();
346 compare += a.Digit();
348 delete Operation(new ParanoidNumber(a), ADD);
351 #ifdef PARANOID_COMPARE_EPSILON
352 CompareForSanity(compare, a.Digit());
358 ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
360 #ifdef PARANOID_COMPARE_EPSILON
361 digit_t compare = Digit();
362 compare -= a.Digit();
364 delete Operation(new ParanoidNumber(a), SUBTRACT);
367 #ifdef PARANOID_COMPARE_EPSILON
368 CompareForSanity(compare, a.Digit());
373 ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
375 #ifdef PARANOID_COMPARE_EPSILON
376 digit_t compare = Digit();
377 compare *= a.Digit();
379 delete Operation(new ParanoidNumber(a), MULTIPLY);
382 #ifdef PARANOID_COMPARE_EPSILON
383 CompareForSanity(compare, a.Digit());
389 ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
391 #ifdef PARANOID_COMPARE_EPSILON
392 digit_t compare = Digit();
393 compare /= a.Digit();
395 delete Operation(new ParanoidNumber(a), DIVIDE);
398 #ifdef PARANOID_COMPARE_EPSILON
399 CompareForSanity(compare, a.Digit());
404 ParanoidNumber & ParanoidNumber::operator=(const digit_t & a)
407 for (int i = 0; i < NOP; ++i)
409 for (auto n : m_next[i])
414 #ifdef PARANOID_CACHE_RESULT
418 #ifdef PARANOID_COMPARE_EPSILON
419 CompareForSanity(a,a);
426 ParanoidNumber * ParanoidNumber::OperationTerm(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
428 ////assert(b->SanityCheck());
429 #ifdef PARANOID_CACHE_RESULTS
430 m_cached_result = NAN;
432 #ifdef PARANOID_SIZE_LIMIT
433 if (m_size >= PARANOID_SIZE_LIMIT)
435 this->operator=(this->Digit());
437 m_value += b->Digit();
439 m_value -= b->Digit();
441 Debug("Cut off %p", this);
444 //Debug("At size limit %d", m_size);
448 if (Floating() && m_value == 0) // 0 + b = b
450 m_value = b->m_value;
454 swap(b->m_next[ADD], b->m_next[SUBTRACT]);
457 for (int i = 0; i < NOP; ++i)
459 m_next[i] = b->m_next[i];
460 b->m_next[i].clear();
463 //assert(SanityCheck());
466 if (b->Floating() && b->m_value == 0) // a + 0 = a
471 if ((NoFactors() && b->NoFactors())
472 || (GetFactors() == b->GetFactors()))
474 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
476 Optype addop = (op == ADD) ? ADD : SUBTRACT;
477 for (auto add : b->m_next[ADD])
479 delete (OperationTerm(add, addop));
481 Optype subop = (op == ADD) ? SUBTRACT : ADD;
482 for (auto sub : b->m_next[SUBTRACT])
483 delete (OperationTerm(sub, subop));
485 b->m_next[ADD].clear();
486 b->m_next[SUBTRACT].clear();
487 //assert(SanityCheck());
494 bool parent = (merge_point == NULL);
495 ParanoidNumber * merge = this;
497 //assert(mop != NOP); // silence compiler warning
500 merge_point = &merge;
505 merge = *merge_point;
509 Optype invop = InverseOp(op); // inverse of p
518 for (auto prev : m_next[invop])
520 if (prev->OperationTerm(b, rev, merge_point, merge_op) == b)
522 //assert(SanityCheck());
527 for (auto next : m_next[op])
529 if (next->OperationTerm(b, fwd, merge_point, merge_op) == b)
531 //assert(SanityCheck());
541 //merge->m_next[*merge_op].push_back(b);
542 m_next[op].push_back(b);
543 #ifdef PARANOID_SIZE_LIMIT
544 m_size += 1+b->m_size;
549 if (m_next[op].size() == 0)
556 //assert(SanityCheck());
561 ParanoidNumber * ParanoidNumber::OperationFactor(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
563 #ifdef PARANOID_CACHE_RESULTS
564 m_cached_result = NAN;
566 #ifdef PARANOID_SIZE_LIMIT
567 if (m_size >= PARANOID_SIZE_LIMIT)
569 this->operator=(this->Digit());
571 m_value *= b->Digit();
573 m_value /= b->Digit();
576 Debug("Cut off %p", this);
582 if (Floating() && m_value == 0)
587 if (Floating() && m_value == 1 && op == MULTIPLY)
589 m_value = b->m_value;
590 for (int i = 0; i < NOP; ++i)
592 for (auto n : m_next[i])
595 swap(m_next[i], b->m_next[i]);
597 //assert(SanityCheck());
600 if (b->Floating() && b->m_value == 1)
602 if (b->Floating() && b->m_value == 0 && op == MULTIPLY)
609 if (NoTerms() && b->NoTerms())
611 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
613 Optype mulop = (op == MULTIPLY) ? MULTIPLY : DIVIDE;
614 for (auto mul : b->m_next[MULTIPLY])
616 delete(OperationFactor(mul, mulop));
618 Optype divop = (op == MULTIPLY) ? DIVIDE : MULTIPLY;
619 for (auto div : b->m_next[DIVIDE])
620 delete(OperationFactor(div, divop));
622 b->m_next[DIVIDE].clear();
623 b->m_next[MULTIPLY].clear();
629 bool parent = (merge_point == NULL);
630 ParanoidNumber * merge = this;
634 merge_point = &merge;
639 merge = *merge_point;
643 Optype invop = InverseOp(op); // inverse of p
652 ParanoidNumber * cpy_b = new ParanoidNumber(*b);
653 for (auto prev : m_next[invop])
655 if (prev->OperationFactor(b, rev, merge_point, merge_op) == b)
657 for (auto add : m_next[ADD])
658 delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
659 for (auto sub : m_next[SUBTRACT])
660 delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
666 for (auto next : m_next[op])
668 if (next->OperationFactor(b, fwd, merge_point, merge_op) == b)
670 for (auto add : m_next[ADD])
672 delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
674 for (auto sub : m_next[SUBTRACT])
676 delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
686 m_next[op].push_back(b);
687 for (auto add : m_next[ADD])
688 delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
689 for (auto sub : m_next[SUBTRACT])
690 delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
692 #ifdef PARANOID_SIZE_LIMIT
693 m_size += 1+b->m_size;
696 //assert(SanityCheck());
706 * Performs the operation on a with argument b (a += b, a -= b, a *= b, a /= b)
707 * @returns b if b can safely be deleted
708 * @returns NULL if b has been merged with a
709 * append indicates that b should be merged
711 ParanoidNumber * ParanoidNumber::Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
718 if (op == SUBTRACT || op == ADD)
719 return OperationTerm(b, op, merge_point, merge_op);
720 if (op == MULTIPLY || op == DIVIDE)
721 return OperationFactor(b, op, merge_point, merge_op);
727 string ParanoidNumber::PStr() const
730 for (int i = 0; i < NOP; ++i)
732 Optype f = Optype(i);
734 for (auto n : m_next[f])
736 s << OpChar(f) << n->PStr();
742 bool ParanoidNumber::Simplify(Optype op)
748 //assert(SanityCheck());
749 vector<ParanoidNumber*> next;
751 swap(m_next[op], next);
753 //assert(m_next[op].size() == 0);
754 //assert(SanityCheck());
758 else if (op == SUBTRACT)
762 vector<ParanoidNumber*> hold[2];
763 if (op == MULTIPLY || op == DIVIDE)
765 swap(m_next[ADD], hold[0]);
766 swap(m_next[SUBTRACT], hold[1]);
769 for (vector<ParanoidNumber*>::iterator n = next.begin(); n != next.end(); ++n)
773 for (vector<ParanoidNumber*>::iterator m = n; m != next.end(); ++m)
780 ParanoidNumber * parent = this;
782 ParanoidNumber * result = (*n)->Operation(*m, fwd, &parent, &mop);
785 #ifdef PARANOID_SIZE_LIMIT
786 m_size -= (1+result->m_size);
800 #ifdef PARANOID_SIZE_LIMIT
801 if (Operation(n, op) == n)
803 m_size -= (1+n->m_size);
807 delete(Operation(n, op));
812 if (op == MULTIPLY || op == DIVIDE)
814 swap(m_next[ADD], hold[0]);
815 swap(m_next[SUBTRACT], hold[1]);
818 set<ParanoidNumber*> s;
819 //if (!SanityCheck(s))
821 // Error("Simplify broke Sanity");
823 return (next.size() > m_next[op].size());
826 bool ParanoidNumber::FullSimplify()
829 result |= Simplify(MULTIPLY);
830 result |= Simplify(DIVIDE);
831 result |= Simplify(ADD);
832 result |= Simplify(SUBTRACT);
836 ParanoidNumber::digit_t ParanoidNumber::Digit() const
839 // Get around the absurd requirement that const correctness be observed.
840 #ifdef PARANOID_CACHE_RESULTS
841 digit_t & result = ((ParanoidNumber*)(this))->m_cached_result;
843 if (!isnan(float(result))) // le sigh ambiguous function compiler warnings
849 for (auto mul : m_next[MULTIPLY])
851 result *= mul->Digit();
853 for (auto div : m_next[DIVIDE])
855 result /= div->Digit();
857 for (auto add : m_next[ADD])
858 result += add->Digit();
859 for (auto sub : m_next[SUBTRACT])
860 result -= sub->Digit();
865 ParanoidNumber::digit_t ParanoidNumber::GetFactors() const
868 for (auto mul : m_next[MULTIPLY])
869 value *= mul->Digit();
870 for (auto div : m_next[DIVIDE])
871 value /= div->Digit();
876 ParanoidNumber::digit_t ParanoidNumber::GetTerms() const
879 for (auto add : m_next[ADD])
880 value += add->Digit();
881 for (auto sub : m_next[SUBTRACT])
882 value -= sub->Digit();
886 bool ParanoidNumber::SanityCheck(set<ParanoidNumber*> & visited) const
890 Error("NULL pointer in tree");
894 if (visited.find((ParanoidNumber*)this) != visited.end())
896 Error("I think I've seen this tree before...");
900 visited.insert((ParanoidNumber*)this);
902 for (auto add : m_next[ADD])
904 if (!add->SanityCheck(visited))
907 for (auto sub : m_next[SUBTRACT])
909 if (!sub->SanityCheck(visited))
912 for (auto mul : m_next[MULTIPLY])
914 if (!mul->SanityCheck(visited))
918 for (auto div : m_next[DIVIDE])
920 if (!div->SanityCheck(visited))
922 if (div->Digit() == 0)
924 Error("Divide by zero");
931 void ParanoidNumber::Negate()
933 swap(m_next[ADD], m_next[SUBTRACT]);
937 #ifdef PARANOID_USE_ARENA
939 void * ParanoidNumber::operator new(size_t s)
941 return g_arena.allocate(s);
944 void ParanoidNumber::operator delete(void * p)
946 g_arena.deallocate(p);
949 ParanoidNumber::Arena::Arena(int64_t block_size) : m_block_size(block_size), m_spare(NULL)
951 m_blocks.push_back({malloc(block_size*sizeof(ParanoidNumber)),0});
954 ParanoidNumber::Arena::~Arena()
956 for (auto block : m_blocks)
963 void * ParanoidNumber::Arena::allocate(size_t s)
967 void * result = m_spare;
972 Block & b = m_blocks.back();
973 void * result = (ParanoidNumber*)(b.memory) + (b.used++);
974 if (b.used >= m_block_size)
977 Debug("Add block of size %d", m_block_size);
978 m_blocks.push_back({malloc(m_block_size*sizeof(ParanoidNumber)), 0});
983 void ParanoidNumber::Arena::deallocate(void * p)
987 #endif //PARANOID_USE_ARENA