From: Sam Moore Date: Sun, 6 Jul 2014 11:44:50 +0000 (+0800) Subject: Improved "realops" tester X-Git-Url: https://git.ucc.asn.au/?p=ipdf%2Fcode.git;a=commitdiff_plain;h=a182cbba4fe9e6a36c5063735dbec1c5340da04c Improved "realops" tester And there are signs of serious problems in Arbint too! --- diff --git a/src/arbint.cpp b/src/arbint.cpp index 6ea4636..59e2252 100644 --- a/src/arbint.cpp +++ b/src/arbint.cpp @@ -17,10 +17,20 @@ #include #include +#include "rational.h" + using namespace std; namespace IPDF { + +/** Absolute value hackery **/ +template <> Arbint Tabs(const Arbint & a) +{ + //Debug("Called"); + return a.Abs(); +} + Arbint::Arbint(int64_t i) : m_digits(1), m_sign(i < 0) { diff --git a/src/arbint.h b/src/arbint.h index 5fefa80..0f047c9 100644 --- a/src/arbint.h +++ b/src/arbint.h @@ -129,6 +129,8 @@ namespace IPDF //inline operator int() const {return int(AsDigit());} unsigned Shrink(); + + inline Arbint Abs() const {Arbint a(*this); a.m_sign = false; return a;} private: Arbint & AddBasic(const Arbint & add); @@ -145,6 +147,8 @@ namespace IPDF }; + + extern "C" { typedef uint64_t digit_t; @@ -154,5 +158,7 @@ extern "C" digit_t div_digits(digit_t * dst, digit_t div, digit_t size, digit_t * rem); } + + } #endif //_ARBINT_H diff --git a/src/rational.h b/src/rational.h index 9dde43a..e74567d 100644 --- a/src/rational.h +++ b/src/rational.h @@ -8,9 +8,17 @@ #include "common.h" #include #include +#include "arbint.h" namespace IPDF { + +template T Tabs(const T & a) +{ + return llabs(a); +} + + /* Recursive version of GCD template @@ -65,7 +73,7 @@ struct Rational Rational(double d=0) : P(d*1e6), Q(1e6) // Possibly the worst thing ever... { Simplify(); - if (!CheckAccuracy(d, "Construct from double")) + //if (!CheckAccuracy(d, "Construct from double")) { //Fatal("Bwah bwah :("); } @@ -93,7 +101,7 @@ struct Rational Q = T(1); return; } - T g = gcd(T(llabs(P)),T(llabs(Q))); + T g = gcd(Tabs(P), Tabs(Q)); //Debug("Got gcd!"); P /= g; Q /= g; @@ -115,10 +123,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); - 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()); - } + //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 @@ -130,19 +138,19 @@ struct Rational Rational operator*(const Rational & r) const { Rational result(P * r.P, Q * r.Q); - 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()); - } + //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 { Rational result(P * r.Q, Q*r.P); - 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()); - } + //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; } @@ -192,6 +200,7 @@ inline Rational pow(const Rational & a, const Rational a(1, 2); + Rational b(3, 9); + Rational c(b); - Rect rect(0,0,1,1); + Rational d(c); + c *= a; + Debug("%s * %s = %s", d.Str().c_str(), a.Str().c_str(), c.Str().c_str()); return 0; - for (unsigned i = 0; i < TEST_CASES; ++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()); - - - - } - printf("Tests done"); + } diff --git a/src/tests/realops.cpp b/src/tests/realops.cpp index 2ccba67..0aeeeb6 100644 --- a/src/tests/realops.cpp +++ b/src/tests/realops.cpp @@ -4,23 +4,91 @@ using namespace std; using namespace IPDF; +#define TEST_CASES 100 + +bool NotEqual(double a, double b, double threshold=1e-4) +{ + return (fabs(a-b) > threshold); +} int main(int argc, char ** argv) { srand(time(NULL)); - Real a = Random(100); - Real b = Random(100); - Debug("a = %f", Float(a)); - Debug("b = %f", Float(b)); - Debug("a + b = %f", Float(a + b)); - Debug("a - b = %f", Float(a - b)); - Debug("a * b = %f", Float(a * b)); - Debug("a / b = %f", Float(a / b)); - Debug("a += b => %f", Float(a += b)); - Debug("a -= b => %f", Float(a -= b)); - Debug("a *= b => %f", Float(a *= b)); - Debug("a /= b => %f", Float(a /= b)); - Debug("a = %f", Float(a)); - Debug("b = %f", Float(b)); + + unsigned failures = 0; + for (unsigned i = 0; i < TEST_CASES; ++i) + { + double da = (double)(rand()) / (double)(rand()); + double db = (double)(rand()) / (double)(rand()); + + Real a(da); + Real b(db); + + + + if (NotEqual(Double(a), da)) + { + failures++; + Warn("a != da; %f vs %f", Double(a), da); + } + if (NotEqual(Double(b), db)) + { + failures++; + Warn("b != db; %f vs %f", Double(b), db); + } + if (NotEqual(Double(a+b), da+db)) + { + failures++; + Warn("a + b = %f should be %f", Double(a+b), da+db); + } + if (NotEqual(Double(a-b), da-db)) + { + failures++; + Warn("a - b = %f should be %f", Double(a-b), da-db); + } + if (NotEqual(Double(a*b), da*db)) + { + failures++; + Warn("a * b = %f should be %f", Double(a*b), da*db); + } + if (NotEqual(Double(a/b), da/db)) + { + failures++; + Warn("a / b = %f should be %f", Double(a/b), da/db); + } + + if (NotEqual(Double(a), da)) + { + failures++; + Warn("a has changed after +-*/ from %f to %f", da, Double(a)); + } + if (NotEqual(Double(b), db)) + { + failures++; + Warn("b has changed after +-*/ from %f to %f", db, Double(b)); + } + + if (NotEqual(Double(a+=b), da+=db)) + { + failures++; + Warn("a += b = %f should be %f", Double(a), da); + } + if (NotEqual(Double(a-=b), da-=db)) + { + failures++; + Warn("a -= b = %f should be %f", Double(a), da); + } + if (NotEqual(Double(a*=b), da*=db)) + { + failures++; + Warn("a *= b = %f should be %f", Double(a), da); + } + if (NotEqual(Double(a/=b), da/=db)) + { + failures++; + Warn("a /= b = %f should be %f", Double(a), da); + } + } + Debug("Completed %u test cases with total of %u operations, %u failures", TEST_CASES, 12*TEST_CASES, failures); }