updated changelog to previous commit

git-svn-id: svn://ttmath.org/publicrep/ttmath/trunk@164 e52654a7-88a9-db11-a3e9-0013d4bc506e

CHANGELOG

@@ -1,3 +1,61 @@
+Version 0.8.5 (2009.06.16):
+    * fixed:   Big::Mod(x) didn't correctly return a carry
+               and the result was sometimes very big (even greater than x)
+    * fixed:   global function Mod(x) didn't set an ErrorCode object
+    * fixed:   global function Round() didn't test a carry
+               now it sets ErrorCode object
+    * changed: function Sin(x) to Sin(x, ErrorCode * err=0)
+               when x was very big the function returns zero
+               now it sets ErrorCode object to err_overflow
+               and the result has a NaN flag set
+               the same is to Cos() function
+    * changed: PrepareSin(x) is using Big::Mod() now when reducing 2PI period
+               should be a little accurate especially on a very big 'x'
+    * changed: uint Mul(const UInt<value_size> & ss2, uint algorithm = 100)
+               void MulBig(const UInt<value_size> & ss2, UInt<value_size*2> & result, uint algorithm = 100)
+               those methods by default use MulFastest() and MulFastestBig()
+    * changed: changed a little Mul2Big() to cooperate with Mul3Big()
+    * added:   uint UInt::Mul3(const UInt<value_size> & ss2)
+               void UInt::Mul3Big(const UInt<value_size> & ss2, UInt<value_size*2> & result)
+               a new multiplication algorithm: Karatsuba multiplication,
+               on a vector UInt<100> with all items different from zero this algorithm is faster
+               about 3 times than Mul2Big(), and on a vector UInt<1000> with all items different from
+               zero this algorithm is faster more than 5 times than Mul2Big()
+               (measured on 32bit platform with GCC 4.3.3 with -O3 and -DTTMATH_RELEASE)
+    * added:   uint MulFastest(const UInt<value_size> & ss2)
+               void MulFastestBig(const UInt<value_size> & ss2, UInt<value_size*2> & result)
+               those methods are trying to select the fastest multiplication algorithm
+    * added:   uint AddVector(const uint * ss1, const uint * ss2, uint ss1_size, uint ss2_size, uint * result)
+               uint SubVector(const uint * ss1, const uint * ss2, uint ss1_size, uint ss2_size, uint * result)
+               three forms: asm x86, asm x86_64, no-asm
+               those methods are used by the Karatsuba multiplication algorithm
+    * added:   to Big<> class: support for NaN flag (Not a Number)
+               bool Big::IsNan() - returns true if the NaN flag is set
+               void Big::SetNan() - sets the NaN flag
+               The NaN flag is set by default after creating an object:
+                 Big<1, 2> a;    // NaN is set (it means the object has not a valid number)
+                 std::cout << a; // cout gives "NaN"
+                 a = 123;        // now NaN is not set
+                 std::cout << a; // cout gives "123"
+               The NaN is set if there was a carry during calculations
+                 a.Mul(very_big_value); // a will have a NaN set
+               The NaN is set if an argument is NaN too
+                 b.SetNan();
+                 a.Add(b);  // a will have NaN because b has NaN too
+               If you try to do something on a NaN object, the result is a NaN too
+                 a.SetNan();
+                 a.Add(2);  // a is still a NaN 
+               The NaN is set if you use incorrect arguments
+                 a.Ln(-10); // a will have the NaN flag
+               The only way to clear the NaN flag is to assign a correct value or other correct object,
+               supposing 'a' has NaN flag, to remove the flag you can either:
+                 a = 10;
+                 a.FromInt(30);
+                 a.SetOne(); 
+                 a.FromBig(other_object_without_nan);
+                 etc.
+               
+               
 Version 0.8.4 (2009.05.08):
     * fixed:   UInt::DivInt() didn't check whether the divisor is zero
                there was a hardware interruption when the divisor was zero

constgen/Makefile

@@ -1,6 +1,6 @@
 o      = main.o
 CC     = g++
-CFLAGS = -s -O2 -DCONSTANTSGENERATOR
+CFLAGS = -s -O2 -DTTMATH_CONSTANTSGENERATOR
 name   = gen
 
 

constgen/main.cpp

@@ -91,7 +91,7 @@ void CalcE()
 	ttmath::Big<1,400> e;
 	ttmath::uint steps;
 
-	// macro CONSTANTSGENERATOR has to be defined	
+	// macro TTMATH_CONSTANTSGENERATOR has to be defined	
 	e.ExpSurrounding0(1, &steps);
 	std::cout << "---------------- e ----------------" << std::endl;
 	e.mantissa.PrintTable(std::cout);
@@ -105,7 +105,7 @@ void CalcLn(int x)
 	ttmath::Big<1,400> ln;
 	ttmath::uint steps;
 
-	// macro CONSTANTSGENERATOR has to be defined	
+	// macro TTMATH_CONSTANTSGENERATOR has to be defined	
 	ln.LnSurrounding1(x, &steps);
 	std::cout << "---------------- ln(" << x << ") ----------------" << std::endl;
 	ln.mantissa.PrintTable(std::cout);

ttmath/ttmath.h

@@ -45,6 +45,12 @@
     \brief Mathematics functions.
 */
 
+#ifdef _MSC_VER
+//warning C4127: conditional expression is constant
+#pragma warning( disable: 4127 )
+#endif
+
+
 #include "ttmathbig.h"
 #include "ttmathobjects.h"
 
@@ -94,10 +100,21 @@ namespace ttmath
 			 -2.7 = -3
 	*/
 	template<class ValueType>
-	ValueType Round(const ValueType & x)
+	ValueType Round(const ValueType & x, ErrorCode * err = 0)
 	{
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return x; // NaN
+		}
+
 		ValueType result( x );
-		result.Round();
+		uint c = result.Round();
+
+		if( err )
+			*err = c ? err_overflow : err_ok;
 
 	return result;
 	}
@@ -118,6 +135,14 @@ namespace ttmath
 	template<class ValueType>
 	ValueType Ceil(const ValueType & x, ErrorCode * err = 0)
 	{
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return x; // NaN
+		}
+
 		ValueType result(x);
 		uint c = 0;
 
@@ -157,6 +182,14 @@ namespace ttmath
 	template<class ValueType>
 	ValueType Floor(const ValueType & x, ErrorCode * err = 0)
 	{
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return x; // NaN
+		}
+
 		ValueType result(x);
 		uint c = 0;
 
@@ -197,8 +230,15 @@ namespace ttmath
 	template<class ValueType>
 	ValueType Ln(const ValueType & x, ErrorCode * err = 0)
 	{
-	ValueType result;
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
 
+		return x; // NaN
+		}
+
+		ValueType result;
 		uint state = result.Ln(x);
 
 		if( err )
@@ -231,8 +271,15 @@ namespace ttmath
 	template<class ValueType>
 	ValueType Log(const ValueType & x, const ValueType & base, ErrorCode * err = 0)
 	{
-	ValueType result;
+		if( x.IsNan() || base.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
 
+		return ValueType(); // default NaN
+		}
+
+		ValueType result;
 		uint state = result.Log(x, base);
 
 		if( err )
@@ -265,8 +312,15 @@ namespace ttmath
 	template<class ValueType>
 	ValueType Exp(const ValueType & x, ErrorCode * err = 0)
 	{
-	ValueType result;
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
 
+		return x; // NaN
+		}
+
+		ValueType result;
 		uint c = result.Exp(x);
 
 		if( err )
@@ -295,7 +349,7 @@ namespace ttmath
 		(you don't have to call this function) 
 	*/
 	template<class ValueType>
-	void PrepareSin(ValueType & x, bool & change_sign)
+	uint PrepareSin(ValueType & x, bool & change_sign)
 	{
 	ValueType temp;
 
@@ -311,12 +365,10 @@ namespace ttmath
 		// we're reducing the period 2*PI
 		// (for big values there'll always be zero)
 		temp.Set2Pi();
-		if( x > temp )
-		{
-			x.Div( temp );
-			x.RemainFraction();
-			x.Mul( temp );
-		}
+		
+		if( x.Mod(temp) )
+			return 1;
+		
 
 		// we're setting 'x' as being in the range of <0, 0.5PI>
 
@@ -337,6 +389,8 @@ namespace ttmath
 			x.Sub( temp );
 			x = temp - x;
 		}
+
+	return 0;
 	}
 
 	
@@ -455,15 +509,38 @@ namespace ttmath
 		this function calculates the Sine
 	*/
 	template<class ValueType>
-	ValueType Sin(ValueType x)
+	ValueType Sin(ValueType x, ErrorCode * err = 0)
 	{
 	using namespace auxiliaryfunctions;
 
-	ValueType one;
+	ValueType one, result;
 	bool change_sign;	
 	
-		PrepareSin( x, change_sign );
-		ValueType result = Sin0pi05( x );
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return result; // NaN is set by default
+		}
+
+		if( err )
+			*err = err_ok;
+
+		if( PrepareSin( x, change_sign ) )
+		{
+			// x is too big, we cannnot reduce the 2*PI period
+			// prior to version 0.8.5 the result was zero
+			
+			// result has NaN flag set by default
+
+			if( err )
+				*err = err_overflow; // maybe another error code? err_improper_argument?
+
+		return result; // NaN is set by default
+		}
+
+		result = Sin0pi05( x );
 	
 		one.SetOne();
 
@@ -488,14 +565,30 @@ namespace ttmath
 		we're using the formula cos(x) = sin(x + PI/2)
 	*/
 	template<class ValueType>
-	ValueType Cos(ValueType x)
+	ValueType Cos(ValueType x, ErrorCode * err = 0)
 	{
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return x; // NaN
+		}
+
 		ValueType pi05;
 		pi05.Set05Pi();
 
-		x.Add( pi05 );
+		uint c = x.Add( pi05 );
 
-	return Sin(x);
+		if( c )
+		{
+			if( err )
+				*err = err_overflow;
+	
+		return ValueType(); // result is undefined (NaN is set by default)
+		}
+
+	return Sin(x, err);
 	}
 	
 
@@ -512,20 +605,22 @@ namespace ttmath
 	template<class ValueType>
 	ValueType Tan(const ValueType & x, ErrorCode * err = 0)
 	{
-		ValueType result = Cos(x);
+		ValueType result = Cos(x, err);
+		
+		if( err && *err != err_ok )
+			return result;
 
 		if( result.IsZero() )
 		{
 			if( err )
 				*err = err_improper_argument;
 
+			result.SetNan();
+
 		return result;
 		}
 
-		if( err )
-			*err = err_ok;
-
-	return Sin(x) / result;
+	return Sin(x, err) / result;
 	}
 
 
@@ -552,20 +647,22 @@ namespace ttmath
 	template<class ValueType>
 	ValueType Cot(const ValueType & x, ErrorCode * err = 0)
 	{
-		ValueType result = Sin(x);
+		ValueType result = Sin(x, err);
+
+		if( err && *err != err_ok )
+			return result;
 
 		if( result.IsZero() )
 		{
 			if( err )
 				*err = err_improper_argument;
 
+			result.SetNan();
+
 		return result;
 		}
 	
-		if( err )
-			*err = err_ok;
-
-	return Cos(x) / result;
+	return Cos(x, err) / result;
 	}
 
 
@@ -746,16 +843,24 @@ namespace ttmath
 	{
 	using namespace auxiliaryfunctions;
 
-		ValueType one;
+		ValueType result, one;
 		one.SetOne();
 		bool change_sign = false;
 
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return result; // NaN is set by default
+		}
+
 		if( x.GreaterWithoutSignThan(one) )
 		{
 			if( err )
 				*err = err_improper_argument;
 
-			return one;
+			return result; // NaN is set by default
 		}
 
 		if( x.IsSign() )
@@ -766,8 +871,6 @@ namespace ttmath
 
 		one.exponent.SubOne(); // =0.5
 
-		ValueType result;
-
 		// asin(-x) = -asin(x)
 		if( x.GreaterWithoutSignThan(one) )
 			result = ASin_1(x);	
