Add quadtree back to the Makefile

[ipdf/code.git] / src / paranoidnumber.cpp

diff --git a/src/paranoidnumber.cpp b/src/paranoidnumber.cpp
index 2510f66..53729a3 100644 (file)
--- a/src/paranoidnumber.cpp
+++ b/src/paranoidnumber.cpp
@@ -3,13 +3,14 @@
 #include <sstream>
 #include <fenv.h>
 #include "log.h"
 #include <sstream>
 #include <fenv.h>
 #include "log.h"

+#include <iostream>

 
 using namespace std;
 namespace IPDF
 {
 
 using namespace std;
 namespace IPDF
 {

-       
+int64_t ParanoidNumber::g_count = 0;

 
 

-ParanoidNumber::ParanoidNumber(const char * str) : m_value(0), m_op(ADD), m_next(NULL)
+ParanoidNumber::ParanoidNumber(const char * str) : m_value(0), m_op(ADD), m_next_term(NULL), m_next_factor(NULL)

 {
        int dp = 0;
        int end = 0;
 {
        int dp = 0;
        int end = 0;


@@ -26,7 +27,10 @@ ParanoidNumber::ParanoidNumber(const char * str) : m_value(0), m_op(ADD), m_next
        {
                ParanoidNumber b(str[i]-'0');
                b*=m;
        {
                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);
                this->operator+=(b);

+               //Debug("Now at %s", Str().c_str());

                m*=10;
        }
        ParanoidNumber n(1);
                m*=10;
        }
        ParanoidNumber n(1);


@@ -34,266 +38,414 @@ ParanoidNumber::ParanoidNumber(const char * str) : m_value(0), m_op(ADD), m_next
        {
                n/=10;
                ParanoidNumber b(str[i]-'0');
        {
                n/=10;
                ParanoidNumber b(str[i]-'0');

+               //Debug("%s * %s", b.Str().c_str(), n.Str().c_str());

                b*=n;
                b*=n;

+               //Debug("b -> %s", b.Str().c_str());
+               //Debug("Add %s", b.Str().c_str());

                this->operator+=(b);
                this->operator+=(b);

-       }
-       
-}
+               //Debug("Now at %s", Str().c_str());

 
 

-ParanoidNumber * ParanoidNumber::InsertAfter(float value, Optype op)
-{
-       return InsertAfter(new ParanoidNumber(value, op));
-}
-       
-ParanoidNumber * ParanoidNumber::InsertAfter(ParanoidNumber * insert)
-{
-       //Debug("Insert {%s} after {%f, %s}",insert->Str().c_str(),
-       //      m_value, g_opstr[m_op]);
-       ParanoidNumber * n = m_next;
-       m_next = insert;
-       
-       ParanoidNumber * p = m_next;
-       while (p->m_next != NULL)
-               p = p->m_next;
-       p->m_next = n;
-       
-       return m_next;
+       }
+       //Debug("Constructed {%s} from %s (%f)", Str().c_str(), str, ToDouble());       

 }
 
 }
 

-

 ParanoidNumber & ParanoidNumber::operator=(const ParanoidNumber & a)
 {
 ParanoidNumber & ParanoidNumber::operator=(const ParanoidNumber & a)
 {

-       const ParanoidNumber * na = &a;
-       ParanoidNumber * nb = this;
-       ParanoidNumber * p = NULL;
-       while (na != NULL && nb != NULL)
+       //TODO: Optimise
+       delete m_next_term;
+       delete m_next_factor;
+       m_op = a.m_op;
+       if (a.m_next_term != NULL)

        {
        {

-               nb->m_value = na->m_value;
-               nb->m_op = na->m_op;
-               na = na->m_next;
-               p = nb;
-               nb = nb->m_next;
-               
-       }       
-       
-       while (na != NULL) // => nb == NULL
-       {
-               InsertAfter(na->m_value, na->m_op);
-               na = na->m_next;
+               m_next_term = new ParanoidNumber(*(a.m_next_term));

        }
        }

-
-       if (nb != NULL)
+       if (a.m_next_factor != NULL)

        {
        {

-               if (p != NULL)
-                       p->m_next = NULL;
-               delete nb;
+               m_next_factor = new ParanoidNumber(*(a.m_next_factor));

        }
        return *this;
 }
 
 ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
 {
        }
        return *this;
 }
 
 ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
 {

-       ParanoidNumber * insert = new ParanoidNumber(a, ADD);
-       if (m_next == NULL || m_next->m_op == ADD || m_next->m_op == SUBTRACT)
+       if (m_next_factor == NULL && a.Floating())

        {
        {

-               InsertAfter(insert);
-               Simplify();
-               return *this;
+               if (ParanoidOp<digit_t>(m_value, a.m_value, ADD))
+               {
+                       Simplify();
+                       return *this;
+               }

        }
        }

+       ParanoidNumber * nt = m_next_term;
+       ParanoidNumber * nf = m_next_factor;

        
        

-       if (m_next->m_op == MULTIPLY) // (a*b) + c == (a+[c/b]) * b
-               insert->operator/=(*m_next);
-       else
-               insert->operator*=(*m_next); // (a/b) + c == (a+[c*b])/b
-       
-       if (insert->m_next != NULL) // neither of the above simplified
+       ParanoidNumber ca(a);
+       if (m_next_factor != NULL)

        {
        {

-               //Debug("{%s} did not simplify, change back to {%s}", insert->Str().c_str(), a.Str().c_str());
-               insert->operator=(a); // Just add as is
-               insert->m_op = ADD;
-               ParanoidNumber * n = this;
-               while (n->m_next != NULL)
-                       n = n->m_next;
-               n->InsertAfter(insert);
-       }
-       else
-       {
-               InsertAfter(insert);
+               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;
+               }
+               

        }
        }

+       
+       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;
 }
        Simplify();
        return *this;
 }

