From b02dcbab39b8c28b9baa41436842ca9fe4ae7ffd Mon Sep 17 00:00:00 2001 From: Sam Moore Date: Sun, 6 Jul 2014 22:41:19 +0800 Subject: [PATCH] More debugging, harder realops test --- src/arbint.cpp | 94 ++++++++++++++++++++++++++++++++++--------- src/arbint.h | 16 ++++++-- src/rational.h | 10 +++-- src/tests/arb.cpp | 5 +++ src/tests/realops.cpp | 7 ++-- 5 files changed, 102 insertions(+), 30 deletions(-) diff --git a/src/arbint.cpp b/src/arbint.cpp index 582d99a..666b10e 100644 --- a/src/arbint.cpp +++ b/src/arbint.cpp @@ -74,13 +74,29 @@ void Arbint::Zero() unsigned Arbint::Shrink() { + 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) @@ -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) { - new_digits.push_back(carry); + Debug("Add carry"); + GrowDigit(carry); } } m_digits.swap(new_digits); m_sign = !(m_sign == mul.m_sign); + Shrink(); return *this; } @@ -119,13 +137,15 @@ void Arbint::Division(const Arbint & div, Arbint & result, Arbint & remainder) c 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; - } + } */ + 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.Shrink(); } Arbint & Arbint::operator+=(const Arbint & add) @@ -162,6 +183,7 @@ Arbint & Arbint::operator+=(const Arbint & add) // a + -b == a - b SubBasic(add); } + Shrink(); return *this; } @@ -174,29 +196,45 @@ Arbint & Arbint::operator-=(const Arbint & sub) 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) - 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) { - 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()); + 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) @@ -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]); - std::vector one_digits(m_digits.size(), 0L); - one_digits[0] = 1; + vector 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()); } + 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; } @@ -303,6 +350,7 @@ Arbint & Arbint::operator>>=(unsigned amount) m_digits[i] |= underflow; underflow = next_underflow; } + Shrink(); 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(); - 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)); @@ -338,7 +390,7 @@ Arbint & Arbint::operator<<=(unsigned amount) } if (overflow != 0L) m_digits.push_back(overflow); - + Shrink(); return *this; } @@ -366,10 +418,12 @@ void Arbint::BitClear(unsigned i) 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); } diff --git a/src/arbint.h b/src/arbint.h index 06008e0..ae0b23e 100644 --- a/src/arbint.h +++ b/src/arbint.h @@ -68,6 +68,13 @@ namespace IPDF 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); @@ -99,6 +106,7 @@ namespace IPDF bool operator==(const Arbint & equ) const; bool operator<(const Arbint & less) const; + inline bool operator!=(const Arbint & equ) const { @@ -130,11 +138,11 @@ namespace IPDF return result; } bool IsZero() const; - - inline operator double() const + + inline operator double() const { 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; @@ -153,6 +161,8 @@ namespace IPDF 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); diff --git a/src/rational.h b/src/rational.h index 2c98e1a..61a38ab 100644 --- a/src/rational.h +++ b/src/rational.h @@ -17,7 +17,6 @@ template T Tabs(const T & a) { return abs(a); } - template <> Arbint Tabs(const Arbint & a); @@ -103,7 +102,6 @@ struct Rational P = Q = T(1); return; } - 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); - result.CheckAccuracy(ToDouble() - r.ToDouble(),"-"); + //result.CheckAccuracy(ToDouble() - r.ToDouble(),"-"); 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(*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;} - 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); diff --git a/src/tests/arb.cpp b/src/tests/arb.cpp index 97e4930..ffbf74e 100644 --- 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..."); diff --git a/src/tests/realops.cpp b/src/tests/realops.cpp index 400ce86..4ac54ff 100644 --- a/src/tests/realops.cpp +++ b/src/tests/realops.cpp @@ -4,7 +4,7 @@ using namespace std; using namespace IPDF; -#define TEST_CASES 100 +#define TEST_CASES 1000 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) { - 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; -- 2.20.1