From: Sam Moore Date: Sun, 6 Jul 2014 09:32:28 +0000 (+0800) Subject: Horribly unreliable Arbint's and Rational's X-Git-Url: https://git.ucc.asn.au/?p=ipdf%2Fcode.git;a=commitdiff_plain;h=39599aa6423d3e0181fbfe2aac8f78f388a3f372;ds=inline Horribly unreliable Arbint's and Rational's Changes from before: Arbint's use a division algorithm based on bitshifting. So, you might get a random segfault (if you scroll you will definitely get a random segfault) and pretty much every single operation with Rational gives you an accuracy warning... And yet somehow I can render a Bezier. Really slowly. --- diff --git a/src/Makefile b/src/Makefile index 94d6f1f..ec8789a 100644 --- a/src/Makefile +++ b/src/Makefile @@ -3,7 +3,7 @@ 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 graphicsbuffer.o framebuffer.o shaderprogram.o stb_truetype.o gl_core44.o add_digits_asm.o sub_digits_asm.o mul_digits_asm.o arbint.o +OBJ = log.o real.o bezier.o document.o objectrenderer.o view.o screen.o vfpu.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 LIB_x86_64 = ../contrib/lib/libSDL2-2.0.so.0 -lGL LIB_i386 = ../contrib/lib32/libSDL2-2.0.so.0 -lGL diff --git a/src/arbint.cpp b/src/arbint.cpp index 748330f..c24eb8b 100644 --- a/src/arbint.cpp +++ b/src/arbint.cpp @@ -22,7 +22,7 @@ using namespace std; namespace IPDF { -Arbint::Arbint(digit_t i) : m_digits(1), m_sign(i < 0) +Arbint::Arbint(int64_t i) : m_digits(1), m_sign(i < 0) { m_digits[0] = llabs(i); } @@ -31,7 +31,7 @@ Arbint::Arbint(unsigned n, digit_t d0, ...) : m_digits(n), m_sign(false) { va_list ap; va_start(ap, d0); - if (n > 1) + if (n >= 1) m_digits[0] = d0; for (unsigned i = 1; i < n; ++i) { @@ -108,32 +108,21 @@ Arbint & Arbint::operator*=(const Arbint & mul) void Arbint::Division(const Arbint & div, Arbint & result, Arbint & remainder) const { - if (div == Arbint(1L)) + remainder = 0; + result = 0; + for (int i = 8*sizeof(digit_t)*m_digits.size(); i >= 0; --i) { - result = *this; - remainder = 0L; - return; - } - if (div == *this) - { - result = 1L; - remainder = 0L; - return; - } - - - result = 0L; - remainder = *this; - Arbint next_rem(remainder); - while ((next_rem -= div) >= Arbint(0L)) - { - //Debug("%li - %li = %li", digit_t(remainder), digit_t(div), digit_t(next_rem)); - //Debug("Sign is %d", next_rem.m_sign); - remainder = next_rem; - result += 1L; + remainder <<= 1; + if (GetBit(i)) + remainder.BitSet(0); + else + remainder.BitClear(0); + if (remainder >= div) + { + remainder -= div; + result.BitSet(i); + } } - - } Arbint & Arbint::operator+=(const Arbint & add) @@ -206,18 +195,8 @@ Arbint & Arbint::SubBasic(const Arbint & sub) string Arbint::Str(const string & base) const { string s(""); - for (unsigned i = 0; i < m_digits.size(); ++i) - { - digit_t w = m_digits[i]; - do - { - digit_t q = w % 10; - w /= 10; - s += ('0' + q); - } - while (w > 0); - if (i+1 < m_digits.size()) s += ","; - } + Arbint cpy(*this); + reverse(s.begin(), s.end()); return s; } @@ -264,13 +243,115 @@ bool Arbint::operator<(const Arbint & less) const string Arbint::DigitStr() const { stringstream ss(""); - ss << std::hex << std::setfill('0'); + //ss << std::hex << std::setfill('0'); for (unsigned i = 0; i < m_digits.size(); ++i) { if (i != 0) ss << ','; - ss << std::setw(2*sizeof(digit_t)) << static_cast(m_digits[i]); + //ss << std::setw(2*sizeof(digit_t)) << static_cast(m_digits[i]); + ss << static_cast(m_digits[i]); } return ss.str(); } +Arbint & Arbint::operator>>=(unsigned amount) +{ + // Shift by whole number of digits + unsigned whole = amount/(8*sizeof(digit_t)); + unsigned old_size = m_digits.size(); + + if (whole >= old_size) + { + m_digits.resize(1); + m_digits[0] = 0L; + return *this; + } + memmove(m_digits.data(), m_digits.data()+whole, sizeof(digit_t)*(old_size-whole)); + m_digits.resize(old_size-whole, 0L); + + // Shift by partial amount + amount = amount %(8*sizeof(digit_t)); + if (amount == 0) + return *this; + + digit_t underflow = 0L; + for (int i = (int)(m_digits.size()-1); i >= 0; --i) + { + unsigned shl = (8*sizeof(digit_t)-amount); + digit_t next_underflow = (m_digits[i] << shl); + //digit_t mask_upper = ~(0L >> amount); + m_digits[i] = (m_digits[i] >> amount);// & mask_upper; + m_digits[i] |= underflow; + underflow = next_underflow; + } + return *this; +} + +Arbint & Arbint::operator<<=(unsigned amount) +{ + // Shift by whole number of digits + unsigned whole = amount/(8*sizeof(digit_t)); + unsigned old_size = m_digits.size(); + m_digits.resize(m_digits.size() + whole); + memmove(m_digits.data()+whole, m_digits.data(), sizeof(digit_t)*old_size); + memset(m_digits.data(), 0L, whole*sizeof(digit_t)); + + + + // Partial shifting + amount = amount % (8*sizeof(digit_t)); + if (amount == 0) + return *this; + + //Debug("Shift by %u from %u", amount, whole); + digit_t overflow = 0L; + for (unsigned i = whole; i < m_digits.size(); ++i) + { + //Debug("Digit is %.16lx", m_digits[i]); + unsigned shr = (8*sizeof(digit_t)-amount); + //Debug("shr is %u", shr); + digit_t next_overflow = (m_digits[i] >> shr); + //Debug("Next overflow %.16lx", next_overflow); + m_digits[i] <<= amount; + //Debug("Before overflow %.16lx", m_digits[i]); + m_digits[i] |= overflow; + overflow = next_overflow; + } + if (overflow != 0L) + m_digits.push_back(overflow); + + return *this; +} + +bool Arbint::GetBit(unsigned i) const +{ + unsigned digit = i/(8*sizeof(digit_t)); + if (digit > m_digits.size()) + return false; + + i = i % (8*sizeof(digit_t)); + + return (m_digits[digit] & (1L << i)); + +} + +void Arbint::BitClear(unsigned i) +{ + unsigned digit = i/(8*sizeof(digit_t)); + if (digit > m_digits.size()) + return; + i = i % (8*sizeof(digit_t)); + m_digits[digit] &= ~(1L << i); +} + +void Arbint::BitSet(unsigned i) +{ + unsigned digit = i/(8*sizeof(digit_t)); + if (digit > m_digits.size()) + { + m_digits.resize(digit+1, 0L); + } + i = i % (8*sizeof(digit_t)); + m_digits[digit] |= (1L << i); +} + } diff --git a/src/arbint.h b/src/arbint.h index 4b819b9..5fefa80 100644 --- a/src/arbint.h +++ b/src/arbint.h @@ -14,18 +14,18 @@ namespace IPDF class Arbint { public: - typedef int64_t digit_t; + typedef uint64_t digit_t; - Arbint(digit_t i); + Arbint(int64_t i); Arbint(const std::vector & digits); Arbint(unsigned n, digit_t d0, ...); Arbint(const std::string & str, const std::string & base="0123456789"); - ~Arbint() {} + virtual ~Arbint() {} Arbint(const Arbint & cpy); - digit_t AsDigit() const + int64_t AsDigit() const { - digit_t digit = (m_digits.size() == 1) ? m_digits[0] : 0xFFFFFFFFFFFFFFFF; + int64_t digit = (m_digits.size() == 1) ? m_digits[0] : 0x7FFFFFFFFFFFFFFF; return (m_sign) ? -digit : digit; } @@ -45,6 +45,9 @@ namespace IPDF Arbint & operator*=(const Arbint & mul); void Division(const Arbint & div, Arbint & result, Arbint & modulo) const; + Arbint & operator<<=(unsigned amount); + Arbint & operator>>=(unsigned amount); + inline Arbint operator+(const Arbint & add) const { Arbint a(*this); @@ -107,11 +110,22 @@ namespace IPDF { return !(this->operator<=(grea)); } - + inline Arbint operator>>(unsigned amount) const + { + Arbint result(*this); + result >>= amount; + return result; + } + inline Arbint operator<<(unsigned amount) const + { + Arbint result(*this); + result <<= amount; + return result; + } bool IsZero() const; //inline operator double() const {return double(AsDigit());} - inline operator digit_t() const {return AsDigit();} + inline operator int64_t() const {return AsDigit();} //inline operator int() const {return int(AsDigit());} unsigned Shrink(); @@ -119,6 +133,11 @@ namespace IPDF Arbint & AddBasic(const Arbint & add); Arbint & SubBasic(const Arbint & sub); + + bool GetBit(unsigned i) const; + void BitClear(unsigned i); + void BitSet(unsigned i); + std::vector m_digits; bool m_sign; @@ -128,10 +147,11 @@ namespace IPDF extern "C" { - typedef int64_t digit_t; + typedef uint64_t digit_t; digit_t add_digits(digit_t * dst, digit_t * add, digit_t size); digit_t sub_digits(digit_t * dst, digit_t * add, digit_t size); digit_t mul_digits(digit_t * dst, digit_t mul, digit_t size); + digit_t div_digits(digit_t * dst, digit_t div, digit_t size, digit_t * rem); } } diff --git a/src/main.cpp b/src/main.cpp index cd74d6d..9af7bbc 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -82,9 +82,9 @@ int main(int argc, char ** argv) else { doc.AddBezierData(Bezier(0,0,0,1,1,0)); - doc.AddBezierData(Bezier(0,0,1,0,0,1)); - doc.AddBezierData(Bezier(0,0,1,1,1,0)); - doc.AddBezierData(Bezier(0,1,1,0,0,1)); + //doc.AddBezierData(Bezier(0,0,1,0,0,1)); + //doc.AddBezierData(Bezier(0,0,1,1,1,0)); + //doc.AddBezierData(Bezier(0,1,1,0,0,1)); @@ -94,14 +94,15 @@ int main(int argc, char ** argv) for (int y = 0; y < 8; ++y) { //doc.Add(static_cast((x^y)%3), Rect(0.2+x-4.0,0.2+y-4.0,0.6,0.6)); - doc.Add(BEZIER, Rect(0.2+x-4.0, 0.2+y-4.0, 0.6,0.6), (x^y)%3); + //doc.Add(BEZIER, Rect(0.2+x-4.0, 0.2+y-4.0, 0.6,0.6), (x^y)%3); } - - //doc.Add(RECT_OUTLINE, Rect(0.1,0.1,0.8,0.8), 0); - //doc.Add(BEZIER, Rect(0,0,1,1), 0); } + Debug("Make rect"); + doc.Add(BEZIER, Rect(0.1,0.1,0.8,0.8), 0); + Debug("Made rect"); } + Debug("Start!"); Rect bounds(b[0],b[1],b[2],b[3]); if (mode == LOOP) diff --git a/src/rational.h b/src/rational.h index 6269ef1..9dde43a 100644 --- a/src/rational.h +++ b/src/rational.h @@ -65,7 +65,10 @@ struct Rational Rational(double d=0) : P(d*1e6), Q(1e6) // Possibly the worst thing ever... { Simplify(); - CheckAccuracy(d, "Construct from double"); + if (!CheckAccuracy(d, "Construct from double")) + { + //Fatal("Bwah bwah :("); + } } Rational(const T & _P, const T & _Q) : P(_P), Q(_Q) @@ -91,7 +94,7 @@ struct Rational return; } T g = gcd(T(llabs(P)),T(llabs(Q))); - Debug("Got gcd!"); + //Debug("Got gcd!"); P /= g; Q /= g; } @@ -112,7 +115,10 @@ struct Rational Rational operator+(const Rational & r) const { Rational result = (r.P == T(0)) ? Rational(P,Q) : Rational(P*r.Q + r.P*Q, Q*r.Q); - result.CheckAccuracy(ToDouble() + r.ToDouble(),"+"); + if (!result.CheckAccuracy(ToDouble() * r.ToDouble(),"+")) + { + Debug("This is %s (%f) and r is %s (%f)", Str().c_str(), ToDouble(), r.Str().c_str(), r.ToDouble()); + } return result; } Rational operator-(const Rational & r) const @@ -155,10 +161,12 @@ struct Rational double ToDouble() const {return (double)(P) / (double)(Q);} bool CheckAccuracy(double d, const char * msg, double threshold = 1e-3) const { - double result = fabs(ToDouble() - d) / d; + double result = fabs(ToDouble() - d); + if (d != 0e0) result /= d; if (result > threshold) { Warn("(%s) : Rational %s (%f) is not close enough at representing %f (%f vs %f)", msg, Str().c_str(), ToDouble(), d, result, threshold); + Backtrace(); return false; } return true; @@ -166,7 +174,7 @@ struct Rational std::string Str() const { std::stringstream s; - s << (int64_t)P << "/" << (int64_t)Q; + s << int64_t(P) << "/" << int64_t(Q); return s.str(); } diff --git a/src/tests/arb.cpp b/src/tests/arb.cpp index 00bfefb..97e4930 100644 --- a/src/tests/arb.cpp +++ b/src/tests/arb.cpp @@ -10,46 +10,94 @@ using namespace IPDF; #define TEST_CASES 100 -int main(int argc, char ** argv) +void TestShifting(unsigned cases) { - + for (unsigned c = 0; c < cases; ++c) + { + Arbint aa(1u, rand()); + Arbint bb(aa); + for (unsigned i = 0; i < 100*sizeof(Arbint::digit_t); ++i) + { + bb <<= i; + //Debug("i is %u bb is %c%s, a was %c%s", i, bb.SignChar(), bb.DigitStr().c_str(), aa.SignChar(), aa.DigitStr().c_str()); + bb >>= i; + if (bb != aa) + { + Fatal("i is %u bb is %c%s, a was %c%s", i, bb.SignChar(), bb.DigitStr().c_str(), aa.SignChar(), aa.DigitStr().c_str()); + } + } + } +} - for (unsigned i = 0; i < TEST_CASES; ++i) +void TestMaths(unsigned cases) +{ + for (unsigned i = 0; i < cases; ++i) { - Arbint::digit_t a = rand(); - Arbint::digit_t b = rand(); + int64_t a = rand(); + int64_t b = rand(); Arbint arb_a(a); Arbint arb_b(b); - Debug("Test %u: a = %li, b = %li", i, a, b); + //Debug("Test %u: a = %li, b = %li", i, a, b); - Debug("a < b = %d", a < b); - Debug("Arb a b = %d", a > b); - Debug("Arb a>b = %d", arb_a > arb_b); + //Debug("a > b = %d", a > b); + //Debug("Arb a>b = %d", arb_a > arb_b); assert((arb_a > arb_b) == (a > b)); - Debug("a + b = %.16lx", a+b); - Debug("Arb a+b = %s", (arb_a+arb_b).DigitStr().c_str()); + //Debug("a + b = %.16lx", a+b); + //Debug("Arb a+b = %s", (arb_a+arb_b).DigitStr().c_str()); assert((arb_a+arb_b).AsDigit() == a + b); - Debug("a - b = %li %.16lx", (a-b), (a-b)); - Debug("Arb a-b = %s %.16lx", (arb_a-arb_b).DigitStr().c_str(), (arb_a-arb_b).AsDigit()); + //Debug("a - b = %li %.16lx", (a-b), (a-b)); + //Debug("Arb a-b = %s %.16lx", (arb_a-arb_b).DigitStr().c_str(), (arb_a-arb_b).AsDigit()); assert((arb_a-arb_b).AsDigit() == a - b); - Debug("a * b = %.16lx", a*b); - Debug("Arb a*b = %s", (arb_a*arb_b).DigitStr().c_str()); + //Debug("a * b = %.16lx", a*b); + //Debug("Arb a*b = %s", (arb_a*arb_b).DigitStr().c_str()); assert((arb_a*arb_b).AsDigit() == a * b); - Debug("a / b = %.16lx", a/b); - Debug("Arb a/b = %s", (arb_a/arb_b).DigitStr().c_str()); + //Debug("a / b = %.16lx", a/b); + //Debug("Arb a/b = %s", (arb_a/arb_b).DigitStr().c_str()); assert((arb_a/arb_b).AsDigit() == a / b); assert((arb_a%arb_b).AsDigit() == a % b); + } +} + +void TestConstructFromDouble(unsigned cases) +{ + for (unsigned c = 0; c < cases; ++c) + { + int64_t test = (int64_t)(rand()*1e6); + double d = (double)(test); + Arbint a(d); + Arbint b(test); + if (a != b) + { + Fatal("test is %li, d is %f, a is %s, b is %s", test, d, a.DigitStr().c_str(), b.DigitStr().c_str()); + } + + } - printf("All single digit tests successful!\n"); +} + +int main(int argc, char ** argv) +{ + Debug("Shift testing..."); + TestShifting(TEST_CASES); + Debug("Left/Right shift testing succeeded"); + + Debug("Testing +-*/%"); + TestMaths(TEST_CASES); + Debug("All (single digit) maths operator tests successful!"); + + Debug("Testing construct from double"); + TestConstructFromDouble(TEST_CASES); + Debug("Construct from double successful"); return 0; } diff --git a/src/tests/arbrationals.cpp b/src/tests/arbrationals.cpp index 3af559f..09e7a78 100644 --- a/src/tests/arbrationals.cpp +++ b/src/tests/arbrationals.cpp @@ -11,20 +11,19 @@ using namespace std; using namespace IPDF; +#include "rect.h" + int main(int argc, char ** argv) { + + Rect rect(0,0,1,1); + return 0; for (unsigned i = 0; i < TEST_CASES; ++i) { - Debug("Cycle %u",i); - Arbint p1 = 1L;//rand(); - Arbint q1 = 2L;//rand(); - Arbint p2 = rand(); - Arbint q2 = rand(); - Debug("Construct rationals"); - Rational a(p1, q1); - Debug("Constructed a %u\n",i); - Rational b(p2, q2); - Debug("Constructed b %u\n",i); + double d = (double)(rand()) / (double)(rand()); + Rational r(d); + Debug("%f -> %s -> %s/%s", d, r.Str().c_str(), r.P.DigitStr().c_str(), r.Q.DigitStr().c_str()); + }