Added really terrible number representation as a series of operations on floats.
It currently gets a few things wrong...
# 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
+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
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
Real d2 = a3*b1 - a1*b3;
Real d3 = a1*b2 - a2*b1;
- if (fabs(d1+d2+d3) < 1e-6)
+ if (Abs(d1+d2+d3) < 1e-6)
{
type = LINE;
//Debug("LINE %s", Str().c_str());
Real delta1 = -(d1*d1);
Real delta2 = d1*d2;
Real delta3 = d1*d3 -(d2*d2);
- if (fabs(delta1+delta2+delta3) < 1e-6)
+ if (Abs(delta1+delta2+delta3) < 1e-6)
{
type = QUADRATIC;
}
Real discriminant = d1*d3*4 -d2*d2;
- if (fabs(discriminant) < 1e-6)
+ if (Abs(discriminant) < 1e-6)
{
type = CUSP;
//Debug("CUSP %s", Str().c_str());
Bezier ReParametrise(const Real& t0, const Real& t1)
{
- Debug("Reparametrise: %f -> %f",t0,t1);
+ Debug("Reparametrise: %f -> %f",Double(t0),Double(t1));
Bezier new_bezier;
// Subdivide to get from [0,t1]
new_bezier = DeCasteljauSubdivideLeft(t1);
// Convert t0 from [0,1] range to [0, t1]
Real new_t0 = t0 / t1;
- Debug("New t0 = %f", new_t0);
+ Debug("New t0 = %f", Double(new_t0));
new_bezier = new_bezier.DeCasteljauSubdivideRight(new_t0);
Debug("%s becomes %s", this->Str().c_str(), new_bezier.Str().c_str());
{
Real t1 = *it;
if (t1 == t0) continue;
- Debug(" -- t0: %f to t1: %f", t0, t1);
+ Debug(" -- t0: %f to t1: %f", Double(t0), Double(t1));
Real ptx, pty;
Evaluate(ptx, pty, ((t1 + t0) / Real(2)));
if (r.PointIn(ptx, pty))
}
else
{
- Debug("Segment removed (point at %f, %f)", ptx, pty);
+ Debug("Segment removed (point at %f, %f)", Double(ptx), Double(pty));
}
t0 = t1;
}
//Real y0(y);
int ascent = 0, descent = 0, line_gap = 0;
stbtt_GetFontVMetrics(&m_font, &ascent, &descent, &line_gap);
- float font_scale = scale / (float)(ascent - descent);
+ Real font_scale = scale;
+ font_scale /= Real(ascent - descent);
Real y_advance = Real(font_scale) * Real(ascent - descent + line_gap);
for (unsigned i = 0; i < text.size(); ++i)
{
{
kerning = stbtt_GetCodepointKernAdvance(&m_font, text[i-1], text[i]);
}
- x += Real(font_scale) * Real(kerning);
+ x += font_scale * Real(kerning);
AddFontGlyphAtPoint(&m_font, text[i], font_scale, x, y);
- x += Real(font_scale) * Real(advance_width);
+ x += font_scale * Real(advance_width);
}
}
{
// hack...
Rect result = view.TransformToViewCoords(Rect(point.x, point.y,1,1));
- int64_t x = result.x*target.w;
- int64_t y = result.y*target.h;
+ int64_t x = Int64(result.x*target.w);
+ int64_t y = Int64(result.y*target.h);
return PixelPoint(x,y);
}
ObjectRenderer::RenderLineOnCPU(pix_bounds.x+pix_bounds.w, pix_bounds.y, pix_bounds.x+pix_bounds.w, pix_bounds.y+pix_bounds.h, target, Colour(0,255,0,0));
}
- unsigned blen = min(max(2U, (unsigned)(target.w/view.GetBounds().w)),
+ unsigned blen = min(max(2U, (unsigned)Int64(Real(target.w)/view.GetBounds().w)),
min((unsigned)(pix_bounds.w+pix_bounds.h)/4 + 1, 100U));
// DeCasteljau Divide the Bezier
while (divisions.size() > 0)
{
Bezier & current = divisions.front();
- RenderLineOnCPU(current.x0, current.y0, current.x3, current.y3, target, c);
+ RenderLineOnCPU(Int64(current.x0), Int64(current.y0), Int64(current.x3), Int64(current.y3), target, c);
divisions.pop();
}
}
--- /dev/null
+#include "paranoidnumber.h"
+
+#include <sstream>
+#include <fenv.h>
+#include "log.h"
+
+using namespace std;
+namespace IPDF
+{
+
+
+ParanoidNumber::ParanoidNumber(const char * str) : m_value(0), m_op(ADD), m_next(NULL)
+{
+ int dp = 0;
+ int end = 0;
+ while (str[dp] != '\0' && str[dp] != '.')
+ {
+ ++dp;
+ ++end;
+ }
+ while (str[end] != '\0')
+ ++end;
+
+ ParanoidNumber m(1);
+ for (int i = dp-1; i >= 0; --i)
+ {
+ ParanoidNumber b(str[i]-'0');
+ b*=m;
+ this->operator+=(b);
+ m*=10;
+ }
+ ParanoidNumber n(1);
+ for (int i = dp+1; i < end; ++i)
+ {
+ n/=10;
+ ParanoidNumber b(str[i]-'0');
+ b*=n;
+ this->operator+=(b);
+ }
+
+}
+
+ParanoidNumber * ParanoidNumber::InsertAfter(float value, Optype op)
+{
+ return InsertAfter(new ParanoidNumber(value, op));
+}
+
+ParanoidNumber * ParanoidNumber::InsertAfter(ParanoidNumber * insert)
+{
+ //Debug("Insert {%s} after {%f, %s}",insert->Str().c_str(),
+ // m_value, g_opstr[m_op]);
+ ParanoidNumber * n = m_next;
+ m_next = insert;
+
+ ParanoidNumber * p = m_next;
+ while (p->m_next != NULL)
+ p = p->m_next;
+ p->m_next = n;
+
+ return m_next;
+}
+
+
+ParanoidNumber & ParanoidNumber::operator=(const ParanoidNumber & a)
+{
+ const ParanoidNumber * na = &a;
+ ParanoidNumber * nb = this;
+ ParanoidNumber * p = NULL;
+ while (na != NULL && nb != NULL)
+ {
+ nb->m_value = na->m_value;
+ nb->m_op = na->m_op;
+ na = na->m_next;
+ p = nb;
+ nb = nb->m_next;
+
+ }
+
+ while (na != NULL) // => nb == NULL
+ {
+ InsertAfter(na->m_value, na->m_op);
+ na = na->m_next;
+ }
+
+ if (nb != NULL)
+ {
+ if (p != NULL)
+ p->m_next = NULL;
+ delete nb;
+ }
+ return *this;
+}
+
+ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
+{
+ ParanoidNumber * insert = new ParanoidNumber(a, ADD);
+ if (m_next == NULL || m_next->m_op == ADD || m_next->m_op == SUBTRACT)
+ {
+ InsertAfter(insert);
+ Simplify();
+ return *this;
+ }
+
+ if (m_next->m_op == MULTIPLY) // (a*b) + c == (a+[c/b]) * b
+ insert->operator/=(*m_next);
+ else
+ insert->operator*=(*m_next); // (a/b) + c == (a+[c*b])/b
+
+ if (insert->m_next != NULL) // neither of the above simplified
+ {
+ //Debug("{%s} did not simplify, change back to {%s}", insert->Str().c_str(), a.Str().c_str());
+ insert->operator=(a); // Just add as is
+ insert->m_op = ADD;
+ ParanoidNumber * n = this;
+ while (n->m_next != NULL)
+ n = n->m_next;
+ n->InsertAfter(insert);
+ }
+ else
+ {
+ InsertAfter(insert);
+ }
+ Simplify();
+ return *this;
+}
+ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
+{
+ ParanoidNumber * insert = new ParanoidNumber(a, SUBTRACT);
+ if (m_next == NULL || m_next->m_op == ADD || m_next->m_op == SUBTRACT)
+ {
+ InsertAfter(insert);
+ Simplify();
+ return *this;
+ }
+
+ if (m_next->m_op == MULTIPLY) // (a*b) - c == (a-[c/b]) * b
+ insert->operator/=(*m_next);
+ else
+ insert->operator*=(*m_next); // (a/b) - c == (a-[c*b])/b
+
+ if (insert->m_next != NULL) // neither of the above simplified
+ {
+ //Debug("{%s} did not simplify, change back to {%s}", insert->Str().c_str(), a.Str().c_str());
+ insert->operator=(a); // Just add as is
+ insert->m_op = SUBTRACT;
+ ParanoidNumber * n = this;
+ while (n->m_next != NULL)
+ n = n->m_next;
+ n->InsertAfter(insert);
+ }
+ else
+ {
+ InsertAfter(insert);
+ }
+ Simplify();
+ return *this;
+}
+
+ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
+{
+ ParanoidNumber * n = m_next;
+ while (n != NULL)
+ {
+ if (n->m_op == ADD || n->m_op == SUBTRACT)
+ {
+ n->operator*=(a);
+ break;
+ }
+ n = n->m_next;
+ }
+
+ InsertAfter(new ParanoidNumber(a, MULTIPLY));
+ Simplify();
+ return *this;
+}
+
+
+ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
+{
+ ParanoidNumber * n = m_next;
+ while (n != NULL)
+ {
+ if (n->m_op == ADD || n->m_op == SUBTRACT)
+ {
+ n->operator/=(a);
+ break;
+ }
+ n = n->m_next;
+ }
+
+ InsertAfter(new ParanoidNumber(a, DIVIDE));
+ Simplify();
+ return *this;
+}
+
+double ParanoidNumber::ToDouble() const
+{
+ double value = (m_op == SUBTRACT) ? -m_value : m_value;
+ const ParanoidNumber * n = m_next;
+ while (n != NULL)
+ {
+ switch (n->m_op)
+ {
+ case ADD:
+ case SUBTRACT:
+ return value + n->ToDouble();
+ break;
+ case MULTIPLY:
+ value *= n->m_value;
+ break;
+ case DIVIDE:
+ value /= n->m_value;
+ break;
+ }
+ n = n->m_next;
+ }
+ return value;
+}
+
+void ParanoidNumber::Simplify()
+{
+ ParanoidNumber * n = m_next;
+ ParanoidNumber * p = this;
+
+ while (n != NULL)
+ {
+
+ float a = p->m_value;
+ switch (n->m_op)
+ {
+ case ADD:
+ n->Simplify();
+ feclearexcept(FE_ALL_EXCEPT);
+ a += n->m_value;
+ break;
+ case SUBTRACT:
+ n->Simplify();
+ feclearexcept(FE_ALL_EXCEPT);
+ a -= n->m_value;
+ break;
+ case MULTIPLY:
+ feclearexcept(FE_ALL_EXCEPT);
+ a *= n->m_value;
+ break;
+ case DIVIDE:
+ feclearexcept(FE_ALL_EXCEPT);
+ a /= n->m_value;
+ break;
+
+ }
+ // can't merge p and n
+ if (fetestexcept(FE_ALL_EXCEPT))
+ {
+ while (n != NULL && n->m_op != ADD && n->m_op != SUBTRACT)
+ n = n->m_next;
+ if (n != NULL)
+ n->Simplify();
+ return;
+ }
+ else
+ {
+ // merge n into p
+ p->m_value = a;
+ p->m_next = n->m_next;
+ n->m_next = NULL;
+ delete n;
+ n = p->m_next;
+ }
+ //Debug(" -> {%s}", Str().c_str());
+ }
+}
+
+string ParanoidNumber::Str() const
+{
+ string result("");
+ switch (m_op)
+ {
+ case ADD:
+ result += " +";
+ break;
+ case SUBTRACT:
+ result += " -";
+ break;
+ case MULTIPLY:
+ result += " *";
+ break;
+ case DIVIDE:
+ result += " /";
+ break;
+ }
+ stringstream s("");
+ s << m_value;
+ result += s.str();
+ if (m_next != NULL)
+ result += m_next->Str();
+ return result;
+}
+
+}
--- /dev/null
+#ifndef _PARANOIDNUMBER_H
+#define _PARANOIDNUMBER_H
+
+#include <list>
+#include <cfloat>
+#include <map>
+#include <string>
+
+namespace IPDF
+{
+ 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(NULL)
+ {
+
+ }
+
+ ParanoidNumber(const ParanoidNumber & cpy) : m_value(cpy.m_value), m_op(cpy.m_op), m_next(NULL)
+ {
+ if (cpy.m_next != NULL)
+ {
+ m_next = new ParanoidNumber(*(cpy.m_next));
+ }
+ }
+
+ ParanoidNumber(const ParanoidNumber & cpy, Optype type) : ParanoidNumber(cpy)
+ {
+ m_op = type;
+ }
+
+ ParanoidNumber(const char * str);
+
+ virtual ~ParanoidNumber()
+ {
+ if (m_next != NULL)
+ delete m_next;
+ }
+
+
+ double ToDouble() const;
+
+ 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);
+
+
+ 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));}
+
+ ParanoidNumber operator+(const ParanoidNumber & a) const
+ {
+ ParanoidNumber result(*this);
+ result += a;
+ return result;
+ }
+ ParanoidNumber operator-(const ParanoidNumber & a) const
+ {
+ ParanoidNumber result(*this);
+ result -= a;
+ return result;
+ }
+ ParanoidNumber operator*(const ParanoidNumber & a) const
+ {
+ ParanoidNumber result(*this);
+ result *= a;
+ return result;
+ }
+ ParanoidNumber operator/(const ParanoidNumber & a) const
+ {
+ ParanoidNumber result(*this);
+ result /= a;
+ return result;
+ }
+
+ std::string Str() const;
+
+ private:
+ void Simplify();
+ ParanoidNumber * InsertAfter(ParanoidNumber * insert);
+ ParanoidNumber * InsertAfter(float value, Optype op);
+
+ float m_value;
+ Optype m_op;
+ ParanoidNumber * m_next;
+
+
+
+
+ };
+
+
+}
+
+#endif //_PARANOIDNUMBER_H
+
for (unsigned i = 0; i < x_ints.size(); ++i)
{
if (debug)
- Debug("X Intersection %u at %f,%f vs %f,%f", i,x_ints[i].x, x_ints[i].y, pt.x, pt.y);
+ Debug("X Intersection %u at %f,%f vs %f,%f", i,Double(x_ints[i].x), Double(x_ints[i].y), Double(pt.x), Double(pt.y));
if (x_ints[i].y >= pt.y)
{
for (unsigned i = 0; i < y_ints.size(); ++i)
{
if (debug)
- Debug("Y Intersection %u at %f,%f vs %f,%f", i,x_ints[i].x, x_ints[i].y, pt.x, pt.y);
+ Debug("Y Intersection %u at %f,%f vs %f,%f", i,Double(y_ints[i].x), Double(y_ints[i].y), Double(pt.x), Double(pt.y));
if (y_ints[i].x >= pt.x)
{
template <> Arbint Tabs(const Arbint & a);
template <> Gmpint Tabs(const Gmpint & a);
+
/* Recursive version of GCD
template <class T>
T gcd(const T & a, const T & b)
Rational & operator-=(const Rational & r) {this->operator=(*this-r); return *this;}
Rational & operator*=(const Rational & r) {this->operator=(*this*r); return *this;}
Rational & operator/=(const Rational & r) {this->operator=(*this/r); return *this;}
+ Rational Sqrt() const
+ {
+ return Rational(sqrt(ToDouble()));
+ }
+
+ int64_t ToInt64() const
+ {
+ return (int64_t)ToDouble();
+ }
double ToDouble() const
{
};
+template <class P>
+Rational<P> Abs(const Rational<P> & a)
+{
+ return Rational<P>(Tabs(a.P), a.Q);
+}
}
typedef Rational<ARBINT> Real;
inline float Float(const Real & r) {return (float)r.ToDouble();}
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();}
#else
#error "Type of Real unspecified."
inline double Double(double f) {return (double)f;}
inline double Double(long double f) {return (double)(f);}
inline double Sqrt(double f) {return sqrt(f);}
+ inline double Abs(double a) {return fabs(a);}
+ inline int64_t Int64(double a){return (int64_t)a;}
inline Real Power(const Real & a, int n)
{
--- /dev/null
+#include "main.h"
+#include "real.h"
+#include <cmath>
+#include <cassert>
+#include <list>
+#include <bitset>
+#include <iostream>
+#include <cfloat>
+#include <fenv.h>
+
+using namespace std;
+using namespace IPDF;
+
+int main(int argc, char ** argv)
+{
+ Debug("FLT_MAX = %.40f", FLT_MAX);
+ Debug("FLT_MIN = %.40f", FLT_MIN);
+ Debug("FLT_EPSILON = %.40f", FLT_EPSILON);
+
+ while (true)
+ {
+ float a; float b;
+ char op;
+ feclearexcept(FE_ALL_EXCEPT);
+ cin >> a >> op >> b;
+
+ float c;
+ feclearexcept(FE_ALL_EXCEPT);
+ switch (op)
+ {
+ case '+':
+ c = a + b;
+ break;
+ case '-':
+ c = a - b;
+ break;
+ case '*':
+ c = a * b;
+ break;
+ case '/':
+ c = a / b;
+ break;
+ }
+ if (fetestexcept(FE_OVERFLOW))
+ Debug("Overflow occured");
+ else if (fetestexcept(FE_UNDERFLOW))
+ Debug("Underflow occured");
+ else if (fetestexcept(FE_INEXACT))
+ Debug("Inexact result");
+
+ printf("%.40f\n", c);
+
+
+ }
+}
--- /dev/null
+#include "main.h"
+#include "real.h"
+#include <cmath>
+#include <cassert>
+#include <list>
+#include <bitset>
+#include <iostream>
+#include <cfloat>
+#include <fenv.h>
+#include "paranoidnumber.h"
+
+using namespace std;
+using namespace IPDF;
+
+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", a.Str().c_str());
+ cout << "0.3 ";
+ float fa = 0.3;
+ double da = 0.3;
+ while (cin.good())
+ {
+ char op;
+ double db;
+ cin >> op >> db;
+ float fb(db);
+ ParanoidNumber b(fb);
+ switch (op)
+ {
+ case '+':
+ a += b;
+ fa += fb;
+ da += db;
+ break;
+ case '-':
+ a -= b;
+ fa -= fb;
+ da -= db;
+ break;
+ case '*':
+ a *= b;
+ fa *= fb;
+ da *= db;
+ break;
+ case '/':
+ a /= b;
+ fa /= fb;
+ da /= db;
+ break;
+ }
+
+ Debug("a is: %s", a.Str().c_str());
+ Debug("a as double: %.40f\n", a.ToDouble());
+ Debug("floats give: %.40f\n", fa);
+ Debug("double gives: %.40f\n", da);
+
+
+ }
+}