More debugging, harder realops test

diff --git a/src/arbint.cpp b/src/arbint.cpp
index 582d99a..666b10e 100644 (file)
--- a/src/arbint.cpp
+++ b/src/arbint.cpp
@@ -74,13 +74,29 @@ void Arbint::Zero()
 
 unsigned Arbint::Shrink()
 {
 
 unsigned Arbint::Shrink()
 {

+       

        if (m_digits.size() <= 1)
                return 0;
        if (m_digits.size() <= 1)
                return 0;

-       unsigned i;
-       for (i = m_digits.size()-1; (i > 0 && m_digits[i] != 0L); --i);
-       unsigned result = m_digits.size() - i;
-       m_digits.resize(i);
-       return result;
+       
+       unsigned shrunk = 0;
+       while (m_digits.size() > 1 && m_digits[m_digits.size()-1] == 0L)
+       {
+               Debug("Shrink 1");
+               m_digits.pop_back();
+               shrunk++;
+       }
+       return shrunk;
+}
+
+void Arbint::GrowDigit(digit_t new_msd)
+{
+       static unsigned total_grows = 0;
+       m_digits.push_back(new_msd);
+       Warn("Arbint grows digit (%.16lx), this->m_digits.size() = %u, total grown = %u", new_msd, m_digits.size(), ++total_grows);
+       if (total_grows++ > 10000)
+       {
+               Fatal("Too many GrowDigit calls!");
+       }

 }
 
 Arbint & Arbint::operator*=(const Arbint & mul)
 }
 
 Arbint & Arbint::operator*=(const Arbint & mul)


@@ -101,12 +117,14 @@ Arbint & Arbint::operator*=(const Arbint & mul)
                digit_t carry = add_digits((digit_t*)new_digits.data(), step.data(), step.size());
                if (carry != 0L)
                {
                digit_t carry = add_digits((digit_t*)new_digits.data(), step.data(), step.size());
                if (carry != 0L)
                {

-                       new_digits.push_back(carry);
+                       Debug("Add carry");
+                       GrowDigit(carry);

                }
        }
        
        m_digits.swap(new_digits);
        m_sign = !(m_sign == mul.m_sign);
                }
        }
        
        m_digits.swap(new_digits);
        m_sign = !(m_sign == mul.m_sign);

+       Shrink();

        return *this;
 }
 
        return *this;
 }
 


@@ -119,13 +137,15 @@ void Arbint::Division(const Arbint & div, Arbint & result, Arbint & remainder) c
                result = *this;
                return;
        }
                result = *this;
                return;
        }

+       /* may break things (even more that is)

        else if (div.m_digits.size() == 1)
        {
                result.m_digits.resize(m_digits.size(), 0L);
                remainder = Arbint(div_digits((digit_t*)&m_digits[0], div.m_digits[0], m_digits.size(), result.m_digits.data()));
                result.m_sign = !(m_sign == div.m_sign);
                return;
        else if (div.m_digits.size() == 1)
        {
                result.m_digits.resize(m_digits.size(), 0L);
                remainder = Arbint(div_digits((digit_t*)&m_digits[0], div.m_digits[0], m_digits.size(), result.m_digits.data()));
                result.m_sign = !(m_sign == div.m_sign);
                return;

-       }
+       } */
+

        for (int i = 8*sizeof(digit_t)*m_digits.size(); i >= 0; --i)
        {
                remainder <<= 1;
        for (int i = 8*sizeof(digit_t)*m_digits.size(); i >= 0; --i)
        {
                remainder <<= 1;


@@ -140,6 +160,7 @@ void Arbint::Division(const Arbint & div, Arbint & result, Arbint & remainder) c
                }
        }
        result.m_sign = !(m_sign == div.m_sign);
                }
        }
        result.m_sign = !(m_sign == div.m_sign);

+       result.Shrink();

 }
 
 Arbint & Arbint::operator+=(const Arbint & add)
 }
 
 Arbint & Arbint::operator+=(const Arbint & add)


