Compile with float/double again.
[ipdf/code.git] / src / paranoidnumber.cpp
index 49cceaf..60802eb 100644 (file)
@@ -6,24 +6,43 @@
 #include <cassert>
 #include <iostream>
 
+// here be many copy paste bugs
+
 using namespace std;
 namespace IPDF
 {
-int64_t ParanoidNumber::g_count = 0;
 
 
+#ifdef PARANOID_USE_ARENA
+ParanoidNumber::Arena ParanoidNumber::g_arena;
+#endif //PARANOID_USE_ARENA
+
 ParanoidNumber::~ParanoidNumber()
 {
-       g_count--;
        for (int i = 0; i < NOP; ++i)
-               delete m_next[i];
+       {
+               for (auto n : m_next[i])
+                       delete n;
+       }
 }
 
-ParanoidNumber::ParanoidNumber(const char * str) : m_value(0)
+ParanoidNumber::ParanoidNumber(const string & str) : m_value(0), m_next()
 {
-       Construct();
+       #ifdef PARANOID_SIZE_LIMIT
+               m_size = 1;
+       #endif
+       #ifdef PARANOID_CACHE_RESULTS
+       m_cache_valid = false;
+       #endif
+       
        int dp = 0;
        int end = 0;
+       bool negate = str[0] == '-';
+       if (negate)
+       {
+               dp++;
+               end++;
+       }
        while (str[dp] != '\0' && str[dp] != '.')
        {
                ++dp;
@@ -32,7 +51,7 @@ ParanoidNumber::ParanoidNumber(const char * str) : m_value(0)
        while (str[end] != '\0')
                ++end;
        ParanoidNumber m(1);
-       for (int i = dp-1; i >= 0; --i)
+       for (int i = dp-1; i >= negate; --i)
        {
                ParanoidNumber b(str[i]-'0');
                b*=m;
@@ -42,80 +61,89 @@ ParanoidNumber::ParanoidNumber(const char * str) : m_value(0)
        ParanoidNumber n(1);
        for (int i = dp+1; i < end; ++i)
        {
-               Debug("{%s} /= 10", n.Str().c_str());
                n/=10;
-               Debug("{%s}", n.Str().c_str());
                ParanoidNumber b(str[i]-'0');
                b*=n;
-               Debug("{%s} += {%s}", Str().c_str(), b.Str().c_str());
                this->operator+=(b);
        }
+       
+       if (negate)
+               Negate();
+       
+       #ifdef PARANOID_COMPARE_EPSILON
+               double d = strtod(str.c_str(), NULL);
+               CompareForSanity(d, d);
+       #endif
 }
 
 ParanoidNumber & ParanoidNumber::operator=(const ParanoidNumber & a)
 {
+       //assert(this != NULL);
+       
+       #ifdef PARANOID_SIZE_LIMIT
+               m_size = a.m_size;
+       #endif
+       
        m_value = a.m_value;
+       #ifdef PARANOID_CACHE_RESULTS
+       m_cached_result = a.m_cached_result;
+       m_cache_valid = a.m_cache_valid;
+       #endif
        for (int i = 0; i < NOP; ++i)
        {
-               if (a.m_next[i] == NULL)
-               {
-                       if (m_next[i] != NULL)
-                               delete m_next[i];
-                       m_next[i] = NULL;
-                       continue;
-               }
-                       
-               if (m_next[i] != NULL)
-               {
-                       m_next[i]->operator=(*(a.m_next[i]));
-               }
-               else
-               {
-                       m_next[i] = new ParanoidNumber(*(a.m_next[i]));
-               }
-       }       
+               for (auto n : m_next[i])
+                       delete n;
+               m_next[i].clear();
+               for (auto n : a.m_next[i])
+                       m_next[i].push_back(new ParanoidNumber(*n));
+       }
+       #ifdef PARANOID_COMPARE_EPSILON
+               CompareForSanity(a.Digit(),a.Digit());
+       #endif
        return *this;
 }
 
 
 string ParanoidNumber::Str() const
 {
+       
+       //assert(this != NULL);
        string result("");
        stringstream s;
        s << (double)m_value;
        result += s.str();
-       if (m_next[MULTIPLY] != NULL)
+       for (auto mul : m_next[MULTIPLY])
        {
                result += "*";
-               if (m_next[MULTIPLY]->m_next[ADD] != NULL || m_next[MULTIPLY]->m_next[SUBTRACT] != NULL)
-                       result += "(" + m_next[MULTIPLY]->Str() + ")";
+               if (!mul->Floating())
+                       result += "(" + mul->Str() + ")";
                else
-                       result += m_next[MULTIPLY]->Str();
+                       result += mul->Str();
        }
-       if (m_next[DIVIDE] != NULL)
+       for (auto div : m_next[DIVIDE])
        {
                result += "/";
-               if (m_next[DIVIDE]->m_next[ADD] != NULL || m_next[DIVIDE]->m_next[SUBTRACT] != NULL)
-                       result += "(" + m_next[DIVIDE]->Str() + ")";
+               if (!div->Floating())
+                       result += "(" + div->Str() + ")";
                else
-                       result += m_next[DIVIDE]->Str();
+                       result += div->Str();
        }       
        
-       if (m_next[ADD] != NULL)
+       for (auto add : m_next[ADD])
        {
                result += "+";
-               if (m_next[ADD]->m_next[MULTIPLY] != NULL || m_next[ADD]->m_next[DIVIDE] != NULL)
-                       result += "(" + m_next[ADD]->Str() + ")";
+               if (!add->Floating())
+                       result += "(" + add->Str() + ")";
                else
-                       result += m_next[ADD]->Str();
+                       result += add->Str();
        }
-       if (m_next[SUBTRACT] != NULL)
+       for (auto sub : m_next[SUBTRACT])
        {
                result += "-";
-               if (m_next[SUBTRACT]->m_next[MULTIPLY] != NULL || m_next[SUBTRACT]->m_next[DIVIDE] != NULL)
-                       result += "(" + m_next[SUBTRACT]->Str() + ")";
+               if (!sub->Floating())
+                       result += "(" + sub->Str() + ")";
                else
-                       result += m_next[SUBTRACT]->Str();
+                       result += sub->Str();
        }
        
 
@@ -125,6 +153,8 @@ string ParanoidNumber::Str() const
 template <>
 bool TrustingOp<float>(float & a, const float & b, Optype op)
 {
+
+       
        feclearexcept(FE_ALL_EXCEPT);
        switch (op)
        {
@@ -138,6 +168,12 @@ bool TrustingOp<float>(float & a, const float & b, Optype op)
                        a *= b;
                        break;
                case DIVIDE:
+                       if (b == 0)
+                       {
+                               a = (a >= 0) ? INFINITY : -INFINITY;
+                               return false;
+                       }
+                       
                        a /= b;
                        break;
                case NOP:
@@ -162,6 +198,45 @@ bool TrustingOp<double>(double & a, const double & b, Optype op)
                        a *= b;
                        break;
                case DIVIDE:
+                       if (b == 0)
+                       {
+                               a = (a >= 0) ? INFINITY : -INFINITY;
+                               return false;
+                       }
+                       a /= b;
+                       break;
+               case NOP:
+                       break;
+       }
+       return !fetestexcept(FE_ALL_EXCEPT);
+}
+
+
+
+template <>
+bool TrustingOp<long double>(long double & a, const long double & b, Optype op)
+{
+
+       
+       feclearexcept(FE_ALL_EXCEPT);
+       switch (op)
+       {
+               case ADD:
+                       a += b;
+                       break;
+               case SUBTRACT:
+                       a -= b;
+                       break;
+               case MULTIPLY:
+                       a *= b;
+                       break;
+               case DIVIDE:
+                       if (b == 0)
+                       {
+                               a = (a >= 0) ? INFINITY : -INFINITY;
+                               return false;
+                       }
+                       
                        a /= b;
                        break;
                case NOP:
@@ -170,6 +245,7 @@ bool TrustingOp<double>(double & a, const double & b, Optype op)
        return !fetestexcept(FE_ALL_EXCEPT);
 }
 
+
 template <>
 bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op)
 {
@@ -201,245 +277,729 @@ bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op)
 }
 
 
+ParanoidNumber & ParanoidNumber::operator+=(const digit_t & a)
+{
+       #ifdef PARANOID_COMPARE_EPSILON
+               digit_t compare = Digit();
+               compare += a;
+       #endif
+       delete Operation(new ParanoidNumber(a), ADD);
+       Simplify(SUBTRACT);
+       Simplify(ADD);
+       #ifdef PARANOID_COMPARE_EPSILON
+               CompareForSanity(compare, a);
+       #endif
+       return *this;
+}
+
+
+ParanoidNumber & ParanoidNumber::operator-=(const digit_t & a)
+{
+       #ifdef PARANOID_COMPARE_EPSILON
+               digit_t compare = Digit();
+               compare -= a;
+       #endif
+       delete Operation(new ParanoidNumber(a), SUBTRACT);
+       Simplify(ADD);
+       Simplify(SUBTRACT);
+       #ifdef PARANOID_COMPARE_EPSILON
+               CompareForSanity(compare, a);
+       #endif
+       return *this;
+}
+
+ParanoidNumber & ParanoidNumber::operator*=(const digit_t & a)
+{
+       #ifdef PARANOID_COMPARE_EPSILON
+               digit_t compare = Digit();
+               compare *= a;
+       #endif
+       delete Operation(new ParanoidNumber(a), MULTIPLY);
+       Simplify(DIVIDE);
+       Simplify(MULTIPLY);
+       #ifdef PARANOID_COMPARE_EPSILON
+               CompareForSanity(compare, a);
+       #endif
+       return *this;
+}
+
+
+ParanoidNumber & ParanoidNumber::operator/=(const digit_t & a)
+{
+       #ifdef PARANOID_COMPARE_EPSILON
+               digit_t compare = Digit();
+               compare /= a;
+       #endif
+       delete Operation(new ParanoidNumber(a), DIVIDE);
+       Simplify(MULTIPLY);
+       Simplify(DIVIDE);
+       #ifdef PARANOID_COMPARE_EPSILON
+               CompareForSanity(compare, a);
+       #endif
+       return *this;
+}
+
+
 ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
 {
+       #ifdef PARANOID_COMPARE_EPSILON
+               digit_t compare = Digit();
+               compare += a.Digit();
+       #endif
        delete Operation(new ParanoidNumber(a), ADD);
+       Simplify(SUBTRACT);
+       Simplify(ADD);
+       #ifdef PARANOID_COMPARE_EPSILON
+               CompareForSanity(compare, a.Digit());
+       #endif
        return *this;
 }
 
 
 ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
 {
+       #ifdef PARANOID_COMPARE_EPSILON
+               digit_t compare = Digit();
+               compare -= a.Digit();
+       #endif
        delete Operation(new ParanoidNumber(a), SUBTRACT);
+       Simplify(ADD);
+       Simplify(SUBTRACT);
+       #ifdef PARANOID_COMPARE_EPSILON
+               CompareForSanity(compare, a.Digit());
+       #endif
        return *this;
 }
 
 ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
 {
+       #ifdef PARANOID_COMPARE_EPSILON
+               digit_t compare = Digit();
+               compare *= a.Digit();
+       #endif
        delete Operation(new ParanoidNumber(a), MULTIPLY);
+       Simplify(DIVIDE);
+       Simplify(MULTIPLY);
+       #ifdef PARANOID_COMPARE_EPSILON
+               CompareForSanity(compare, a.Digit());
+       #endif
        return *this;
 }
 
 
 ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
 {
+       #ifdef PARANOID_COMPARE_EPSILON
+               digit_t compare = Digit();
+               compare /= a.Digit();
+       #endif
        delete Operation(new ParanoidNumber(a), DIVIDE);
+       Simplify(MULTIPLY);
+       Simplify(DIVIDE);
+       #ifdef PARANOID_COMPARE_EPSILON
+               CompareForSanity(compare, a.Digit());
+       #endif
        return *this;
 }
 
-/**
- * Performs the operation on a with argument b (a += b, a -= b, a *= b, a /= b)
- * @returns b if b can safely be deleted
- * @returns NULL if b has been merged with a
- * append indicates that b should be merged
- */
-ParanoidNumber * ParanoidNumber::Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** parent)
+ParanoidNumber & ParanoidNumber::operator=(const digit_t & a)
 {
-       if (b == NULL)
-               return NULL;
+
+       for (int i = 0; i < NOP; ++i)
+       {
+               for (auto n : m_next[i])
+                       delete n;
+               m_next[i].clear();
+       }
+       m_value = a;
+       #ifdef PARANOID_CACHE_RESULTS
+       m_cached_result = a;
+       m_cache_valid = true;
+       #endif
+
+       #ifdef PARANOID_COMPARE_EPSILON
+               CompareForSanity(a,a);
+       #endif
+       
+       return *this;
+}
+
+// a + b
+ParanoidNumber * ParanoidNumber::OperationTerm(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
+{
+       ////assert(b->SanityCheck());
+       #ifdef PARANOID_CACHE_RESULTS
+       m_cache_valid = false;
+       #endif
+       #ifdef PARANOID_SIZE_LIMIT
+               if (m_size + b->m_size >= PARANOID_SIZE_LIMIT)
+               {
+                       this->operator=(this->Digit());
+                       if (op == ADD)
+                               m_value += b->Digit();
+                       else
+                               m_value -= b->Digit();
+                       m_size = 1;
+                       #ifdef PARANOID_CACHE_RESULTS
+                               m_cached_result = m_value;
+                               m_cache_valid = true;
+                       #endif
+                       return b;
+               }
+               //Debug("At size limit %d", m_size);
+       #endif 
+       
+
+       if (Floating() && m_value == 0) // 0 + b = b
+       {
+               m_value = b->m_value;
+               if (op == SUBTRACT)
+               {
+                       m_value = -m_value;
+                       swap(b->m_next[ADD], b->m_next[SUBTRACT]);
+               }
                
-       Optype invop = InverseOp(op); // inverse of p
-       ParanoidNumber * append_at = this;
+               for (int i = 0; i < NOP; ++i)
+               {
+                       m_next[i] = b->m_next[i];
+                       b->m_next[i].clear();
+               }
+               
+               //assert(SanityCheck());
+               return b;
+       }
+       if (b->Floating() && b->m_value == 0) // a + 0 = a
+               return b;
+               
+
        
-       if (Floating())
+       if ((NoFactors() && b->NoFactors())
+               || (GetFactors() == b->GetFactors()))
        {
-               if ((op == ADD || op == SUBTRACT) && (m_value == 0))
+               if (ParanoidOp<digit_t>(m_value, b->m_value, op))
                {
-                       m_value = b->m_value;
-                       for (int i = 0; i < NOP; ++i)
+                       Optype addop = (op == ADD) ? ADD : SUBTRACT;
+                       for (auto add : b->m_next[ADD])
                        {
-                               m_next[i] = b->m_next[i];
-                               b->m_next[i] = NULL;
+                               delete (OperationTerm(add, addop));
                        }
+                       Optype subop = (op == ADD) ? SUBTRACT : ADD;
+                       for (auto sub : b->m_next[SUBTRACT])
+                               delete (OperationTerm(sub, subop));
+                               
+                       b->m_next[ADD].clear();
+                       b->m_next[SUBTRACT].clear();
+                       //assert(SanityCheck());
                        return b;
                }
-               if ((op == MULTIPLY) && (m_value == 1))
+       }
+
+       
+       
+       bool parent = (merge_point == NULL);
+       ParanoidNumber * merge = this;
+       Optype mop = op;
+       //assert(mop != NOP); // silence compiler warning
+       if (parent)
+       {
+               merge_point = &merge;
+               merge_op = &mop;
+       }
+       else
+       {
+               merge = *merge_point;
+               mop = *merge_op;
+       }
+               
+       Optype invop = InverseOp(op); // inverse of p
+       Optype fwd = op;
+       Optype rev = invop;
+       if (op == SUBTRACT)
+       {
+               fwd = ADD;
+               rev = SUBTRACT;
+       }
+       
+       for (auto prev : m_next[invop])
+       {
+               if (prev->OperationTerm(b, rev, merge_point, merge_op) == b)
                {
-                       m_value = b->m_value;
-                       for (int i = 0; i < NOP; ++i)
-                       {
-                               m_next[i] = b->m_next[i];
-                               b->m_next[i] = NULL;
-                       }
-                       return b;
+                       //assert(SanityCheck());
                        return b;
                }
                
        }
-       
-       if (b->Floating())
+       for (auto next : m_next[op])
        {
-               if ((op == ADD || op == SUBTRACT) && (b->m_value == 0))
+               if (next->OperationTerm(b, fwd, merge_point, merge_op) == b)
+               {
+                       //assert(SanityCheck());
                        return b;
-               if ((op == MULTIPLY || op == DIVIDE) && (b->m_value == 1))
+               }
+       }
+       
+
+       
+       
+       if (parent)
+       {
+               //merge->m_next[*merge_op].push_back(b);
+               m_next[op].push_back(b);
+               #ifdef PARANOID_SIZE_LIMIT
+                       m_size += b->m_size;
+               #endif  
+       }
+       else
+       {
+               if (m_next[op].size() == 0)
+               {
+                       *merge_point = this;
+                       *merge_op = op;
+               }
+       }
+
+       //assert(SanityCheck());
+
+       return NULL;
+}
+
+ParanoidNumber * ParanoidNumber::OperationFactor(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
+{
+       #ifdef PARANOID_CACHE_RESULTS
+       m_cache_valid = false;
+       #endif
+       #ifdef PARANOID_SIZE_LIMIT
+               if (m_size + b->m_size >= PARANOID_SIZE_LIMIT)
+               {
+                       this->operator=(this->Digit());
+                       if (op == MULTIPLY)
+                               m_value *= b->Digit();
+                       else
+                               m_value /= b->Digit();
+                       m_size = 1;
+                       #ifdef PARANOID_CACHE_RESULTS
+                               m_cached_result = m_value;
+                               m_cache_valid = true;
+                       #endif
+                       //Debug("Cut off %p", this);
                        return b;
+                       
+               }
+       #endif  
+
+       if (Floating() && m_value == 0)
+       {
+               return b;
        }
        
-       // Operation can be applied directly to the m_value of this and b
-       // ie: op is + or - and this and b have no * or / children
-       // or: op is * or / and this and b have no + or - children
-       if (Pure(op) && (b->Pure(op))) 
+       if (Floating() && m_value == 1 && op == MULTIPLY)
        {
-               if (ParanoidOp<digit_t>(m_value, b->m_value, op)) // op applied successfully...
-               {       
-                       Simplify(op);
-                       Simplify(invop);
-                       for (int i = 0; i < NOP; ++i) // Try applying b's children to this
+               m_value = b->m_value;
+               for (int i = 0; i < NOP; ++i)
+               {
+                       for (auto n : m_next[i])
+                               delete (n);
+                       m_next[i].clear();
+                       swap(m_next[i], b->m_next[i]);
+               }
+               //assert(SanityCheck());
+               return b;
+       }
+       if (b->Floating() && b->m_value == 1)
+               return b;
+       if (b->Floating() && b->m_value == 0 && op == MULTIPLY)
+       {
+               operator=(*b);
+               return b;
+       }
+
+               
+       if (NoTerms() && b->NoTerms())
+       {
+               if (ParanoidOp<digit_t>(m_value, b->m_value, op))
+               {
+                       Optype mulop = (op == MULTIPLY) ? MULTIPLY : DIVIDE;
+                       for (auto mul : b->m_next[MULTIPLY])
                        {
-                               delete Operation(b->m_next[i], Optype(i));
-                               b->m_next[i] = NULL;
+                               delete(OperationFactor(mul, mulop));
                        }
-                       return b; // can delete b
+                       Optype divop = (op == MULTIPLY) ? DIVIDE : MULTIPLY;
+                       for (auto div : b->m_next[DIVIDE])
+                               delete(OperationFactor(div, divop));
+                               
+                       b->m_next[DIVIDE].clear();
+                       b->m_next[MULTIPLY].clear();
+                       return b;               
                }
        }
        
-       // Try to simplify the cases:
-       // a + b*c == (a/c + b)*c
-       // a + b/c == (a*c + b)/c
-       else if ((op == ADD || op == SUBTRACT) &&
-                       (Pure(op) || b->Pure(op)))
+               
+       bool parent = (merge_point == NULL);
+       ParanoidNumber * merge = this;
+       Optype mop = op;
+       if (parent)
        {
+               merge_point = &merge;
+               merge_op = &mop;        
+       }
+       else
+       {
+               merge = *merge_point;
+               mop = *merge_op;
+       }
                
-               Debug("Simplify: {%s} %c {%s}", Str().c_str(), OpChar(op), b->Str().c_str());
-               Optype adj[] = {MULTIPLY, DIVIDE};
-               for (int i = 0; i < 2; ++i)
-               {
-
-                       Optype f = adj[i];
-                       Optype invf = InverseOp(f);
-                       
-                       Debug("Try %c", OpChar(f));
-                       
-                       if (m_next[f] == NULL && b->m_next[f] == NULL)
-                               continue;
+       Optype invop = InverseOp(op); // inverse of p
+       Optype fwd = op;
+       Optype rev = invop;
+       if (op == DIVIDE)
+       {
+               fwd = MULTIPLY;
+               rev = DIVIDE;
+       }
 
-                       ParanoidNumber * tmp_a = new ParanoidNumber(*this);
-                       ParanoidNumber * tmp_b = new ParanoidNumber(*b);
+       ParanoidNumber * cpy_b = new ParanoidNumber(*b);
+       for (auto prev : m_next[invop])
+       {
+               if (prev->OperationFactor(b, rev, merge_point, merge_op) == b)
+               {
+                       for (auto add : m_next[ADD])
+                               delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
+                       for (auto sub : m_next[SUBTRACT])
+                               delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
                                
-               
-                       ParanoidNumber * af = (tmp_a->m_next[f] != NULL) ? new ParanoidNumber(*(tmp_a->m_next[f])) : NULL;
-                       ParanoidNumber * bf = (tmp_b->m_next[f] != NULL) ? new ParanoidNumber(*(tmp_b->m_next[f])) : NULL;
-                       
-                       Debug("{%s} %c {%s}", tmp_a->Str().c_str(), OpChar(op), tmp_b->Str().c_str());
-                       Debug("{%s} %c {%s}", tmp_a->Str().c_str(), OpChar(op), tmp_b->Str().c_str());
-                       if (tmp_a->Operation(af, invf) != af || tmp_b->Operation(bf, invf) != bf)
+                       delete(cpy_b);
+                       return b;
+               }
+       }
+       for (auto next : m_next[op])
+       {
+               if (next->OperationFactor(b, fwd, merge_point, merge_op) == b)
+               {
+                       for (auto add : m_next[ADD])
                        {
-                               delete af;
-                               delete bf;
-                               delete tmp_a;
-                               delete tmp_b;
-                               continue;
+                               delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
                        }
-                       Debug("{%s} %c {%s}", tmp_a->Str().c_str(), OpChar(op), tmp_b->Str().c_str());
-                       
-                       if (tmp_a->Operation(bf, invf) == bf && tmp_b->Operation(af, invf) == af) // a / c simplifies
-                       {  
-                               if (tmp_a->Operation(tmp_b, op) != NULL) // (a/c) + b simplifies
-                               {
-                                       this->operator=(*tmp_a);
-                                       if (bf != NULL)
-                                               delete Operation(bf, f);
-                                       if (af != NULL)
-                                               delete Operation(af, f);
-                                       delete tmp_a;
-                                       delete tmp_b;
-                                       return b; // It simplified after all!
-                               }
-                               else
-                               {
-                                       tmp_b = NULL;
-                                       delete af;
-                                       delete bf;
-                               }       
+                       for (auto sub : m_next[SUBTRACT])
+                       {
+                               delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
                        }
-                       //Debug("tmp_a : %s", tmp_a->PStr().c_str());
-                       //Debug("tmp_b : %s", tmp_b->PStr().c_str());
-                       delete tmp_a;
-                       delete tmp_b;
+                       delete(cpy_b);
+                       return b;
                }
        }
        
-               // See if operation can be applied to children of this in the same dimension
+       if (parent)
+       {
+               //assert(b != NULL);
+               m_next[op].push_back(b);
+               for (auto add : m_next[ADD])
+                       delete(add->OperationFactor(new ParanoidNumber(*cpy_b), op));
+               for (auto sub : m_next[SUBTRACT])
+                       delete(sub->OperationFactor(new ParanoidNumber(*cpy_b), op));
+                       
+               #ifdef PARANOID_SIZE_LIMIT
+                       m_size += b->m_size;
+               #endif  
+       }
+       //assert(SanityCheck());
+
+
+       
+       return NULL;    
+}
+
+
+
+/**
+ * Performs the operation on a with argument b (a += b, a -= b, a *= b, a /= b)
+ * @returns b if b can safely be deleted
+ * @returns NULL if b has been merged with a
+ * append indicates that b should be merged
+ */
+ParanoidNumber * ParanoidNumber::Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
+{
+
+       if (b == NULL)
+               return NULL;
+
+       
+       if (op == SUBTRACT || op == ADD)
+               return OperationTerm(b, op, merge_point, merge_op);
+       if (op == MULTIPLY || op == DIVIDE)
+               return OperationFactor(b, op, merge_point, merge_op);
+       return b;
+}
+
+
+
+string ParanoidNumber::PStr() const
+{
+       stringstream s;
+       for (int i = 0; i < NOP; ++i)
        {
-               // (a / b) / c = a / (b*c)
-               // (a * b) / c = a * (b/c)
-               // (a / b) * c = a / (b/c)
-               // (a * b) * c = a * (b*c)
-               // (a + b) + c = a + (b+c)
-               // (a - b) + c = a - (b-c)
-               // (a + b) - c = a + (b-c)
-               // (a - b) - c = a - (b+c)
-               Optype fwd(op);
-               Optype rev(invop);
-               if (op == DIVIDE || op == SUBTRACT)
+               Optype f = Optype(i);
+               s << this;
+               for (auto n : m_next[f])
                {
-                       fwd = invop;
-                       rev = op;
+                       s << OpChar(f) << n->PStr();
                }
-               // opposite direction first (because ideally things will cancel each other out...)
-               if (m_next[invop] != NULL && m_next[invop]->Operation(b, rev, &append_at) != NULL)
-                       return b;
-               // forward direction
-               if (m_next[op] != NULL && m_next[op]->Operation(b, fwd, &append_at) != NULL) 
-                       return b;
        }
+       return s.str();
+}
+
+bool ParanoidNumber::Simplify(Optype op)
+{
        
-       // At this point, we have no choice but to merge 'b' with this ParanoidNumber
+       if (Floating())
+               return false;
+               
+       //assert(SanityCheck());
+       vector<ParanoidNumber*> next;
+       next.clear();
+       swap(m_next[op], next);
+       m_next[op].clear();
+       //assert(m_next[op].size() == 0);
+       //assert(SanityCheck());
+       Optype fwd = op;
+       if (op == DIVIDE)
+               fwd = MULTIPLY;
+       else if (op == SUBTRACT)
+               fwd = ADD;
+               
+               
+       vector<ParanoidNumber*> hold[2];
+       if (op == MULTIPLY || op == DIVIDE)
+       {
+               swap(m_next[ADD], hold[0]);
+               swap(m_next[SUBTRACT], hold[1]);
+       }
        
-       // we are a child; the merge operation needs to be applied by the root, so leave
-       if (parent != NULL) 
+       for (vector<ParanoidNumber*>::iterator n = next.begin(); n != next.end(); ++n)
        {
-               if (m_next[op] == NULL)
-                       *parent = this; // last element in list
-               return NULL;
+               if (*n == NULL)
+                       continue;
+               for (vector<ParanoidNumber*>::iterator m = n; m != next.end(); ++m)
+               {
+                       if ((*m) == (*n))
+                               continue;
+                       if (*m == NULL)
+                               continue;
+                               
+                       ParanoidNumber * parent = this;
+                       Optype mop = op;
+                       ParanoidNumber * result = (*n)->Operation(*m, fwd, &parent, &mop);
+                       if (result != NULL)
+                       {
+                               #ifdef PARANOID_SIZE_LIMIT
+                                       m_size -= result->m_size;
+                               #endif
+                               *m = NULL;
+                               delete(result);
+                       }
+               }
        }
        
-       append_at->m_next[op] = b; // Merge with b
        
-       // MULTIPLY and DIVIDE operations need to be performed on each term in the ADD/SUBTRACT dimension
-       if (op == DIVIDE || op == MULTIPLY)
+       
+       for (auto n : next)
        {
-               // apply the operation to each term
-               if (m_next[ADD] != NULL) delete m_next[ADD]->Operation(new ParanoidNumber(*b), op);
-               if (m_next[SUBTRACT] != NULL) delete m_next[SUBTRACT]->Operation(new ParanoidNumber(*b), op);
-               
-               // try and simplify this by adding the terms (you never know...)
-               Simplify(ADD);
-               Simplify(SUBTRACT);
+               if (n != NULL)
+               {               
+                       #ifdef PARANOID_SIZE_LIMIT
+                               if (Operation(n, op) == n)
+                               {
+                                       m_size -= n->m_size;
+                                       delete n;
+                               }
+                       #else   
+                               delete(Operation(n, op));
+                       #endif 
+               }
        }
-       // failed to simplify
-       return NULL;
+       
+       if (op == MULTIPLY || op == DIVIDE)
+       {
+               swap(m_next[ADD], hold[0]);
+               swap(m_next[SUBTRACT], hold[1]);
+       }
+       
+       set<ParanoidNumber*> s;
+       //if (!SanityCheck(s))
+       //{
+       //      Error("Simplify broke Sanity");
+       //}
+       return (next.size() > m_next[op].size());
 }
 
-bool ParanoidNumber::Simplify(Optype op)
+bool ParanoidNumber::FullSimplify()
+{
+       bool result = false;
+       result |= Simplify(MULTIPLY);
+       result |= Simplify(DIVIDE);
+       result |= Simplify(ADD);
+       result |= Simplify(SUBTRACT);
+       return result;
+}
+
+ParanoidNumber::digit_t ParanoidNumber::Digit() const
 {
-       ParanoidNumber * n = m_next[op];
-       m_next[op] = NULL;
-       if (Operation(n, Optype(op)))
+
+       // Get around the absurd requirement that const correctness be observed.
+       #ifdef PARANOID_CACHE_RESULTS
+       if (m_cache_valid) // le sigh ambiguous function compiler warnings
+               return m_cached_result;
+               
+       digit_t & result = ((ParanoidNumber*)(this))->m_cached_result;
+       
+       #else
+               digit_t result;
+       #endif
+       result = m_value;
+       for (auto mul : m_next[MULTIPLY])
        {
-               delete n;
-               return true;
+               result *= mul->Digit();
        }
-       else
+       for (auto div : m_next[DIVIDE])
        {
-               m_next[op] = n;
-               return false;
+               result /= div->Digit();
        }
+       for (auto add : m_next[ADD])
+               result += add->Digit();
+       for (auto sub : m_next[SUBTRACT])
+               result -= sub->Digit();
+               
+       #ifdef PARANOID_CACHE_RESULTS
+               ((ParanoidNumber*)(this))->m_cache_valid = true;
+       #endif
+       return result;
+               
 }
 
-string ParanoidNumber::PStr() const
+ParanoidNumber::digit_t ParanoidNumber::GetFactors() const
 {
-       stringstream s;
-       for (int i = 0; i < NOP; ++i)
+       digit_t value = 1;
+       for (auto mul : m_next[MULTIPLY])
+               value *= mul->Digit();
+       for (auto div : m_next[DIVIDE])
+               value /= div->Digit();
+       return value;
+}
+
+
+ParanoidNumber::digit_t ParanoidNumber::GetTerms() const
+{
+       digit_t value = 0;
+       for (auto add : m_next[ADD])
+               value += add->Digit();
+       for (auto sub : m_next[SUBTRACT])
+               value -= sub->Digit();
+       return value;
+}
+
+bool ParanoidNumber::SanityCheck(set<ParanoidNumber*> & visited) const
+{
+       if (this == NULL)
        {
-               Optype f = Optype(i);
-               s << this << OpChar(f) << m_next[f] << "\n";
+               Error("NULL pointer in tree");
+               return false;
        }
-       return s.str();
+               
+       if (visited.find((ParanoidNumber*)this) != visited.end())
+       {
+               Error("I think I've seen this tree before...");
+               return false;
+       }
+       
+       visited.insert((ParanoidNumber*)this);
+               
+       for (auto add : m_next[ADD])
+       {
+               if (!add->SanityCheck(visited))
+                       return false;
+       }
+       for (auto sub : m_next[SUBTRACT])
+       {
+               if (!sub->SanityCheck(visited))
+                       return false;
+       }
+       for (auto mul : m_next[MULTIPLY])
+       {
+               if (!mul->SanityCheck(visited))
+                       return false;
+       }
+       
+       for (auto div : m_next[DIVIDE])
+       {
+               if (!div->SanityCheck(visited))
+                       return false;
+               if (div->Digit() == 0)
+               {
+                       Error("Divide by zero");
+                       return false;
+               }
+       }
+       return true;
 }
 
+void ParanoidNumber::Negate()
+{
+       swap(m_next[ADD], m_next[SUBTRACT]);
+       m_value = -m_value;
+       #ifdef PARANOID_CACHE_RESULTS
+               m_cached_result = -m_cached_result;
+       #endif
+}
 
+#ifdef PARANOID_USE_ARENA
 
+void * ParanoidNumber::operator new(size_t s)
+{
+       return g_arena.allocate(s);
+}
 
+void ParanoidNumber::operator delete(void * p)
+{
+       g_arena.deallocate(p);
+}
+
+ParanoidNumber::Arena::Arena(int64_t block_size) : m_block_size(block_size), m_spare(NULL)
+{
+       m_blocks.push_back({malloc(block_size*sizeof(ParanoidNumber)),0});
+}
+
+ParanoidNumber::Arena::~Arena()
+{
+       for (auto block : m_blocks)
+       {
+               free(block.memory);
+       }
+       m_blocks.clear();
+}
+
+void * ParanoidNumber::Arena::allocate(size_t s)
+{
+       if (m_spare != NULL)
+       {
+               void * result = m_spare;
+               m_spare = NULL;
+               return result;
+       }
+               
+       Block & b = m_blocks.back();
+       void * result = (ParanoidNumber*)(b.memory) + (b.used++);
+       if (b.used >= m_block_size)
+       {
+               m_block_size *= 2;
+               Debug("Add block of size %d", m_block_size);
+               m_blocks.push_back({malloc(m_block_size*sizeof(ParanoidNumber)), 0});
+       }
+       return result;
+}
+
+void ParanoidNumber::Arena::deallocate(void * p)
+{
+       m_spare = p;
+}
+#endif //PARANOID_USE_ARENA
 
 }

UCC git Repository :: git.ucc.asn.au