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_cache_valid = false;
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_RESULTS
89 m_cached_result = a.m_cached_result;
90 m_cache_valid = a.m_cache_valid;
92 for (int i = 0; i < NOP; ++i)
94 for (auto n : m_next[i])
97 for (auto n : a.m_next[i])
98 m_next[i].push_back(new ParanoidNumber(*n));
100 #ifdef PARANOID_COMPARE_EPSILON
101 CompareForSanity(a.Digit(),a.Digit());
107 string ParanoidNumber::Str() const
110 //assert(this != NULL);
113 s << (double)m_value;
115 for (auto mul : m_next[MULTIPLY])
118 if (!mul->Floating())
119 result += "(" + mul->Str() + ")";
121 result += mul->Str();
123 for (auto div : m_next[DIVIDE])
126 if (!div->Floating())
127 result += "(" + div->Str() + ")";
129 result += div->Str();
132 for (auto add : m_next[ADD])
135 if (!add->Floating())
136 result += "(" + add->Str() + ")";
138 result += add->Str();
140 for (auto sub : m_next[SUBTRACT])
143 if (!sub->Floating())
144 result += "(" + sub->Str() + ")";
146 result += sub->Str();
154 bool TrustingOp<float>(float & a, const float & b, Optype op)
158 feclearexcept(FE_ALL_EXCEPT);
173 a = (a >= 0) ? INFINITY : -INFINITY;
182 return !fetestexcept(FE_ALL_EXCEPT);
186 bool TrustingOp<double>(double & a, const double & b, Optype op)
188 feclearexcept(FE_ALL_EXCEPT);
203 a = (a >= 0) ? INFINITY : -INFINITY;
211 return !fetestexcept(FE_ALL_EXCEPT);
217 bool TrustingOp<long double>(long double & a, const long double & b, Optype op)
221 feclearexcept(FE_ALL_EXCEPT);
236 a = (a >= 0) ? INFINITY : -INFINITY;
245 return !fetestexcept(FE_ALL_EXCEPT);
250 bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op)
258 exact = (abs(sa) <= 127);
262 exact = (abs(sa) <= 127);
266 exact = (abs(sa) <= 127);
269 exact = (b != 0 && sa > b && sa % b == 0);
280 ParanoidNumber & ParanoidNumber::operator+=(const digit_t & a)
282 #ifdef PARANOID_COMPARE_EPSILON
283 digit_t compare = Digit();
286 delete Operation(new ParanoidNumber(a), ADD);
289 #ifdef PARANOID_COMPARE_EPSILON
290 CompareForSanity(compare, a);
296 ParanoidNumber & ParanoidNumber::operator-=(const digit_t & a)
298 #ifdef PARANOID_COMPARE_EPSILON
299 digit_t compare = Digit();
302 delete Operation(new ParanoidNumber(a), SUBTRACT);
305 #ifdef PARANOID_COMPARE_EPSILON
306 CompareForSanity(compare, a);
311 ParanoidNumber & ParanoidNumber::operator*=(const digit_t & a)
313 #ifdef PARANOID_COMPARE_EPSILON
314 digit_t compare = Digit();
317 delete Operation(new ParanoidNumber(a), MULTIPLY);
320 #ifdef PARANOID_COMPARE_EPSILON
321 CompareForSanity(compare, a);
327 ParanoidNumber & ParanoidNumber::operator/=(const digit_t & a)
329 #ifdef PARANOID_COMPARE_EPSILON
330 digit_t compare = Digit();
333 delete Operation(new ParanoidNumber(a), DIVIDE);
336 #ifdef PARANOID_COMPARE_EPSILON
337 CompareForSanity(compare, a);
343 ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
345 #ifdef PARANOID_COMPARE_EPSILON
346 digit_t compare = Digit();
347 compare += a.Digit();
349 delete Operation(new ParanoidNumber(a), ADD);
352 #ifdef PARANOID_COMPARE_EPSILON
353 CompareForSanity(compare, a.Digit());
359 ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
361 #ifdef PARANOID_COMPARE_EPSILON
362 digit_t compare = Digit();
363 compare -= a.Digit();
365 delete Operation(new ParanoidNumber(a), SUBTRACT);
368 #ifdef PARANOID_COMPARE_EPSILON
369 CompareForSanity(compare, a.Digit());
374 ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
376 #ifdef PARANOID_COMPARE_EPSILON
377 digit_t compare = Digit();
378 compare *= a.Digit();
380 delete Operation(new ParanoidNumber(a), MULTIPLY);
383 #ifdef PARANOID_COMPARE_EPSILON
384 CompareForSanity(compare, a.Digit());
390 ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
392 #ifdef PARANOID_COMPARE_EPSILON
393 digit_t compare = Digit();
394 compare /= a.Digit();
396 delete Operation(new ParanoidNumber(a), DIVIDE);
399 #ifdef PARANOID_COMPARE_EPSILON
400 CompareForSanity(compare, a.Digit());
405 ParanoidNumber & ParanoidNumber::operator=(const digit_t & a)
408 for (int i = 0; i < NOP; ++i)
410 for (auto n : m_next[i])
415 #ifdef PARANOID_CACHE_RESULTS
417 m_cache_valid = true;
420 #ifdef PARANOID_COMPARE_EPSILON
421 CompareForSanity(a,a);
428 ParanoidNumber * ParanoidNumber::OperationTerm(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
430 ////assert(b->SanityCheck());
431 #ifdef PARANOID_CACHE_RESULTS
432 m_cache_valid = false;
434 #ifdef PARANOID_SIZE_LIMIT
435 if (m_size + b->m_size >= PARANOID_SIZE_LIMIT)
437 this->operator=(this->Digit());
439 m_value += b->Digit();
441 m_value -= b->Digit();
443 #ifdef PARANOID_CACHE_RESULTS
444 m_cached_result = m_value;
445 m_cache_valid = true;
449 //Debug("At size limit %d", m_size);
453 if (Floating() && m_value == 0) // 0 + b = b
455 m_value = b->m_value;
459 swap(b->m_next[ADD], b->m_next[SUBTRACT]);
462 for (int i = 0; i < NOP; ++i)
464 m_next[i] = b->m_next[i];
465 b->m_next[i].clear();
468 //assert(SanityCheck());
471 if (b->Floating() && b->m_value == 0) // a + 0 = a
476 if ((NoFactors() && b->NoFactors())
477 || (GetFactors() == b->GetFactors()))
479 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
481 Optype addop = (op == ADD) ? ADD : SUBTRACT;
482 for (auto add : b->m_next[ADD])
484 delete (OperationTerm(add, addop));
486 Optype subop = (op == ADD) ? SUBTRACT : ADD;
487 for (auto sub : b->m_next[SUBTRACT])
488 delete (OperationTerm(sub, subop));
490 b->m_next[ADD].clear();
491 b->m_next[SUBTRACT].clear();
492 //assert(SanityCheck());
499 bool parent = (merge_point == NULL);
500 ParanoidNumber * merge = this;
502 //assert(mop != NOP); // silence compiler warning
505 merge_point = &merge;
510 merge = *merge_point;
514 Optype invop = InverseOp(op); // inverse of p
523 for (auto prev : m_next[invop])
525 if (prev->OperationTerm(b, rev, merge_point, merge_op) == b)
527 //assert(SanityCheck());
532 for (auto next : m_next[op])
534 if (next->OperationTerm(b, fwd, merge_point, merge_op) == b)
536 //assert(SanityCheck());
546 //merge->m_next[*merge_op].push_back(b);
547 m_next[op].push_back(b);
548 #ifdef PARANOID_SIZE_LIMIT
554 if (m_next[op].size() == 0)
561 //assert(SanityCheck());
566 ParanoidNumber * ParanoidNumber::OperationFactor(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
568 #ifdef PARANOID_CACHE_RESULTS
569 m_cache_valid = false;
571 #ifdef PARANOID_SIZE_LIMIT
572 if (m_size + b->m_size >= PARANOID_SIZE_LIMIT)
574 this->operator=(this->Digit());
576 m_value *= b->Digit();
578 m_value /= b->Digit();
580 #ifdef PARANOID_CACHE_RESULTS
581 m_cached_result = m_value;
582 m_cache_valid = true;
584 //Debug("Cut off %p", this);
590 if (Floating() && m_value == 0)
595 if (Floating() && m_value == 1 && op == MULTIPLY)
597 m_value = b->m_value;
598 for (int i = 0; i < NOP; ++i)
600 for (auto n : m_next[i])
603 swap(m_next[i], b->m_next[i]);
605 //assert(SanityCheck());
608 if (b->Floating() && b->m_value == 1)
610 if (b->Floating() && b->m_value == 0 && op == MULTIPLY)
617 if (NoTerms() && b->NoTerms())
619 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
621 Optype mulop = (op == MULTIPLY) ? MULTIPLY : DIVIDE;
622 for (auto mul : b->m_next[MULTIPLY])
624 delete(OperationFactor(mul, mulop));
626 Optype divop = (op == MULTIPLY) ? DIVIDE : MULTIPLY;
627 for (auto div : b->m_next[DIVIDE])
628 delete(OperationFactor(div, divop));
630 b->m_next[DIVIDE].clear();
631 b->m_next[MULTIPLY].clear();
637 bool parent = (merge_point == NULL);
638 ParanoidNumber * merge = this;
642 merge_point = &merge;
647 merge = *merge_point;
651 Optype invop = InverseOp(op); // inverse of p
660 ParanoidNumber * cpy_b = new ParanoidNumber(*b);
661 for (auto prev : m_next[invop])
663 if (prev->OperationFactor(b, rev, merge_point, merge_op) == b)
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));
674 for (auto next : m_next[op])
676 if (next->OperationFactor(b, fwd, merge_point, merge_op) == b)
678 for (auto add : m_next[ADD])
680 delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
682 for (auto sub : m_next[SUBTRACT])
684 delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
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));
700 #ifdef PARANOID_SIZE_LIMIT
704 //assert(SanityCheck());
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
719 ParanoidNumber * ParanoidNumber::Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
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);
735 string ParanoidNumber::PStr() const
738 for (int i = 0; i < NOP; ++i)
740 Optype f = Optype(i);
742 for (auto n : m_next[f])
744 s << OpChar(f) << n->PStr();
750 bool ParanoidNumber::Simplify(Optype op)
756 //assert(SanityCheck());
757 vector<ParanoidNumber*> next;
759 swap(m_next[op], next);
761 //assert(m_next[op].size() == 0);
762 //assert(SanityCheck());
766 else if (op == SUBTRACT)
770 vector<ParanoidNumber*> hold[2];
771 if (op == MULTIPLY || op == DIVIDE)
773 swap(m_next[ADD], hold[0]);
774 swap(m_next[SUBTRACT], hold[1]);
777 for (vector<ParanoidNumber*>::iterator n = next.begin(); n != next.end(); ++n)
781 for (vector<ParanoidNumber*>::iterator m = n; m != next.end(); ++m)
788 ParanoidNumber * parent = this;
790 ParanoidNumber * result = (*n)->Operation(*m, fwd, &parent, &mop);
793 #ifdef PARANOID_SIZE_LIMIT
794 m_size -= result->m_size;
808 #ifdef PARANOID_SIZE_LIMIT
809 if (Operation(n, op) == n)
815 delete(Operation(n, op));
820 if (op == MULTIPLY || op == DIVIDE)
822 swap(m_next[ADD], hold[0]);
823 swap(m_next[SUBTRACT], hold[1]);
826 set<ParanoidNumber*> s;
827 //if (!SanityCheck(s))
829 // Error("Simplify broke Sanity");
831 return (next.size() > m_next[op].size());
834 bool ParanoidNumber::FullSimplify()
837 result |= Simplify(MULTIPLY);
838 result |= Simplify(DIVIDE);
839 result |= Simplify(ADD);
840 result |= Simplify(SUBTRACT);
844 ParanoidNumber::digit_t ParanoidNumber::Digit() const
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;
852 digit_t & result = ((ParanoidNumber*)(this))->m_cached_result;
858 for (auto mul : m_next[MULTIPLY])
860 result *= mul->Digit();
862 for (auto div : m_next[DIVIDE])
864 result /= div->Digit();
866 for (auto add : m_next[ADD])
867 result += add->Digit();
868 for (auto sub : m_next[SUBTRACT])
869 result -= sub->Digit();
871 #ifdef PARANOID_CACHE_RESULTS
872 ((ParanoidNumber*)(this))->m_cache_valid = true;
878 ParanoidNumber::digit_t ParanoidNumber::GetFactors() const
881 for (auto mul : m_next[MULTIPLY])
882 value *= mul->Digit();
883 for (auto div : m_next[DIVIDE])
884 value /= div->Digit();
889 ParanoidNumber::digit_t ParanoidNumber::GetTerms() const
892 for (auto add : m_next[ADD])
893 value += add->Digit();
894 for (auto sub : m_next[SUBTRACT])
895 value -= sub->Digit();
899 bool ParanoidNumber::SanityCheck(set<ParanoidNumber*> & visited) const
903 Error("NULL pointer in tree");
907 if (visited.find((ParanoidNumber*)this) != visited.end())
909 Error("I think I've seen this tree before...");
913 visited.insert((ParanoidNumber*)this);
915 for (auto add : m_next[ADD])
917 if (!add->SanityCheck(visited))
920 for (auto sub : m_next[SUBTRACT])
922 if (!sub->SanityCheck(visited))
925 for (auto mul : m_next[MULTIPLY])
927 if (!mul->SanityCheck(visited))
931 for (auto div : m_next[DIVIDE])
933 if (!div->SanityCheck(visited))
935 if (div->Digit() == 0)
937 Error("Divide by zero");
944 void ParanoidNumber::Negate()
946 swap(m_next[ADD], m_next[SUBTRACT]);
948 #ifdef PARANOID_CACHE_RESULTS
949 m_cached_result = -m_cached_result;
953 #ifdef PARANOID_USE_ARENA
955 void * ParanoidNumber::operator new(size_t s)
957 return g_arena.allocate(s);
960 void ParanoidNumber::operator delete(void * p)
962 g_arena.deallocate(p);
965 ParanoidNumber::Arena::Arena(int64_t block_size) : m_block_size(block_size), m_spare(NULL)
967 m_blocks.push_back({malloc(block_size*sizeof(ParanoidNumber)),0});
970 ParanoidNumber::Arena::~Arena()
972 for (auto block : m_blocks)
979 void * ParanoidNumber::Arena::allocate(size_t s)
983 void * result = m_spare;
988 Block & b = m_blocks.back();
989 void * result = (ParanoidNumber*)(b.memory) + (b.used++);
990 if (b.used >= m_block_size)
993 Debug("Add block of size %d", m_block_size);
994 m_blocks.push_back({malloc(m_block_size*sizeof(ParanoidNumber)), 0});
999 void ParanoidNumber::Arena::deallocate(void * p)
1003 #endif //PARANOID_USE_ARENA