@@ -162,6 +183,7 @@ Arbint & Arbint::operator+=(const Arbint & add)
                // a + -b == a - b
                SubBasic(add);
        }
                // a + -b == a - b
                SubBasic(add);
        }

+       Shrink();

        return *this;
 }
 
        return *this;
 }
 


@@ -174,29 +196,45 @@ Arbint & Arbint::operator-=(const Arbint & sub)
 
 Arbint & Arbint::AddBasic(const Arbint & add)
 {
 
 Arbint & Arbint::AddBasic(const Arbint & add)
 {

-       if (add.m_digits.size() >= m_digits.size())
+       while (m_digits.size() < add.m_digits.size())

        {
        {

-               m_digits.resize(add.m_digits.size()+1,0L);
+               
+               Debug("Size is %u, add's size is %u", m_digits.size(), add.m_digits.size());
+               GrowDigit(0L);

        }
        }

+       //m_digits.resize(add.m_digits.size()+1,0L);

        
        digit_t carry = add_digits((digit_t*)m_digits.data(), 
                        (digit_t*)add.m_digits.data(), add.m_digits.size());
        if (carry != 0L)
        
        digit_t carry = add_digits((digit_t*)m_digits.data(), 
                        (digit_t*)add.m_digits.data(), add.m_digits.size());
        if (carry != 0L)

-               m_digits[m_digits.size()-1] = carry;
-       else if (m_digits.back() == 0L)
-               m_digits.resize(m_digits.size()-1);
+       {
+               Debug("Grow carry %lu", carry);
+               GrowDigit(carry);
+       }
+       Shrink();

        return *this;
 }
 
 Arbint & Arbint::SubBasic(const Arbint & sub)
 {
        return *this;
 }
 
 Arbint & Arbint::SubBasic(const Arbint & sub)
 {

-       if (sub.m_digits.size() >= m_digits.size())
+       bool fith = false;
+       while (sub.m_digits.size() > m_digits.size())

        {
        {

-               m_digits.resize(sub.m_digits.size(),0L);
+               fith = true;
+               GrowDigit(0L);

        }
        }

+       if (fith)
+       {
+               Debug("START sub was %c%s, I am %c%s", SignChar(), sub.DigitStr().c_str(), SignChar(), DigitStr().c_str());
+       }
+       

        digit_t borrow = sub_digits((digit_t*)m_digits.data(), 
                        (digit_t*)sub.m_digits.data(), sub.m_digits.size());
                
        digit_t borrow = sub_digits((digit_t*)m_digits.data(), 
                        (digit_t*)sub.m_digits.data(), sub.m_digits.size());
                

+       if (fith)
+       {
+               Debug("SUB_DIGITS -> sub was %c%s, I am %c%s", SignChar(), sub.DigitStr().c_str(), SignChar(), DigitStr().c_str());
+       }

                
        //TODO: Write ASM to do this bit?
        if (borrow != 0L)
                
        //TODO: Write ASM to do this bit?
        if (borrow != 0L)


@@ -204,10 +242,19 @@ Arbint & Arbint::SubBasic(const Arbint & sub)
                m_sign = !m_sign;
                for (unsigned i = 0; i < m_digits.size(); ++i)
                        m_digits[i] = (~m_digits[i]);
                m_sign = !m_sign;
                for (unsigned i = 0; i < m_digits.size(); ++i)
                        m_digits[i] = (~m_digits[i]);

-               std::vector<digit_t> one_digits(m_digits.size(), 0L);
-               one_digits[0] = 1;
+               vector<digit_t> one_digits(m_digits.size(), 0L);
+               one_digits[0] = 1L;

                add_digits((digit_t*)m_digits.data(), (digit_t*)one_digits.data(), m_digits.size());
        }
                add_digits((digit_t*)m_digits.data(), (digit_t*)one_digits.data(), m_digits.size());
        }

