Slightly better results
[ipdf/code.git] / src / paranoidnumber.cpp
index 2bb44da..b45c71e 100644 (file)
@@ -3,15 +3,28 @@
 #include <sstream>
 #include <fenv.h>
 #include "log.h"
+#include <cassert>
 #include <iostream>
 
 using namespace std;
 namespace IPDF
 {
-       
+int64_t ParanoidNumber::g_count = 0;
+
+
+ParanoidNumber::~ParanoidNumber()
+{
+       g_count--;
+       for (int i = 0; i < NOP; ++i)
+       {
+               for (auto n : m_next[i])
+                       delete n;
+       }
+}
 
-ParanoidNumber::ParanoidNumber(const char * str) : m_value(0), m_op(ADD), m_next_term(NULL), m_next_factor(NULL)
+ParanoidNumber::ParanoidNumber(const char * str) : m_value(0), m_cached_result(0)
 {
+       Construct();
        int dp = 0;
        int end = 0;
        while (str[dp] != '\0' && str[dp] != '.')
@@ -21,16 +34,12 @@ ParanoidNumber::ParanoidNumber(const char * str) : m_value(0), m_op(ADD), m_next
        }
        while (str[end] != '\0')
                ++end;
-               
        ParanoidNumber m(1);
        for (int i = dp-1; i >= 0; --i)
        {
                ParanoidNumber b(str[i]-'0');
                b*=m;
-               //Debug("m is %s", m.Str().c_str());
-               //Debug("Add %s", b.Str().c_str());
                this->operator+=(b);
-               //Debug("Now at %s", Str().c_str());
                m*=10;
        }
        ParanoidNumber n(1);
@@ -38,297 +47,504 @@ ParanoidNumber::ParanoidNumber(const char * str) : m_value(0), m_op(ADD), m_next
        {
                n/=10;
                ParanoidNumber b(str[i]-'0');
-               //Debug("%s * %s", b.Str().c_str(), n.Str().c_str());
                b*=n;
-               //Debug("b -> %s", b.Str().c_str());
-               //Debug("Add %s", b.Str().c_str());
                this->operator+=(b);
-               //Debug("Now at %s", Str().c_str());
-
        }
-       Debug("Constructed {%s} from %s (%f)", Str().c_str(), str, ToDouble()); 
 }
 
 ParanoidNumber & ParanoidNumber::operator=(const ParanoidNumber & a)
 {
-       //TODO: Optimise
-       delete m_next_term;
-       delete m_next_factor;
-       m_op = a.m_op;
-       if (a.m_next_term != NULL)
+       m_value = a.m_value;
+       m_cached_result = a.m_cached_result;
+       for (int i = 0; i < NOP; ++i)
        {
-               m_next_term = new ParanoidNumber(*(a.m_next_term));
-       }
-       if (a.m_next_factor != NULL)
-       {
-               m_next_factor = new ParanoidNumber(*(a.m_next_factor));
-       }
+               for (unsigned j = 0; j < m_next[i].size() && j < a.m_next[i].size(); ++j)
+               {
+                       m_next[i][j]->operator=(*(a.m_next[i][j]));
+               }
+               
+               for (unsigned j = a.m_next[i].size(); j < m_next[i].size(); ++j)
+               {
+                       delete m_next[i][j];
+               }
+               m_next[i].resize(a.m_next[i].size());
+       }       
        return *this;
 }
 
-ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
+
+string ParanoidNumber::Str() const
 {
-       // this = v + t + (a)
-       // -> v + (a) + t
-       
-       ParanoidNumber * nt = m_next_term;
-       ParanoidNumber * nf = m_next_factor;
+       string result("");
+       stringstream s;
+       s << (double)m_value;
+       result += s.str();
+       for (auto mul : m_next[MULTIPLY])
+       {
+               result += "*";
+               if (!mul->Floating())
+                       result += "(" + mul->Str() + ")";
+               else
+                       result += mul->Str();
+       }
+       for (auto div : m_next[DIVIDE])
+       {
+               result += "/";
+               if (!div->Floating())
+                       result += "(" + div->Str() + ")";
+               else
+                       result += div->Str();
+       }       
        
-       ParanoidNumber ca(a);
-       if (m_next_factor != NULL)
+       for (auto add : m_next[ADD])
        {
-               if (m_next_factor->m_op == MULTIPLY)
-                       ca /= (*m_next_factor);
+               result += "+";
+               if (!add->Floating())
+                       result += "(" + add->Str() + ")";
                else
-                       ca *= (*m_next_factor);
-                       
-               if (ca.Floating())
-               {
-                       m_next_factor = NULL;
-                       m_next_term = NULL;
-                       operator+=(ca);
-                       m_next_factor = nf;
-                       m_next_term = nt;
-                       Simplify();
-                       return *this;
-               }
-               
+                       result += add->Str();
+       }
+       for (auto sub : m_next[SUBTRACT])
+       {
+               result += "-";
+               if (!sub->Floating())
+                       result += "(" + sub->Str() + ")";
+               else
+                       result += sub->Str();
        }
        
-       m_next_term = new ParanoidNumber(a, ADD);
-       ParanoidNumber * t = m_next_term;
-       while (t->m_next_term != NULL)
-               t = t->m_next_term;
-       t->m_next_term = nt;
-       //Debug("Simplify {%s} after add", Str().c_str());
-       Simplify();
-       return *this;
+
+       return result;
 }
 
-ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
+template <>
+bool TrustingOp<float>(float & a, const float & b, Optype op)
 {
-       // this = v + t + (a)
-       // -> v + (a) + t
-       
-       ParanoidNumber * nt = m_next_term;
-       ParanoidNumber * nf = m_next_factor;
-       
-       ParanoidNumber ca(a, SUBTRACT);
-       if (m_next_factor != NULL)
+       feclearexcept(FE_ALL_EXCEPT);
+       switch (op)
        {
-               if (m_next_factor->m_op == MULTIPLY)
-                       ca /= (*m_next_factor);
-               else
-                       ca *= (*m_next_factor);
-                       
-               if (ca.Floating())
-               {
-                       m_next_factor = NULL;
-                       m_next_term = NULL;
-                       operator-=(ca);
-                       m_next_factor = nf;
-                       m_next_term = nt;
-                       Simplify();
-                       return *this;
-               }
-               
+               case ADD:
+                       a += b;
+                       break;
+               case SUBTRACT:
+                       a -= b;
+                       break;
+               case MULTIPLY:
+                       a *= b;
+                       break;
+               case DIVIDE:
+                       a /= b;
+                       break;
+               case NOP:
+                       break;
        }
-       
-       m_next_term = new ParanoidNumber(a,SUBTRACT);
-       ParanoidNumber * t = m_next_term;
-       while (t->m_next_term != NULL)
-       {
-               t->m_op = SUBTRACT;
-               t = t->m_next_term;
-       }
-       t->m_op = SUBTRACT;
-       //Debug("next term {%s}", m_next_term->Str().c_str());
-       t->m_next_term = nt;
-       //Debug("Simplify {%s} after sub", Str().c_str());
-       Simplify();
+       return !fetestexcept(FE_ALL_EXCEPT);
+}
+
+template <>
+bool TrustingOp<double>(double & a, const 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:
+                       a /= b;
+                       break;
+               case NOP:
+                       break;
+       }
+       return !fetestexcept(FE_ALL_EXCEPT);
+}
+
+template <>
+bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op)
+{
+       int16_t sa(a);
+       bool exact = true;
+       switch (op)
+       {
+               case ADD:
+                       sa += b;
+                       exact = (abs(sa) <= 127);
+                       break;
+               case SUBTRACT:
+                       sa -= b;
+                       exact = (abs(sa) <= 127);
+                       break;
+               case MULTIPLY:
+                       sa *= b;
+                       exact = (abs(sa) <= 127);
+                       break;
+               case DIVIDE:
+                       exact = (b != 0 && sa > b && sa % b == 0);
+                       sa /= b;
+                       break;
+               case NOP:
+                       break;
+       }
+       a = (int8_t)(sa);
+       return exact;
+}
+
+
+ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
+{
+       delete Operation(new ParanoidNumber(a), ADD);
+       Simplify(ADD);
+       Simplify(SUBTRACT);
+       return *this;
+}
+
+
+ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
+{
+       delete Operation(new ParanoidNumber(a), SUBTRACT);
+       //Simplify(SUBTRACT);
+       //Simplify(ADD);
        return *this;
 }
 
 ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
 {
-       Debug("{%s} *= {%s}", Str().c_str(), a.Str().c_str());
-       // this = (vf + t) * (a)
-       ParanoidNumber * nf = m_next_factor;
-       m_next_factor = new ParanoidNumber(a, MULTIPLY);
-       ParanoidNumber * t = m_next_factor;
-       while (t->m_next_factor != NULL)
-               t = t->m_next_factor;
-       t->m_next_factor = nf;
-       if (m_next_term != NULL)
-               m_next_term->operator*=(a);
-       //Debug("Simplify after mul");
-       Debug("Simplify {%s}", Str().c_str());
-       Simplify();
+       delete Operation(new ParanoidNumber(a), MULTIPLY);
        return *this;
 }
 
 
 ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
 {
-       //Debug("Called %s /= %s", Str().c_str(), a.Str().c_str());
-       // this = (vf + t) * (a)
-       ParanoidNumber * nf = m_next_factor;
-       m_next_factor = new ParanoidNumber(a, DIVIDE);
-       ParanoidNumber * t = m_next_factor;
-       while (t->m_next_factor != NULL)
-               t = t->m_next_factor;
-       t->m_next_factor = nf;
-       if (m_next_term != NULL)
-               m_next_term->operator/=(a);
-       Simplify();
+       delete Operation(new ParanoidNumber(a), DIVIDE);
        return *this;
 }
 
