#include <sstream>
#include <fenv.h>
#include "log.h"
+#include <iostream>
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;
{
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);
{
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());
-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)
{
- 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)
{
- 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;
}
+
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)
{
- 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();
+ //Debug("Simplified to {%s}", Str().c_str());
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;
}
-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("");
- 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:
- result += " +";
+ a += b;
break;
case SUBTRACT:
- result += " -";
+ a -= b;
break;
case MULTIPLY:
- result += " *";
+ a *= b;
break;
case DIVIDE:
- result += " /";
+ a /= b;
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;
}
}