+       if (fith)
+       {
+               Debug("END -> sub was %c%s, I am %c%s", sub.SignChar(), sub.DigitStr().c_str(), SignChar(), DigitStr().c_str());
+       }
+       Shrink();
+       if (fith)
+       {
+               Debug("SHRUNK -> sub was %c%s, I am %c%s", sub.SignChar(), sub.DigitStr().c_str(), SignChar(), DigitStr().c_str());
+       }

        return *this;
 }
 
        return *this;
 }
 


@@ -303,6 +350,7 @@ Arbint & Arbint::operator>>=(unsigned amount)
                m_digits[i] |= underflow;
                underflow = next_underflow;
        }
                m_digits[i] |= underflow;
                underflow = next_underflow;
        }

+       Shrink();

        return *this;
 }
 
        return *this;
 }
 


@@ -311,7 +359,11 @@ Arbint & Arbint::operator<<=(unsigned amount)
        // Shift by whole number of digits
        unsigned whole = amount/(8*sizeof(digit_t));
        unsigned old_size = m_digits.size();
        // 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);
+       for (unsigned i = 0; i < whole; ++i)
+       {
+               
+               GrowDigit(0L);//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));
        
        memmove(m_digits.data()+whole, m_digits.data(), sizeof(digit_t)*old_size);
        memset(m_digits.data(), 0L, whole*sizeof(digit_t));
        


@@ -338,7 +390,7 @@ Arbint & Arbint::operator<<=(unsigned amount)
        }
        if (overflow != 0L)
                m_digits.push_back(overflow);
        }
        if (overflow != 0L)
                m_digits.push_back(overflow);

-               
+       Shrink();

        return *this;
 }
 
        return *this;
 }
 


@@ -366,10 +418,12 @@ void Arbint::BitClear(unsigned i)
 void Arbint::BitSet(unsigned i)
 {
        unsigned digit = i/(8*sizeof(digit_t));
 void Arbint::BitSet(unsigned i)
 {
        unsigned digit = i/(8*sizeof(digit_t));

-       if (digit >= m_digits.size())
+       while (m_digits.size() < digit+1)

        {
        {

-               m_digits.resize(digit+1, 0L);
+               Debug("Grow BitSet Size %u, digit %u", m_digits.size(), digit);
+               GrowDigit(0L);

        }
        }

+               

        i = i % (8*sizeof(digit_t));
        m_digits[digit] |= (1L << i);           
 }
        i = i % (8*sizeof(digit_t));
        m_digits[digit] |= (1L << i);           
 }

diff --git a/src/arbint.h b/src/arbint.h
index 06008e0..ae0b23e 100644 (file)
--- a/src/arbint.h
+++ b/src/arbint.h
@@ -68,6 +68,13 @@ namespace IPDF
                                return a;
                        }
                        
                                return a;
                        }
                        

