From 4b4b30392a02187ae18aa76bdd48b2cf54994f67 Mon Sep 17 00:00:00 2001 From: Tomasz Sowa Date: Sun, 1 Nov 2009 20:26:01 +0000 Subject: [PATCH] changed: algorithms in Big::Sqrt() and ttmath::Root(x ; n) they were not too much accurate for some integers e.g. Root(16;4) returned a value very closed to 2 (not exactly 2) git-svn-id: svn://ttmath.org/publicrep/ttmath/trunk@231 e52654a7-88a9-db11-a3e9-0013d4bc506e --- CHANGELOG | 6 +- ttmath/ttmath.h | 304 ++++++++++++++++++++++++++------------------- ttmath/ttmathbig.h | 33 ++++- 3 files changed, 207 insertions(+), 136 deletions(-) diff --git a/CHANGELOG b/CHANGELOG index 10c2f3c..0f6906b 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -45,12 +45,16 @@ Version 0.9.0 prerelease (2009.11.01): uint FromString(const std::string & string, const Conv & conv, const wchar_t **, bool *) uint FromString(const std::wstring & string, const Conv & conv, const wchar_t **, bool *) * added: UInt::Sqrt() - a new algorithm for calculating the square root - * added: to the parser: function frac() - remains fraction + * added: to the parser: function frac() - returns a value without the integer part + (only fraction remains) * changed: Factorial() is using the Gamma() function now * changed: Big::Div(ss2) Big::Mod(ss2) they return 2 when ss2 is zero previously returned 1 + * changed: algorithms in Big::Sqrt() and ttmath::Root(x ; n) + they were not too much accurate for some integers + e.g. Root(16;4) returned a value very closed to 2 (not exactly 2) * removed: Parser<>::SetFactorialMax() method the factorial() is such a fast now that we don't need the method longer * removed: ErrorCode::err_too_big_factorial diff --git a/ttmath/ttmath.h b/ttmath/ttmath.h index 94a0eb4..2d64dc2 100644 --- a/ttmath/ttmath.h +++ b/ttmath/ttmath.h @@ -1856,146 +1856,189 @@ namespace ttmath namespace auxiliaryfunctions { - template - bool RootCheckIndexSign(ValueType & x, const ValueType & index, ErrorCode * err) + template + bool RootCheckIndexSign(ValueType & x, const ValueType & index, ErrorCode * err) + { + if( index.IsSign() ) { - if( index.IsSign() ) - { - // index cannot be negative - if( err ) - *err = err_improper_argument; + // index cannot be negative + if( err ) + *err = err_improper_argument; - x.SetNan(); + x.SetNan(); - return true; - } - - return false; + return true; } - - template - bool RootCheckIndexZero(ValueType & x, const ValueType & index, ErrorCode * err) - { - if( index.IsZero() ) - { - if( x.IsZero() ) - { - // there isn't root(0;0) - we assume it's not defined - if( err ) - *err = err_improper_argument; - - x.SetNan(); - - return true; - } - - // root(x;0) is 1 (if x!=0) - x.SetOne(); - - if( err ) - *err = err_ok; - - return true; - } - - return false; - } + return false; + } - template - bool RootCheckIndexOne(const ValueType & index, ErrorCode * err) - { - ValueType one; - one.SetOne(); - - if( index == one ) - { - //root(x;1) is x - // we do it because if we used the PowFrac function - // we would lose the precision - if( err ) - *err = err_ok; - - return true; - } - - return false; - } - - - template - bool RootCheckIndexFrac(ValueType & x, const ValueType & index, ErrorCode * err) - { - if( !index.IsInteger() ) - { - // index must be integer - if( err ) - *err = err_improper_argument; - - x.SetNan(); - - return true; - } - - return false; - } - - - template - bool RootCheckXZero(ValueType & x, ErrorCode * err) + template + bool RootCheckIndexZero(ValueType & x, const ValueType & index, ErrorCode * err) + { + if( index.IsZero() ) { if( x.IsZero() ) { - // root(0;index) is zero (if index!=0) - // RootCheckIndexZero() must be called beforehand - x.SetZero(); - + // there isn't root(0;0) - we assume it's not defined if( err ) - *err = err_ok; + *err = err_improper_argument; + + x.SetNan(); return true; } + + // root(x;0) is 1 (if x!=0) + x.SetOne(); - return false; - } - - - template - bool RootCheckIndex(ValueType & x, const ValueType & index, ErrorCode * err, bool * change_sign) - { - *change_sign = false; - - if( index.Mod2() ) - { - // index is odd (1,3,5...) - if( x.IsSign() ) - { - *change_sign = true; - x.Abs(); - } - } - else - { - // index is even - // x cannot be negative - if( x.IsSign() ) - { - if( err ) - *err = err_improper_argument; - - x.SetNan(); - - return true; - } - } - - return false; + if( err ) + *err = err_ok; + + return true; } + return false; } + template + bool RootCheckIndexOne(const ValueType & index, ErrorCode * err) + { + ValueType one; + one.SetOne(); + + if( index == one ) + { + //root(x;1) is x + // we do it because if we used the PowFrac function + // we would lose the precision + if( err ) + *err = err_ok; + + return true; + } + + return false; + } + + + template + bool RootCheckIndexTwo(ValueType & x, const ValueType & index, ErrorCode * err) + { + if( index == 2 ) + { + x = Sqrt(x, err); + + return true; + } + + return false; + } + + + template + bool RootCheckIndexFrac(ValueType & x, const ValueType & index, ErrorCode * err) + { + if( !index.IsInteger() ) + { + // index must be integer + if( err ) + *err = err_improper_argument; + + x.SetNan(); + + return true; + } + + return false; + } + + + template + bool RootCheckXZero(ValueType & x, ErrorCode * err) + { + if( x.IsZero() ) + { + // root(0;index) is zero (if index!=0) + // RootCheckIndexZero() must be called beforehand + x.SetZero(); + + if( err ) + *err = err_ok; + + return true; + } + + return false; + } + + + template + bool RootCheckIndex(ValueType & x, const ValueType & index, ErrorCode * err, bool * change_sign) + { + *change_sign = false; + + if( index.Mod2() ) + { + // index is odd (1,3,5...) + if( x.IsSign() ) + { + *change_sign = true; + x.Abs(); + } + } + else + { + // index is even + // x cannot be negative + if( x.IsSign() ) + { + if( err ) + *err = err_improper_argument; + + x.SetNan(); + + return true; + } + } + + return false; + } + + + template + uint RootCorrectInteger(ValueType & old_x, ValueType & x, const ValueType & index) + { + if( !old_x.IsInteger() || x.IsInteger() || !index.exponent.IsSign() ) + return 0; + + // old_x is integer, + // x is not integer, + // index is relatively small (index.exponent<0 or index.exponent<=0) + // (because we're using a special powering algorithm Big::PowUInt()) + + uint c = 0; + + ValueType temp(x); + c += temp.Round(); + + ValueType temp_round(temp); + c += temp.PowUInt(index); + + if( temp == old_x ) + x = temp_round; + + return (c==0)? 0 : 1; + } + + + + } // namespace auxiliaryfunctions + + + /*! indexth Root of x index must be integer and not negative <0;1;2;3....) @@ -2026,30 +2069,33 @@ namespace ttmath if( RootCheckIndexSign(x, index, err) ) return x; if( RootCheckIndexZero(x, index, err) ) return x; if( RootCheckIndexOne ( index, err) ) return x; + if( RootCheckIndexTwo (x, index, err) ) return x; if( RootCheckIndexFrac(x, index, err) ) return x; if( RootCheckXZero (x, err) ) return x; // index integer and index!=0 // x!=0 - uint c = 0; + ValueType old_x(x); bool change_sign; + if( RootCheckIndex(x, index, err, &change_sign ) ) return x; - ValueType newindex; - newindex.SetOne(); - c += newindex.Div(index); - c += x.PowFrac(newindex); // here can only be a carry + ValueType temp; + uint c = 0; + + // we're using the formula: root(x ; n) = exp( ln(x) / n ) + c += temp.Ln(x); + c += temp.Div(index); + c += x.Exp(temp); if( change_sign ) { - // the value of x should be different from zero - // (x is actually tested by RootCheckXZero) - TTMATH_ASSERT( x.IsZero() == false ) - + // x is different from zero x.SetSign(); } + c += RootCorrectInteger(old_x, x, index); if( err ) *err = c ? err_overflow : err_ok; diff --git a/ttmath/ttmathbig.h b/ttmath/ttmathbig.h index 46b1b8e..f2653f4 100644 --- a/ttmath/ttmathbig.h +++ b/ttmath/ttmathbig.h @@ -1511,8 +1511,7 @@ public: /*! this function calculates the square root - - Sqrt(9) = 3 + e.g. let this=9 then this.Sqrt() gives 3 return: 0 - ok 1 - carry @@ -1529,11 +1528,33 @@ public: if( IsZero() ) return 0; - Big pow; - pow.Set05(); + Big old(*this); + Big ln; + uint c = 0; - // PowFrac can return only a carry because x is greater than zero - return PowFrac(pow); + // we're using the formula: sqrt(x) = e ^ (ln(x) / 2) + c += ln.Ln(*this); + c += ln.exponent.SubOne(); // ln = ln / 2 + c += Exp(ln); + + // above formula doesn't give accurate results for some integers + // e.g. Sqrt(81) would not be 9 but a value very closed to 9 + // we're rounding the result, calculating result*result and comparing + // with the old value, if they are equal then the result is an integer too + + if( !c && old.IsInteger() && !IsInteger() ) + { + Big temp(*this); + c += temp.Round(); + + Big temp2(temp); + c += temp.Mul(temp2); + + if( temp == old ) + *this = temp2; + } + + return CheckCarry(c); }