Merge branch 'master' of git.ucc.asn.au:/ipdf/code

Dammit I forgot to `git stash; git pull; git stash apply`

diff --git a/src/arbint.cpp b/src/arbint.cpp
index c24eb8b..59e2252 100644 (file)
--- a/src/arbint.cpp
+++ b/src/arbint.cpp
@@ -17,10 +17,20 @@
 #include <sstream>
 #include <cstdarg>
 
+#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)
 {
@@ -52,20 +62,14 @@ Arbint::Arbint(const vector<digit_t> & digits) : m_digits(digits), m_sign(false)
 
 Arbint & Arbint::operator=(const Arbint & cpy)
 {
-       memmove(m_digits.data(), cpy.m_digits.data(), 
-               sizeof(digit_t)*min(m_digits.size(), cpy.m_digits.size()));
-       if (cpy.m_digits.size() > m_digits.size())
-       {
-               unsigned old_size = m_digits.size();
-               m_digits.resize(cpy.m_digits.size());
-               memset(m_digits.data()+old_size, 0, sizeof(digit_t)*m_digits.size()-old_size);
-       }
+       m_digits = cpy.m_digits;
+       m_sign = cpy.m_sign;
        return *this;
 }
 
 void Arbint::Zero()
 {
-       memset(m_digits.data(), 0, sizeof(digit_t)*m_digits.size());
+       m_digits.resize(1, 0L);
 }
 
 unsigned Arbint::Shrink()
@@ -325,7 +329,7 @@ Arbint & Arbint::operator<<=(unsigned amount)
 bool Arbint::GetBit(unsigned i) const
 {
        unsigned digit = i/(8*sizeof(digit_t));
-       if (digit > m_digits.size())
+       if (digit >= m_digits.size())
                return false;
                
        i = i % (8*sizeof(digit_t));
@@ -337,7 +341,7 @@ bool Arbint::GetBit(unsigned i) const
 void Arbint::BitClear(unsigned i)
 {
        unsigned digit = i/(8*sizeof(digit_t));
-       if (digit > m_digits.size())
+       if (digit >= m_digits.size())
                return;
        i = i % (8*sizeof(digit_t));
        m_digits[digit] &= ~(1L << i);  
@@ -346,7 +350,7 @@ void Arbint::BitClear(unsigned i)
 void Arbint::BitSet(unsigned i)
 {
        unsigned digit = i/(8*sizeof(digit_t));
-       if (digit > m_digits.size())
+       if (digit >= m_digits.size())
        {
                m_digits.resize(digit+1, 0L);
        }

diff --git a/src/arbint.h b/src/arbint.h
index 5fefa80..0f047c9 100644 (file)
--- 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 (file)
--- a/src/rational.h
+++ b/src/rational.h
@@ -8,9 +8,17 @@
 #include "common.h"
 #include <cmath>
 #include <cassert>
+#include "arbint.h"
 
 namespace IPDF
 {
+       
+template <class T> T Tabs(const T & a)
+{
+       return llabs(a);
+}
+
+
 
 /* Recursive version  of GCD
 template <class T>
@@ -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<int64_t> pow(const Rational<int64_t> & a, const Rational<int64_t
 
 
 
+
 }
 
 #endif //_RATIONAL_H

diff --git a/src/tests/arbrationals.cpp b/src/tests/arbrationals.cpp
index 09e7a78..f26314c 100644 (file)
--- a/src/tests/arbrationals.cpp
+++ b/src/tests/arbrationals.cpp
@@ -15,17 +15,14 @@ using namespace IPDF;
 
 int main(int argc, char ** argv)
 {
+
+       Rational<Arbint> a(1, 2);
+       Rational<Arbint> b(3, 9);
+       Rational<Arbint> c(b);
        
-       Rect rect(0,0,1,1);
+       Rational<Arbint> 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<Arbint> 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 (file)
--- 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);
 
 }