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)
437 m_value += b->Digit() / GetFactors();
441 m_value -= b->Digit() / GetFactors();
445 //Debug("At size limit %d", m_size);
449 if (Floating() && m_value == 0) // 0 + b = b
451 m_value = b->m_value;
455 swap(b->m_next[ADD], b->m_next[SUBTRACT]);
458 for (int i = 0; i < NOP; ++i)
460 m_next[i] = b->m_next[i];
461 b->m_next[i].clear();
464 //assert(SanityCheck());
467 if (b->Floating() && b->m_value == 0) // a + 0 = a
472 if ((NoFactors() && b->NoFactors())
473 || (GetFactors() == b->GetFactors()))
475 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
477 Optype addop = (op == ADD) ? ADD : SUBTRACT;
478 for (auto add : b->m_next[ADD])
480 delete (OperationTerm(add, addop));
482 Optype subop = (op == ADD) ? SUBTRACT : ADD;
483 for (auto sub : b->m_next[SUBTRACT])
484 delete (OperationTerm(sub, subop));
486 b->m_next[ADD].clear();
487 b->m_next[SUBTRACT].clear();
488 //assert(SanityCheck());
495 bool parent = (merge_point == NULL);
496 ParanoidNumber * merge = this;
498 //assert(mop != NOP); // silence compiler warning
501 merge_point = &merge;
506 merge = *merge_point;
510 Optype invop = InverseOp(op); // inverse of p
519 for (auto prev : m_next[invop])
521 if (prev->OperationTerm(b, rev, merge_point, merge_op) == b)
523 //assert(SanityCheck());
528 for (auto next : m_next[op])
530 if (next->OperationTerm(b, fwd, merge_point, merge_op) == b)
532 //assert(SanityCheck());
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;
550 if (m_next[op].size() == 0)
557 //assert(SanityCheck());
562 ParanoidNumber * ParanoidNumber::OperationFactor(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
564 #ifdef PARANOID_CACHE_RESULTS
565 m_cached_result = NAN;
567 #ifdef PARANOID_SIZE_LIMIT
568 if (m_size >= PARANOID_SIZE_LIMIT)
571 m_value *= b->Digit();
573 m_value /= b->Digit();
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);
584 if (Floating() && m_value == 0)
589 if (Floating() && m_value == 1 && op == MULTIPLY)
591 m_value = b->m_value;
592 for (int i = 0; i < NOP; ++i)
594 for (auto n : m_next[i])
597 swap(m_next[i], b->m_next[i]);
599 //assert(SanityCheck());
602 if (b->Floating() && b->m_value == 1)
604 if (b->Floating() && b->m_value == 0 && op == MULTIPLY)
611 if (NoTerms() && b->NoTerms())
613 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
615 Optype mulop = (op == MULTIPLY) ? MULTIPLY : DIVIDE;
616 for (auto mul : b->m_next[MULTIPLY])
618 delete(OperationFactor(mul, mulop));
620 Optype divop = (op == MULTIPLY) ? DIVIDE : MULTIPLY;
621 for (auto div : b->m_next[DIVIDE])
622 delete(OperationFactor(div, divop));
624 b->m_next[DIVIDE].clear();
625 b->m_next[MULTIPLY].clear();
631 bool parent = (merge_point == NULL);
632 ParanoidNumber * merge = this;
636 merge_point = &merge;
641 merge = *merge_point;
645 Optype invop = InverseOp(op); // inverse of p
654 ParanoidNumber * cpy_b = new ParanoidNumber(*b);
655 for (auto prev : m_next[invop])
657 if (prev->OperationFactor(b, rev, merge_point, merge_op) == b)
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));
668 for (auto next : m_next[op])
670 if (next->OperationFactor(b, fwd, merge_point, merge_op) == b)
672 for (auto add : m_next[ADD])
674 delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
676 for (auto sub : m_next[SUBTRACT])
678 delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
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));
694 #ifdef PARANOID_SIZE_LIMIT
695 m_size += 1+b->m_size;
698 //assert(SanityCheck());
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
713 ParanoidNumber * ParanoidNumber::Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
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);
729 string ParanoidNumber::PStr() const
732 for (int i = 0; i < NOP; ++i)
734 Optype f = Optype(i);
736 for (auto n : m_next[f])
738 s << OpChar(f) << n->PStr();
744 bool ParanoidNumber::Simplify(Optype op)
750 //assert(SanityCheck());
751 vector<ParanoidNumber*> next;
753 swap(m_next[op], next);
755 //assert(m_next[op].size() == 0);
756 //assert(SanityCheck());
760 else if (op == SUBTRACT)
764 vector<ParanoidNumber*> hold[2];
765 if (op == MULTIPLY || op == DIVIDE)
767 swap(m_next[ADD], hold[0]);
768 swap(m_next[SUBTRACT], hold[1]);
771 for (vector<ParanoidNumber*>::iterator n = next.begin(); n != next.end(); ++n)
775 for (vector<ParanoidNumber*>::iterator m = n; m != next.end(); ++m)
782 ParanoidNumber * parent = this;
784 ParanoidNumber * result = (*n)->Operation(*m, fwd, &parent, &mop);
787 #ifdef PARANOID_SIZE_LIMIT
788 m_size -= (1+result->m_size);
802 #ifdef PARANOID_SIZE_LIMIT
803 if (Operation(n, op) == n)
805 m_size -= (1+n->m_size);
809 delete(Operation(n, op));
814 if (op == MULTIPLY || op == DIVIDE)
816 swap(m_next[ADD], hold[0]);
817 swap(m_next[SUBTRACT], hold[1]);
820 set<ParanoidNumber*> s;
821 //if (!SanityCheck(s))
823 // Error("Simplify broke Sanity");
825 return (next.size() > m_next[op].size());
828 bool ParanoidNumber::FullSimplify()
831 result |= Simplify(MULTIPLY);
832 result |= Simplify(DIVIDE);
833 result |= Simplify(ADD);
834 result |= Simplify(SUBTRACT);
838 ParanoidNumber::digit_t ParanoidNumber::Digit() const
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;
845 if (!isnan(float(result))) // le sigh ambiguous function compiler warnings
851 for (auto mul : m_next[MULTIPLY])
853 result *= mul->Digit();
855 for (auto div : m_next[DIVIDE])
857 result /= div->Digit();
859 for (auto add : m_next[ADD])
860 result += add->Digit();
861 for (auto sub : m_next[SUBTRACT])
862 result -= sub->Digit();
867 ParanoidNumber::digit_t ParanoidNumber::GetFactors() const
870 for (auto mul : m_next[MULTIPLY])
871 value *= mul->Digit();
872 for (auto div : m_next[DIVIDE])
873 value /= div->Digit();
878 ParanoidNumber::digit_t ParanoidNumber::GetTerms() const
881 for (auto add : m_next[ADD])
882 value += add->Digit();
883 for (auto sub : m_next[SUBTRACT])
884 value -= sub->Digit();
888 bool ParanoidNumber::SanityCheck(set<ParanoidNumber*> & visited) const
892 Error("NULL pointer in tree");
896 if (visited.find((ParanoidNumber*)this) != visited.end())
898 Error("I think I've seen this tree before...");
902 visited.insert((ParanoidNumber*)this);
904 for (auto add : m_next[ADD])
906 if (!add->SanityCheck(visited))
909 for (auto sub : m_next[SUBTRACT])
911 if (!sub->SanityCheck(visited))
914 for (auto mul : m_next[MULTIPLY])
916 if (!mul->SanityCheck(visited))
920 for (auto div : m_next[DIVIDE])
922 if (!div->SanityCheck(visited))
924 if (div->Digit() == 0)
926 Error("Divide by zero");
933 void ParanoidNumber::Negate()
935 swap(m_next[ADD], m_next[SUBTRACT]);
939 #ifdef PARANOID_USE_ARENA
941 void * ParanoidNumber::operator new(size_t s)
943 return g_arena.allocate(s);
946 void ParanoidNumber::operator delete(void * p)
948 g_arena.deallocate(p);
951 ParanoidNumber::Arena::Arena(int64_t block_size) : m_block_size(block_size), m_spare(NULL)
953 m_blocks.push_back({malloc(block_size*sizeof(ParanoidNumber)),0});
956 ParanoidNumber::Arena::~Arena()
958 for (auto block : m_blocks)
965 void * ParanoidNumber::Arena::allocate(size_t s)
969 void * result = m_spare;
974 Block & b = m_blocks.back();
975 void * result = (ParanoidNumber*)(b.memory) + (b.used++);
976 if (b.used >= m_block_size)
979 Debug("Add block of size %d", m_block_size);
980 m_blocks.push_back({malloc(m_block_size*sizeof(ParanoidNumber)), 0});
985 void ParanoidNumber::Arena::deallocate(void * p)
989 #endif //PARANOID_USE_ARENA