@@ -974,6 +1077,9 @@ namespace ttmath
 		one.SetOne();
 		bool change_sign = false;
 
+		if( x.IsNan() )
+			return result; // NaN is set by default
+
 		// if x is negative we're using the formula:
 		// atan(-x) = -atan(x)
 		if( x.IsSign() )
@@ -1054,6 +1160,14 @@ namespace ttmath
 	template<class ValueType>
 	ValueType Sinh(const ValueType & x, ErrorCode * err = 0)
 	{
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return x; // NaN
+		}
+
 		ValueType ex, emx;
 		uint c = 0;
 
@@ -1078,6 +1192,14 @@ namespace ttmath
 	template<class ValueType>
 	ValueType Cosh(const ValueType & x, ErrorCode * err = 0)
 	{
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return x; // NaN
+		}
+
 		ValueType ex, emx;
 		uint c = 0;
 
@@ -1102,6 +1224,14 @@ namespace ttmath
 	template<class ValueType>
 	ValueType Tanh(const ValueType & x, ErrorCode * err = 0)
 	{
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return x; // NaN
+		}
+
 		ValueType ex, emx, nominator, denominator;
 		uint c = 0;
 
@@ -1142,12 +1272,20 @@ namespace ttmath
 	template<class ValueType>
 	ValueType Coth(const ValueType & x, ErrorCode * err = 0)
 	{
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return x; // NaN
+		}
+
 		if( x.IsZero() )
 		{
 			if( err )
 				*err = err_improper_argument;
 
-			return x;
+			return ValueType(); // NaN is set by default
 		}
 
 		ValueType ex, emx, nominator, denominator;
@@ -1199,6 +1337,14 @@ namespace ttmath
 	template<class ValueType>
 	ValueType ASinh(const ValueType & x, ErrorCode * err = 0)
 	{
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return x; // NaN
+		}
+
 		ValueType xx(x), one, result;
 		uint c = 0;
 		one.SetOne();
@@ -1227,6 +1373,14 @@ namespace ttmath
 	template<class ValueType>
 	ValueType ACosh(const ValueType & x, ErrorCode * err = 0)
 	{
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return x; // NaN
+		}
+
 		ValueType xx(x), one, result;
 		uint c = 0;
 		one.SetOne();
@@ -1236,7 +1390,7 @@ namespace ttmath
 			if( err )
 				*err = err_improper_argument;
 
-		return result;
+		return result; // NaN is set by default
 		}
 
 		c += xx.Mul(x);
@@ -1268,6 +1422,14 @@ namespace ttmath
 	template<class ValueType>
 	ValueType ATanh(const ValueType & x, ErrorCode * err = 0)
 	{
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return x; // NaN
+		}
+
 		ValueType nominator(x), denominator, one, result;
 		uint c = 0;
 		one.SetOne();
@@ -1277,7 +1439,7 @@ namespace ttmath
 			if( err )
 				*err = err_improper_argument;
 
-		return result;
+		return result; // NaN is set by default
 		}
 
 		c += nominator.Add(one);
@@ -1313,6 +1475,14 @@ namespace ttmath
 	template<class ValueType>
 	ValueType ACoth(const ValueType & x, ErrorCode * err = 0)
 	{
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return x; // NaN
+		}
+
 		ValueType nominator(x), denominator(x), one, result;
 		uint c = 0;
 		one.SetOne();
@@ -1322,7 +1492,7 @@ namespace ttmath
 			if( err )
 				*err = err_improper_argument;
 
-		return result;
+		return result; // NaN is set by default
 		}
 
 		c += nominator.Add(one);
@@ -1371,6 +1541,14 @@ namespace ttmath
 	ValueType result, temp;
 	uint c = 0;
 
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return result; // NaN is set by default
+		}
+
 		result = x;
 
 		// it is better to make division first and then multiplication
@@ -1399,6 +1577,14 @@ namespace ttmath
 	ValueType result, delimiter;
 	uint c = 0;
 
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return result; // NaN is set by default
+		}
+
 		result = 180;
 		c += result.Mul(x);
 
@@ -1436,12 +1622,12 @@ namespace ttmath
 	ValueType delimiter, multipler;
 	uint c = 0;
 
-		if( m.IsSign() || s.IsSign() )
+		if( d.IsNan() || m.IsNan() || s.IsNan() || m.IsSign() || s.IsSign() )
 		{
 			if( err )
 				*err = err_improper_argument;
 
-			return delimiter;
+		return delimiter ; // NaN is set by default
 		}
 
 		multipler = 60;
@@ -1490,6 +1676,14 @@ namespace ttmath
 	ValueType result, temp;
 	uint c = 0;
 
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return result; // NaN is set by default
+		}
+
 		result = x;
 
 		// it is better to make division first and then multiplication
@@ -1518,6 +1712,14 @@ namespace ttmath
 	ValueType result, delimiter;
 	uint c = 0;
 
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return result; // NaN is set by default
+		}
+
 		result = 200;
 		c += result.Mul(x);
 
@@ -1542,6 +1744,14 @@ namespace ttmath
 	ValueType result, temp;
 	uint c = 0;
 
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return result; // NaN is set by default
+		}
+
 		result = x;
 
 		temp = 200;
@@ -1584,6 +1794,14 @@ namespace ttmath
 	ValueType result, temp;
 	uint c = 0;
 
+		if( x.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return result; // NaN is set by default
+		}
+
 		result = x;
 
 		temp = 180;
@@ -1617,12 +1835,12 @@ namespace ttmath
 	template<class ValueType>
 	ValueType Sqrt(ValueType x, ErrorCode * err = 0)
 	{
-		if( x.IsSign() )
+		if( x.IsNan() || x.IsSign() )
 		{
 			if( err )
 				*err = err_improper_argument;
 
-		return x;
+		return ValueType(); // NaN is set by default
 		}
 
 		if( x.IsZero() )
@@ -1660,6 +1878,8 @@ namespace ttmath
 				if( err )
 					*err = err_improper_argument;
 
+				x.SetNan();
+
 			return true;
 			}
 
@@ -1678,6 +1898,8 @@ namespace ttmath
 					if( err )
 						*err = err_improper_argument;
 
+					x.SetNan();
+
 				return true;
 				}
 		
@@ -1695,7 +1917,7 @@ namespace ttmath
 
 
 		template<class ValueType>
-		bool RootCheckIndexOne(ValueType & x, const ValueType & index, ErrorCode * err)
+		bool RootCheckIndexOne(const ValueType & index, ErrorCode * err)
 		{
 			ValueType one;
 			one.SetOne();
@@ -1727,6 +1949,8 @@ namespace ttmath
 				if( err )
 					*err = err_improper_argument;
 
+				x.SetNan();
+
 			return true;
 			}
 
@@ -1735,11 +1959,12 @@ namespace ttmath
 
 
 		template<class ValueType>
-		bool RootCheckXZero(ValueType & x, const ValueType & index, ErrorCode * err)
+		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 )
@@ -1775,6 +2000,8 @@ namespace ttmath
 					if( err )
 						*err = err_improper_argument;
 
+					x.SetNan();
+
 					return true;
 				}
 			}