-
-
-void ParanoidNumber::SimplifyTerms()
-{ 
-       //Debug("Simplify {%s}", Str().c_str()); 
-       if (m_next_term == NULL)
+// a + b
+ParanoidNumber * ParanoidNumber::OperationTerm(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
+{
+       m_cached_result = nan("");
+       if (Floating() && m_value == 0) // 0 + b = b
        {
-               //Debug("No terms!");
-               return;
+               m_value = b->m_value;
+               if (op == SUBTRACT)
+               {
+                       m_value = -m_value;
+                       swap(b->m_next[ADD], b->m_next[SUBTRACT]);
+               }
+               
+               for (int i = 0; i < NOP; ++i)
+               {
+                       m_next[i] = b->m_next[i];
+                       b->m_next[i].clear();
+               }
+               return b;
        }
+       if (b->Floating() && b->m_value == 0) // a + 0 = a
+               return b;
+               
 
-       for (ParanoidNumber * a = this; a != NULL; a = a->m_next_term)
+       
+       if ((NoFactors() && b->NoFactors())
+               || (GetFactors() == b->GetFactors()))
        {
-               ParanoidNumber * bprev = a;
-               ParanoidNumber * b = a->m_next_term;
-               while (b != NULL)
+               if (ParanoidOp<digit_t>(m_value, b->m_value, op))
                {
-                       //Debug("Simplify factors of %s", b->Str().c_str());
-                       b->SimplifyFactors();
-                       if (b->m_next_factor != NULL)
+                       Optype addop = (op == ADD) ? ADD : SUBTRACT;
+                       for (auto add : b->m_next[ADD])
                        {
-                               bprev = b;
-                               b = b->m_next_term;
-                               continue;
+                               delete OperationTerm(add, addop);
                        }
-                       float f = a->m_value;
-                       feclearexcept(FE_ALL_EXCEPT);
-                       switch (b->m_op)
-                       {
-                               case ADD:
-                                       f += b->m_value;
-                                       break;
-                               case SUBTRACT:
-                                       f -= b->m_value;
-                                       break;
-                               default:
-                                       Fatal("Unexpected %c in term list...", OpChar(b->m_op));
-                                       break;
-                       }
-                       if (!fetestexcept(FE_ALL_EXCEPT))
-                       {
-                               a->m_value = f;
-                               bprev->m_next_term = b->m_next_term;
-                               b->m_next_term = NULL;
-                               delete b;
-                               b = bprev;
-                       }
-                       bprev = b;
-                       b = b->m_next_term;
+                       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();
+                       return b;
                }
        }
+
+       
+       
+       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)
+                       return b;
+               
+       }
+       for (auto next : m_next[op])
+       {
+               if (next->OperationTerm(b, fwd, merge_point, merge_op) == b)
+                       return b;
+       }
+       
+
+       
+       
+       if (parent)
+       {
+               //merge->m_next[*merge_op].push_back(b);
+               m_next[op].push_back(b);
+       }
+       else
+       {
+               if (m_next[op].size() == 0)
+               {
+                       *merge_point = this;
+                       *merge_op = op;
+               }
+       }
+       return NULL;
 }
 
-void ParanoidNumber::SimplifyFactors()
-{ 
-       //Debug("Simplify {%s}", Str().c_str()); 
-       if (m_next_factor == NULL)
+ParanoidNumber * ParanoidNumber::OperationFactor(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
+{
+       m_cached_result = nan("");
+       if (Floating() && m_value == 0)
+       {
+               return b;
+       }
+       
+       if (Floating() && m_value == 1 && op == MULTIPLY)
        {
-               //Debug("No factors!");
-               return;
+               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]);
+               }
+               return b;
        }
+       if (b->Floating() && b->m_value == 1)
+               return b;
+       
 
