Paranoia is setting in
authorSam Moore <matches@ucc.asn.au>
Thu, 11 Sep 2014 16:37:35 +0000 (00:37 +0800)
committerSam Moore <matches@ucc.asn.au>
Thu, 11 Sep 2014 16:37:35 +0000 (00:37 +0800)
Paranoid Numbers are hard to simplify properly. But they work. Kind of...

12 files changed:
src/Makefile
src/main.cpp
src/main.h
src/objectrenderer.cpp
src/paranoidnumber.cpp
src/paranoidnumber.h
src/progressbar.h [new file with mode: 0644]
src/real.cpp
src/real.h
src/tests/calculator.cpp
src/tests/paranoidcalculator.cpp
src/tests/paranoidtester.cpp

index 1e4d83f..2b5fc73 100644 (file)
@@ -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/<target>` 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/<target>` 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)
index dc20d75..d65dc60 100644 (file)
 
 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));
index 3de819c..d20ea52 100644 (file)
@@ -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())
                {
index 8f4c55f..fe97afd 100644 (file)
@@ -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);
 }
        
index 2bb44da..b11543e 100644 (file)
@@ -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<digit_t>(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<digit_t>(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<digit_t>(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<digit_t>(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<digit_t>(a->m_value, b->Head<digit_t>(), 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<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;
                        }
-                       //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>(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>(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;
+}
+
 }
index a2e8542..f767ffa 100644 (file)
@@ -6,19 +6,95 @@
 #include <map>
 #include <string>
 #include "log.h"
+#include <fenv.h>
 
-
+#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 <class T> bool TrustingOp(T & a, const T & b, Optype op);
+
+       /** Performs an operation _only_ if the result would be exact **/
+       template <class T> bool ParanoidOp(T & a, const T & b, Optype op)
+       {
+               T cpy(a);
+               if (TrustingOp<T>(cpy, b, op))
+               {
+                       a = cpy;
+                       return true;
+               }
+               return false;
+       }
+
+
+       template <> bool TrustingOp<float>(float & a, const float & b, Optype op);
+       template <> bool TrustingOp<double>(double & a, const double & b, Optype op);
+       template <> bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op);
+       
+       // Attempt to comine two terms: a*b + c*d or a/b + c/d
+       template <class T> 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<T>(c, b, DIVIDE) || ParanoidOp(d, b, DIVIDE))
+                               && TrustingOp<T>(c, d, MULTIPLY) && TrustingOp<T>(a,c,ADD)
+                               && TrustingOp<T>(a, b, MULTIPLY)) // (a + (cd)/b) * b
+                       {
+                               aa = a;
+                               bb = 1;
+                               cc = 1;
+                               dd = 1;
+                               return true;
+                       }
+                       if ((ParanoidOp<T>(a, d, DIVIDE) || ParanoidOp(b, d, DIVIDE))
+                               && TrustingOp<T>(a, b, MULTIPLY) && TrustingOp<T>(a,c,ADD)
+                               && TrustingOp<T>(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<T>(a, d, MULTIPLY) && TrustingOp<T>(c, b, MULTIPLY)
+                               && TrustingOp<T>(a, c, ADD) && TrustingOp<T>(b, d, MULTIPLY))
+                       {
+                               cc = 1;
+                               dd = 1;
+                               if (ParanoidOp<T>(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 <class T> T Convert() const;
+                       template <class T> T AddTerms() const;
+                       template <class T> T MultiplyFactors() const;
+                       template <class T> T Head() const {return (m_op == SUBTRACT) ? T(-m_value) : T(m_value);}
+                       
+
+                       
+                       
                        double ToDouble() const {return Convert<double>();}
                        float ToFloat() const {return Convert<float>();}
+                       digit_t Digit() const {return Convert<digit_t>();}
                        
                        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 <class T>
-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<T>();
-                               break;
-                       case DIVIDE:
-                               value /= n->Convert<T>();
-                               break;
-                       default:
-                               Fatal("Shouldn't happen");
-                               break;
-               }
+               value += a->Head<T>() * a->MultiplyFactors<T>();
        }
-       n = m_next_term;
-       if (n != NULL)
+       return value;
+}
+
+template <class T>
+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<T>();
-                               break;
-                       default:
-                               Fatal("Shouldn't happen");
-               }
+               if (a->m_op == DIVIDE)
+                       value /= (a->Head<T>() + a->AddTerms<T>());     
+               else
+                       value *= (a->Head<T>() + a->AddTerms<T>());     
        }
        return value;
 }
 
 
+
+template <class T>
+T ParanoidNumber::Convert() const
+{
+       return Head<T>() * MultiplyFactors<T>() + AddTerms<T>();
+}
+
+
+
 }
 
 #endif //_PARANOIDNUMBER_H
diff --git a/src/progressbar.h b/src/progressbar.h
new file mode 100644 (file)
index 0000000..f00027a
--- /dev/null
@@ -0,0 +1,49 @@
+#include <sys/time.h>
+#include <unistd.h>
+#include <cstdio>
+
+/**
+ * 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;
+}
index aa0cfc8..bed613f 100644 (file)
@@ -1,8 +1,9 @@
 #include "real.h"
+#include <fenv.h>
 
 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<int64_t>", 
                "Rational<Arbint>",
-               "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
 
+
+
+
 }
index b107649..9ea302d 100644 (file)
 #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 <mpreal.h>
-#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<int64_t> 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<ARBINT> 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<ARBINT> Sqrt(const Rational<ARBINT> & 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<Real, Real> & p) : x(p.first), y(p.second) {}
+               #if REALTYPE != REAL_IRRAM
                Vec2(const std::pair<int64_t, int64_t> & p) : x(p.first), y(p.second) {}
-       
+               #else
+               Vec2(const std::pair<int64_t, int64_t> & 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
index ab100b9..13903d6 100644 (file)
@@ -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;
index 9c5a88a..156247f 100644 (file)
@@ -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 == '/')
index bf1dc0b..98c056d 100644 (file)
@@ -8,14 +8,16 @@
 #include <cfloat>
 #include <fenv.h>
 #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);
 }

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