From 20788af97c06b76040ea2de5ab3ddc683a261365 Mon Sep 17 00:00:00 2001 From: Sam Moore Date: Fri, 12 Sep 2014 00:37:35 +0800 Subject: [PATCH] Paranoia is setting in Paranoid Numbers are hard to simplify properly. But they work. Kind of... --- src/Makefile | 64 +++++++-- src/main.cpp | 6 +- src/main.h | 9 +- src/objectrenderer.cpp | 4 +- src/paranoidnumber.cpp | 220 ++++++++++++++++++++++--------- src/paranoidnumber.h | 176 ++++++++++++++++++------- src/progressbar.h | 49 +++++++ src/real.cpp | 11 +- src/real.h | 54 +++++--- src/tests/calculator.cpp | 3 + src/tests/paranoidcalculator.cpp | 16 ++- src/tests/paranoidtester.cpp | 47 +++++-- 12 files changed, 489 insertions(+), 170 deletions(-) create mode 100644 src/progressbar.h diff --git a/src/Makefile b/src/Makefile index 1e4d83f..2b5fc73 100644 --- a/src/Makefile +++ b/src/Makefile @@ -3,14 +3,14 @@ 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 MAIN = main.o -OBJ = log.o real.o bezier.o document.o objectrenderer.o view.o screen.o vfpu.o quadtree.o graphicsbuffer.o framebuffer.o shaderprogram.o stb_truetype.o gl_core44.o add_digits_asm.o sub_digits_asm.o mul_digits_asm.o div_digits_asm.o arbint.o moc_controlpanel.o controlpanel.o path.o paranoidnumber.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 QT_INCLUDE := -I/usr/share/qt4/mkspecs/linux-g++-64 -I. -I/usr/include/qt4/QtCore -I/usr/include/qt4/QtGui -I/usr/include/qt4 -I. -Itests -I. QT_DEF := -DQT_NO_DEBUG -DQT_GUI_LIB -DQT_CORE_LIB QT_LIB := -L/usr/lib/x86_64-linux-gnu -lQtGui -lQtCore -lpthread -LIB_x86_64 = ../contrib/lib/libSDL2-2.0.so.0 -lGL -lmpfr -lgmp $(QT_LIB) -LIB_i386 = ../contrib/lib32/libSDL2-2.0.so.0 -lGL -lgmp +LIB_x86_64 = ../contrib/lib/libSDL2-2.0.so.0 -lGL +LIB_i386 = ../contrib/lib32/libSDL2-2.0.so.0 -lGL LIB_i686 = $(LIB_i386) MAINRPATH_x86_64 = -Wl,-rpath,'$$ORIGIN/../contrib/lib' @@ -21,7 +21,7 @@ TESTRPATH_i386 = -Wl,-rpath,'$$ORIGIN/../../contrib/lib32' TESTRPATH_i686 = $(TESTRPATH_i386) OBJPATHS = $(OBJ:%=../obj/%) DEPS := $(OBJPATHS:%.o=%.d) -CFLAGS_x86_64 := -I../contrib/include/SDL2 -I`pwd` $(QT_INCLUDE) +CFLAGS_x86_64 := -I../contrib/include/SDL2 -I`pwd` CFLAGS_i386 := -I../contrib/include32/SDL2 -I`pwd` CFLAGS_i686 := $(CFLAGS_i386) @@ -31,23 +31,61 @@ MAINRPATH := $(MAINRPATH_$(ARCH)) TESTRPATH := $(TESTRPATH_$(ARCH)) CFLAGS := $(CFLAGS_$(ARCH)) -REALTYPE=1 -LINKOBJ = $(OBJPATHS) - RM = rm -f BIN = ../bin/ipdf +REALTYPE=1 +CONTROLPANEL=enabled +QUADTREE=disabled +DEF = -DREALTYPE=$(REALTYPE) + +## Only link with things we care about + +ifeq ($(QUADTREE),enabled) + OBJ := $(OBJ) quadtree.o +else + DEF := $(DEF) -DQUADTREE_DISABLED +endif + +ifeq ($(CONTROLPANEL),enabled) + LIB := $(LIB) $(QT_LIB) + DEF := $(DEF) $(QT_DEF) + CFLAGS := $(CFLAGS) $(QT_INCLUDE) + OBJ := $(OBJ) moc_controlpanel.o controlpanel.o +else + DEF := $(DEF) -DCONTROLPANEL_DISABLED +endif + +ifeq ($REALTYPE),3) + OBJ := $(OBJ) vfpu.o +endif + +ifeq ($(REALTYPE),5) + OBJ := $(OBJ) add_digits_asm.o sub_digits_asm.o mul_digits_asm.o div_digits_asm.o arbint.o + LIB := $(LIB) -lgmp +endif + +ifeq ($(REALTYPE),6) + LIB := $(LIB) -lgmp -lmpfr +endif + +ifeq ($(REALTYPE),7) + LIB := $(LIB) -L../contrib/iRRAM/lib -liRRAM -lgmp -lmpfr + CFLAGS := $(CFLAGS) -I../contrib/iRRAM/include +endif + +LINKOBJ = $(OBJPATHS) -all : REAL = 1 +all : REALTYPE = 1 all : $(BIN) -single : REAL = 0 +single : REALTYPE = 0 single : $(BIN) -double : REAL = 1 +double : REALTYPE = 1 double : $(BIN) -DEF = -DREAL=$(REALTYPE) -DQUADTREE_DISABLED + CFLAGS := $(CFLAGS) $(QT_DEF) @@ -62,8 +100,8 @@ movie : $(BIN) ../tools/stream_plot.py $(RM) ../data/movie.ogv ./ipdf | tee ../data/performance.dat | ../tools/stream_plot.py 2>/dev/null & recordmydesktop --fps 10 --on-the-fly-encoding -o ../data/movie.ogv -# The tests will compile with the default REAL definition -# To change that you can run as `make DEFS="REAL=X" tests/` where X is your chosen type +# The tests will compile with the default REALTYPE definition +# To change that you can run as `make DEFS="REALTYPE=X" tests/` where X is your chosen type # But remember to make clean first. tests/% : tests/%.cpp ../obj/tests/%.o $(LINKOBJ) $(CXX) $(CFLAGS) -o $@.test $(LINKOBJ) ../obj/$@.o $(LIB) $(TESTRPATH) diff --git a/src/main.cpp b/src/main.cpp index dc20d75..d65dc60 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -11,6 +11,10 @@ int main(int argc, char ** argv) { + #if REALTYPE == REAL_IRRAM + iRRAM_initialize(argc,argv); + #endif + #ifndef __STDC_IEC_559__ Warn("__STDC_IEC_559__ not defined. IEEE 754 floating point not fully supported.\n"); #endif @@ -18,7 +22,7 @@ int main(int argc, char ** argv) // We want to crash if we ever get a NaN. feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); - Debug("Compiled with REAL = %d => \"%s\" sizeof(Real) == %d bytes", REAL, g_real_name[REAL], sizeof(Real)); + Debug("Compiled with REAL = %d => \"%s\" sizeof(Real) == %d bytes", REALTYPE, g_real_name[REALTYPE], sizeof(Real)); Document doc("","fonts/ComicSans.ttf"); srand(time(NULL)); diff --git a/src/main.h b/src/main.h index 3de819c..d20ea52 100644 --- a/src/main.h +++ b/src/main.h @@ -53,7 +53,7 @@ void RatCatcher(int x, int y, int buttons, int wheel, Screen * scr, View * view) } if (buttons) { - #if REAL == REAL_RATIONAL + #if REALTYPE == REAL_RATIONAL view->Translate(Real(oldx, scr->ViewportWidth()) -Real(x,scr->ViewportWidth()), Real(oldy, scr->ViewportHeight()) - Real(y,scr->ViewportHeight())); #else view->Translate(Real(oldx-x)/Real(scr->ViewportWidth()), Real(oldy-y)/Real(scr->ViewportHeight())); @@ -69,7 +69,7 @@ void RatCatcher(int x, int y, int buttons, int wheel, Screen * scr, View * view) if (wheel) { - #if REAL == REAL_RATIONAL + #if REALTYPE == REAL_RATIONAL view->ScaleAroundPoint(Real(x,scr->ViewportWidth()), Real(y,scr->ViewportHeight()), Real(20-wheel, 20)); #else view->ScaleAroundPoint(Real(x)/Real(scr->ViewportWidth()),Real(y)/Real(scr->ViewportHeight()), Real(expf(-wheel/20.f))); @@ -126,7 +126,10 @@ inline void MainLoop(Document & doc, Screen & scr, View & view) scr.DebugFontPrintF("[GPU] Render took %lf ms (%lf FPS) (total %lf s, avg FPS %lf)\n", gpu_frame*1e3, 1.0/gpu_frame, total_gpu_time, frames/total_gpu_time); scr.DebugFontPrintF("[REALTIME] Render+Present+Cruft took %lf ms (%lf FPS) (total %lf s, avg FPS %lf)\n", real_frame*1e3, 1.0/real_frame, total_real_time,frames/total_real_time); scr.DebugFontPrintF("View bounds: %s\n", view.GetBounds().Str().c_str()); - scr.DebugFontPrintF("type of Real == %s\n", g_real_name[REAL]); + scr.DebugFontPrintF("type of Real == %s\n", g_real_name[REALTYPE]); + #if REALTYPE == REAL_MPFRCPP + scr.DebugFontPrintf("Precision: %s\nRounding: %s\n"); + #endif if (view.UsingGPUTransform()) { diff --git a/src/objectrenderer.cpp b/src/objectrenderer.cpp index 8f4c55f..fe97afd 100644 --- a/src/objectrenderer.cpp +++ b/src/objectrenderer.cpp @@ -219,8 +219,8 @@ ObjectRenderer::PixelPoint ObjectRenderer::CPUPointLocation(const Vec2 & point, { // hack... Rect result = view.TransformToViewCoords(Rect(point.x, point.y,1,1)); - int64_t x = Int64(result.x*target.w); - int64_t y = Int64(result.y*target.h); + int64_t x = Int64(result.x)*target.w; + int64_t y = Int64(result.y)*target.h; return PixelPoint(x,y); } diff --git a/src/paranoidnumber.cpp b/src/paranoidnumber.cpp index 2bb44da..b11543e 100644 --- a/src/paranoidnumber.cpp +++ b/src/paranoidnumber.cpp @@ -8,7 +8,7 @@ using namespace std; 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) { @@ -46,7 +46,7 @@ ParanoidNumber::ParanoidNumber(const char * str) : m_value(0), m_op(ADD), m_next //Debug("Now at %s", Str().c_str()); } - Debug("Constructed {%s} from %s (%f)", Str().c_str(), str, ToDouble()); + //Debug("Constructed {%s} from %s (%f)", Str().c_str(), str, ToDouble()); } ParanoidNumber & ParanoidNumber::operator=(const ParanoidNumber & a) @@ -68,9 +68,15 @@ ParanoidNumber & ParanoidNumber::operator=(const ParanoidNumber & a) 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)) + { + Simplify(); + return *this; + } + } ParanoidNumber * nt = m_next_term; ParanoidNumber * nf = m_next_factor; @@ -109,7 +115,15 @@ 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)) + { + Simplify(); + return *this; + } + } + ParanoidNumber * nt = m_next_term; ParanoidNumber * nf = m_next_factor; @@ -151,35 +165,59 @@ ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a) ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a) { - Debug("{%s} *= {%s}", Str().c_str(), a.Str().c_str()); + + //if (m_value == 0) + // return *this; + //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; + 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 = nf; + t->m_next_factor = new ParanoidNumber(a, MULTIPLY); + if (m_next_term != NULL) m_next_term->operator*=(a); - //Debug("Simplify after mul"); - Debug("Simplify {%s}", Str().c_str()); + + //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 * nf = m_next_factor; - m_next_factor = new ParanoidNumber(a, DIVIDE); - ParanoidNumber * t = m_next_factor; + ParanoidNumber * t = this; while (t->m_next_factor != NULL) + { t = t->m_next_factor; - t->m_next_factor = nf; + } + t->m_next_factor = new ParanoidNumber(a, DIVIDE); + if (m_next_term != NULL) m_next_term->operator/=(a); + Simplify(); return *this; } @@ -188,6 +226,7 @@ ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a) void ParanoidNumber::SimplifyTerms() { + //Debug("Simplify {%s}", Str().c_str()); if (m_next_term == NULL) { @@ -197,8 +236,13 @@ void ParanoidNumber::SimplifyTerms() for (ParanoidNumber * a = this; a != NULL; a = a->m_next_term) { - ParanoidNumber * bprev = a; ParanoidNumber * b = a->m_next_term; + if (a->m_next_factor != NULL) + { + continue; + } + + ParanoidNumber * bprev = a; while (b != NULL) { //Debug("Simplify factors of %s", b->Str().c_str()); @@ -209,28 +253,16 @@ void ParanoidNumber::SimplifyTerms() b = b->m_next_term; continue; } - 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)) + bool simplify = false; + simplify = ParanoidOp(a->m_value, b->Head(), ADD); + if (simplify) { - 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; } @@ -239,6 +271,7 @@ void ParanoidNumber::SimplifyTerms() void ParanoidNumber::SimplifyFactors() { + //Debug("Simplify {%s}", Str().c_str()); if (m_next_factor == NULL) { @@ -248,6 +281,9 @@ void ParanoidNumber::SimplifyFactors() 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) @@ -259,37 +295,21 @@ void ParanoidNumber::SimplifyFactors() b = b->m_next_factor; continue; } - float f = a->m_value; - feclearexcept(FE_ALL_EXCEPT); - switch (b->m_op) + + Optype op = b->m_op; + if (a->m_op == DIVIDE) { - 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; + op = (b->m_op == DIVIDE) ? MULTIPLY : DIVIDE; } - if (!fetestexcept(FE_ALL_EXCEPT)) - { - - a->m_value = f; + + 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; } - //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; } @@ -306,29 +326,101 @@ string ParanoidNumber::Str() const { string result(""); stringstream s; - - s << m_value; + s << (double)m_value; if (m_next_factor != NULL) { - result += OpChar(m_op); - result += "("; result += s.str(); - result += m_next_factor->Str(); - result += ")"; + 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 += OpChar(m_op); 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 & a, const float & 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(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 & 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; +} + } diff --git a/src/paranoidnumber.h b/src/paranoidnumber.h index a2e8542..f767ffa 100644 --- a/src/paranoidnumber.h +++ b/src/paranoidnumber.h @@ -6,19 +6,95 @@ #include #include #include "log.h" +#include - +#define PARANOID_DIGIT_T int8_t // we could theoretically replace this with a template + // but let's not do that... namespace IPDF { + typedef enum {ADD, SUBTRACT, MULTIPLY, DIVIDE} Optype; + + /** Performs an operation, returning if the result was exact **/ + // NOTE: DIFFERENT to ParanoidOp (although that wraps to this...) + template bool TrustingOp(T & a, const T & b, Optype op); + + /** Performs an operation _only_ if the result would be exact **/ + template bool ParanoidOp(T & a, const T & b, Optype op) + { + T cpy(a); + if (TrustingOp(cpy, b, op)) + { + a = cpy; + return true; + } + 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; + } + class ParanoidNumber { + public: - typedef enum {ADD, SUBTRACT, MULTIPLY, DIVIDE} Optype; - - ParanoidNumber(float value=0, Optype type = ADD) : m_value(value), m_op(type), m_next_term(NULL), m_next_factor(NULL) + 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) { - + Construct(); } ParanoidNumber(const ParanoidNumber & cpy) : m_value(cpy.m_value), m_op(cpy.m_op), m_next_term(NULL), m_next_factor(NULL) @@ -31,22 +107,16 @@ namespace IPDF { m_next_factor = new ParanoidNumber(*(cpy.m_next_factor)); } + Construct(); } - ParanoidNumber(const ParanoidNumber & cpy, Optype type) : m_value(cpy.m_value), m_op(type), m_next_term(NULL), m_next_factor(NULL) + ParanoidNumber(const ParanoidNumber & cpy, Optype type) : ParanoidNumber(cpy) { - if (cpy.m_next_term != NULL) - { - m_next_term = new ParanoidNumber(*(cpy.m_next_term)); - } - if (cpy.m_next_factor != NULL) - { - m_next_factor = new ParanoidNumber(*(cpy.m_next_factor)); - } + m_op = type; } ParanoidNumber(const char * str); - ParanoidNumber(const std::string & str) : ParanoidNumber(str.c_str()) {} + ParanoidNumber(const std::string & str) : ParanoidNumber(str.c_str()) {Construct();} virtual ~ParanoidNumber() { @@ -54,11 +124,23 @@ namespace IPDF delete m_next_term; if (m_next_factor != NULL) delete m_next_factor; + 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);} + + + + 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 Sunken() const {return !Floating();} // I could not resist... @@ -109,62 +191,56 @@ namespace IPDF return opch[(int)op]; } + static int64_t Paranoia() {return g_count;} + private: + static int64_t g_count; void Simplify(); void SimplifyTerms(); void SimplifyFactors(); - float m_value; + digit_t m_value; Optype m_op; ParanoidNumber * m_next_term; ParanoidNumber * m_next_factor; - - - - - - }; - template -T ParanoidNumber::Convert() const +T ParanoidNumber::AddTerms() const { - T value = (m_op == SUBTRACT) ? -m_value : m_value; - const ParanoidNumber * n = m_next_factor; - if (n != NULL) + T value(0); + for (ParanoidNumber * a = m_next_term; a != NULL; a = a->m_next_term) { - switch (n->m_op) - { - case MULTIPLY: - value *= n->Convert(); - break; - case DIVIDE: - value /= n->Convert(); - break; - default: - Fatal("Shouldn't happen"); - break; - } + value += a->Head() * a->MultiplyFactors(); } - n = m_next_term; - if (n != NULL) + return value; +} + +template +T ParanoidNumber::MultiplyFactors() const +{ + T value(1); + for (ParanoidNumber * a = m_next_factor; a != NULL; a = a->m_next_factor) { - switch (n->m_op) - { - case ADD: - case SUBTRACT: - value += n->Convert(); - break; - default: - Fatal("Shouldn't happen"); - } + if (a->m_op == DIVIDE) + value /= (a->Head() + a->AddTerms()); + else + value *= (a->Head() + a->AddTerms()); } return value; } + +template +T ParanoidNumber::Convert() const +{ + return Head() * MultiplyFactors() + AddTerms(); +} + + + } #endif //_PARANOIDNUMBER_H diff --git a/src/progressbar.h b/src/progressbar.h new file mode 100644 index 0000000..f00027a --- /dev/null +++ b/src/progressbar.h @@ -0,0 +1,49 @@ +#include +#include +#include + +/** + * Print a progress bar to stderr + elapsed time since execution + * Stolen from some of my Codejam stuff, but modified to actually print the correct times... + */ +inline void ProgressBar(int64_t c, int64_t nCases, int len=20) +{ + static timeval startTime; + static timeval thisTime; + static timeval lastTime; + static int64_t total_usec = 0; + gettimeofday(&thisTime, NULL); + if (c <= 0) + { + startTime = thisTime; + lastTime = thisTime; + total_usec = 0; + } + + + fprintf(stderr, "\33[2K\r"); + //fprintf(stderr, "\n"); + + fprintf(stderr, "["); + int i = 0; + for (i = 0; i < ((len*c)/nCases); ++i) + fprintf(stderr, "="); + fprintf(stderr, ">"); + for (; i < len; ++i) + fprintf(stderr, " "); + + + int64_t delta_usec = 1000000*(thisTime.tv_sec-lastTime.tv_sec) + (thisTime.tv_usec-lastTime.tv_usec); + total_usec += delta_usec; + + int64_t ds = delta_usec / 1000000; + int64_t dus = delta_usec - 1000000*ds; + int64_t ts = total_usec / 1000000; + int64_t tus = total_usec - 1000000*ts; + + fprintf(stderr, "] Case #%li: %02li:%02li.%06li Total: %02li:%02li.%06li| ", c, + ds/60, ds%60, dus, + ts/60, ts%60, tus); + + lastTime = thisTime; +} diff --git a/src/real.cpp b/src/real.cpp index aa0cfc8..bed613f 100644 --- a/src/real.cpp +++ b/src/real.cpp @@ -1,8 +1,9 @@ #include "real.h" +#include namespace IPDF { - // Maps the REAL to a string + // Maps the REALTYPE to a string const char * g_real_name[] = { "single", "double", @@ -10,14 +11,18 @@ namespace IPDF "VFPU", "Rational", "Rational", - "mpfrc++ real" + "mpfrc++ real", + "iRRAM REAL" }; -#if REAL == REAL_RATIONAL_ARBINT +#if REALTYPE == REAL_RATIONAL_ARBINT template <> Gmpint Tabs(const Gmpint & a) { return a.Abs(); } #endif + + + } diff --git a/src/real.h b/src/real.h index b107649..9ea302d 100644 --- a/src/real.h +++ b/src/real.h @@ -12,48 +12,53 @@ #define REAL_RATIONAL 4 #define REAL_RATIONAL_ARBINT 5 #define REAL_MPFRCPP 6 +#define REAL_IRRAM 7 -#ifndef REAL - #error "REAL was not defined!" +#ifndef REALTYPE + #error "REALTYPE was not defined!" #endif -#if REAL == REAL_VFPU +#if REALTYPE == REAL_VFPU #include "vfpu.h" #endif -#if REAL == REAL_RATIONAL +#if REALTYPE == REAL_RATIONAL #include "rational.h" -#endif //REAL +#endif //REALTYPE -#if REAL == REAL_RATIONAL_ARBINT +#if REALTYPE == REAL_RATIONAL_ARBINT #include "rational.h" #include "arbint.h" #include "gmpint.h" -#endif //REAL +#endif //REALTYPE -#if REAL == REAL_MPFRCPP +#if REALTYPE == REAL_MPFRCPP #include -#endif //REAL +#endif //REALTYPE + +#if REALTYPE == REAL_IRRAM + #include "../contrib/iRRAM/include/iRRAM/lib.h" +#endif namespace IPDF { extern const char * g_real_name[]; -#if REAL == REAL_SINGLE +#if REALTYPE == REAL_SINGLE typedef float Real; -#elif REAL == REAL_DOUBLE +#elif REALTYPE == REAL_DOUBLE typedef double Real; -#elif REAL == REAL_LONG_DOUBLE +#elif REALTYPE == REAL_LONG_DOUBLE typedef long double Real; -#elif REAL == REAL_VFPU +#elif REALTYPE == REAL_VFPU typedef VFPU::VFloat Real; inline float Float(const Real & r) {return r.m_value;} inline double Double(const Real & r) {return r.m_value;} -#elif REAL == REAL_RATIONAL +#elif REALTYPE == REAL_RATIONAL typedef Rational Real; inline float Float(const Real & r) {return (float)r.ToDouble();} inline double Double(const Real & r) {return r.ToDouble();} -#elif REAL == REAL_RATIONAL_ARBINT +#elif REALTYPE == REAL_RATIONAL_ARBINT #define ARBINT Gmpint // Set to Gmpint or Arbint here typedef Rational Real; @@ -61,21 +66,26 @@ namespace IPDF inline double Double(const Real & r) {return r.ToDouble();} inline int64_t Int64(const Real & r) {return r.ToInt64();} inline Rational Sqrt(const Rational & r) {return r.Sqrt();} -#elif REAL == REAL_MPFRCPP +#elif REALTYPE == REAL_MPFRCPP typedef mpfr::mpreal Real; inline double Double(const Real & r) {return r.toDouble();} inline float Float(const Real & r) {return r.toDouble();} inline int64_t Int64(const Real & r) {return r.toLong();} inline Real Sqrt(const Real & r) {return mpfr::sqrt(r, mpfr::mpreal::get_default_rnd());} inline Real Abs(const Real & r) {return mpfr::abs(r, mpfr::mpreal::get_default_rnd());} +#elif REALTYPE == REAL_IRRAM + typedef iRRAM::REAL Real; + inline double Double(const Real & r) {return r.as_double(53);} + inline float Float(const Real & r) {return r.as_double(53);} + inline int64_t Int64(const Real & r) {return (int64_t)r.as_double(53);} + inline Real Sqrt(const Real & r) {return iRRAM::sqrt(r);} #else #error "Type of Real unspecified." -#endif //REAL +#endif //REALTYPE // Allow us to call Float on the primative types // Useful so I can template some things that could be either (a more complicated) Real or a primitive type // Mostly in the testers. - inline float Float(float f) {return (float)f;} inline float Float(double f) {return (float)f;} inline float Float(long double f) {return (float)(f);} inline double Double(float f) {return (double)f;} @@ -103,8 +113,12 @@ namespace IPDF Vec2() : x(0), y(0) {} Vec2(Real _x, Real _y) : x(_x), y(_y) {} Vec2(const std::pair & p) : x(p.first), y(p.second) {} + #if REALTYPE != REAL_IRRAM Vec2(const std::pair & p) : x(p.first), y(p.second) {} - + #else + Vec2(const std::pair & p) : x((int)p.first), y((int)p.second) {} + // Apparently iRRAM didn't implement -= for a constant argument + #endif bool operator==(const Vec2& other) const { return (x == other.x) && (y == other.y); } bool operator!=(const Vec2& other) const { return !(*this == other); } @@ -124,6 +138,8 @@ namespace IPDF }; + + } #endif //_REAL_H diff --git a/src/tests/calculator.cpp b/src/tests/calculator.cpp index ab100b9..13903d6 100644 --- a/src/tests/calculator.cpp +++ b/src/tests/calculator.cpp @@ -11,6 +11,9 @@ using namespace IPDF; int main(int argc, char ** argv) { + #if REALTYPE == REAL_IRRAM + iRRAM_initialize(argc,argv); + #endif while (cin.good()) { double da; double db; diff --git a/src/tests/paranoidcalculator.cpp b/src/tests/paranoidcalculator.cpp index 9c5a88a..156247f 100644 --- a/src/tests/paranoidcalculator.cpp +++ b/src/tests/paranoidcalculator.cpp @@ -17,17 +17,19 @@ int main(int argc, char ** argv) Debug("FLT_MAX = %.40f", FLT_MAX); Debug("FLT_MIN = %.40f", FLT_MIN); Debug("FLT_EPSILON = %.40f", FLT_EPSILON); - - ParanoidNumber a("0.3"); - Debug("start at {%s} = %lf", a.Str().c_str(), a.ToDouble()); - cout << "0.3 "; - float fa = 0.3; - double da = 0.3; + Debug("Sizeof ParanoidNumber::digit_t is %u", sizeof(ParanoidNumber::digit_t)); + Debug("Sizeof ParanoidNumber is %u", sizeof(ParanoidNumber)); + + string token(""); + cin >> token; + ParanoidNumber a(token.c_str()); + double da = a.ToDouble(); + float fa = da; while (cin.good()) { char op; cin >> op; - string token(""); + token = ""; for (char c = cin.peek(); cin.good() && !iswspace(c); c = cin.peek()) { if (c == '+' || c == '-' || c == '*' || c == '/') diff --git a/src/tests/paranoidtester.cpp b/src/tests/paranoidtester.cpp index bf1dc0b..98c056d 100644 --- a/src/tests/paranoidtester.cpp +++ b/src/tests/paranoidtester.cpp @@ -8,14 +8,16 @@ #include #include #include "paranoidnumber.h" +#include "progressbar.h" using namespace std; using namespace IPDF; -string RandomNumberAsString(int digits = 6) +string RandomNumberAsString(int max_digits = 12) { string result(""); - int dp = 1+(rand() % 3); + int digits = 1+(rand() % max_digits); + int dp = (rand() % digits)+1; for (int i = 0; i < digits; ++i) { if (i == dp) @@ -28,15 +30,25 @@ string RandomNumberAsString(int digits = 6) return result; } -#define TEST_CASES 10000 +bool CloseEnough(double d, ParanoidNumber & p) +{ + double pd = p.ToDouble(); + + if (d == 0) + return fabs(pd) <= 1e-6; + return fabs((fabs(pd - d) / d)) <= 1e-6; +} + +#define TEST_CASES 1000 int main(int argc, char ** argv) { - + srand(time(NULL)); 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); if (fabs(a.ToDouble() - da) > 1e-6) @@ -46,15 +58,27 @@ int main(int argc, char ** argv) } char opch[] = {'+','-','*','/'}; - + int opcount[] = {0,0,0,0}; for (int i = 0; i < TEST_CASES; ++i) { + ProgressBar(i, TEST_CASES); fprintf(stderr, "%.30lf (%d)+ (%d)- (%d)* (%d)/ [Paranoia %ld]",diff, opcount[0],opcount[1],opcount[2],opcount[3], ParanoidNumber::Paranoia()); number = RandomNumberAsString(); ParanoidNumber b(number); float fb = strtof(number.c_str(), NULL); double db = strtod(number.c_str(), NULL); + if (db == 0) + { + --i; + continue; + } long double ldb = strtold(number.c_str(), NULL); - int op = (rand() % 4); + int op = rand() % 4;//= 2*(rand() % 2); + while (op == 1) + op = rand() % 4; + //if (rand() % 2 == 0) + //op = 2+(rand() % 2); + + opcount[op]++; ParanoidNumber olda(a); double oldda(da); switch (op) @@ -84,16 +108,23 @@ int main(int argc, char ** argv) lda /= ldb; break; } - if (fabs(a.ToDouble() - da) > 1.0 ) + diff = 100.0*(fabs(a.ToDouble() - da) / da); + if (!CloseEnough(da, 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()); - Fatal("%lf, expected aboout %lf", a.ToDouble(), da); + Error("PN String: %s", a.Str().c_str()); + Error("Diff is %.40lf", diff); + Fatal("%.40lf, expected aboout %.40lf", a.ToDouble(), da); + + } + } printf("ParanoidNumber: {%s} = %.40lf\n", a.Str().c_str(), a.ToDouble()); printf("float: %.40f\n", fa); printf("double: %.40lf\n", da); printf("long double: %.40Lf\n", lda); + printf("diff %.40lf\n", diff); } -- 2.20.1