X-Git-Url: https://git.ucc.asn.au/?p=ipdf%2Fcode.git;a=blobdiff_plain;f=src%2Frational.h;h=b01950533ad6810653c724f8924f68d4f982690c;hp=b3709ad7f47b067f30986c90ab09f00bb17d8d4f;hb=da646c739f87bf28c5a7af2bc180b93b3444321b;hpb=33356addacfe4296ecb613c6c4696f082e351159 diff --git a/src/rational.h b/src/rational.h index b3709ad..b019505 100644 --- a/src/rational.h +++ b/src/rational.h @@ -8,28 +8,46 @@ #include "common.h" #include #include +#include "arbint.h" +#include "gmpint.h" +#include +#include namespace IPDF { + +template T Tabs(const T & a) +{ + return abs(a); +} +template <> Arbint Tabs(const Arbint & a); +template <> Gmpint Tabs(const Gmpint & a); -/** Greatest Common Divisor - Euclid's algorithm **/ - +/* Recursive version of GCD template T gcd(const T & a, const T & b) { - if (a == 1 || a == 0) return 1; - if (b == 0) return a; - if (b == a) return a; + Debug("Called on %li/%li", int64_t(a), int64_t(b)); + if (a == T(1) || a == T(0)) return T(1); + if (b == T(0)) return a; + if (b == a) + { + Debug("Equal!"); + return a; + } + Debug("Not equal!"); if (a > b) return gcd(a-b,b); return gcd(a, b-a); } +*/ -/* +/** Greatest Common Divisor of p and q **/ template T gcd(const T & p, const T & q) { - Debug("p/q = % + + T g(1); T big(p); T small(q); @@ -38,24 +56,26 @@ T gcd(const T & p, const T & q) big = q; small = p; } - if (small == 0) + if (small == T(0)) return g; - while ((g = big % small) > 0) + while ((g = big % small) > T(0)) { + //Debug("big = %li, small = %li", int64_t(big), int64_t(small)); big = small; small = g; + //Debug("Loop %u", ++count); } return small; -} -*/ +} + + template struct Rational { /** Construct from a double.**/ - Rational(double d = 0) : P(d*1e3), Q(1e3) // Possibly the worst thing ever... + Rational(double d=0) : P(d*1e6), Q(1e6) // Possibly the worst thing ever... { Simplify(); - //CheckAccuracy(d, "Construct from double"); } Rational(const T & _P, const T & _Q) : P(_P), Q(_Q) @@ -70,13 +90,23 @@ struct Rational void Simplify() { - if (Q < 0) + if (Q < T(0)) { - P = (P < 0) ? -P : P; + P = -P; Q = -Q; } - - T g = gcd(llabs(P),llabs(Q)); + if (P == T(0)) + { + Q = T(1); + return; + } + if (P == Q) + { + P = Q = T(1); + return; + } + T g = gcd(Tabs(P), Tabs(Q)); + //Debug("Got gcd!"); P /= g; Q /= g; } @@ -94,59 +124,71 @@ struct Rational bool operator>=(const Rational & r) const {return *this == r || *this > r;} bool operator!=(const Rational & r) const {return !(*this == r);} - - - /* Rational operator+(const Rational & r) const { - Rational result = (r.P == 0) ? Rational(P,Q) : Rational(P*r.Q + r.P*Q, Q*r.Q); - result.CheckAccuracy(ToDouble() + r.ToDouble(),"+"); + 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()); + //} return result; } Rational operator-(const Rational & r) const { - Rational result = (r.P == 0) ? Rational(P,Q) : Rational(P*r.Q - r.P*Q, Q*r.Q); - result.CheckAccuracy(ToDouble() - r.ToDouble(),"-"); + Rational result = (r.P == T(0)) ? Rational(P,Q) : Rational(P*r.Q - r.P*Q, Q*r.Q); + //result.CheckAccuracy(ToDouble() - r.ToDouble(),"-"); return result; } - */ 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 = (r.P == 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()); - } + 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()); + //} return result; } - Rational operator+(const Rational & r) const {return Rational(ToDouble()+r.ToDouble());} - Rational operator-(const Rational & r) const {return Rational(ToDouble()-r.ToDouble());} + /** To cheat, use these **/ + //Rational operator+(const Rational & r) const {return Rational(ToDouble()+r.ToDouble());} + //Rational operator-(const Rational & r) const {return Rational(ToDouble()-r.ToDouble());} //Rational operator*(const Rational & r) const {return Rational(ToDouble()*r.ToDouble());} //Rational operator/(const Rational & r) const {return Rational(ToDouble()/r.ToDouble());} - Rational & operator=(const Rational & r) {P = r.P; Q = r.Q; return *this;} + Rational operator-() const {Rational r(*this); r.P = -r.P; return r;} + Rational & operator=(const Rational & r) {P = r.P; Q = r.Q; Simplify(); 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 & operator*=(const Rational & r) {this->operator=(*this*r); return *this;} Rational & operator/=(const Rational & r) {this->operator=(*this/r); return *this;} - double ToDouble() const {return (double)(P) / (double)(Q);} - bool CheckAccuracy(double d, const char * msg, double threshold = 1e-6) const + double ToDouble() const { - double result = fabs(ToDouble() - d) / d; + T num = P, denom = Q; + while (Tabs(num) > T(DBL_MAX)) + { + num /= T(16); + denom /= T(16); + } + return ((double)(num))/((double)(denom)); + } + bool CheckAccuracy(double d, const char * msg, double threshold = 1e-3) const + { + 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; @@ -154,7 +196,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(); } @@ -162,13 +204,8 @@ struct Rational T Q; }; -inline Rational pow(const Rational & a, const Rational & b) -{ - //TODO:Implement properly - int64_t P = std::pow((double)a.P, b.ToDouble()); - int64_t Q = std::pow((double)a.Q, b.ToDouble()); - return Rational(P, Q); -} + + }