1 #include "paranoidnumber.h"
14 #ifdef PARANOID_USE_ARENA
15 ParanoidNumber::Arena ParanoidNumber::g_arena;
16 #endif //PARANOID_USE_ARENA
18 ParanoidNumber::~ParanoidNumber()
20 for (int i = 0; i < NOP; ++i)
22 for (auto n : m_next[i])
27 ParanoidNumber::ParanoidNumber(const string & str) : m_value(0), m_next()
29 #ifdef PARANOID_SIZE_LIMIT
32 #ifdef PARANOID_CACHE_RESULTS
33 m_cached_result = NAN;
38 while (str[dp] != '\0' && str[dp] != '.')
43 while (str[end] != '\0')
46 for (int i = dp-1; i >= 0; --i)
48 ParanoidNumber b(str[i]-'0');
54 for (int i = dp+1; i < end; ++i)
57 ParanoidNumber b(str[i]-'0');
63 ParanoidNumber & ParanoidNumber::operator=(const ParanoidNumber & a)
65 //assert(this != NULL);
67 #ifdef PARANOID_SIZE_LIMIT
72 #ifdef PARANOID_CACHE_RESULT
73 m_cached_result = a.m_cached_result;
75 for (int i = 0; i < NOP; ++i)
77 for (auto n : m_next[i])
80 for (auto n : a.m_next[i])
81 m_next[i].push_back(new ParanoidNumber(*n));
84 for (unsigned j = 0; j < m_next[i].size() && j < a.m_next[i].size(); ++j)
86 if (a.m_next[i][j] != NULL)
87 m_next[i][j]->operator=(*(a.m_next[i][j]));
90 for (unsigned j = a.m_next[i].size(); j < m_next[i].size(); ++j)
94 m_next[i].resize(a.m_next[i].size());
101 string ParanoidNumber::Str() const
104 //assert(this != NULL);
107 s << (double)m_value;
109 for (auto mul : m_next[MULTIPLY])
112 if (!mul->Floating())
113 result += "(" + mul->Str() + ")";
115 result += mul->Str();
117 for (auto div : m_next[DIVIDE])
120 if (!div->Floating())
121 result += "(" + div->Str() + ")";
123 result += div->Str();
126 for (auto add : m_next[ADD])
129 if (!add->Floating())
130 result += "(" + add->Str() + ")";
132 result += add->Str();
134 for (auto sub : m_next[SUBTRACT])
137 if (!sub->Floating())
138 result += "(" + sub->Str() + ")";
140 result += sub->Str();
148 bool TrustingOp<float>(float & a, const float & b, Optype op)
152 feclearexcept(FE_ALL_EXCEPT);
167 a = (a >= 0) ? INFINITY : -INFINITY;
176 return !fetestexcept(FE_ALL_EXCEPT);
180 bool TrustingOp<double>(double & a, const double & b, Optype op)
182 feclearexcept(FE_ALL_EXCEPT);
197 a = (a >= 0) ? INFINITY : -INFINITY;
205 return !fetestexcept(FE_ALL_EXCEPT);
211 bool TrustingOp<long double>(long double & a, const long double & b, Optype op)
215 feclearexcept(FE_ALL_EXCEPT);
230 a = (a >= 0) ? INFINITY : -INFINITY;
239 return !fetestexcept(FE_ALL_EXCEPT);
244 bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op)
252 exact = (abs(sa) <= 127);
256 exact = (abs(sa) <= 127);
260 exact = (abs(sa) <= 127);
263 exact = (b != 0 && sa > b && sa % b == 0);
274 ParanoidNumber & ParanoidNumber::operator+=(const digit_t & a)
277 //assert(this != NULL);
278 delete Operation(new ParanoidNumber(a), ADD);
285 ParanoidNumber & ParanoidNumber::operator-=(const digit_t & a)
287 delete Operation(new ParanoidNumber(a), SUBTRACT);
293 ParanoidNumber & ParanoidNumber::operator*=(const digit_t & a)
295 delete Operation(new ParanoidNumber(a), MULTIPLY);
302 ParanoidNumber & ParanoidNumber::operator/=(const digit_t & a)
304 delete Operation(new ParanoidNumber(a), DIVIDE);
311 ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
314 //assert(this != NULL);
315 delete Operation(new ParanoidNumber(a), ADD);
322 ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
324 delete Operation(new ParanoidNumber(a), SUBTRACT);
330 ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
332 delete Operation(new ParanoidNumber(a), MULTIPLY);
339 ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
341 delete Operation(new ParanoidNumber(a), DIVIDE);
347 ParanoidNumber & ParanoidNumber::operator=(const digit_t & a)
349 for (int i = 0; i < NOP; ++i)
351 for (auto n : m_next[i])
355 #ifdef PARANOID_CACHE_RESULT
362 ParanoidNumber * ParanoidNumber::OperationTerm(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
364 ////assert(b->SanityCheck());
365 #ifdef PARANOID_CACHE_RESULTS
366 m_cached_result = NAN;
368 #ifdef PARANOID_SIZE_LIMIT
369 if (m_size > PARANOID_SIZE_LIMIT)
372 m_value += b->Digit();
374 m_value -= b->Digit();
377 //Debug("At size limit %d", m_size);
381 if (Floating() && m_value == 0) // 0 + b = b
383 m_value = b->m_value;
387 swap(b->m_next[ADD], b->m_next[SUBTRACT]);
390 for (int i = 0; i < NOP; ++i)
392 m_next[i] = b->m_next[i];
393 b->m_next[i].clear();
396 //assert(SanityCheck());
399 if (b->Floating() && b->m_value == 0) // a + 0 = a
404 if ((NoFactors() && b->NoFactors())
405 || (GetFactors() == b->GetFactors()))
407 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
409 Optype addop = (op == ADD) ? ADD : SUBTRACT;
410 for (auto add : b->m_next[ADD])
412 delete (OperationTerm(add, addop));
414 Optype subop = (op == ADD) ? SUBTRACT : ADD;
415 for (auto sub : b->m_next[SUBTRACT])
416 delete (OperationTerm(sub, subop));
418 b->m_next[ADD].clear();
419 b->m_next[SUBTRACT].clear();
420 //assert(SanityCheck());
427 bool parent = (merge_point == NULL);
428 ParanoidNumber * merge = this;
430 //assert(mop != NOP); // silence compiler warning
433 merge_point = &merge;
438 merge = *merge_point;
442 Optype invop = InverseOp(op); // inverse of p
451 for (auto prev : m_next[invop])
453 if (prev->OperationTerm(b, rev, merge_point, merge_op) == b)
455 //assert(SanityCheck());
460 for (auto next : m_next[op])
462 if (next->OperationTerm(b, fwd, merge_point, merge_op) == b)
464 //assert(SanityCheck());
474 //merge->m_next[*merge_op].push_back(b);
475 m_next[op].push_back(b);
476 #ifdef PARANOID_SIZE_LIMIT
477 m_size += 1+b->m_size;
482 if (m_next[op].size() == 0)
489 //assert(SanityCheck());
494 ParanoidNumber * ParanoidNumber::OperationFactor(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
496 //assert(SanityCheck());
497 //assert(b->SanityCheck());
498 #ifdef PARANOID_CACHE_RESULTS
499 m_cached_result = NAN;
501 #ifdef PARANOID_SIZE_LIMIT
502 if (m_size > PARANOID_SIZE_LIMIT)
505 m_value *= b->Digit();
507 m_value /= b->Digit();
508 //Debug("At size limit %d", m_size);
513 if (Floating() && m_value == 0)
518 if (Floating() && m_value == 1 && op == MULTIPLY)
520 m_value = b->m_value;
521 for (int i = 0; i < NOP; ++i)
523 for (auto n : m_next[i])
526 swap(m_next[i], b->m_next[i]);
528 //assert(SanityCheck());
531 if (b->Floating() && b->m_value == 1)
536 if (NoTerms() && b->NoTerms())
538 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
540 Optype mulop = (op == MULTIPLY) ? MULTIPLY : DIVIDE;
541 for (auto mul : b->m_next[MULTIPLY])
543 delete(OperationFactor(mul, mulop));
545 Optype divop = (op == MULTIPLY) ? DIVIDE : MULTIPLY;
546 for (auto div : b->m_next[DIVIDE])
547 delete(OperationFactor(div, divop));
549 b->m_next[DIVIDE].clear();
550 b->m_next[MULTIPLY].clear();
556 bool parent = (merge_point == NULL);
557 ParanoidNumber * merge = this;
561 merge_point = &merge;
566 merge = *merge_point;
570 Optype invop = InverseOp(op); // inverse of p
579 ParanoidNumber * cpy_b = new ParanoidNumber(*b);
580 for (auto prev : m_next[invop])
582 if (prev->OperationFactor(b, rev, merge_point, merge_op) == b)
584 for (auto add : m_next[ADD])
585 delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
586 for (auto sub : m_next[SUBTRACT])
587 delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
593 for (auto next : m_next[op])
595 if (next->OperationFactor(b, fwd, merge_point, merge_op) == b)
597 for (auto add : m_next[ADD])
599 delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
601 for (auto sub : m_next[SUBTRACT])
603 delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
613 m_next[op].push_back(b);
614 for (auto add : m_next[ADD])
615 delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
616 for (auto sub : m_next[SUBTRACT])
617 delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
619 #ifdef PARANOID_SIZE_LIMIT
620 m_size += 1+b->m_size;
623 //assert(SanityCheck());
633 * Performs the operation on a with argument b (a += b, a -= b, a *= b, a /= b)
634 * @returns b if b can safely be deleted
635 * @returns NULL if b has been merged with a
636 * append indicates that b should be merged
638 ParanoidNumber * ParanoidNumber::Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
645 if (op == SUBTRACT || op == ADD)
646 return OperationTerm(b, op, merge_point, merge_op);
647 if (op == MULTIPLY || op == DIVIDE)
648 return OperationFactor(b, op, merge_point, merge_op);
654 string ParanoidNumber::PStr() const
657 for (int i = 0; i < NOP; ++i)
659 Optype f = Optype(i);
661 for (auto n : m_next[f])
663 s << OpChar(f) << n->PStr();
669 bool ParanoidNumber::Simplify(Optype op)
675 //assert(SanityCheck());
676 vector<ParanoidNumber*> next;
678 swap(m_next[op], next);
680 //assert(m_next[op].size() == 0);
681 //assert(SanityCheck());
685 else if (op == SUBTRACT)
689 vector<ParanoidNumber*> hold[2];
690 if (op == MULTIPLY || op == DIVIDE)
692 swap(m_next[ADD], hold[0]);
693 swap(m_next[SUBTRACT], hold[1]);
696 for (vector<ParanoidNumber*>::iterator n = next.begin(); n != next.end(); ++n)
700 for (vector<ParanoidNumber*>::iterator m = n; m != next.end(); ++m)
707 ParanoidNumber * parent = this;
709 ParanoidNumber * result = (*n)->Operation(*m, fwd, &parent, &mop);
712 #ifdef PARANOID_SIZE_LIMIT
713 m_size -= (1+result->m_size);
727 #ifdef PARANOID_SIZE_LIMIT
728 if (Operation(n, op) == n)
730 m_size -= (1+n->m_size);
734 delete(Operation(n, op));
739 if (op == MULTIPLY || op == DIVIDE)
741 swap(m_next[ADD], hold[0]);
742 swap(m_next[SUBTRACT], hold[1]);
745 set<ParanoidNumber*> s;
746 //if (!SanityCheck(s))
748 // Error("Simplify broke Sanity");
750 return (next.size() > m_next[op].size());
753 bool ParanoidNumber::FullSimplify()
756 result |= Simplify(MULTIPLY);
757 result |= Simplify(DIVIDE);
758 result |= Simplify(ADD);
759 result |= Simplify(SUBTRACT);
763 ParanoidNumber::digit_t ParanoidNumber::Digit() const
766 // Get around the absurd requirement that const correctness be observed.
767 #ifdef PARANOID_CACHE_RESULTS
768 digit_t & result = ((ParanoidNumber*)(this))->m_cached_result;
770 if (!isnan(float(result))) // le sigh ambiguous function compiler warnings
776 for (auto mul : m_next[MULTIPLY])
778 result *= mul->Digit();
780 for (auto div : m_next[DIVIDE])
782 result /= div->Digit();
784 for (auto add : m_next[ADD])
785 result += add->Digit();
786 for (auto sub : m_next[SUBTRACT])
787 result -= sub->Digit();
792 ParanoidNumber::digit_t ParanoidNumber::GetFactors() const
795 for (auto mul : m_next[MULTIPLY])
796 value *= mul->Digit();
797 for (auto div : m_next[DIVIDE])
798 value /= div->Digit();
803 ParanoidNumber::digit_t ParanoidNumber::GetTerms() const
806 for (auto add : m_next[ADD])
807 value += add->Digit();
808 for (auto sub : m_next[SUBTRACT])
809 value -= sub->Digit();
813 bool ParanoidNumber::SanityCheck(set<ParanoidNumber*> & visited) const
817 Error("NULL pointer in tree");
821 if (visited.find((ParanoidNumber*)this) != visited.end())
823 Error("I think I've seen this tree before...");
827 visited.insert((ParanoidNumber*)this);
829 for (auto add : m_next[ADD])
831 if (!add->SanityCheck(visited))
834 for (auto sub : m_next[SUBTRACT])
836 if (!sub->SanityCheck(visited))
839 for (auto mul : m_next[MULTIPLY])
841 if (!mul->SanityCheck(visited))
845 for (auto div : m_next[DIVIDE])
847 if (!div->SanityCheck(visited))
853 #ifdef PARANOID_USE_ARENA
855 void * ParanoidNumber::operator new(size_t s)
857 return g_arena.allocate(s);
860 void ParanoidNumber::operator delete(void * p)
862 g_arena.deallocate(p);
865 ParanoidNumber::Arena::Arena(int64_t block_size) : m_block_size(block_size), m_spare(NULL)
867 m_blocks.push_back({malloc(block_size*sizeof(ParanoidNumber)),0});
870 ParanoidNumber::Arena::~Arena()
872 for (auto block : m_blocks)
879 void * ParanoidNumber::Arena::allocate(size_t s)
883 void * result = m_spare;
888 Block & b = m_blocks.back();
889 void * result = (ParanoidNumber*)(b.memory) + (b.used++);
890 if (b.used >= m_block_size)
893 Debug("Add block of size %d", m_block_size);
894 m_blocks.push_back({malloc(m_block_size*sizeof(ParanoidNumber)), 0});
899 void ParanoidNumber::Arena::deallocate(void * p)
903 #endif //PARANOID_USE_ARENA