From: Sam Moore Date: Tue, 16 Sep 2014 10:46:43 +0000 (+0800) Subject: Merge branch 'master' of git.ucc.asn.au:/ipdf/code X-Git-Url: https://git.ucc.asn.au/?p=ipdf%2Fcode.git;a=commitdiff_plain;h=77137590512d969da2d54d9ba53d76836a290c6a;hp=888817a67a9d840be66b52811b01eb77f10ff3e6 Merge branch 'master' of git.ucc.asn.au:/ipdf/code sulix fixed the quad trees and I forgot about `git stash` again. --- diff --git a/src/Makefile b/src/Makefile index 668f0a9..38f69fd 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,7 +1,7 @@ #Makefile ARCH := $(shell uname -m) # TODO: stb_truetype doesn't compile with some of these warnings. -CXX = g++ -std=gnu++0x -g -Wall -Werror -Wshadow -pedantic -rdynamic +CXX = g++ -std=c++11 -g -Wall -Werror -Wshadow -pedantic -rdynamic MAIN = main.o OBJ = log.o real.o bezier.o document.o objectrenderer.o view.o screen.o graphicsbuffer.o framebuffer.o shaderprogram.o stb_truetype.o gl_core44.o path.o paranoidnumber.o quadtree.o diff --git a/src/paranoidnumber.cpp b/src/paranoidnumber.cpp index 53729a3..b45c71e 100644 --- a/src/paranoidnumber.cpp +++ b/src/paranoidnumber.cpp @@ -3,6 +3,7 @@ #include #include #include "log.h" +#include #include using namespace std; @@ -10,8 +11,20 @@ namespace IPDF { int64_t ParanoidNumber::g_count = 0; -ParanoidNumber::ParanoidNumber(const char * str) : m_value(0), m_op(ADD), m_next_term(NULL), m_next_factor(NULL) + +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_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,341 +47,73 @@ 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_next_term = new ParanoidNumber(*(a.m_next_term)); - } - if (a.m_next_factor != NULL) + m_value = a.m_value; + m_cached_result = a.m_cached_result; + for (int i = 0; i < NOP; ++i) { - m_next_factor = new ParanoidNumber(*(a.m_next_factor)); - } - return *this; -} - -ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a) -{ - if (m_next_factor == NULL && a.Floating()) - { - if (ParanoidOp(m_value, a.m_value, ADD)) + for (unsigned j = 0; j < m_next[i].size() && j < a.m_next[i].size(); ++j) { - Simplify(); - return *this; - } - } - ParanoidNumber * nt = m_next_term; - ParanoidNumber * nf = m_next_factor; - - ParanoidNumber ca(a); - if (m_next_factor != NULL) - { - 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[i][j]->operator=(*(a.m_next[i][j])); } - } - - 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) -{ - // this = v + t + (a) - // -> v + (a) + t - if (m_next_factor == NULL && a.Floating()) - { - if (ParanoidOp(m_value, a.m_value, ADD)) + for (unsigned j = a.m_next[i].size(); j < m_next[i].size(); ++j) { - Simplify(); - return *this; + delete m_next[i][j]; } - } - - ParanoidNumber * nt = m_next_term; - ParanoidNumber * nf = m_next_factor; - - ParanoidNumber ca(a, SUBTRACT); - if (m_next_factor != NULL) - { - 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) - { - 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) -{ - - //if (m_value == 0) - // return *this; - //Debug("{%s} *= {%s}", Str().c_str(), a.Str().c_str()); - // this = (vf + t) * (a) - if (a.Floating() && ParanoidOp(m_value, a.m_value, MULTIPLY)) - { - if (m_next_term != NULL) - m_next_term->operator*=(a); - Simplify(); - return *this; - } - - 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) -{ - - - - if (a.Floating() && ParanoidOp(m_value, a.m_value, DIVIDE)) - { - if (m_next_term != NULL) - m_next_term->operator/=(a); - Simplify(); - return *this; - } - - //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(); + m_next[i].resize(a.m_next[i].size()); + } return *this; } - -void ParanoidNumber::SimplifyTerms() -{ - - //Debug("Simplify {%s}", Str().c_str()); - if (m_next_term == NULL) - { - //Debug("No terms!"); - return; - } - - for (ParanoidNumber * a = this; a != NULL; a = a->m_next_term) - { - ParanoidNumber * b = a->m_next_term; - if (a->m_next_factor != NULL && !a->m_next_factor->Floating()) - { - continue; - } - - ParanoidNumber * bprev = a; - while (b != NULL) - { - //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 ab = (a->m_next_factor != NULL) ? a->m_next_factor->Head() : 1; - digit_t bc(b->Head()); - digit_t bd = (b->m_next_factor != NULL) ? b->m_next_factor->Head() : 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(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(a->m_value, b->Head(), 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; - } - } -} - -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) - { - 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(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; - } - } -} - -void ParanoidNumber::Simplify() -{ - SimplifyFactors(); - SimplifyTerms(); -} - string ParanoidNumber::Str() const { string result(""); stringstream s; s << (double)m_value; - - if (m_next_factor != NULL) + result += s.str(); + for (auto mul : m_next[MULTIPLY]) { - result += s.str(); - result += OpChar(m_next_factor->m_op); - if (m_next_factor->m_next_term != NULL) - result += "(" + m_next_factor->Str() + ")"; + result += "*"; + if (!mul->Floating()) + result += "(" + mul->Str() + ")"; else - result += m_next_factor->Str(); + result += mul->Str(); } - else + for (auto div : m_next[DIVIDE]) + { + result += "/"; + if (!div->Floating()) + result += "(" + div->Str() + ")"; + else + result += div->Str(); + } + + for (auto add : m_next[ADD]) { - result += s.str(); + result += "+"; + if (!add->Floating()) + result += "(" + add->Str() + ")"; + else + result += add->Str(); } - - if (m_next_term != NULL) + for (auto sub : m_next[SUBTRACT]) { - result += " "; - result += OpChar(m_next_term->m_op); - result += m_next_term->Str(); + result += "-"; + if (!sub->Floating()) + result += "(" + sub->Str() + ")"; + else + result += sub->Str(); } + + return result; } @@ -394,6 +135,8 @@ bool TrustingOp(float & a, const float & b, Optype op) case DIVIDE: a /= b; break; + case NOP: + break; } return !fetestexcept(FE_ALL_EXCEPT); } @@ -416,6 +159,8 @@ bool TrustingOp(double & a, const double & b, Optype op) case DIVIDE: a /= b; break; + case NOP: + break; } return !fetestexcept(FE_ALL_EXCEPT); } @@ -443,9 +188,363 @@ bool TrustingOp(int8_t & a, const int8_t & b, Optype op) 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) +{ + delete Operation(new ParanoidNumber(a), MULTIPLY); + return *this; +} + + +ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a) +{ + delete Operation(new ParanoidNumber(a), DIVIDE); + return *this; +} + +// 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 + { + 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; + + + + if ((NoFactors() && b->NoFactors()) + || (GetFactors() == b->GetFactors())) + { + if (ParanoidOp(m_value, b->m_value, op)) + { + Optype addop = (op == ADD) ? ADD : SUBTRACT; + for (auto add : b->m_next[ADD]) + { + 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(); + 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; +} + +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) + { + 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; + + + + if (NoTerms() && b->NoTerms()) + { + if (ParanoidOp(m_value, b->m_value, op)) + { + Optype mulop = (op == MULTIPLY) ? MULTIPLY : DIVIDE; + for (auto mul : b->m_next[MULTIPLY]) + { + delete OperationFactor(mul, mulop); + } + 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; + } + } + + + 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; +} + + + +/** + * 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) + { + Optype f = Optype(i); + s << this; + for (auto n : m_next[f]) + { + s << OpChar(f) << n->PStr(); + } + } + return s.str(); +} + +bool ParanoidNumber::Simplify(Optype op) +{ + vector next(0); + swap(m_next[op], next); + for (auto n : next) + { + delete Operation(n, op); + } + 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]) + { + m_cached_result *= mul->Digit(); + } + 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; +} + + } diff --git a/src/paranoidnumber.h b/src/paranoidnumber.h index 9005a15..9e63a23 100644 --- a/src/paranoidnumber.h +++ b/src/paranoidnumber.h @@ -7,13 +7,31 @@ #include #include "log.h" #include +#include +#include #define PARANOID_DIGIT_T float // we could theoretically replace this with a template // but let's not do that... namespace IPDF { - typedef enum {ADD, SUBTRACT, MULTIPLY, DIVIDE} Optype; + typedef enum {ADD, SUBTRACT, MULTIPLY, DIVIDE, NOP} Optype; + inline Optype InverseOp(Optype op) + { + return ((op == ADD) ? SUBTRACT : + (op == SUBTRACT) ? ADD : + (op == MULTIPLY) ? DIVIDE : + (op == DIVIDE) ? MULTIPLY : + (op == NOP) ? NOP : NOP); + } + + + inline char OpChar(int op) + { + static char opch[] = {'+','-','*','/'}; + return (op < NOP && op >= 0) ? opch[op] : '?'; + } + /** Performs an operation, returning if the result was exact **/ // NOTE: DIFFERENT to ParanoidOp (although that wraps to this...) @@ -30,137 +48,91 @@ namespace IPDF } return false; } - - template <> bool TrustingOp(float & a, const float & b, Optype op); template <> bool TrustingOp(double & a, const double & b, Optype op); template <> bool TrustingOp(int8_t & a, const int8_t & b, Optype op); - // Attempt to comine two terms: a*b + c*d or a/b + c/d - template bool CombineTerms(T & aa, Optype aop, T & bb, T & cc, Optype cop, T & dd) - { - T a(aa); T b(bb); T c(cc); T d(dd); - if (aop == MULTIPLY && cop == MULTIPLY) // a*b + c*d - { - - if ((ParanoidOp(c, b, DIVIDE) || ParanoidOp(d, b, DIVIDE)) - && TrustingOp(c, d, MULTIPLY) && TrustingOp(a,c,ADD) - && TrustingOp(a, b, MULTIPLY)) // (a + (cd)/b) * b - { - aa = a; - bb = 1; - cc = 1; - dd = 1; - return true; - } - if ((ParanoidOp(a, d, DIVIDE) || ParanoidOp(b, d, DIVIDE)) - && TrustingOp(a, b, MULTIPLY) && TrustingOp(a,c,ADD) - && TrustingOp(a, d, MULTIPLY)) // ((ab)/d + c)*d - { - aa = a; - bb = 1; - cc = 1; - dd = 1; - return true; - } - return false; - } - else if (aop == DIVIDE && cop == DIVIDE) - { - - - if (TrustingOp(a, d, MULTIPLY) && TrustingOp(c, b, MULTIPLY) - && TrustingOp(a, c, ADD) && TrustingOp(b, d, MULTIPLY)) - { - cc = 1; - dd = 1; - if (ParanoidOp(a, b, DIVIDE)) - { - aa = a; - bb = 1; - return true; - } - aa = a; - bb = b; - return true; - } - return false; - } - return false; - } - + /** + * A ParanoidNumber + * Idea: Perform regular floating point arithmetic but rearrange operations to only ever use exact results + * Memory Usage: O(all of it) + * CPU Usage: O(all of it) + * Accuracy: O(gives better result for 0.3+0.3+0.3, gives same result for everything else, or worse result) + * + * The ParanoidNumber basically stores 4 linked lists which can be split into two "dimensions" + * 1. Terms to ADD and terms to SUBTRACT + * 2. Factors to MULTIPLY and DIVIDE + * Because ADD and SUBTRACT are inverse operations and MULTIPLY and DIVIDE are inverse operations + * See paranoidnumber.cpp and the ParanoidNumber::Operation function + */ class ParanoidNumber { public: typedef PARANOID_DIGIT_T digit_t; - ParanoidNumber(digit_t value=0, Optype type = ADD) : m_value(value), m_op(type), m_next_term(NULL), m_next_factor(NULL) + ParanoidNumber(digit_t value=0) : m_value(value), m_cached_result(value) { Construct(); } - ParanoidNumber(const ParanoidNumber & cpy) : m_value(cpy.m_value), m_op(cpy.m_op), m_next_term(NULL), m_next_factor(NULL) + ParanoidNumber(const ParanoidNumber & cpy) : m_value(cpy.m_value), m_cached_result(cpy.m_cached_result) { - if (cpy.m_next_term != NULL) - { - m_next_term = new ParanoidNumber(*(cpy.m_next_term)); - } - if (cpy.m_next_factor != NULL) + Construct(); + for (int i = 0; i < NOP; ++i) { - m_next_factor = new ParanoidNumber(*(cpy.m_next_factor)); + for (auto next : cpy.m_next[i]) + m_next[i].push_back(new ParanoidNumber(*next)); } - Construct(); - } - - ParanoidNumber(const ParanoidNumber & cpy, Optype type) : ParanoidNumber(cpy) - { - m_op = type; } ParanoidNumber(const char * str); - ParanoidNumber(const std::string & str) : ParanoidNumber(str.c_str()) {Construct();} + ParanoidNumber(const std::string & str) : ParanoidNumber(str.c_str()) {} + + virtual ~ParanoidNumber(); - virtual ~ParanoidNumber() + inline void Construct() { - if (m_next_term != NULL) - delete m_next_term; - if (m_next_factor != NULL) - delete m_next_factor; - g_count--; + g_count++; } - inline void Construct() {g_count++;} - template T Convert() const; - template T AddTerms() const; - template T MultiplyFactors() const; - template T Head() const {return (m_op == SUBTRACT) ? T(-m_value) : T(m_value);} - + digit_t GetFactors(); + digit_t GetTerms(); + + double ToDouble() {return (double)Digit();} + digit_t Digit(); - - double ToDouble() const {return Convert();} - float ToFloat() const {return Convert();} - digit_t Digit() const {return Convert();} - - bool Floating() const {return (m_next_term == NULL && m_next_factor == NULL);} + bool Floating() const + { + return NoFactors() && NoTerms(); + } bool Sunken() const {return !Floating();} // I could not resist... + bool NoFactors() const {return (m_next[MULTIPLY].size() == 0 && m_next[DIVIDE].size() == 0);} + bool NoTerms() const {return (m_next[ADD].size() == 0 && m_next[SUBTRACT].size() == 0);} + ParanoidNumber & operator+=(const ParanoidNumber & a); ParanoidNumber & operator-=(const ParanoidNumber & a); ParanoidNumber & operator*=(const ParanoidNumber & a); ParanoidNumber & operator/=(const ParanoidNumber & a); ParanoidNumber & operator=(const ParanoidNumber & a); + ParanoidNumber * OperationTerm(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point = NULL, Optype * mop = NULL); + ParanoidNumber * OperationFactor(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point = NULL, Optype * mop = NULL); + ParanoidNumber * TrivialOp(ParanoidNumber * b, Optype op); + ParanoidNumber * Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point = NULL, Optype * mop = NULL); + bool Simplify(Optype op); + bool FullSimplify(); - bool operator<(const ParanoidNumber & a) const {return ToDouble() < a.ToDouble();} - bool operator<=(const ParanoidNumber & a) const {return this->operator<(a) || this->operator==(a);} - bool operator>(const ParanoidNumber & a) const {return !(this->operator<=(a));} - bool operator>=(const ParanoidNumber & a) const {return !(this->operator<(a));} - bool operator==(const ParanoidNumber & a) const {return ToDouble() == a.ToDouble();} - bool operator!=(const ParanoidNumber & a) const {return !(this->operator==(a));} + bool operator<(ParanoidNumber & a) {return ToDouble() < a.ToDouble();} + bool operator<=(ParanoidNumber & a) {return this->operator<(a) || this->operator==(a);} + bool operator>(ParanoidNumber & a) {return !(this->operator<=(a));} + bool operator>=(ParanoidNumber & a) {return !(this->operator<(a));} + bool operator==(ParanoidNumber & a) {return ToDouble() == a.ToDouble();} + bool operator!=(ParanoidNumber & a) {return !(this->operator==(a));} ParanoidNumber operator+(const ParanoidNumber & a) const { @@ -188,13 +160,29 @@ namespace IPDF } std::string Str() const; - static char OpChar(Optype op) + + ParanoidNumber * CopyTerms() { - static char opch[] = {'+','-','*','/'}; - return opch[(int)op]; + ParanoidNumber * copy = new ParanoidNumber(*this); + copy->m_value = 0; + copy->Simplify(ADD); + copy->Simplify(SUBTRACT); + return copy; } - + + ParanoidNumber * CopyFactors() + { + ParanoidNumber * copy = new ParanoidNumber(*this); + copy->m_value = 1; + copy->Simplify(MULTIPLY); + copy->Simplify(DIVIDE); + return copy; + } + + static int64_t Paranoia() {return g_count;} + + std::string PStr() const; private: static int64_t g_count; @@ -205,45 +193,37 @@ namespace IPDF digit_t m_value; Optype m_op; - ParanoidNumber * m_next_term; - ParanoidNumber * m_next_factor; + std::vector m_next[4]; + digit_t m_cached_result; + bool m_cache_valid; }; + + + template -T ParanoidNumber::AddTerms() const +T ParanoidNumber::Convert() const { - T value(0); - for (ParanoidNumber * a = m_next_term; a != NULL; a = a->m_next_term) + if (!isnan(m_cached_result)) + return (T)m_cached_result; + T value(m_value); + for (auto mul : m_next[MULTIPLY]) { - value += a->Head() * a->MultiplyFactors(); + value *= mul->Convert(); } - return value; -} - -template -T ParanoidNumber::MultiplyFactors() const -{ - T value(1); - for (ParanoidNumber * a = m_next_factor; a != NULL; a = a->m_next_factor) + for (auto div : m_next[DIVIDE]) { - if (a->m_op == DIVIDE) - value /= (a->Head() + a->AddTerms()); - else - value *= (a->Head() + a->AddTerms()); + value /= div->Convert(); } + for (auto add : m_next[ADD]) + value += add->Convert(); + for (auto sub : m_next[SUBTRACT]) + value -= sub->Convert(); return value; } -template -T ParanoidNumber::Convert() const -{ - return Head() * MultiplyFactors() + AddTerms(); -} - - - } #endif //_PARANOIDNUMBER_H diff --git a/src/tests/paranoidcalculator.cpp b/src/tests/paranoidcalculator.cpp index eb30860..c8f7524 100644 --- a/src/tests/paranoidcalculator.cpp +++ b/src/tests/paranoidcalculator.cpp @@ -29,6 +29,7 @@ int main(int argc, char ** argv) float fa = da; while (cin.good()) { + Debug("a is {%s} \"%.40lf\"", a.Str().c_str(), a.ToDouble()); char op; cin >> op; token = ""; @@ -41,11 +42,14 @@ int main(int argc, char ** argv) token += c; c = cin.get(); } - Debug("String is %s", token.c_str()); + + //Debug("String is \"%s\"", token.c_str()); float fb = strtof(token.c_str(), NULL); double db = strtod(token.c_str(), NULL); ParanoidNumber b(token.c_str()); - Debug("b is {%s} %lf", b.Str().c_str(), b.ToDouble()); + + Debug("b is {%s} \"%lf\"", b.Str().c_str(), b.ToDouble()); + Debug("db is %lf", db); switch (op) { case '+': @@ -70,12 +74,12 @@ int main(int argc, char ** argv) break; } - Debug("a is: %s", a.Str().c_str()); - Debug("a as double: %.40f\n", a.ToDouble()); - Debug("a as float: %.40f\n", a.ToFloat()); - Debug("a as int64_t: %ld\n", a.Convert()); - Debug("floats give: %.40f\n", fa); - Debug("double gives: %.40f\n", da); + Debug("a is: {%s}", a.Str().c_str()); + Debug("a as double: %.40lf", a.ToDouble()); + //Debug("a as float: %.40f", a.ToFloat()); + //Debug("a as int64_t: %ld", a.Convert()); + //Debug("floats give: %.40f", fa); + Debug("double gives: %.40lf", da); } diff --git a/src/tests/paranoidtester.cpp b/src/tests/paranoidtester.cpp index 84c4bb0..47969b9 100644 --- a/src/tests/paranoidtester.cpp +++ b/src/tests/paranoidtester.cpp @@ -13,7 +13,7 @@ using namespace std; using namespace IPDF; -string RandomNumberAsString(int max_digits = 12) +string RandomNumberAsString(int max_digits = 3) { string result(""); int digits = 1+(rand() % max_digits); @@ -30,31 +30,219 @@ string RandomNumberAsString(int max_digits = 12) return result; } -bool CloseEnough(double d, ParanoidNumber & p) +bool CloseEnough(long double d, ParanoidNumber & p, long double eps = 1e-6) { - double pd = p.ToDouble(); + long double pd = p.Convert(); if (d == 0) - return fabs(pd) <= 1e-6; - return fabs((fabs(pd - d) / d)) <= 1e-6; + return fabs(pd) <= eps; + return fabs((fabs(pd - d) / d)) <= eps; +} + +void TestOp(ParanoidNumber & p, double & d, Optype op, const double amount) +{ + string p0str(p.Str()); + double p0 = p.ToDouble(); + switch (op) + { + case ADD: + p += amount; + d += amount; + break; + case SUBTRACT: + p -= amount; + d -= amount; + break; + case MULTIPLY: + p *= amount; + d *= amount; + break; + case DIVIDE: + p /= amount; + d /= amount; + break; + default: + break; + } + if (!CloseEnough(d, p)) + { + Debug("%lf %c= %lf failed", p0, OpChar(op), amount); + Debug("%lf vs %lf", p.ToDouble(), d); + Debug("Before: {%s}\n", p0str.c_str()); + Debug("After: {%s}\n", p.Str().c_str()); + Fatal(":-("); + } + +} + +void TestAddSubIntegers(int max=100) +{ + Debug("Test add/sub integers 0 -> %i", max); + ParanoidNumber p; + double d(0); + for (int a = 0; a < max; ++a) + { + TestOp(p, d, ADD, a); + for (int b = 0; b < max; ++b) + { + TestOp(p, d, SUBTRACT, b); + } + for (int b = 0; b < max; ++b) + { + TestOp(p, d, ADD, b); + } + } + for (int a = 0; a < max; ++a) + { + TestOp(p, d, SUBTRACT, a); + for (int b = 0; b < max; ++b) + { + TestOp(p, d, ADD, b); + } + for (int b = 0; b < max; ++b) + { + TestOp(p, d, SUBTRACT, b); + } + } + Debug("PN Yields: %.40lf", p.ToDouble()); + Debug("Doubles Yield: %.40lf", d); + Debug("Complete!"); + +} + +void TestMulDivIntegers(int max=50) +{ + Debug("Test mul/div integers 1 -> %i", max); + ParanoidNumber p(1.0); + double d(1.0); + for (int a = 1; a < max; ++a) + { + TestOp(p, d, MULTIPLY, a); + for (int b = 1; b < max; ++b) + { + TestOp(p, d, DIVIDE, b); + } + for (int b = 1; b < max; ++b) + { + TestOp(p, d, MULTIPLY, b); + } + } + for (int a = 1; a < max; ++a) + { + TestOp(p, d, DIVIDE, a); + for (int b = 1; b < max; ++b) + { + TestOp(p, d, MULTIPLY, b); + } + for (int b = 1; b < max; ++b) + { + TestOp(p, d, DIVIDE, b); + } + } + Debug("PN Yields: %.40lf", p.ToDouble()); + Debug("Doubles Yield: %.40lf", d); + Debug("Complete!"); + +} + +void TestRandomisedOps(int test_cases = 1000, int ops_per_case = 1, int max_digits = 4) +{ + Debug("Test %i*%i randomised ops (max digits = %i)", test_cases, ops_per_case, max_digits); + long double eps = 1e-2; //* (1e4*ops_per_case); + for (int i = 0; i < test_cases; ++i) + { + string s = RandomNumberAsString(max_digits); + ParanoidNumber a(s); + + double da(a.ToDouble()); + for (int j = 1; j <= ops_per_case; ++j) + { + double da2(a.ToDouble()); + s = RandomNumberAsString(max_digits); + ParanoidNumber b(s); + double db(b.ToDouble()); + + + + Optype op = Optype(rand() % 4); + + ParanoidNumber a_before(a); + + + switch (op) + { + case ADD: + a += b; + da += db; + da2 += db; + break; + case SUBTRACT: + a -= b; + da -= db; + da2 -= db; + break; + case MULTIPLY: + a *= b; + da *= db; + da2 *= db; + break; + case DIVIDE: + if (db == 0) + { + --i; + } + else + { + a /= b; + da /= db; + da2 /= db; + } + break; + case NOP: + break; + } + if (!CloseEnough(da2, a, eps)) + { + Error("{%s} %c= {%s}", a_before.Str().c_str(), OpChar(op), b.Str().c_str()); + Error("{%s}", a.Str().c_str()); + Error("double Yields: %.40lf", da); + Error("PN Yields: %.40lf", a.ToDouble()); + Fatal("Failed on case %i", i*ops_per_case + j-1); + } + } + if (!CloseEnough(da, a, eps)) + { + Warn("double Yields: %.40lf", da); + Warn("PN Yields: %.40lf", a.ToDouble()); + } + } + Debug("Complete!"); + } #define TEST_CASES 1000 int main(int argc, char ** argv) { - srand(time(NULL)); + TestAddSubIntegers(); + TestMulDivIntegers(); + for (int i = 1; i <= 100; ++i) + TestRandomisedOps(1000, i); + return 0; + srand(0);//time(NULL)); //always test off same set string number(RandomNumberAsString()); ParanoidNumber a(number); + float fa = strtof(number.c_str(), NULL); double da = strtod(number.c_str(), NULL); double diff = 0; long double lda = strtold(number.c_str(), NULL); - + Debug("a is %s", a.Str().c_str()); if (fabs(a.ToDouble() - da) > 1e-6) { Error("double %lf, pn %lf {%s}", da, a.ToDouble(), a.Str().c_str()); Fatal("Didn't construct correctly off %s", number.c_str()); + } char opch[] = {'+','-','*','/'}; @@ -109,15 +297,15 @@ int main(int argc, char ** argv) break; } diff = 100.0*(fabs(a.ToDouble() - da) / da); - if (!CloseEnough(da, a)) + if (!CloseEnough(lda, a)) { Error("Op %i: ParanoidNumber probably doesn't work", i); Error("Operation: %lf %c %lf", oldda, opch[op], db); Error("As PN: %lf %c %lf", olda.ToDouble(), opch[op], b.ToDouble()); + Error("PN String before: %s", olda.Str().c_str()); Error("PN String: %s", a.Str().c_str()); - Error("Diff is %.40lf", diff); Error("LONG double gives %.40llf", lda); - Fatal("%.40lf, expected aboout %.40lf", a.ToDouble(), da); + Fatal("%.40llf, expected aboout %.40llf", a.Convert(), lda); }