-       for (ParanoidNumber * a = this; a != NULL; a = a->m_next_factor)
+               
+       if (NoTerms() && b->NoTerms())
        {
-               ParanoidNumber * bprev = a;
-               ParanoidNumber * b = a->m_next_factor;
-               while (b != NULL)
+               if (ParanoidOp<digit_t>(m_value, b->m_value, op))
                {
-                       b->SimplifyTerms();
-                       if (b->m_next_term != NULL)
-                       {
-                               bprev = b;
-                               b = b->m_next_factor;
-                               continue;
-                       }
-                       float f = a->m_value;
-                       feclearexcept(FE_ALL_EXCEPT);
-                       switch (b->m_op)
+                       Optype mulop = (op == MULTIPLY) ? MULTIPLY : DIVIDE;
+                       for (auto mul : b->m_next[MULTIPLY])
                        {
-                               case MULTIPLY:
-                                       if (a->m_op != DIVIDE) 
-                                               f *= b->m_value;
-                                       else
-                                               f /= b->m_value;
-                                       break;
-                               case DIVIDE:
-                                       if (a->m_op != DIVIDE)
-                                               f /= b->m_value;
-                                       else
-                                               f *= b->m_value;
-                                       break;
-                               default:
-                                       Fatal("Unexpected %c in factor list...",OpChar(b->m_op));
-                                       break;
+                               delete OperationFactor(mul, mulop);
                        }
-                       if (!fetestexcept(FE_ALL_EXCEPT))
-                       {
+                       Optype divop = (op == MULTIPLY) ? DIVIDE : MULTIPLY;
+                       for (auto div : b->m_next[DIVIDE])
+                               delete OperationFactor(div, divop);
                                
-                               a->m_value = f;
-                               bprev->m_next_factor = b->m_next_factor;
-                               b->m_next_factor = NULL;
-                               delete b;
-                               b = bprev;
-                       }
-                       //else
-                               //Debug("Failed to simplify %f %c %f", a->m_value, OpChar(b->m_op), b->m_value);
-                       bprev = b;
-                       b = b->m_next_factor;
+                       b->m_next[DIVIDE].clear();
+                       b->m_next[MULTIPLY].clear();
+                       
+                       
+                       
+                       return b;               
                }
        }
+       
+               
+       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;
+       }
+               
+       Optype invop = InverseOp(op); // inverse of p
+       Optype fwd = op;
+       Optype rev = invop;
+       if (op == DIVIDE)
+       {
+               fwd = MULTIPLY;
+               rev = DIVIDE;
+       }
+
+       ParanoidNumber * cpy_b = NULL;
+       
+       if (m_next[ADD].size() > 0 || m_next[SUBTRACT].size() > 0)
+       {
+               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);
+                               
+                       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 add->OperationFactor(new ParanoidNumber(*cpy_b), op);
+                       for (auto sub : m_next[SUBTRACT])
+                               delete sub->OperationFactor(new ParanoidNumber(*cpy_b), op);
+                       delete cpy_b;
+                       return b;
+               }
+       }
+       
+       if (parent)
+       {
+               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);
+       }
+       return NULL;    
 }
 
-void ParanoidNumber::Simplify()
+
+
+/**
+ * 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)
 {
-       SimplifyFactors();
-       SimplifyTerms();
+
+       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::Str() const
+
+
+string ParanoidNumber::PStr() const
 {
-       string result("");
        stringstream s;
-       
-       s << m_value;
-       
-       if (m_next_factor != NULL)
+       for (int i = 0; i < NOP; ++i)
        {
-               result += OpChar(m_op);
-               result += "(";
-               result += s.str();
-               result += m_next_factor->Str();
-               result += ")";
+               Optype f = Optype(i);
+               s << this;
+               for (auto n : m_next[f])
+               {
+                       s << OpChar(f) << n->PStr();
+               }
        }
-       else
+       return s.str();
+}
+
+bool ParanoidNumber::Simplify(Optype op)
+{
+       vector<ParanoidNumber*> next(0);
+       swap(m_next[op], next);
+       for (auto n : next)
        {
-               result += OpChar(m_op);
-               result += s.str();
+               delete Operation(n, op);
        }
-               
-       if (m_next_term != NULL)
+       return (next.size() > m_next[op].size());
+}
+
+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()
+{
+       if (!isnan(m_cached_result))
+               return m_cached_result;
+       m_cached_result = m_value;
+       for (auto mul : m_next[MULTIPLY])
        {
-               result += " ";
-               result += m_next_term->Str();
+               m_cached_result *= mul->Digit();
        }
-       return result;
+       for (auto div : m_next[DIVIDE])
+       {
+               m_cached_result /= div->Digit();
+       }
+       for (auto add : m_next[ADD])
+               m_cached_result += add->Digit();
+       for (auto sub : m_next[SUBTRACT])
+               m_cached_result -= sub->Digit();
+       return m_cached_result;
+               
+}
+
+ParanoidNumber::digit_t ParanoidNumber::GetFactors()
+{
+       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()
+{
+       digit_t value = 0;
+       for (auto add : m_next[ADD])
+               value += add->Digit();
+       for (auto sub : m_next[SUBTRACT])
+               value -= sub->Digit();
+       return value;
+}
+
+
 }

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