@@ -1804,11 +2031,19 @@ namespace ttmath
 	{
 		using namespace auxiliaryfunctions;
 
+		if( x.IsNan() || index.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return ValueType(); // NaN is set by default
+		}
+
 		if( RootCheckIndexSign(x, index, err) ) return x;
 		if( RootCheckIndexZero(x, index, err) ) return x;
-		if( RootCheckIndexOne (x, index, err) ) return x;
+		if( RootCheckIndexOne (   index, err) ) return x;
 		if( RootCheckIndexFrac(x, index, err) ) return x;
-		if( RootCheckXZero(x, index, err) ) return x;
+		if( RootCheckXZero    (x,        err) ) return x;
 
 		// index integer and index!=0
 		// x!=0
@@ -1859,13 +2094,17 @@ namespace ttmath
 
 		while( !carry && multipler<maxvalue  )
 		{
-			if( stop && stop->WasStopSignal() )
+			if( stop && (multipler & 127)==0 ) // it means 'stop && (multipler % 128)==0'
+			{
+				// after each 128 iterations we make a test
+				if( stop->WasStopSignal() )
 				{
 					if( err )
 						*err = err_interrupt;
 
 					return 2;
 				}
+			}
 			
 			++multipler;
 			carry += result.MulUInt(multipler);
@@ -1888,19 +2127,25 @@ namespace ttmath
 
 		one.SetOne();
 		uint carry = 0;
+		uint iter = 1; // only for testing the stop object
 
 		while( !carry && multipler < x )
 		{
-			if( stop && stop->WasStopSignal() )
+			if( stop && (iter & 31)==0 ) // it means 'stop && (iter % 32)==0'
+			{
+				// after each 32 iterations we make a test
+				if( stop->WasStopSignal() )
 				{
 					if( err )
 						*err = err_interrupt;
 
 					return 2;
 				}
+			}
 			
 			carry += multipler.Add(one);
 			carry += result.Mul(multipler);
+			++iter;
 		}
 
 		if( err )
@@ -1927,16 +2172,16 @@ namespace ttmath
 	static History<ValueType> history;
 	ValueType result;
 
-		result.SetOne();
-		
-		if( x.IsSign() )
+		if( x.IsNan() || x.IsSign() )
 		{
 			if( err )
 				*err = err_improper_argument;
 
-		return result;
+		return result; // NaN set by default
 		}
 
+		result.SetOne();
+
 		if( !x.exponent.IsSign() && !x.exponent.IsZero() )
 		{
 			// when x.exponent>0 there's no sense to calculate the formula
@@ -1945,6 +2190,8 @@ namespace ttmath
 			if( err )
 				*err = err_overflow;
 
+			result.SetNan();
+
 		return result;
 		}
 
@@ -1963,8 +2210,11 @@ namespace ttmath
 			status = FactorialMore(x, err, stop, result);
 
 		if( status == 2 )
+		{
 			// the calculation has been interrupted
+			result.SetNan();
 			return result;
+		}
 
 		err_tmp = status==1 ? err_overflow : err_ok;
 		history.Add(x, result, err_tmp);
@@ -2008,14 +2258,25 @@ namespace ttmath
 
 		e.g.
 		mod( 12.6 ;  3) =  0.6   because 12.6  = 3*4 + 0.6
-		mod(-12.6 ;  3) = -0.6
+		mod(-12.6 ;  3) = -0.6   bacause -12.6 = 3*(-4) + (-0.6)
 		mod( 12.6 ; -3) =  0.6
 		mod(-12.6 ; -3) = -0.6
 	*/
 	template<class ValueType>
-	ValueType Mod(ValueType a, const ValueType & b)
+	ValueType Mod(ValueType a, const ValueType & b, ErrorCode * err = 0)
 	{
-		a.Mod(b);
+		if( a.IsNan() || b.IsNan() )
+		{
+			if( err )
+				*err = err_improper_argument;
+
+		return ValueType(); // NaN is set by default
+		}
+
+		uint c = a.Mod(b);
+
+		if( err )
+			*err = c ? err_overflow : err_ok;
 
 	return a;
 	}
@@ -2031,4 +2292,10 @@ namespace ttmath
 */
 #include "ttmathparser.h"
 
+
+#ifdef _MSC_VER
+#pragma warning( default: 4127 )
+//warning C4127: conditional expression is constant
+#endif
+
 #endif

ttmath/ttmathbig.h

@@ -87,14 +87,35 @@ unsigned char info;
 /*!
 	the number of a bit from 'info' which means that a value is with a sign
 	(when the bit is set)
-
-	/at the moment the rest bits from 'info' are not used/
 */
 #define TTMATH_BIG_SIGN 128
 
 
+/*!
+	Not a number
+	if this bit is set that there is no a valid number
+*/
+#define TTMATH_BIG_NAN  64
 
 
+
+
+	/*!
+		this method sets NaN if there was a carry (and returns 1 in such a case)
+
+		c can be 0, 1 or other value different from zero
+	*/
+	uint CheckCarry(uint c)
+	{
+		if( c != 0 )
+		{
+			SetNan();
+			return 1;
+		}
+
+	return 0;
+	}
+
 public:
 
 
@@ -184,6 +205,15 @@ public:
 	}
 
 
+	/*!
+		this method sets NaN flag (Not a Number)
+		when this flag is set that means there is no a valid number
+	*/
+	void SetNan()
+	{
+		info |= TTMATH_BIG_NAN;	
+	}
+
 
 private:
 
@@ -504,6 +534,8 @@ public:
 	{
 		/*
 			we only have to test the mantissa
+			also we don't check the NaN flag
+			(maybe this method should return false if there is NaN flag set?)
 		*/
 		return mantissa.IsZero();
 	}
@@ -511,6 +543,7 @@ public:
 
 	/*!
 		this method returns true when there's the sign set
+		also we don't check the NaN flag
 	*/
 	bool IsSign() const
 	{
@@ -518,6 +551,16 @@ public:
 	}
 
 
+	/*!
+		this method returns true when there is not a valid number
+	*/
+	bool IsNan() const
+	{
+		return (info & TTMATH_BIG_NAN) == TTMATH_BIG_NAN;
+	}
+
+
+
 	/*!
 		this method clears the sign
 		(there'll be an absolute value)
@@ -540,6 +583,10 @@ public:
 	*/
 	void Sgn()
 	{
+		// we have to check the NaN flag, because the next SetOne() method would clear it
+		if( IsNan() )
+			return;
+
 		if( IsSign() )
 		{
 			SetOne();
@@ -572,6 +619,7 @@ public:
 
 	/*!
 		this method changes the sign
+		when there is a value of zero then the sign is not changed
 
 			e.g.
 			-1 -> 1
@@ -579,6 +627,8 @@ public:
 	*/
 	void ChangeSign()
 	{
+		// we don't have to check the NaN flag here
+
 		if( info & TTMATH_BIG_SIGN )
 		{
 			info &= ~TTMATH_BIG_SIGN;
@@ -614,6 +664,9 @@ public:
 
 	uint c = 0;
 
+		if( IsNan() || ss2.IsNan() )
+			return CheckCarry(1);
+
 		exp_offset.Sub( ss2.exponent );
 		exp_offset.Abs();
 
@@ -655,14 +708,19 @@ public:
 			// there shouldn't be a carry here because
 			// (1) (2) guarantee that the mantissa of this
 			// is greater than or equal to the mantissa of the ss2
-			uint c_temp = mantissa.Sub(ss2.mantissa);
+			
+			#ifdef TTMATH_DEBUG
+			// this is to get rid of a warning saying that c_temp is not used
+			uint c_temp = /* mantissa.Sub(ss2.mantissa); */
+			#endif
+			mantissa.Sub(ss2.mantissa);
 
 			TTMATH_ASSERT( c_temp == 0 )
 		}
 
 		c += Standardizing();
 
-	return (c==0)? 0 : 1;
+	return CheckCarry(c);
 	}
 
 
@@ -690,6 +748,9 @@ public:
 	*/
 	uint BitAnd(Big<exp, man> ss2)
 	{
+		if( IsNan() || ss2.IsNan() )
+			return CheckCarry(1);
+
 		if( IsSign() || ss2.IsSign() )
 			return 2;
 
@@ -723,7 +784,7 @@ public:
 
 		c += Standardizing();
 
-	return (c==0)? 0 : 1;
+	return CheckCarry(c);
 	}
 
 
@@ -738,6 +799,9 @@ public:
 	*/
 	uint BitOr(Big<exp, man> ss2)
 	{
+		if( IsNan() || ss2.IsNan() )
+			return CheckCarry(1);
+
 		if( IsSign() || ss2.IsSign() )
 			return 2;
 		
@@ -768,7 +832,7 @@ public:
 
 		c += Standardizing();
 
-	return (c==0)? 0 : 1;
+	return CheckCarry(c);
 	}
 
 
@@ -783,6 +847,9 @@ public:
 	*/
 	uint BitXor(Big<exp, man> ss2)
 	{
+		if( IsNan() || ss2.IsNan() )
+			return CheckCarry(1);
+
 		if( IsSign() || ss2.IsSign() )
 			return 2;
 		
@@ -813,7 +880,7 @@ public:
 
 		c += Standardizing();
 
-	return (c==0)? 0 : 1;
+	return CheckCarry(c);
 	}
 
 
@@ -828,6 +895,9 @@ public:
 	UInt<man+1> man_result;
 	uint i,c = 0;
 
+		if( IsNan() )
+			return 1;
+
 		// man_result = mantissa * ss2.mantissa
 		mantissa.MulInt(ss2, man_result);
 
@@ -856,7 +926,7 @@ public:
 
 		c += Standardizing();
 
-	return (c==0)? 0 : 1;
+	return CheckCarry(c);
 	}
 
 
@@ -867,6 +937,9 @@ public:
 	*/
 	uint MulInt(sint ss2)
 	{
+		if( IsNan() )
+			return 1;
+
 		if( ss2 == 0 )
 		{
 			SetZero();
@@ -906,6 +979,9 @@ public:
 	UInt<man*2> man_result;
 	uint i,c;
 
+		if( IsNan() || ss2.IsNan() )
+			return CheckCarry(1);
+
 		// man_result = mantissa * ss2.mantissa
 		mantissa.MulBig(ss2.mantissa, man_result);
 
@@ -936,7 +1012,7 @@ public:
 
 		c += Standardizing();
 
-	return (c==0)? 0 : 1;
+	return CheckCarry(c);
 	}
 	
 
@@ -952,13 +1028,10 @@ public:
 
 	UInt<man*2> man1;
 	UInt<man*2> man2;
-	uint i,c;
+	uint i,c = 0;
 		
-		if( ss2.IsZero() )
-		{
-			// we don't divide by zero
-			return 1;
-		}
+		if( IsNan() || ss2.IsNan() || ss2.IsZero() )
+			return CheckCarry(1);
 
 		for(i=0 ; i<man ; ++i)
 		{
@@ -976,7 +1049,9 @@ public:
 
 		i = man1.CompensationToLeft();
 
-		c  = exponent.Sub(i);
+		if( i )
+			c += exponent.Sub(i);
+
 		c += exponent.Sub(ss2.exponent);
 		
 		for(i=0 ; i<man ; ++i)
@@ -989,7 +1064,7 @@ public:
 
 		c += Standardizing();
 
-	return (c==0)? 0 : 1;
+	return CheckCarry(c);
 	}
 
 
@@ -998,12 +1073,12 @@ public:
 
 		e.g.
 		 12.6 mod  3 =  0.6   because  12.6 = 3*4 + 0.6
-		-12.6 mod  3 = -0.6
+		-12.6 mod  3 = -0.6   bacause -12.6 = 3*(-4) + (-0.6)
 		 12.6 mod -3 =  0.6
 		-12.6 mod -3 = -0.6
 
 		it means:
-		in other words: this(old) = ss2 * q + this(new -- result)
+		in other words: this(old) = ss2 * q + this(new)
 	*/
 	uint Mod(const Big<exp, man> & ss2)
 	{
@@ -1011,16 +1086,26 @@ public:
 
 	uint c = 0;
 
+		if( IsNan() || ss2.IsNan() )
+			return CheckCarry(1);
+
+		if( !SmallerWithoutSignThan(ss2) )
+		{
 			Big<exp, man> temp(*this);
 
-		c += temp.Div(ss2);
+			c = temp.Div(ss2);
 			temp.SkipFraction();
 			c += temp.Mul(ss2);
 			c += Sub(temp);
 
-	return (c==0)? 0 : 1;
+			if( !SmallerWithoutSignThan( ss2 ) )
+				c += 1;
 		}
 
+	return CheckCarry(c);
+	}
+
+
 
 
 	/*!
@@ -1037,30 +1122,35 @@ public:
 	template<uint pow_size>
 	uint Pow(UInt<pow_size> pow)
 	{
+		if( IsNan() )
+			return 1;
+
 		if( pow.IsZero() && IsZero() )
+		{
 			// we don't define zero^zero
+			SetNan();
 			return 2;
+		}
 
 		Big<exp, man> start(*this), start_temp;
 		Big<exp, man> result;
 		result.SetOne();
+		uint c = 0;
 
-		while( !pow.IsZero() )
+		while( !c && !pow.IsZero() )
 		{
 			if( pow.table[0] & 1 )
-				if( result.Mul(start) )
-					return 1;
+				c += result.Mul(start);
 
 			start_temp = start;
-			if( start.Mul(start_temp) )
-				return 1;
+			c += start.Mul(start_temp);
 
 			pow.Rcr(1);
 		}
 
 		*this = result;
 
-	return 0;
+	return CheckCarry(c);
 	}
 
 
@@ -1076,27 +1166,29 @@ public:
 	template<uint pow_size>
 	uint Pow(Int<pow_size> pow)
 	{
+		if( IsNan() )
+			return 1;
+
 		if( !pow.IsSign() )
 			return Pow( UInt<pow_size>(pow) );
 
 		if( IsZero() )
+		{
 			// if 'p' is negative then
 			// 'this' must be different from zero
+			SetNan();
 			return 2;
+		}
 
-		if( pow.ChangeSign() )
-			return 1;
+		uint c = pow.ChangeSign();
 
 		Big<exp, man> t(*this);
-		uint c_temp = t.Pow( UInt<pow_size>(pow) );
-		if( c_temp > 0 )
-			return c_temp;
+		c += t.Pow( UInt<pow_size>(pow) ); // here can only be a carry (return:1)
 
 		SetOne();
-		if( Div(t) )
-			return 1;
+		c += Div(t);
 
-	return 0;
+	return CheckCarry(c);
 	}
 
 
@@ -1133,8 +1225,14 @@ public:
 	*/
 	uint PowUInt(Big<exp, man> pow)
 	{
+		if( IsNan() || pow.IsNan() )
+			return CheckCarry(1);
+
 		if( pow.IsZero() && IsZero() )
+		{
+			SetNan();
 			return 2;
+		}
 
 		if( pow.IsSign() )
 			pow.Abs();
@@ -1143,27 +1241,26 @@ public:
 		Big<exp, man> result;
 		Big<exp, man> one;
 		Int<exp> e_one;
+		uint c = 0;
 
 		e_one.SetOne();
 		one.SetOne();
 		result = one;
 
-		while( pow >= one )
+		while( !c && pow >= one )
 		{
 			if( pow.Mod2() )
-				if( result.Mul(start) )
-					return 1;
+				c += result.Mul(start);
 
 			start_temp = start;
-			if( start.Mul(start_temp) )
-				return 1;
+			c += start.Mul(start_temp);
 
-			pow.exponent.Sub( e_one );
+			c += pow.exponent.Sub( e_one );
 		}
 
 		*this = result;
 
-	return 0;
+	return CheckCarry(c);
 	}
 
 
@@ -1181,24 +1278,27 @@ public:
 	{
 		TTMATH_REFERENCE_ASSERT( pow )
 	
+		if( IsNan() || pow.IsNan() )
+			return CheckCarry(1);
+
 		if( !pow.IsSign() )
 			return PowUInt(pow);
 
 		if( IsZero() )
+		{
 			// if 'pow' is negative then
 			// 'this' must be different from zero
+			SetNan();
 			return 2;
+		}
 
 		Big<exp, man> temp(*this);
-		uint c_temp = temp.PowUInt(pow);
-		if( c_temp > 0 )
-			return c_temp;
+		uint c = temp.PowUInt(pow); // here can only be a carry (result:1)
 
 		SetOne();
-		if( Div(temp) )
-			return 1;
+		c += Div(temp);
 
-	return 0;
+	return CheckCarry(c);
 	}
 
 
@@ -1216,16 +1316,22 @@ public:
 	{
 		TTMATH_REFERENCE_ASSERT( pow )
 
+		if( IsNan() || pow.IsNan() )
+			return CheckCarry(1);
+
 		Big<exp, man> temp;
 		uint c = temp.Ln(*this);
 
-		if( c!= 0 )
+		if( c != 0 ) // can be 2 from Ln()
+		{
+			SetNan();
 			return c;
+		}
 
 		c += temp.Mul(pow);
 		c += Exp(temp);
 
-	return (c==0)? 0 : 1;
+	return CheckCarry(c);
 	}
 
 
@@ -1243,11 +1349,17 @@ public:
 	{
 		TTMATH_REFERENCE_ASSERT( pow )
 
+		if( IsNan() || pow.IsNan() )
+			return CheckCarry(1);
+
 		if( IsZero() )
 		{
 			// 0^pow will be 0 only for pow>0
 			if( pow.IsSign() || pow.IsZero() )
+			{
+				SetNan();
 				return 2;
+			}
 
 			SetZero();
 
@@ -1269,7 +1381,7 @@ public:
 
 private:
 
-#ifdef CONSTANTSGENERATOR
+#ifdef TTMATH_CONSTANTSGENERATOR
 public:
 #endif
 
@@ -1291,19 +1403,17 @@ public:
 		denominator.SetOne();
 		denominator_i.SetOne();
 
-		// every 'step_test' times we make a test
-		const uint step_test = 5;
 		uint i;
 		old_value = *this;
 
-		// we begin from 1 in order to not testing at the beginning
-	#ifdef CONSTANTSGENERATOR
+		// we begin from 1 in order to not test at the beginning
+	#ifdef TTMATH_CONSTANTSGENERATOR
 		for(i=1 ; true ; ++i)
 	#else
 		for(i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
 	#endif
 		{
-			bool testing = ((i % step_test) == 0);
+			bool testing = ((i & 3) == 0); // it means '(i % 4) == 0'
 
 			next_part = numerator;
 
@@ -1316,13 +1426,15 @@ public:
 			// there shouldn't be a carry here
 			Add( next_part );
 
-			if( testing && old_value==*this )
+			if( testing )
+			{
+				if( old_value == *this )
 					// we've added next few parts of the formula but the result
 					// is still the same then we break the loop
 					break;
 				else
 					old_value = *this;
-
+			}
 
 			// we set the denominator and the numerator for a next part of the formula
 			if( denominator_i.Add(one) )
@@ -1358,6 +1470,9 @@ public:
 	{
 	uint c = 0;
 		
+		if( x.IsNan() )
+			return CheckCarry(1);
+
 		if( x.IsZero() )
 		{
 			SetOne();
@@ -1410,7 +1525,7 @@ public:
 			c += PowUInt(e_);
 		}
 	
-	return (c==0)? 0 : 1;
+	return CheckCarry(c);
 	}
 
 
@@ -1418,7 +1533,7 @@ public:
 
 private:
 
-#ifdef CONSTANTSGENERATOR
+#ifdef TTMATH_CONSTANTSGENERATOR
 public:
 #endif
 
@@ -1454,20 +1569,17 @@ public:
 		SetZero();
 
 		old_value = *this;
-		
-		// every 'step_test' times we make a test
-		const uint step_test = 5;
 		uint i;
 
 
-	#ifdef CONSTANTSGENERATOR
+	#ifdef TTMATH_CONSTANTSGENERATOR
 		for(i=1 ; true ; ++i)
 	#else
-		// we begin from 1 in order to not testing at the beginning
+		// we begin from 1 in order to not test at the beginning
 		for(i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
 	#endif
 		{
-			bool testing = ((i % step_test) == 0);
+			bool testing = ((i & 3) == 0); // it means '(i % 4) == 0'
 
 			next_part = x1;
 
@@ -1480,12 +1592,15 @@ public:
 			// there shouldn't be a carry here
 			Add(next_part);
 
-			if( testing && old_value == *this )
+			if( testing )
+			{
+				if( old_value == *this )
 					// we've added next (step_test) parts of the formula but the result
 					// is still the same then we break the loop
 					break;
 				else
 					old_value = *this;
+			}
 
 			if( x1.Mul(x2) )
 				// if there is a carry here the result we return as good
@@ -1523,15 +1638,21 @@ public:
 
 		return values:
 			0 - ok
-			1 - overflow
+			1 - overflow (carry)
 			2 - incorrect argument (x<=0)
 	*/
 	uint Ln(const Big<exp,man> & x)
 	{
 		TTMATH_REFERENCE_ASSERT( x )
 
+		if( x.IsNan() )
+			return CheckCarry(1);
+
 		if( x.IsSign() || x.IsZero() )
+		{
+			SetNan();
 			return 2;
+		}
 
 		// m will be the value of the mantissa in range <1,2)
 		Big<exp,man> m(x);
@@ -1550,7 +1671,7 @@ public:
 		c += exponent_temp.Mul(ln2);
 		c += Add(exponent_temp);
 
-	return (c==0)? 0 : 1;
+	return CheckCarry(c);
 	}
 
 
@@ -1572,14 +1693,23 @@ public:
 		TTMATH_REFERENCE_ASSERT( base )
 		TTMATH_REFERENCE_ASSERT( x )
 
+		if( x.IsNan() || base.IsNan() )
+			return CheckCarry(1);
+
 		if( x.IsSign() || x.IsZero() )
+		{
+			SetNan();
 			return 2;
+		}
 
 		Big<exp,man> denominator;;
 		denominator.SetOne();
 
 		if( base.IsSign() || base.IsZero() || base==denominator )
+		{
+			SetNan();
 			return 3;
+		}
 		
 		if( x == denominator ) // (this is: if x == 1)
 		{
@@ -1588,14 +1718,14 @@ public:
 			return 0;
 		}
 
-		// another error values we've tested at the start
-		// there can be only a carry
+		// another error values we've tested at the beginning
+		// there can only be a carry
 		uint c = Ln(x);
 
 		c += denominator.Ln(base);
 		c += Div(denominator);
 
-	return (c==0)? 0 : 1;
+	return CheckCarry(c);
 	}
 
 
@@ -1616,9 +1746,15 @@ public:
 	{
 		info = another.info;
 
-		if( exponent.FromInt(another.exponent) )
+		if( IsNan() )
 			return 1;
 
+		if( exponent.FromInt(another.exponent) )
+		{
+			SetNan();
+			return 1;
+		}
+
 		uint man_len_min = (man < another_man)? man : another_man;
 		uint i;
 		uint c = 0;
@@ -1632,8 +1768,8 @@ public:
 
 		// MS Visual Express 2005 reports a warning (in the lines with 'uint man_diff = ...'):
 		// warning C4307: '*' : integral constant overflow
-		// but we're using 'if( man > another_man )' and 'if( man < another_man )' and there'll be no such a situation here
-		#ifndef __GNUC__
+		// but we're using 'if( man > another_man )' and 'if( man < another_man )' and there'll be no such situation here
+		#ifdef _MSC_VER
 		#pragma warning( disable: 4307 )
 		#endif
 
@@ -1649,14 +1785,14 @@ public:
 			c += exponent.AddInt(man_diff, 0);
 		}
 		
-		#ifndef __GNUC__
+		#ifdef _MSC_VER
 		#pragma warning( default: 4307 )
 		#endif
 
 		// mantissa doesn't have to be standardized (either the highest bit is set or all bits are equal zero)
 		CorrectZero();
 
-	return (c == 0 )? 0 : 1;
+	return CheckCarry(c);
 	}
 
 
@@ -1895,9 +2031,10 @@ public:
 			// If E=2047 and F is zero and S is 1, then V=-Infinity
 			// If E=2047 and F is zero and S is 0, then V=Infinity
 
-			// at the moment we do not support NaN, -Infinity and +Infinity
+			// we do not support -Infinity and +Infinity
+			// we assume that there is always NaN 
 
-			SetZero();
+			SetNan();
 		}
 		else
 		if( e > 0 )
@@ -2008,9 +2145,10 @@ public:
 			// If E=2047 and F is zero and S is 1, then V=-Infinity
 			// If E=2047 and F is zero and S is 0, then V=Infinity
 
-			// at the moment we do not support NaN, -Infinity and +Infinity
+			// we do not support -Infinity and +Infinity
+			// we assume that there is always NaN 
 
-			SetZero();
+			SetNan();
 		}
 		else
 		if( e > 0 )
@@ -2100,12 +2238,19 @@ public:
 			return 0;
 		}
 
+		if( IsNan() )
+		{
+			result = ToDouble_SetDouble( false, 2047, 0, false, true);
+
+		return 0;
+		}
+
 		sint e_correction = sint(man*TTMATH_BITS_PER_UINT) - 1;
 
 		if( exponent >= 1024 - e_correction )
 		{
 			// +/- infinity
-			result = ToDouble_SetDouble( IsSign(), 2047, 0, true);
+			result = ToDouble_SetDouble( 0, 2047, 0, true);
 
 		return 1;
 		}
@@ -2140,7 +2285,7 @@ private:
 #ifdef TTMATH_PLATFORM32
 
 	// 32bit platforms
-	double ToDouble_SetDouble(bool is_sign, uint e, sint move, bool infinity = false) const
+	double ToDouble_SetDouble(bool is_sign, uint e, sint move, bool infinity = false, bool nan = false) const
 	{
 		union 
 		{
@@ -2155,6 +2300,12 @@ private:
 
 		temp.u[1] |= (e << 20) & 0x7FF00000u;
 
+		if( nan )
+		{
+			temp.u[0] |= 1;
+			return temp.d;
+		}
+
 		if( infinity )
 			return temp.d;
 
@@ -2177,7 +2328,7 @@ private:
 #else
 
 	// 64bit platforms
-	double ToDouble_SetDouble(bool is_sign, uint e, sint move, bool infinity = false) const
+	double ToDouble_SetDouble(bool is_sign, uint e, sint move, bool infinity = false, bool nan = false) const
 	{
 		union 
 		{
@@ -2192,6 +2343,12 @@ private:
 		                
 		temp.u |= (e << 52) & 0x7FF0000000000000ul;
 
+		if( nan )
+		{
+			temp.u |= 1;
+			return temp.d;
+		}
+
 		if( infinity )
 			return temp.d;
 
@@ -2465,13 +2622,21 @@ public:
 		FromBig(value);
 	}
 
+
 	/*!
 		a default constructor
 
-		warning: we don't set any of the members to zero etc.
+		we don't set any of the members to zero
+		only NaN flag is set
 	*/
 	Big()
 	{
+		info = TTMATH_BIG_NAN;
+
+		/*
+			we're directly setting 'info' (instead of calling SetNan())
+			in order to get rid of a warning saying that 'info' is uninitialized
+		*/
 	}
 
 
@@ -2541,7 +2706,8 @@ public:
 		output:
 			return value:
 			0 - ok and 'result' will be an object of type std::string which holds the value
-			1 - if there was a carry 
+			1 - if there was a carry (shoudn't be in a normal situation - if is that means there
+			    is somewhere an error in the library)
 	*/
 	uint ToString(	std::string & result,
 					uint base                  = 10,
@@ -2552,8 +2718,15 @@ public:
 					char decimal_point = TTMATH_COMMA_CHARACTER_1 ) const
 	{
 		static char error_overflow_msg[] = "overflow";
+		static char error_nan_msg[]      = "NaN";
 		result.erase();
 
+		if( IsNan() )
+		{
+			result = error_nan_msg;
+			return 0;
+		}
+
 		if( base<2 || base>16 )
 		{
 			result = error_overflow_msg;
@@ -2951,7 +3124,7 @@ private:
 			else
 				was_carry = false;
 
-			new_man[i] = UInt<man>::DigitToChar( digit );
+			new_man[i] = static_cast<char>( UInt<man>::DigitToChar(digit) );
 		}
 
 		if( i<0 && was_carry )
@@ -3248,6 +3421,8 @@ public:
 
 		if( base<2 || base>16 )
 		{
+			SetNan();
+
 			if( after_source )
 				*after_source = source;
 
@@ -3277,7 +3452,7 @@ public:
 		if( value_read )
 			*value_read = value_read_temp;
 
-	return (c==0)? 0 : 1;
+	return CheckCarry(c);
 	}
 
 
@@ -3417,9 +3592,9 @@ private:
 
 		it is called when the base is 10 and some digits were read before
 	*/
-	int FromString_ReadScientificIfExists(const char * & source)
+	uint FromString_ReadScientificIfExists(const char * & source)
 	{
-	int c = 0;
+	uint c = 0;
 
 		bool scientific_read = false;
 		const char * before_scientific = source;
@@ -3578,6 +3753,7 @@ public:
 		and returns the result
 
 		(in other words it treats 'this' and 'ss2' as values without a sign)
+		we don't check the NaN flag
 	*/
 	bool SmallerWithoutSignThan(const Big<exp,man> & ss2) const
 	{
@@ -3612,6 +3788,7 @@ public:
 		and returns the result
 
 		(in other words it treats 'this' and 'ss2' as values without a sign)
+		we don't check the NaN flag
 	*/
 	bool GreaterWithoutSignThan(const Big<exp,man> & ss2) const
 	{
@@ -3646,6 +3823,7 @@ public:
 		and returns the result
 
 		(in other words it treats 'this' and 'ss2' as values without a sign)
+		we don't check the NaN flag
 	*/
 	bool EqualWithoutSign(const Big<exp,man> & ss2) const
 	{
@@ -3852,7 +4030,7 @@ public:
 	*/
 	void SkipFraction()
 	{
-		if( IsZero() )
+		if( IsNan() || IsZero() )
 			return;
 
 		if( !exponent.IsSign() )
@@ -3886,7 +4064,7 @@ public:
 	*/
 	void RemainFraction()
 	{
-		if( IsZero() )
+		if( IsNan() || IsZero() )
 			return;
 
 		if( !exponent.IsSign() )
@@ -3936,6 +4114,9 @@ public:
 	Big<exp,man> half;
 	uint c;
 
+		if( IsNan() )
+			return 1;
+
 		half.Set05();
 
 		if( IsSign() )
@@ -3951,7 +4132,7 @@ public:
 
 		SkipFraction();
 
-	return c;
+	return CheckCarry(c);
 	}
 
 

ttmath/ttmathint.h

@@ -1316,5 +1316,4 @@ public:
 
 } // namespace
 
-
 #endif

ttmath/ttmathparser.h

@@ -1,5 +1,5 @@
 /*
- * This file is a part of TTMath Mathematical Library
+ * This file is a part of TTMath Bignum Library
  * and is distributed under the (new) BSD licence.
  * Author: Tomasz Sowa <t.sowa@slimaczek.pl>
  */
@@ -708,7 +708,11 @@ void Sin(int sindex, int amount_of_args, ValueType & result)
 	if( amount_of_args != 1 )
 		Error( err_improper_amount_of_arguments );
 
-	result = ttmath::Sin( ConvertAngleToRad(stack[sindex].value) );
+	ErrorCode err;
+	result = ttmath::Sin( ConvertAngleToRad(stack[sindex].value), &err );
+
+	if(err != err_ok)
+		Error( err );
 }
 
 void Cos(int sindex, int amount_of_args, ValueType & result)
@@ -716,7 +720,11 @@ void Cos(int sindex, int amount_of_args, ValueType & result)
 	if( amount_of_args != 1 )
 		Error( err_improper_amount_of_arguments );
 
-	result = ttmath::Cos( ConvertAngleToRad(stack[sindex].value) );
+	ErrorCode err;
+	result = ttmath::Cos( ConvertAngleToRad(stack[sindex].value), &err );
+
+	if(err != err_ok)
+		Error( err );
 }
 
 void Tan(int sindex, int amount_of_args, ValueType & result)
@@ -757,7 +765,10 @@ void Round(int sindex, int amount_of_args, ValueType & result)
 	if( amount_of_args != 1 )
 		Error( err_improper_amount_of_arguments );
 
-	result = ttmath::Round(stack[sindex].value);
+	result = stack[sindex].value;
+
+	if( result.Round() )
+		Error( err_overflow );
 }
 
 
@@ -973,7 +984,7 @@ void Not(int sindex, int amount_of_args, ValueType & result)
 
 void DegToRad(int sindex, int amount_of_args, ValueType & result)
 {
-	ErrorCode err;
+	ErrorCode err = err_ok;
 
 	if( amount_of_args == 1 )
 	{
@@ -1052,7 +1063,7 @@ void RadToGrad(int sindex, int amount_of_args, ValueType & result)
 
 void DegToGrad(int sindex, int amount_of_args, ValueType & result)
 {
-	ErrorCode err;
+	ErrorCode err = err_ok;
 
 	if( amount_of_args == 1 )
 	{
@@ -1555,7 +1566,7 @@ int character;
 
 	do
 	{
-		result   += character;
+		result   += static_cast<char>( character );
 		character = * ++pstring;
 	}
 	while(	(character>='a' && character<='z') ||

ttmath/ttmathtypes.h

@@ -64,7 +64,7 @@
 */
 #define TTMATH_MAJOR_VER		0
 #define TTMATH_MINOR_VER		8
-#define TTMATH_REVISION_VER		4
+#define TTMATH_REVISION_VER		5
 #define TTMATH_PRERELEASE_VER	0
 
 
@@ -237,6 +237,17 @@ namespace ttmath
 
 
 
+/*!
+	this is a limit when calculating Karatsuba multiplication
+	if the size of a vector is smaller than TTMATH_USE_KARATSUBA_MULTIPLICATION_FROM_SIZE
+	the Karatsuba algorithm will use standard schoolbook multiplication
+*/
+#ifdef __GNUC__
+#define TTMATH_USE_KARATSUBA_MULTIPLICATION_FROM_SIZE 3
+#else
+#define TTMATH_USE_KARATSUBA_MULTIPLICATION_FROM_SIZE 5
+#endif
+
 
 namespace ttmath
 {
@@ -349,7 +360,7 @@ namespace ttmath
 			foo.Add(foo);
 		but there are only few methods which can do that
 	*/
-	class ReferenceError : public std::logic_error, ExceptionInfo
+	class ReferenceError : public std::logic_error, public ExceptionInfo
 	{
 	public:
 
@@ -381,7 +392,7 @@ namespace ttmath
 		the name and the line of a file where the macro TTMATH_ASSERT
 		was used)
 	*/
-	class RuntimeError : public std::runtime_error, ExceptionInfo
+	class RuntimeError : public std::runtime_error, public ExceptionInfo
 	{
 	public:
 

ttmath/ttmathuint.h

@@ -40,6 +40,7 @@
 #ifndef headerfilettmathuint
 #define headerfilettmathuint
 
+
 /*!
 	\file ttmathuint.h
     \brief template class UInt<uint>
@@ -48,6 +49,7 @@
 #include <iostream>
 #include <iomanip>
 
+
 #include "ttmathtypes.h"
 
 
@@ -840,8 +842,10 @@ public:
 
 	/*!
 		the multiplication 'this' = 'this' * ss2
+
+		algorithm: 100 - means automatically choose the fastest algorithm
 	*/
-	uint Mul(const UInt<value_size> & ss2, uint algorithm = 2)
+	uint Mul(const UInt<value_size> & ss2, uint algorithm = 100)
 	{
 		switch( algorithm )
 		{
@@ -849,8 +853,14 @@ public:
 			return Mul1(ss2);
 
 		case 2:
-		default:
 			return Mul2(ss2);
+
+		case 3:
+			return Mul3(ss2);
+
+		case 100:
+		default:
+			return MulFastest(ss2);
 		}
 	}
 
@@ -860,10 +870,12 @@ public:
 
 		since the 'result' is twice bigger than 'this' and 'ss2' 
 		this method never returns a carry
+
+		algorithm: 100 - means automatically choose the fastest algorithm
 	*/
 	void MulBig(const UInt<value_size> & ss2,
 				UInt<value_size*2> & result, 
-				uint algorithm = 2)
+				uint algorithm = 100)
 	{
 		switch( algorithm )
 		{
@@ -871,8 +883,14 @@ public:
 			return Mul1Big(ss2, result);
 
 		case 2:
-		default:
 			return Mul2Big(ss2, result);
+
+		case 3:
+			return Mul3Big(ss2, result);
+
+		case 100:
+		default:
+			return MulFastestBig(ss2, result);
 		}
 	}
 
@@ -964,7 +982,7 @@ public:
 	uint Mul2(const UInt<value_size> & ss2)
 	{
 	UInt<value_size*2> result;
-	uint i;
+	uint i, c = 0;
 
 		Mul2Big(ss2, result);
 	
@@ -975,11 +993,14 @@ public:
 		// testing carry
 		for( ; i<value_size*2 ; ++i)
 			if( result.table[i] != 0 )
-				return 1;
+			{
+				c = 1;
+				break;
+			}
 
 		TTMATH_LOG("UInt::Mul2")
 
-	return 0;
+	return c;
 	}
 
 
@@ -991,44 +1012,385 @@ public:
 	*/
 	void Mul2Big(const UInt<value_size> & ss2, UInt<value_size*2> & result)
 	{
-	uint r2,r1;
-	uint x1size=value_size, x2size=value_size;
+		Mul2Big2<value_size>(table, ss2.table, result);
+
+		TTMATH_LOG("UInt::Mul2Big")
+	}
+
+
+private:
+
+	/*!
+		an auxiliary method for calculating the multiplication 
+
+		arguments we're taking as pointers (this is to improve the Mul3Big2()- avoiding
+		unnecessary copying objects), the result should be taken as a pointer too,
+		but at the moment there is no method AddTwoInts() which can operate on pointers
+	*/
+	template<uint ss_size>
+	void Mul2Big2(const uint * ss1, const uint * ss2, UInt<ss_size*2> & result)
+	{
+	uint x1size  = ss_size, x2size  = ss_size;
 	uint x1start = 0,       x2start = 0;
 
+		if( ss_size > 2 )
+		{	
+			// if the ss_size is smaller than or equal to 2
+			// there is no sense to set x1size (and others) to another values
+
+			for(x1size=ss_size ; x1size>0 && ss1[x1size-1]==0 ; --x1size);
+			for(x2size=ss_size ; x2size>0 && ss2[x2size-1]==0 ; --x2size);
+
+			for(x1start=0 ; x1start<x1size && ss1[x1start]==0 ; ++x1start);
+			for(x2start=0 ; x2start<x2size && ss2[x2start]==0 ; ++x2start);
+		}
+
+		Mul2Big3<ss_size>(ss1, ss2, result, x1start, x1size, x2start, x2size);
+	}
+
+
+
+	/*!
+		an auxiliary method for calculating the multiplication 
+	*/
+	template<uint ss_size>
+	void Mul2Big3(const uint * ss1, const uint * ss2, UInt<ss_size*2> & result, uint x1start, uint x1size, uint x2start, uint x2size)
+	{
+	uint r2, r1;
+
 		result.SetZero();
 
-		if( value_size > 2 )
+		if( x1size==0 || x2size==0 )
+			return;
+
+		for(uint x1=x1start ; x1<x1size ; ++x1)
 		{
-			// if the value_size is smaller than or equal to 2
-			// there is no sense to set x1size (and others) to another values
+			for(uint x2=x2start ; x2<x2size ; ++x2)
+			{
+				MulTwoWords(ss1[x1], ss2[x2], &r2, &r1);
+				result.AddTwoInts(r2, r1, x2+x1);
+				// here will never be a carry
+			}
+		}
+	}
+
+
+public:
+
+
+	/*!
+		multiplication: this = this * ss2
+
+		This is Karatsuba Multiplication algorithm, we're using it when value_size is greater than
+		or equal to TTMATH_USE_KARATSUBA_MULTIPLICATION_FROM_SIZE macro (defined in ttmathuint.h).
+		If value_size is smaller then we're using Mul2Big() instead.
+
+		Karatsuba multiplication:
+		Assume we have:
+			this = x = x1*B^m + x0
+			ss2  = y = y1*B^m + y0
+		where x0 and y0 are less than B^m
+		the product from multiplication we can show as:
+	    x*y = (x1*B^m + x0)(y1*B^m + y0) = z2*B^(2m) + z1*B^m + z0
+		where
+		    z2 = x1*y1
+			z1 = x1*y0 + x0*y1
+			z0 = x0*y0 
+		this is standard schoolbook algorithm with O(n^2), Karatsuba observed that z1 can be given in other form:
+			z1 = (x1 + x0)*(y1 + y0) - z2 - z0    / z1 = (x1*y1 + x1*y0 + x0*y1 + x0*y0) - x1*y1 - x0*y0 = x1*y0 + x0*y1 /
+		and to calculate the multiplication we need only three multiplications (with some additions and subtractions)			
+
+		Our objects 'this' and 'ss2' we divide into two parts and by using recurrence we calculate the multiplication.
+		Karatsuba multiplication has O( n^(ln(3)/ln(2)) )
+	*/
+	uint Mul3(const UInt<value_size> & ss2)
+	{
+	UInt<value_size*2> result;
+	uint i, c = 0;
+
+		Mul3Big(ss2, result);
+	
+		// copying result
+		for(i=0 ; i<value_size ; ++i)
+			table[i] = result.table[i];
+
+		// testing carry
+		for( ; i<value_size*2 ; ++i)
+			if( result.table[i] != 0 )
+			{
+				c = 1;
+				break;
+			}
+
+		TTMATH_LOG("UInt::Mul3")
+
+	return c;
+	}
+
+
+
+	/*!
+		multiplication: result = this * ss2
+
+		result is twice bigger than this and ss2,
+		this method never returns carry,
+		(Karatsuba multiplication)
+	*/
+	void Mul3Big(const UInt<value_size> & ss2, UInt<value_size*2> & result)
+	{
+		Mul3Big2<value_size>(table, ss2.table, result.table);
+
+		TTMATH_LOG("UInt::Mul3Big")
+	}
+
+
+
+private:
+
+	/*!
+		an auxiliary method for calculating the Karatsuba multiplication
+
+		result_size is equal ss_size*2
+	*/
+	template<uint ss_size>
+	void Mul3Big2(const uint * ss1, const uint * ss2, uint * result)
+	{
+	const uint * x1, * x0, * y1, * y0;
+
+
+		if( ss_size>1 && ss_size<TTMATH_USE_KARATSUBA_MULTIPLICATION_FROM_SIZE )
+		{
+			UInt<ss_size*2> res;
+			Mul2Big2<ss_size>(ss1, ss2, res);
+
+			for(uint i=0 ; i<ss_size*2 ; ++i)
+				result[i] = res.table[i];
+
+		return;
+		}
+		else
+		if( ss_size == 1 )
+		{
+			return MulTwoWords(*ss1, *ss2, &result[1], &result[0]);
+		}
+
+
+		if( (ss_size & 1) == 1 )
+		{
+			// ss_size is odd
+			x0 = ss1;
+			y0 = ss2;
+			x1 = ss1 + ss_size / 2 + 1;
+			y1 = ss2 + ss_size / 2 + 1;
+
+			// the second vectors (x1 and y1) are smaller about one from the first ones (x0 and y0)
+			Mul3Big3<ss_size/2 + 1, ss_size/2, ss_size*2>(x1, x0, y1, y0, result);
+		}
+		else
+		{
+			// ss_size is even
+			x0 = ss1;
+			y0 = ss2;
+			x1 = ss1 + ss_size / 2;
+			y1 = ss2 + ss_size / 2;
+			
+			// all four vectors (x0 x1 y0 y1) are equal in size
+			Mul3Big3<ss_size/2, ss_size/2, ss_size*2>(x1, x0, y1, y0, result);
+		}
+	}
+
+
+
+#ifdef _MSC_VER
+#pragma warning (disable : 4717)
+//warning C4717: recursive on all control paths, function will cause runtime stack overflow
+//we have the stop point in Mul3Big2() method
+#endif
+
+
+	/*!
+		an auxiliary method for calculating the Karatsuba multiplication
+
+			x = x1*B^m + x0
+			y = y1*B^m + y0
+
+			first_size  - is the size of vectors: x0 and y0
+			second_size - is the size of vectors: x1 and y1 (can be either equal first_size or smaller about one from first_size)
+
+			x*y = (x1*B^m + x0)(y1*B^m + y0) = z2*B^(2m) + z1*B^m + z0
+		      where
+			   z0 = x0*y0 
+			   z2 = x1*y1
+			   z1 = (x1 + x0)*(y1 + y0) - z2 - z0
+	*/
+	template<uint first_size, uint second_size, uint result_size>
+	void Mul3Big3(const uint * x1, const uint * x0, const uint * y1, const uint * y0, uint * result)
+	{
+	uint i, c, xc, yc;
+
+		UInt<first_size>   temp, temp2;
+		UInt<first_size*3> z1;
+
+		// z0 and z2 we store directly in the result (we don't use any temporary variables)
+		Mul3Big2<first_size>(x0, y0, result);                  // z0
+		Mul3Big2<second_size>(x1, y1, result+first_size*2);    // z2
+
+		// now we calculate z1
+		// temp  = (x0 + x1)
+		// temp2 = (y0 + y1)
+		// we're using temp and temp2 with UInt<first_size>, although there can be a carry but 
+		// we simple remember it in xc and yc (xc and yc can be either 0 or 1),
+		// and (x0 + x1)*(y0 + y1) we calculate in this way (schoolbook algorithm):
+		// 
+		//                 xc     |     temp
+		//                 yc     |     temp2
+		//               --------------------
+		//               (temp    *   temp2)
+		//               xc*temp2 |
+		//               yc*temp  |
+		//       xc*yc |                     
+		//       ----------     z1     --------
+		//
+		// and the result is never larger in size than 3*first_size
+
+		xc = AddVector(x0, x1, first_size, second_size, temp.table);
+		yc = AddVector(y0, y1, first_size, second_size, temp2.table);
+
+		Mul3Big2<first_size>(temp.table, temp2.table, z1.table);
+
+		// clearing the rest of z1
+		for(i=first_size*2 ; i<first_size*3 ; ++i)
+			z1.table[i] = 0;
+
+		
+		if( xc )
+		{
+			c = AddVector(z1.table+first_size, temp2.table, first_size*3-first_size, first_size, z1.table+first_size);
+			TTMATH_ASSERT( c==0 )
+		}
+
+		if( yc )
+		{
+			c = AddVector(z1.table+first_size, temp.table, first_size*3-first_size, first_size, z1.table+first_size);
+			TTMATH_ASSERT( c==0 )
+		}
+
+
+		if( xc && yc )
+		{
+			for( i=first_size*2 ; i<first_size*3 ; ++i )
+				if( ++z1.table[i] != 0 )
+					break;  // break if there was no carry 
+		}
+
+		// z1 = z1 - z2
+		c = SubVector(z1.table, result+first_size*2, first_size*3, second_size*2, z1.table);
+		TTMATH_ASSERT(c==0)
+
+		// z1 = z1 - z0
+		c = SubVector(z1.table, result, first_size*3, first_size*2, z1.table);
+		TTMATH_ASSERT(c==0)
+
+		// here we've calculated the z1
+		// now we're adding it to the result
+
+		if( first_size > second_size )
+		{
+			uint z1_size = result_size - first_size;
+			TTMATH_ASSERT( z1_size <= first_size*3 )
+
+			for(i=z1_size ; i<first_size*3 ; ++i)
+				TTMATH_ASSERT( z1.table[i] == 0 )
+				;
+			
+			c = AddVector(result+first_size, z1.table, result_size-first_size, z1_size, result+first_size);
+			TTMATH_ASSERT(c==0)
+		}
+		else
+		{
+			c = AddVector(result+first_size, z1.table, result_size-first_size, first_size*3, result+first_size);
+			TTMATH_ASSERT(c==0)
+		}
+	}
+
+
+#ifdef _MSC_VER
+#pragma warning (default : 4717)
+#endif
+
+
+public:
+
+
+	/*!
+		multiplication this = this * ss2
+	*/
+	uint MulFastest(const UInt<value_size> & ss2)
+	{
+	UInt<value_size*2> result;
+	uint i, c = 0;
+
+		MulFastestBig(ss2, result);
+	
+		// copying result
+		for(i=0 ; i<value_size ; ++i)
+			table[i] = result.table[i];
+
+		// testing carry
+		for( ; i<value_size*2 ; ++i)
+			if( result.table[i] != 0 )
+			{
+				c = 1;
+				break;
+			}
+
+		TTMATH_LOG("UInt::MulFastest")
+
+	return c;
+	}
+
+
+	/*!
+		multiplication result = this * ss2
+
+		this method is trying to select the fastest algorithm
+		(in the future this method can be improved)
+	*/
+	void MulFastestBig(const UInt<value_size> & ss2, UInt<value_size*2> & result)
+	{
+		if( value_size < TTMATH_USE_KARATSUBA_MULTIPLICATION_FROM_SIZE )
+			return Mul2Big(ss2, result);
+
+		uint x1size  = value_size, x2size  = value_size;
+		uint x1start = 0,          x2start = 0;
 
 		for(x1size=value_size ; x1size>0 && table[x1size-1]==0 ; --x1size);
 		for(x2size=value_size ; x2size>0 && ss2.table[x2size-1]==0 ; --x2size);
 
 		if( x1size==0 || x2size==0 )
 		{
-				TTMATH_LOG("UInt::Mul2Big")
+			// either 'this' or 'ss2' is equal zero - the result is zero too
+			result.SetZero();
 			return;
 		}
 
 		for(x1start=0 ; x1start<x1size && table[x1start]==0 ; ++x1start);
 		for(x2start=0 ; x2start<x2size && ss2.table[x2start]==0 ; ++x2start);
-		}
 
-		for(uint x1=x1start ; x1<x1size ; ++x1)
-		{
-			for(uint x2=x2start ; x2<x2size ; ++x2)
-			{
-				MulTwoWords(table[x1], ss2.table[x2], &r2, &r1);
-				result.AddTwoInts(r2,r1,x2+x1);
-				// here will never be a carry
-			}
-		}
+		uint distancex1 = x1size - x1start;
+		uint distancex2 = x2size - x2start;
 
-		TTMATH_LOG("UInt::Mul2Big")
-	}
+		if( distancex1 < 3 || distancex2 < 3 )
+			// either 'this' or 'ss2' have only 2 (or 1) item different from zero (side by side)
+			// (this condition in the future can be improved)
+			return Mul2Big3<value_size>(table, ss2.table, result, x1start, x1size, x2start, x2size);
 
 
+		// Karatsuba multiplication
+		Mul3Big(ss2, result);
+
+		TTMATH_LOG("UInt::MulFastestBig")
+	}
 
 
 	/*!
@@ -1453,7 +1815,7 @@ private:
 
 		if( Div2_DivisorGreaterOrEqual(	divisor, remainder,
 										table_id, index,
-										divisor_table_id, divisor_index) )
+										divisor_index) )
 		{
 			TTMATH_LOG("UInt::Div2_FindLeadingBitsAndCheck")
 			return 0;
@@ -1473,7 +1835,7 @@ private:
 	bool Div2_DivisorGreaterOrEqual(	const UInt<value_size> & divisor,
 										UInt<value_size> * remainder, 
 										uint table_id, uint index,
-										uint divisor_table_id, uint divisor_index  )
+										uint divisor_index  )
 	{
 		if( divisor_index > index )
 		{
@@ -2346,7 +2708,7 @@ public:
 		do
 		{
 			temp.DivInt(b, &rem);
-			character = DigitToChar( rem );
+			character = static_cast<char>( DigitToChar(rem) );
 			result.insert(result.begin(), character);
 		}
 		while( !temp.IsZero() );
@@ -2909,12 +3271,13 @@ private:
 	uint Rcr2(uint bits, uint c);
 
 public:
-	
 	uint Add(const UInt<value_size> & ss2, uint c=0);
 	uint AddInt(uint value, uint index = 0);
 	uint AddTwoInts(uint x2, uint x1, uint index);
+	static uint AddVector(const uint * ss1, const uint * ss2, uint ss1_size, uint ss2_size, uint * result);
 	uint Sub(const UInt<value_size> & ss2, uint c=0);
 	uint SubInt(uint value, uint index = 0);
+	static uint SubVector(const uint * ss1, const uint * ss2, uint ss1_size, uint ss2_size, uint * result);
 	static sint FindLeadingBitInWord(uint x);
 	static uint SetBitInWord(uint & value, uint bit);
 	static void MulTwoWords(uint a, uint b, uint * result_high, uint * result_low);
@@ -2922,6 +3285,23 @@ public:
 };
 
 
+
+/*!
+	this specialization is needed in order to not confused the compiler "error: ISO C++ forbids zero-size array"
+	when compiling Mul3Big2() method
+*/
+template<>
+class UInt<0>
+{
+public:
+	uint table[1];
+
+	void Mul2Big(const UInt<0> &, UInt<0> &) { TTMATH_ASSERT(false) };
+	void SetZero() { TTMATH_ASSERT(false) };
+	uint AddTwoInts(uint, uint, uint) { TTMATH_ASSERT(false) return 0; };
+};
+
+
 } //namespace
 
 

ttmath/ttmathuint_noasm.h

@@ -95,7 +95,7 @@ namespace ttmath
 		for(i=0 ; i<value_size ; ++i)
 			c = AddTwoWords(table[i], ss2.table[i], c, &table[i]);
 
-		TTMATH_LOG("UInt_noasm::Add")
+		TTMATH_LOG("UInt::Add")
 	
 	return c;
 	}
@@ -131,7 +131,7 @@ namespace ttmath
 		for(i=index+1 ; i<value_size && c ; ++i)
 			c = AddTwoWords(table[i], 0, c, &table[i]);
 
-		TTMATH_LOG("UInt_noasm::AddInt")
+		TTMATH_LOG("UInt::AddInt")
 	
 	return c;
 	}
@@ -184,13 +184,54 @@ namespace ttmath
 		for(i=index+2 ; i<value_size && c ; ++i)
 			c = AddTwoWords(table[i], 0, c, &table[i]);
 
-		TTMATH_LOG("UInt64::AddTwoInts")
+		TTMATH_LOG("UInt::AddTwoInts")
 	
 	return c;
 	}
 
 
 
+	/*!
+		this static method addes one vector to the other
+		'ss1' is larger in size or equal to 'ss2'
+
+		ss1 points to the first (larger) vector
+		ss2 points to the second vector
+		ss1_size - size of the ss1 (and size of the result too)
+		ss2_size - size of the ss2
+		result - is the result vector (which has size the same as ss1: ss1_size)
+
+		Example:  ss1_size is 5, ss2_size is 3
+		ss1:      ss2:   result (output):
+		  5        1         5+1
+		  4        3         4+3
+		  2        7         2+7
+		  6                  6
+		  9                  9
+	  of course the carry is propagated and will be returned from the last item
+	  (this method is used by the Karatsuba multiplication algorithm)
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::AddVector(const uint * ss1, const uint * ss2, uint ss1_size, uint ss2_size, uint * result)
+	{
+	uint i, c = 0;
+
+		TTMATH_ASSERT( ss1_size >= ss2_size )
+		
+		for(i=0 ; i<ss2_size ; ++i)
+			c = AddTwoWords(ss1[i], ss2[i], c, &result[i]);
+
+		for( ; i<ss1_size ; ++i)
+			c = AddTwoWords(ss1[i], 0, c, &result[i]);
+
+		TTMATH_LOG("UInt::AddVector")
+
+	return c;
+	}
+
+
+
+
 	template<uint value_size>
 	uint UInt<value_size>::SubTwoWords(uint a, uint b, uint carry, uint * result)
 	{
@@ -232,7 +273,7 @@ namespace ttmath
 		for(i=0 ; i<value_size ; ++i)
 			c = SubTwoWords(table[i], ss2.table[i], c, &table[i]);
 
-		TTMATH_LOG("UInt_noasm::Sub")
+		TTMATH_LOG("UInt::Sub")
 
 	return c;
 	}
@@ -270,12 +311,51 @@ namespace ttmath
 		for(i=index+1 ; i<value_size && c ; ++i)
 			c = SubTwoWords(table[i], 0, c, &table[i]);
 
-		TTMATH_LOG("UInt_noasm::SubInt")
+		TTMATH_LOG("UInt::SubInt")
 	
 	return c;
 	}
 
 
+	/*!
+		this static method subtractes one vector from the other
+		'ss1' is larger in size or equal to 'ss2'
+
+		ss1 points to the first (larger) vector
+		ss2 points to the second vector
+		ss1_size - size of the ss1 (and size of the result too)
+		ss2_size - size of the ss2
+		result - is the result vector (which has size the same as ss1: ss1_size)
+
+		Example:  ss1_size is 5, ss2_size is 3
+		ss1:      ss2:   result (output):
+		  5        1         5-1
+		  4        3         4-3
+		  2        7         2-7
+		  6                  6-1  (the borrow from previous item)
+		  9                  9
+		                 return (carry): 0
+	  of course the carry (borrow) is propagated and will be returned from the last item
+	  (this method is used by the Karatsuba multiplication algorithm)
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::SubVector(const uint * ss1, const uint * ss2, uint ss1_size, uint ss2_size, uint * result)
+	{
+	uint i, c = 0;
+
+		TTMATH_ASSERT( ss1_size >= ss2_size )
+		
+		for(i=0 ; i<ss2_size ; ++i)
+			c = SubTwoWords(ss1[i], ss2[i], c, &result[i]);
+
+		for( ; i<ss1_size ; ++i)
+			c = SubTwoWords(ss1[i], 0, c, &result[i]);
+
+		TTMATH_LOG("UInt::SubVector")
+
+	return c;
+	}
+
 
 
 	/*!
@@ -305,7 +385,7 @@ namespace ttmath
 			c        = new_c;
 		}
 
-		TTMATH_LOG("UInt64::Rcl2_one")
+		TTMATH_LOG("UInt::Rcl2_one")
 
 	return c;
 	}
@@ -344,7 +424,7 @@ namespace ttmath
 			c        = new_c;
 		}
 
-		TTMATH_LOG("UInt64::Rcr2_one")
+		TTMATH_LOG("UInt::Rcr2_one")
 
 	return c;
 	}
@@ -421,7 +501,7 @@ namespace ttmath
 			c        = new_c;
 		}
 
-		TTMATH_LOG("UInt64::Rcr2")
+		TTMATH_LOG("UInt::Rcr2")
 
 	return (c & TTMATH_UINT_HIGHEST_BIT) ? 1 : 0;
 	}

ttmath/ttmathuint_x86_64.h

@@ -78,7 +78,6 @@ namespace ttmath
 	uint b = value_size;
 	uint * p1 = table;
 	const uint * p2 = ss2.table;
-	uint dummy, dummy2;
 
 		// we don't have to use TTMATH_REFERENCE_ASSERT here
 		// this algorithm doesn't require it
@@ -88,13 +87,15 @@ namespace ttmath
 		#endif
 
 		#ifdef __GNUC__
+		uint dummy, dummy2;
+
 			/*
 				this part should be compiled with gcc
 			*/
 			__asm__ __volatile__(
 	
 				"xorq %%rdx, %%rdx				\n"
-				"neg %%rax						\n"     // CF=1 if rax!=0 , CF=0 if rax==0
+				"negq %%rax						\n"     // CF=1 if rax!=0 , CF=0 if rax==0
 
 			"1:									\n"
 				"movq (%%rsi,%%rdx,8), %%rax	\n"
@@ -112,7 +113,7 @@ namespace ttmath
 
 		#endif
 
-		TTMATH_LOG("UInt64::Add")
+		TTMATH_LOG("UInt::Add")
 	
 	return c;
 	}
@@ -145,7 +146,6 @@ namespace ttmath
 	uint b = value_size;
 	uint * p1 = table;
 	uint c;
-	uint dummy, dummy2;
 
 		TTMATH_ASSERT( index < value_size )
 
@@ -154,6 +154,7 @@ namespace ttmath
 		#endif
 
 		#ifdef __GNUC__
+		uint dummy, dummy2;
 	
 			__asm__ __volatile__(
 
@@ -173,12 +174,12 @@ namespace ttmath
 				"movzx %%al, %%rdx				\n"
 
 				: "=d" (c),    "=a" (dummy), "=c" (dummy2)
-				: "a" (value), "c" (b), "0" (index), "b" (p1)
+				: "0" (index), "1" (value),  "2" (b), "b" (p1)
 				: "cc", "memory" );
 
 		#endif
 
-		TTMATH_LOG("UInt64::AddInt")
+		TTMATH_LOG("UInt::AddInt")
 	
 	return c;
 	}
@@ -223,7 +224,6 @@ namespace ttmath
 	uint b = value_size;
 	uint * p1 = table;
 	uint c;
-	uint dummy, dummy2;
 
 		TTMATH_ASSERT( index < value_size - 1 )
 
@@ -232,6 +232,8 @@ namespace ttmath
 		#endif
 
 		#ifdef __GNUC__
+		uint dummy, dummy2;
+
 			__asm__ __volatile__(
 			
 				"subq %%rdx, %%rcx 				\n"
@@ -254,12 +256,94 @@ namespace ttmath
 				"movzx %%al, %%rax				\n"
 
 				: "=a" (c), "=c" (dummy), "=d" (dummy2)
-				: "1" (b), "2" (index), "b" (p1), "S" (x1), "0" (x2)
+				: "0" (x2), "1" (b),      "2" (index), "b" (p1), "S" (x1)
 				: "cc", "memory" );
 
 		#endif
 
-		TTMATH_LOG("UInt64::AddTwoInts")
+		TTMATH_LOG("UInt::AddTwoInts")
+
+	return c;
+	}
+
+
+
+	/*!
+		this static method addes one vector to the other
+		'ss1' is larger in size or equal to 'ss2'
+
+		ss1 points to the first (larger) vector
+		ss2 points to the second vector
+		ss1_size - size of the ss1 (and size of the result too)
+		ss2_size - size of the ss2
+		result - is the result vector (which has size the same as ss1: ss1_size)
+
+		Example:  ss1_size is 5, ss2_size is 3
+		ss1:      ss2:   result (output):
+		  5        1         5+1
+		  4        3         4+3
+		  2        7         2+7
+		  6                  6
+		  9                  9
+	  of course the carry is propagated and will be returned from the last item
+	  (this method is used by the Karatsuba multiplication algorithm)
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::AddVector(const uint * ss1, const uint * ss2, uint ss1_size, uint ss2_size, uint * result)
+	{
+		TTMATH_ASSERT( ss1_size >= ss2_size )
+
+		uint rest = ss1_size - ss2_size;
+		uint c;
+
+		#ifndef __GNUC__
+			#error "another compiler than GCC is currently not supported in 64bit mode"
+		#endif
+
+		#ifdef __GNUC__
+		uint dummy1, dummy2, dummy3;	
+			
+			//	this part should be compiled with gcc
+		
+			__asm__ __volatile__(
+				"mov %%rdx, %%r8					\n"
+				"xor %%rdx, %%rdx					\n"   // rdx = 0, cf = 0
+			"1:										\n"
+				"mov (%%rsi,%%rdx,8), %%rax			\n"
+				"adc (%%rbx,%%rdx,8), %%rax			\n"
+				"mov %%rax, (%%rdi,%%rdx,8)			\n"
+
+				"inc %%rdx							\n"
+				"dec %%rcx							\n"
+			"jnz 1b									\n"
+
+				"adc %%rcx, %%rcx					\n"   // rcx has the cf state
+
+				"or %%r8, %%r8						\n"
+				"jz 3f								\n"
+				
+				"xor %%rbx, %%rbx					\n"   // ebx = 0
+				"neg %%rcx							\n"   // setting cf from rcx
+				"mov %%r8, %%rcx					\n"   // rcx=rest and is != 0
+			"2:										\n"
+				"mov (%%rsi, %%rdx, 8), %%rax		\n"
+				"adc %%rbx, %%rax 					\n"
+				"mov %%rax, (%%rdi, %%rdx, 8)		\n"
+
+				"inc %%rdx							\n"
+				"dec %%rcx							\n"
+			"jnz 2b									\n"
+
+				"adc %%rcx, %%rcx					\n"
+			"3:										\n"
+
+				: "=a" (dummy1), "=b" (dummy2), "=c" (c),       "=d" (dummy3)
+				:                "1" (ss2),     "2" (ss2_size), "3" (rest),   "S" (ss1),  "D" (result)
+				: "%r8", "cc", "memory" );
+
+		#endif
+
+		TTMATH_LOG("UInt::AddVector")
 
 	return c;
 	}
@@ -284,7 +368,7 @@ namespace ttmath
 	uint b = value_size;
 	uint * p1 = table;
 	const uint * p2 = ss2.table;
-	uint dummy, dummy2;
+	
 
 		// we don't have to use TTMATH_REFERENCE_ASSERT here
 		// this algorithm doesn't require it
@@ -294,10 +378,12 @@ namespace ttmath
 		#endif
 
 		#ifdef __GNUC__
+		uint dummy, dummy2;
+
 			__asm__  __volatile__(
 	
 				"xorq %%rdx, %%rdx				\n"
-				"neg %%rax						\n"     // CF=1 if rax!=0 , CF=0 if rax==0
+				"negq %%rax						\n"     // CF=1 if rax!=0 , CF=0 if rax==0
 
 			"1:									\n"
 				"movq (%%rsi,%%rdx,8), %%rax	\n"
@@ -313,10 +399,9 @@ namespace ttmath
 				: "0" (b),  "1" (c), "b" (p1), "S" (p2)
 				: "cc", "memory" );
 
-
 		#endif
 
-		TTMATH_LOG("UInt64::Sub")
+		TTMATH_LOG("UInt::Sub")
 
 	return c;
 	}
@@ -374,17 +459,105 @@ namespace ttmath
 				"movzx %%al, %%rdx				\n"
 
 				: "=d" (c),    "=a" (dummy), "=c" (dummy2)
-				: "1" (value), "2" (b), "0" (index), "b" (p1)
+				: "0" (index), "1" (value),  "2" (b), "b" (p1)
 				: "cc", "memory" );
 
 		#endif
 
-		TTMATH_LOG("UInt64::SubInt")
+		TTMATH_LOG("UInt::SubInt")
 
 	return c;
 	}
 
 
+
+
+	/*!
+		this static method subtractes one vector from the other
+		'ss1' is larger in size or equal to 'ss2'
+
+		ss1 points to the first (larger) vector
+		ss2 points to the second vector
+		ss1_size - size of the ss1 (and size of the result too)
+		ss2_size - size of the ss2
+		result - is the result vector (which has size the same as ss1: ss1_size)
+
+		Example:  ss1_size is 5, ss2_size is 3
+		ss1:      ss2:   result (output):
+		  5        1         5-1
+		  4        3         4-3
+		  2        7         2-7
+		  6                  6-1  (the borrow from previous item)
+		  9                  9
+		               return (carry): 0
+	  of course the carry (borrow) is propagated and will be returned from the last item
+	  (this method is used by the Karatsuba multiplication algorithm)
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::SubVector(const uint * ss1, const uint * ss2, uint ss1_size, uint ss2_size, uint * result)
+	{
+		TTMATH_ASSERT( ss1_size >= ss2_size )
+
+		uint rest = ss1_size - ss2_size;
+		uint c;
+
+		#ifndef __GNUC__
+			#error "another compiler than GCC is currently not supported in 64bit mode"
+		#endif
+
+		#ifdef __GNUC__
+
+		/*
+			the asm code is nearly the same as in AddVector
+			only two instructions 'adc' are changed to 'sbb'
+		*/
+		uint dummy1, dummy2, dummy3;
+
+			__asm__ __volatile__(
+				"mov %%rdx, %%r8					\n"
+				"xor %%rdx, %%rdx					\n"   // rdx = 0, cf = 0
+			"1:										\n"
+				"mov (%%rsi,%%rdx,8), %%rax			\n"
+				"sbb (%%rbx,%%rdx,8), %%rax			\n"
+				"mov %%rax, (%%rdi,%%rdx,8)			\n"
+
+				"inc %%rdx							\n"
+				"dec %%rcx							\n"
+			"jnz 1b									\n"
+
+				"adc %%rcx, %%rcx					\n"   // rcx has the cf state
+
+				"or %%r8, %%r8						\n"
+				"jz 3f								\n"
+				
+				"xor %%rbx, %%rbx					\n"   // ebx = 0
+				"neg %%rcx							\n"   // setting cf from rcx
+				"mov %%r8, %%rcx					\n"   // rcx=rest and is != 0
+			"2:										\n"
+				"mov (%%rsi, %%rdx, 8), %%rax		\n"
+				"sbb %%rbx, %%rax 					\n"
+				"mov %%rax, (%%rdi, %%rdx, 8)		\n"
+
+				"inc %%rdx							\n"
+				"dec %%rcx							\n"
+			"jnz 2b									\n"
+
+				"adc %%rcx, %%rcx					\n"
+			"3:										\n"
+
+				: "=a" (dummy1), "=b" (dummy2), "=c" (c),       "=d" (dummy3)
+				:                "1" (ss2),     "2" (ss2_size), "3" (rest),   "S" (ss1),  "D" (result)
+				: "%r8", "cc", "memory" );
+
+		#endif
+
+		TTMATH_LOG("UInt::SubVector")
+
+	return c;
+	}
+
+
+
 	/*!
 		this method moves all bits into the left hand side
 		return value <- this <- c
@@ -404,17 +577,19 @@ namespace ttmath
 	{
 	sint b = value_size;
 	uint * p1 = table;
-	uint dummy, dummy2;
+	
 
 		#ifndef __GNUC__
 			#error "another compiler than GCC is currently not supported in 64bit mode"
 		#endif
 
 		#ifdef __GNUC__
+		uint dummy, dummy2;
+
 		__asm__  __volatile__(
 		
 			"xorq %%rdx, %%rdx			\n"   // rdx=0
-			"neg %%rax					\n"   // CF=1 if rax!=0 , CF=0 if rax==0
+			"negq %%rax					\n"   // CF=1 if rax!=0 , CF=0 if rax==0
 
 		"1:								\n"
 			"rclq $1, (%%rbx, %%rdx, 8)	\n"
@@ -426,12 +601,12 @@ namespace ttmath
 			"adcq %%rcx, %%rcx			\n"
 
 			: "=c" (c), "=a" (dummy), "=d" (dummy2)
-			: "1" (c), "0" (b), "b" (p1)
+			: "0" (b),  "1" (c), "b" (p1)
 			: "cc", "memory" );
 	
 		#endif
 
-		TTMATH_LOG("UInt64::Rcl2_one")
+		TTMATH_LOG("UInt::Rcl2_one")
 
 	return c;
 	}
@@ -456,16 +631,18 @@ namespace ttmath
 	{
 	sint b = value_size;
 	uint * p1 = table;
-	uint dummy;
+	
 
 		#ifndef __GNUC__
 			#error "another compiler than GCC is currently not supported in 64bit mode"
 		#endif
 
 		#ifdef __GNUC__
+		uint dummy;
+
 		__asm__  __volatile__(
 
-			"neg %%rax						\n"   // CF=1 if rax!=0 , CF=0 if rax==0
+			"negq %%rax						\n"   // CF=1 if rax!=0 , CF=0 if rax==0
 
 		"1:									\n"
 			"rcrq $1, -8(%%rbx, %%rcx, 8)	\n"
@@ -476,12 +653,12 @@ namespace ttmath
 			"adcq %%rcx, %%rcx				\n"
 
 			: "=c" (c), "=a" (dummy)
-			: "1" (c), "0" (b), "b" (p1)
+			: "0" (b),  "1" (c), "b" (p1)
 			: "cc", "memory" );
 
 		#endif
 
-		TTMATH_LOG("UInt64::Rcr2_one")
+		TTMATH_LOG("UInt::Rcr2_one")
 
 	return c;
 	}
@@ -509,13 +686,15 @@ namespace ttmath
 
 	uint b = value_size;
 	uint * p1 = table;
-	uint dummy, dummy2, dummy3;
+
 
 		#ifndef __GNUC__
 			#error "another compiler than GCC is currently not supported in 64bit mode"
 		#endif
 
 		#ifdef __GNUC__
+		uint dummy, dummy2, dummy3;
+
 		__asm__  __volatile__(
 		
 			"movq %%rcx, %%rsi				\n"
@@ -528,7 +707,6 @@ namespace ttmath
 
 			"xorq %%rdx, %%rdx				\n"
 			"movq %%rdx, %%rsi				\n"
-
 			"orq %%rax, %%rax				\n"
 			"cmovnz %%r8, %%rsi				\n"
 
@@ -553,7 +731,7 @@ namespace ttmath
 
 		#endif
 
-		TTMATH_LOG("UInt64::Rcl2")
+		TTMATH_LOG("UInt::Rcl2")
 
 	return c;
 	}
@@ -602,7 +780,6 @@ namespace ttmath
 			"movq %%rdx, %%rsi				\n"
 			"addq %%rdi, %%rdx				\n"
 			"decq %%rdx						\n"
-
 			"orq %%rax, %%rax				\n"
 			"cmovnz %%R8, %%rsi				\n"
 
@@ -628,7 +805,7 @@ namespace ttmath
 
 		#endif
 
-		TTMATH_LOG("UInt64::Rcr2")
+		TTMATH_LOG("UInt::Rcr2")
 
 	return c;
 	}
@@ -643,22 +820,24 @@ namespace ttmath
 	template<uint value_size>
 	sint UInt<value_size>::FindLeadingBitInWord(uint x)
 	{
-	register sint result;
+	sint result;
+
 
 		#ifndef __GNUC__
 			#error "another compiler than GCC is currently not supported in 64bit mode"
 		#endif
 
 		#ifdef __GNUC__
-			__asm__  __volatile__(
+		uint dummy;
 
-			"bsrq %1, %0		\n"
-			"jnz 1f				\n"
-			"movq $-1, %0		\n"
-			"1:					\n"
+				__asm__ (
 
-			: "=R" (result)
-			: "R" (x)
+				"movq $-1, %1          \n"
+				"bsrq %2, %0           \n"
+				"cmovz %1, %0          \n"
+
+				: "=r" (result), "=&r" (dummy)
+				: "r" (x)
 				: "cc" );
 
 		#endif
@@ -695,10 +874,10 @@ namespace ttmath
 		#endif
 
 		#ifdef __GNUC__
-			__asm__  __volatile__(
+
+			__asm__ (
 
 			"btsq %%rbx, %%rax		\n"
-
 			"setc %%bl				\n"
 			"movzx %%bl, %%rbx		\n"
 			
@@ -742,8 +921,8 @@ namespace ttmath
 		this has no effect in visual studio but it's usefull when
 		using gcc and options like -O
 	*/
-	register uint result1_;
-	register uint result2_;
+	uint result1_;
+	uint result2_;
 
 		#ifndef __GNUC__
 			#error "another compiler than GCC is currently not supported in 64bit mode"
@@ -751,7 +930,7 @@ namespace ttmath
 
 		#ifdef __GNUC__
 
-		__asm__ __volatile__(
+		__asm__ (
 		
 			"mulq %%rdx			\n"
 
@@ -793,8 +972,8 @@ namespace ttmath
 	template<uint value_size>
 	void UInt<value_size>::DivTwoWords(uint a,uint b, uint c, uint * r, uint * rest)
 	{
-		register uint r_;
-		register uint rest_;
+		uint r_;
+		uint rest_;
 		/*
 			these variables have similar meaning like those in
 			the multiplication algorithm MulTwoWords
@@ -808,7 +987,7 @@ namespace ttmath
 
 		#ifdef __GNUC__
 		
-			__asm__ __volatile__(
+			__asm__ (
 
 			"divq %%rcx				\n"