+

 ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
 {
 ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
 {

-       ParanoidNumber * insert = new ParanoidNumber(a, SUBTRACT);
-       if (m_next == NULL || m_next->m_op == ADD || m_next->m_op == SUBTRACT)
+       // this = v + t + (a)
+       // -> v + (a) + t
+       if (m_next_factor == NULL && a.Floating())

        {
        {

-               InsertAfter(insert);
-               Simplify();
-               return *this;
+               if (ParanoidOp<digit_t>(m_value, a.m_value, ADD))
+               {
+                       Simplify();
+                       return *this;
+               }

        }
        }

+
+       ParanoidNumber * nt = m_next_term;
+       ParanoidNumber * nf = m_next_factor;

        
        

-       if (m_next->m_op == MULTIPLY) // (a*b) - c == (a-[c/b]) * b
-               insert->operator/=(*m_next);
-       else
-               insert->operator*=(*m_next); // (a/b) - c == (a-[c*b])/b
-       
-       if (insert->m_next != NULL) // neither of the above simplified
+       ParanoidNumber ca(a, SUBTRACT);
+       if (m_next_factor != NULL)

        {
        {

-               //Debug("{%s} did not simplify, change back to {%s}", insert->Str().c_str(), a.Str().c_str());
-               insert->operator=(a); // Just add as is
-               insert->m_op = SUBTRACT;
-               ParanoidNumber * n = this;
-               while (n->m_next != NULL)
-                       n = n->m_next;
-               n->InsertAfter(insert);
-       }       
-       else
+               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;
+               }
+               
+       }
+       
+       m_next_term = new ParanoidNumber(a,SUBTRACT);
+       ParanoidNumber * t = m_next_term;
+       while (t->m_next_term != NULL)

        {
        {

-               InsertAfter(insert);
+               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 *this;
 }
 
 ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
 {
        Simplify();
        return *this;
 }
 
 ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
 {

-       ParanoidNumber * n = m_next;
-       while (n != NULL)
+
+       //if (m_value == 0)
+       //              return *this;
+       //Debug("{%s} *= {%s}", Str().c_str(), a.Str().c_str());
+       // this = (vf + t) * (a)
+       if (a.Floating() && ParanoidOp<digit_t>(m_value, a.m_value, MULTIPLY))

        {
        {

-               if (n->m_op == ADD || n->m_op == SUBTRACT)
-               {
-                       n->operator*=(a);
-                       break;
-               }
-               n = n->m_next;
+               if (m_next_term != NULL)
+                       m_next_term->operator*=(a);
+               Simplify();
+               return *this;

        }
        }

-               
-       InsertAfter(new ParanoidNumber(a, MULTIPLY));
+       
+       ParanoidNumber * t = this;
+       while (t->m_next_factor != NULL)
+               t = t->m_next_factor;
+       t->m_next_factor = new ParanoidNumber(a, MULTIPLY);
+
+       if (m_next_term != NULL)
+               m_next_term->operator*=(a);
+
+       //Debug("Simplify {%s}", Str().c_str());

        Simplify();
        Simplify();

+       //Debug("Simplified to {%s}", Str().c_str());

        return *this;
 }
 
 
 ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
 {
        return *this;
 }
 
 
 ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
 {

-       ParanoidNumber * n = m_next;
-       while (n != NULL)
+               
+
+               
+       if (a.Floating() && ParanoidOp<digit_t>(m_value, a.m_value, DIVIDE))

        {
        {

-               if (n->m_op == ADD || n->m_op == SUBTRACT)
-               {
-                       n->operator/=(a);
-                       break;
-               }
-               n = n->m_next;
+               if (m_next_term != NULL)
+                       m_next_term->operator/=(a);
+               Simplify();
+               return *this;

        }
        }

-               
-       InsertAfter(new ParanoidNumber(a, DIVIDE));
+       
+       //Debug("Called %s /= %s", Str().c_str(), a.Str().c_str());
+       // this = (vf + t) * (a)
+       ParanoidNumber * t = this;
+       while (t->m_next_factor != NULL)
+       {
+               t = t->m_next_factor;
+       }
+       t->m_next_factor = new ParanoidNumber(a, DIVIDE);
+
+       if (m_next_term != NULL)
+               m_next_term->operator/=(a);
+

        Simplify();
        return *this;
 }
 
        Simplify();
        return *this;
 }
 

-double ParanoidNumber::ToDouble() const
-{
-       double value = (m_op == SUBTRACT) ? -m_value : m_value;
-       const ParanoidNumber * n = m_next;
-       while (n != NULL)
+
+
+void ParanoidNumber::SimplifyTerms()
+{ 
+
+       //Debug("Simplify {%s}", Str().c_str()); 
+       if (m_next_term == NULL)

        {
        {

-               switch (n->m_op)
-               {
-                       case ADD:
-                       case SUBTRACT:
-                               return value + n->ToDouble();
-                               break;
-                       case MULTIPLY:
-                               value *= n->m_value;
-                               break;
-                       case DIVIDE:
-                               value /= n->m_value;
-                               break;
-               }
-               n = n->m_next;
+               //Debug("No terms!");
+               return;

        }
        }

-       return value;
-}

 
 

-void ParanoidNumber::Simplify()
-{
-       ParanoidNumber * n = m_next;
-       ParanoidNumber * p = this;
-       
-       while (n != NULL)
+       for (ParanoidNumber * a = this; a != NULL; a = a->m_next_term)

        {
        {

-               
-               float a = p->m_value;
-               switch (n->m_op)
+               ParanoidNumber * b = a->m_next_term;
+               if (a->m_next_factor != NULL && !a->m_next_factor->Floating())

                {
                {

-                       case ADD:
-                               n->Simplify();
-                               feclearexcept(FE_ALL_EXCEPT);
-                               a += n->m_value;
-                               break;
-                       case SUBTRACT:
-                               n->Simplify();
-                               feclearexcept(FE_ALL_EXCEPT);
-                               a -= n->m_value;
-                               break;
-                       case MULTIPLY:
-                               feclearexcept(FE_ALL_EXCEPT);
-                               a *= n->m_value;
-                               break;
-                       case DIVIDE:
-                               feclearexcept(FE_ALL_EXCEPT);
-                               a /= n->m_value;
-                               break;
-                               
+                       continue;

                }
                }

-               // can't merge p and n
-               if (fetestexcept(FE_ALL_EXCEPT))
+               
+               ParanoidNumber * bprev = a;
+               while (b != NULL)

                {
                {

-                       while (n != NULL && n->m_op != ADD && n->m_op != SUBTRACT)
-                               n = n->m_next;
-                       if (n != NULL)
-                               n->Simplify();
-                       return;
+                       //Debug("Simplify factors of %s", b->Str().c_str());
+                       b->SimplifyFactors();
+                       if (b->m_next_factor != NULL && !b->m_next_factor->Floating())
+                       {
+                               bprev = b;
+                               b = b->m_next_term;
+                               continue;
+                       }
+                       
+                       bool simplify = false;
+                       if (a->m_next_factor != NULL || b->m_next_factor != NULL)
+                       {
+                               digit_t aa(a->Head<digit_t>());
+                               digit_t ab = (a->m_next_factor != NULL) ? a->m_next_factor->Head<digit_t>() : 1;
+                               digit_t bc(b->Head<digit_t>());
+                               digit_t bd = (b->m_next_factor != NULL) ? b->m_next_factor->Head<digit_t>() : 1;
+                               Optype aop = (a->m_next_factor != NULL) ? a->m_next_factor->m_op : DIVIDE;
+                               Optype cop = (b->m_next_factor != NULL) ? b->m_next_factor->m_op : DIVIDE;
+                               simplify = CombineTerms<digit_t>(aa, aop, ab, bc, cop, bd);
+                               if (simplify)
+                               {
+                                       a->m_value = aa;
+                                       if (a->m_next_factor != NULL)
+                                               a->m_next_factor->m_value = ab;
+                                       else if (ab != 1)
+                                       {
+                                               a->m_next_factor = b->m_next_factor;
+                                               b->m_next_factor = NULL;
+                                               a->m_next_factor->m_value = ab;
+                                       }
+                               }
+                       }
+                       else
+                       {
+                               simplify = ParanoidOp<digit_t>(a->m_value, b->Head<digit_t>(), ADD);
+                       }
+                       if (simplify)
+                       {
+                               bprev->m_next_term = b->m_next_term;
+                               b->m_next_term = NULL;
+                               delete b;
+                               b = bprev;
+                       }
+                       
+                       bprev = b;
+                       b = b->m_next_term;

                }
                }

-               else
+       }
+}
+
+void ParanoidNumber::SimplifyFactors()
+{ 
+
+       //Debug("Simplify {%s}", Str().c_str()); 
+       if (m_next_factor == NULL)
+       {
+               //Debug("No factors!");
+               return;
+       }
+
+       for (ParanoidNumber * a = this; a != NULL; a = a->m_next_factor)
+       {
+               if ((a->m_op != ADD || a->m_op != SUBTRACT) && a->m_next_term != NULL)
+                       continue;
+                       
+               ParanoidNumber * bprev = a;
+               ParanoidNumber * b = a->m_next_factor;
+               while (b != NULL)

                {
                {

-                       //  merge n into p
-                       p->m_value = a;
-                       p->m_next = n->m_next;
-                       n->m_next = NULL;
-                       delete n;
-                       n = p->m_next;          
+                       b->SimplifyTerms();
+                       if (b->m_next_term != NULL)
+                       {
+                               bprev = b;
+                               b = b->m_next_factor;
+                               continue;
+                       }
+               
+                       Optype op = b->m_op;
+                       if (a->m_op == DIVIDE)
+                       {
+                               op = (b->m_op == DIVIDE) ? MULTIPLY : DIVIDE;
+                       }
+                       
+                       if (ParanoidOp<digit_t>(a->m_value, b->m_value, op))
+                       {       
+
+                               bprev->m_next_factor = b->m_next_factor;
+                               b->m_next_factor = NULL;
+                               delete b;
+                               b = bprev;
+                       }
+                       bprev = b;
+                       b = b->m_next_factor;

                }
                }

-               //Debug("  -> {%s}", Str().c_str());

        }
 }
 
        }
 }
 

+void ParanoidNumber::Simplify()
+{
+       SimplifyFactors();
+       SimplifyTerms();
+}
+

 string ParanoidNumber::Str() const
 {
        string result("");
 string ParanoidNumber::Str() const
 {
        string result("");

-       switch (m_op)
+       stringstream s;
+       s << (double)m_value;
+       
+       if (m_next_factor != NULL)
+       {
+               result += s.str();
+               result += OpChar(m_next_factor->m_op);
+               if (m_next_factor->m_next_term != NULL)
+                       result += "(" + m_next_factor->Str() + ")";
+               else
+                       result += m_next_factor->Str();
+       }
+       else
+       {
+               result += s.str();
+       }
+               
+       if (m_next_term != NULL)
+       {
+               result += " ";
+               result += OpChar(m_next_term->m_op);
+               result += m_next_term->Str();
+       }
+       return result;
+}
+
+template <>
+bool TrustingOp<float>(float & a, const float & b, Optype op)
+{
+       feclearexcept(FE_ALL_EXCEPT);
+       switch (op)

        {
                case ADD:
        {
                case ADD:

-                       result += " +";
+                       a += b;

                        break;
                case SUBTRACT:
                        break;
                case SUBTRACT:

-                       result += " -";
+                       a -= b;

                        break;
                case MULTIPLY:
                        break;
                case MULTIPLY:

-                       result += " *";
+                       a *= b;

                        break;
                case DIVIDE:
                        break;
                case DIVIDE:

-                       result += " /";
+                       a /= b;

                        break;
        }
                        break;
        }

-       stringstream s("");
-       s << m_value;
-       result += s.str();
-       if (m_next != NULL)
-               result += m_next->Str();
-       return result;
+       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;
+       }
+       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;
+       }
+       a = (int8_t)(sa);
+       return exact;

 }
 
 }
 }
 
 }