+                       inline Arbint operator-() const
+                       {
+                               Arbint a(*this);
+                               a.m_sign = !a.m_sign;
+                               return a;
+                       }
+                       

                        inline Arbint operator*(const Arbint & mul) const
                        {
                                Arbint a(*this);
                        inline Arbint operator*(const Arbint & mul) const
                        {
                                Arbint a(*this);


@@ -99,6 +106,7 @@ namespace IPDF
                        
                        bool operator==(const Arbint & equ) const;
                        bool operator<(const Arbint & less) const;
                        
                        bool operator==(const Arbint & equ) const;
                        bool operator<(const Arbint & less) const;

+                       

 
                        inline bool operator!=(const Arbint & equ) const 
                        {
 
                        inline bool operator!=(const Arbint & equ) const 
                        {


@@ -130,11 +138,11 @@ namespace IPDF
                                return result;
                        }
                        bool IsZero() const;
                                return result;
                        }
                        bool IsZero() const;

-                       
-                       inline operator double() const 
+
+                       inline operator double() const

                        {
                                double acc = 0;
                        {
                                double acc = 0;

-                               for(int i = m_digits.size()-1; i >= 0; --i)
+                               for (int i = m_digits.size()-1; i >= 0; --i)

                                {
                                        acc += (double)m_digits[i];
                                        acc *= (double)UINT64_MAX + 1.0;
                                {
                                        acc += (double)m_digits[i];
                                        acc *= (double)UINT64_MAX + 1.0;


@@ -153,6 +161,8 @@ namespace IPDF
                        Arbint & AddBasic(const Arbint & add);
                        Arbint & SubBasic(const Arbint & sub);
                        
                        Arbint & AddBasic(const Arbint & add);
                        Arbint & SubBasic(const Arbint & sub);
                        

+                       void GrowDigit(digit_t new_msd); // add a new most significant digit
+                       

                        bool GetBit(unsigned i) const;
                        void BitClear(unsigned i);
                        void BitSet(unsigned i);
                        bool GetBit(unsigned i) const;
                        void BitClear(unsigned i);
                        void BitSet(unsigned i);

diff --git a/src/rational.h b/src/rational.h
index 2c98e1a..61a38ab 100644 (file)
--- a/src/rational.h
+++ b/src/rational.h
@@ -17,7 +17,6 @@ template <class T> T Tabs(const T & a)
 {
        return abs(a);
 }
 {
        return abs(a);
 }

-

 template <> Arbint Tabs(const Arbint & a);
 
 
 template <> Arbint Tabs(const Arbint & a);
 
 


@@ -103,7 +102,6 @@ struct Rational
                        P = Q = T(1);
                        return;
                }
                        P = Q = T(1);
                        return;
                }

-               

                T g = gcd(Tabs(P), Tabs(Q));
                //Debug("Got gcd!");
                P /= g;
                T g = gcd(Tabs(P), Tabs(Q));
                //Debug("Got gcd!");
                P /= g;


@@ -135,7 +133,7 @@ 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);
        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(),"-");
+               //result.CheckAccuracy(ToDouble() - r.ToDouble(),"-");

                return result;
        }
        Rational operator*(const Rational & r) const 
                return result;
        }
        Rational operator*(const Rational & r) const 


@@ -163,13 +161,17 @@ struct Rational
        //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(*this); r.P = -r.P;}

        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;}
 
        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);}
+       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);
        bool CheckAccuracy(double d, const char * msg, double threshold = 1e-3) const
        {
                double result = fabs(ToDouble() - d);

diff --git a/src/tests/arb.cpp b/src/tests/arb.cpp
index 97e4930..ffbf74e 100644 (file)
--- a/src/tests/arb.cpp
+++ b/src/tests/arb.cpp
@@ -85,6 +85,11 @@ void TestConstructFromDouble(unsigned cases)
        }
 }
 
        }
 }
 

+void TestSubtraction()
+{
+       
+}
+

 int main(int argc, char ** argv)
 {
        Debug("Shift testing...");
 int main(int argc, char ** argv)
 {
        Debug("Shift testing...");

diff --git a/src/tests/realops.cpp b/src/tests/realops.cpp
index 400ce86..4ac54ff 100644 (file)
--- a/src/tests/realops.cpp
+++ b/src/tests/realops.cpp
@@ -4,7 +4,7 @@
 using namespace std;
 using namespace IPDF;
 
 using namespace std;
 using namespace IPDF;
 

-#define TEST_CASES 100
+#define TEST_CASES 1000

 
 static double g_totalerror = 0;
 
 
 static double g_totalerror = 0;
 


@@ -22,8 +22,9 @@ int main(int argc, char ** argv)
        unsigned failures = 0;
        for (unsigned i = 0; i < TEST_CASES; ++i)
        {
        unsigned failures = 0;
        for (unsigned i = 0; i < TEST_CASES; ++i)
        {

-               double da = (double)(rand()%100 + 1) / (double)(rand()%100 + 1);
-               double db = (double)(rand()%100 + 1) / (double)(rand()%100 + 1);
+               Debug("Test %u of %u", i, TEST_CASES);
+               double da = (double)(rand() + 1) / (double)(rand() + 1);
+               double db = (double)(rand() + 1) / (double)(rand() + 1);

                
                if (rand() % 2 == 0)
                        da = -da;
                
                if (rand() % 2 == 0)
                        da = -da;