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

@@ -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

ttmath/ttmath.h

@@ -1924,6 +1924,20 @@ namespace ttmath
 	}
 
 
+	template<class ValueType>
+	bool RootCheckIndexTwo(ValueType & x, const ValueType & index, ErrorCode * err)
+	{
+		if( index == 2 )
+		{
+			x = Sqrt(x, err);
+
+		return true;
+		}
+
+	return false;
+	}
+
+
 	template<class ValueType>
 	bool RootCheckIndexFrac(ValueType & x, const ValueType & index, ErrorCode * err)
 	{
@@ -1993,9 +2007,38 @@ namespace ttmath
 	return false;
 	}
 
+
+	template<class ValueType>
+	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;

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<exp, man> pow;
-		pow.Set05();
+		Big<exp, man> old(*this);
+		Big<exp, man> 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<exp, man> temp(*this);
+			c += temp.Round();
+
+			Big<exp, man> temp2(temp);
+			c += temp.Mul(temp2);
+
+			if( temp == old )
+				*this = temp2;
+		}
+
+	return CheckCarry(c);
 	}