updated changelog to previous commit

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

CHANGELOG

@@ -1,3 +1,90 @@
+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
+               (now the method returns one)
+    * fixed:   the problem with GCC optimization on x86_64
+               sometimes when using -O2 or -O3 GCC doesn't set correctly
+               the stack pointer (actually the stack is used for other things)
+               and you can't use instructions like push/pop in assembler code.
+               All the asm code in x86_64 have been rewritten, now instructions
+               push/pop are not used, other thing which have access to stack 
+               (like "m" (mask) constraints in Rcl2 and Rcr2) have also gone away,
+               now the library works well with -O2 and -O3 and the asm code
+               is a little faster
+    * added:   UInt::PrintLog(const char * msg, std::ostream & output)
+               used (for debugging purposes) by macro TTMATH_LOG(msg)
+               (it is used in nearly all methods in UInt class)
+    * added:   macro TTMATH_DEBUG_LOG: when defined then TTMATH_LOG() 
+               put some debug information (to std::cout)
+    * added:   ttmathuint_x86.h, ttmathuint_x86_64.h, ttmathuint_noasm.h, 
+               all the methods which are using assembler code have been 
+               rewritten to no-asm forms, now we have:
+               1. asm for x86      file: ttmathuint_x86.h
+               2. asm for x86_64   file: ttmathuint_x86_64.h
+               3. no asm           file: ttmathuint_noasm.h
+                  (it's used when macro TTMATH_NOASM is defined)
+               The third form can be used on x86 and x86_64 as well and
+               on other platforms with a little effort.
+                              
+
 Version 0.8.3 (2009.04.06):
     * fixed:   RclMoveAllWords() and RcrMoveAllWords() sometimes didn't return
                the proper carry, (when 'bits' was greater than or equal to 'value_size')

README

@@ -21,6 +21,6 @@ This means only C++ developers can use this library and one thing they have
 to do is to use 'include' directive of the preprocessor. How big the
 values can be is set directly in the source code by the programmer.
 
-Author: Tomasz Sowa <t.sowa@slimaczek.pl>
-Project pages: http://ttmath.slimaczek.pl
+Author: Tomasz Sowa <t.sowa@ttmath.org>
+Project pages: http://www.ttmath.org
 http://sourceforge.net/projects/ttmath

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"
 
@@ -53,6 +59,13 @@
 
 namespace ttmath
 {
+	/*
+ 	 *
+	 *  functions defined here are used only with Big<> types
+	 *
+	 *
+	 */
+
 
 	/*
  	 *
@@ -87,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;
 	}
@@ -111,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;
 
@@ -150,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;
 
@@ -190,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 )
@@ -224,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 )
@@ -258,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 )
@@ -288,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;
 
@@ -304,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>
 
@@ -330,6 +389,8 @@ namespace ttmath
 			x.Sub( temp );
 			x = temp - x;
 		}
+
+	return 0;
 	}
 
 	
@@ -416,7 +477,7 @@ namespace ttmath
 			if( c )
 				// Sin is from <-1,1> and cannot make an overflow
 				// but the carry can be from the Taylor series
-				// (then we only breaks our calculations)
+				// (then we only break our calculations)
 				break;
 
 			if( addition )
@@ -448,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();
 
@@ -481,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);
 	}
 	
 
@@ -505,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;
 	}
 
 
@@ -545,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;
 	}
 
 
@@ -739,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() )
@@ -759,8 +871,6 @@ namespace ttmath
 
 		one.exponent.SubOne(); // =0.5
 
-		ValueType result;
-
 		// asin(-x) = -asin(x)
 		if( x.GreaterWithoutSignThan(one) )
 			result = ASin_1(x);	
@@ -967,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() )
@@ -1047,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;
 
@@ -1071,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;
 
@@ -1095,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;
 
@@ -1135,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;
@@ -1192,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();
@@ -1220,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();
@@ -1229,7 +1390,7 @@ namespace ttmath
 			if( err )
 				*err = err_improper_argument;
 
-		return result;
+		return result; // NaN is set by default
 		}
 
 		c += xx.Mul(x);
@@ -1261,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();
@@ -1270,7 +1439,7 @@ namespace ttmath
 			if( err )
 				*err = err_improper_argument;
 
-		return result;
+		return result; // NaN is set by default
 		}
 
 		c += nominator.Add(one);
@@ -1306,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();
@@ -1315,7 +1492,7 @@ namespace ttmath
 			if( err )
 				*err = err_improper_argument;
 
-		return result;
+		return result; // NaN is set by default
 		}
 
 		c += nominator.Add(one);
@@ -1364,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
@@ -1392,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);
 
@@ -1429,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;
@@ -1483,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
@@ -1511,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);
 
@@ -1535,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;
@@ -1577,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;
@@ -1610,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() )
@@ -1653,6 +1878,8 @@ namespace ttmath
 				if( err )
 					*err = err_improper_argument;
 
+				x.SetNan();
+
 			return true;
 			}
 
@@ -1671,6 +1898,8 @@ namespace ttmath
 					if( err )
 						*err = err_improper_argument;
 
+					x.SetNan();
+
 				return true;
 				}
 		
@@ -1688,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();
@@ -1720,6 +1949,8 @@ namespace ttmath
 				if( err )
 					*err = err_improper_argument;
 
+				x.SetNan();
+
 			return true;
 			}
 
@@ -1728,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 )
@@ -1768,6 +2000,8 @@ namespace ttmath
 					if( err )
 						*err = err_improper_argument;
 
+					x.SetNan();
+
 					return true;
 				}
 			}
@@ -1797,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
@@ -1852,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);
@@ -1881,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 )
@@ -1920,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
@@ -1938,6 +2190,8 @@ namespace ttmath
 			if( err )
 				*err = err_overflow;
 
+			result.SetNan();
+
 		return result;
 		}
 
@@ -1956,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);
@@ -2001,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;
 	}
@@ -2024,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);
 	}
 
 
@@ -1874,7 +2010,7 @@ public:
 		// error but I leave it at the moment as is
 		TTMATH_ASSERT( sizeof(double) == 8 )
 
-		// I am not sure what will be on a plaltform which has 
+		// I am not sure what will be on a platform which has 
 		// a different endianness... but we use this library only
 		// on x86 and amd (intel) 64 bits (as there's a lot of assembler code)
 		union 
@@ -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);
 	}
 
 
@@ -3367,10 +3542,9 @@ private:
 
 		// we don't remove any white characters here
 
-		// this is only to avoid getting a warning about an uninitialized object
-		// gcc 4.1.2 reports: 'old_value.info' may be used uninitialized in this function
+		// this is only to avoid getting a warning about an uninitialized object 'old_value' which GCC reports
 		// (in fact we will initialize it later when the condition 'testing' is fulfilled)
-		old_value.info = 0;
+		old_value.SetZero();
 
 		power.SetOne();
 
@@ -3418,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;
@@ -3579,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
 	{
@@ -3613,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
 	{
@@ -3647,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
 	{
@@ -3853,7 +4030,7 @@ public:
 	*/
 	void SkipFraction()
 	{
-		if( IsZero() )
+		if( IsNan() || IsZero() )
 			return;
 
 		if( !exponent.IsSign() )
@@ -3887,7 +4064,7 @@ public:
 	*/
 	void RemainFraction()
 	{
-		if( IsZero() )
+		if( IsNan() || IsZero() )
 			return;
 
 		if( !exponent.IsSign() )
@@ -3937,6 +4114,9 @@ public:
 	Big<exp,man> half;
 	uint c;
 
+		if( IsNan() )
+			return 1;
+
 		half.Set05();
 
 		if( IsSign() )
@@ -3952,7 +4132,7 @@ public:
 
 		SkipFraction();
 
-	return c;
+	return CheckCarry(c);
 	}
 
 

ttmath/ttmathint.h

@@ -647,7 +647,7 @@ public:
 	return 0;
 	}
 
-	// converting from Int
+
 
 	/*!
 		the default assignment operator
@@ -714,7 +714,6 @@ public:
 	}
 
 
-	// converting from UInt
 
 	/*!
 		this operator converts an UInt<another_size> type to this class
@@ -760,7 +759,7 @@ public:
 		FromUInt(u);
 	}
  
-	// 
+
 
 
 #ifdef TTMATH_PLATFORM64
@@ -1317,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		3
+#define TTMATH_REVISION_VER		5
 #define TTMATH_PRERELEASE_VER	0
 
 
@@ -83,7 +83,7 @@
 	 gcc -DTTMATH_RELEASE -o myprogram myprogram.cpp 
 	or by defining this macro in your code before using any header files of this library
 
-	if TTMATH_RELEASE is not set then TTMATH_DEBUG is set	
+	if TTMATH_RELEASE is not set then TTMATH_DEBUG is set automatically
 */
 #ifndef TTMATH_RELEASE
 	#define TTMATH_DEBUG
@@ -120,6 +120,20 @@ namespace ttmath
 	typedef unsigned int uint;
 	typedef signed   int sint;
 
+	/*!
+		this type is twice bigger than uint
+		(64bit on a 32bit platforms)
+
+		although C++ Standard - ANSI ISO IEC 14882:2003 doesn't define such a type (long long) 
+		but it is defined in C99 and in upcoming C++0x /3.9.1 (2)/ and many compilers support it
+
+		this type is used in UInt::MulTwoWords and UInt::DivTwoWords when macro TTMATH_NOASM is defined
+		but only on a 32bit platform
+	*/
+	#ifdef TTMATH_NOASM
+		typedef unsigned long long int ulint;
+	#endif
+
 	/*!
 		how many bits there are in the uint type
 	*/
@@ -151,6 +165,15 @@ namespace ttmath
 	typedef unsigned long uint;
 	typedef signed   long sint;
 
+	/*!
+		on 64bit platform we do not define ulint
+		sizeof(long long) is 8 (64bit) but we need 128bit
+
+		on 64 bit platform (when there is defined TTMATH_NOASM macro)
+		methods UInt::MulTwoWords and UInt::DivTwoWords are using other algorithms than those on 32 bit
+	*/
+	//typedef unsigned long long int ulint;
+
 	/*!
 		how many bits there are in the uint type
 	*/
@@ -214,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
 {
@@ -326,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:
 
@@ -358,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:
 
@@ -407,6 +441,19 @@ namespace ttmath
 	#endif
 
 
+
+	#ifdef TTMATH_DEBUG_LOG
+
+		#define TTMATH_LOG(msg) \
+			PrintLog(msg, std::cout);		
+
+	#else
+
+		#define TTMATH_LOG(msg)
+
+	#endif
+
+
 } // namespace
 
 

ttmath/ttmathuint_noasm.h

@@ -0,0 +1,937 @@
+/*
+ * This file is a part of TTMath Bignum Library
+ * and is distributed under the (new) BSD licence.
+ * Author: Tomasz Sowa <t.sowa@slimaczek.pl>
+ */
+
+/* 
+ * Copyright (c) 2006-2009, Tomasz Sowa
+ * All rights reserved.
+ * 
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ * 
+ *  * Redistributions of source code must retain the above copyright notice,
+ *    this list of conditions and the following disclaimer.
+ *    
+ *  * Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *    
+ *  * Neither the name Tomasz Sowa nor the names of contributors to this
+ *    project may be used to endorse or promote products derived
+ *    from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+ * THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef headerfilettmathuint_noasm
+#define headerfilettmathuint_noasm
+
+
+#ifdef TTMATH_NOASM
+
+/*!
+	\file ttmathuint_noasm.h
+    \brief template class UInt<uint> with methods without any assembler code
+
+	this file is included at the end of ttmathuint.h
+*/
+
+
+namespace ttmath
+{
+	template<uint value_size>
+	uint UInt<value_size>::AddTwoWords(uint a, uint b, uint carry, uint * result)
+	{
+	uint temp;
+
+		if( carry == 0 )
+		{
+			temp = a + b;
+
+			if( temp < a )
+				carry = 1;
+		}
+		else
+		{
+			carry = 1;
+			temp  = a + b + carry;
+
+			if( temp > a ) // !(temp<=a)
+				carry = 0;
+		}
+
+		*result = temp;
+
+	return carry;
+	}
+
+
+
+	/*!
+		this method adding ss2 to the this and adding carry if it's defined
+		(this = this + ss2 + c)
+
+		c must be zero or one (might be a bigger value than 1)
+		function returns carry (1) (if it was)
+	*/
+	
+	template<uint value_size>
+	uint UInt<value_size>::Add(const UInt<value_size> & ss2, uint c)
+	{
+	uint i;
+
+		for(i=0 ; i<value_size ; ++i)
+			c = AddTwoWords(table[i], ss2.table[i], c, &table[i]);
+
+		TTMATH_LOG("UInt::Add")
+	
+	return c;
+	}
+
+
+	/*!
+		this method adds one word (at a specific position)
+		and returns a carry (if it was)
+
+		if we've got (value_size=3):
+			table[0] = 10;
+			table[1] = 30;
+			table[2] = 5;	
+		and we call:
+			AddInt(2,1)
+		then it'll be:
+			table[0] = 10;
+			table[1] = 30 + 2;
+			table[2] = 5;
+
+		of course if there was a carry from table[2] it would be returned
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::AddInt(uint value, uint index)
+	{
+	uint i, c;
+
+		TTMATH_ASSERT( index < value_size )
+
+
+		c = AddTwoWords(table[index], value, 0, &table[index]);
+
+		for(i=index+1 ; i<value_size && c ; ++i)
+			c = AddTwoWords(table[i], 0, c, &table[i]);
+
+		TTMATH_LOG("UInt::AddInt")
+	
+	return c;
+	}
+
+
+
+
+
+	/*!
+		this method adds only two unsigned words to the existing value
+		and these words begin on the 'index' position
+		(it's used in the multiplication algorithm 2)
+
+		index should be equal or smaller than value_size-2 (index <= value_size-2)
+		x1 - lower word, x2 - higher word
+
+		for example if we've got value_size equal 4 and:
+			table[0] = 3
+			table[1] = 4
+			table[2] = 5
+			table[3] = 6
+		then let
+			x1 = 10
+			x2 = 20
+		and
+			index = 1
+
+		the result of this method will be:
+			table[0] = 3
+			table[1] = 4 + x1 = 14
+			table[2] = 5 + x2 = 25
+			table[3] = 6
+		
+		and no carry at the end of table[3]
+
+		(of course if there was a carry in table[2](5+20) then 
+		this carry would be passed to the table[3] etc.)
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::AddTwoInts(uint x2, uint x1, uint index)
+	{
+	uint i, c;
+
+		TTMATH_ASSERT( index < value_size - 1 )
+
+
+		c = AddTwoWords(table[index],   x1, 0, &table[index]);
+		c = AddTwoWords(table[index+1], x2, c, &table[index+1]);
+
+		for(i=index+2 ; i<value_size && c ; ++i)
+			c = AddTwoWords(table[i], 0, c, &table[i]);
+
+		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)
+	{
+		if( carry == 0 )
+		{
+			*result = a - b;
+
+			if( a < b )
+				carry = 1;
+		}
+		else
+		{
+			carry   = 1;
+			*result = a - b - carry;
+
+			if( a > b ) // !(a <= b )
+				carry = 0;
+		}
+
+	return carry;
+	}
+
+
+
+
+	/*!
+		this method's subtracting ss2 from the 'this' and subtracting
+		carry if it has been defined
+		(this = this - ss2 - c)
+
+		c must be zero or one (might be a bigger value than 1)
+		function returns carry (1) (if it was)
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::Sub(const UInt<value_size> & ss2, uint c)
+	{
+	uint i;
+
+		for(i=0 ; i<value_size ; ++i)
+			c = SubTwoWords(table[i], ss2.table[i], c, &table[i]);
+
+		TTMATH_LOG("UInt::Sub")
+
+	return c;
+	}
+
+
+
+
+	/*!
+		this method subtracts one word (at a specific position)
+		and returns a carry (if it was)
+
+		if we've got (value_size=3):
+			table[0] = 10;
+			table[1] = 30;
+			table[2] = 5;	
+		and we call:
+			SubInt(2,1)
+		then it'll be:
+			table[0] = 10;
+			table[1] = 30 - 2;
+			table[2] = 5;
+
+		of course if there was a carry from table[2] it would be returned
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::SubInt(uint value, uint index)
+	{
+	uint i, c;
+
+		TTMATH_ASSERT( index < value_size )
+
+
+		c = SubTwoWords(table[index], value, 0, &table[index]);
+
+		for(i=index+1 ; i<value_size && c ; ++i)
+			c = SubTwoWords(table[i], 0, c, &table[i]);
+
+		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;
+	}
+
+
+
+	/*!
+		this method moves all bits into the left hand side
+		return value <- this <- c
+
+		the lowest *bit* will be held the 'c' and
+		the state of one additional bit (on the left hand side)
+		will be returned
+
+		for example:
+		let this is 001010000
+		after Rcl2_one(1) there'll be 010100001 and Rcl2_one returns 0
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::Rcl2_one(uint c)
+	{
+	uint i, new_c;
+
+		if( c != 0 )
+			c = 1;
+
+		for(i=0 ; i<value_size ; ++i)
+		{
+			new_c    = (table[i] & TTMATH_UINT_HIGHEST_BIT) ? 1 : 0;
+			table[i] = (table[i] << 1) | c;
+			c        = new_c;
+		}
+
+		TTMATH_LOG("UInt::Rcl2_one")
+
+	return c;
+	}
+
+
+
+
+
+
+
+	/*!
+		this method moves all bits into the right hand side
+		c -> this -> return value
+
+		the highest *bit* will be held the 'c' and
+		the state of one additional bit (on the right hand side)
+		will be returned
+
+		for example:
+		let this is 000000010
+		after Rcr2_one(1) there'll be 100000001 and Rcr2_one returns 0
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::Rcr2_one(uint c)
+	{
+	sint i; // signed i
+	uint new_c;
+
+		if( c != 0 )
+			c = TTMATH_UINT_HIGHEST_BIT;
+
+		for(i=sint(value_size)-1 ; i>=0 ; --i)
+		{
+			new_c    = (table[i] & 1) ? TTMATH_UINT_HIGHEST_BIT : 0;
+			table[i] = (table[i] >> 1) | c;
+			c        = new_c;
+		}
+
+		TTMATH_LOG("UInt::Rcr2_one")
+
+	return c;
+	}
+
+
+
+
+	/*!
+		this method moves all bits into the left hand side
+		return value <- this <- c
+
+		the lowest *bits* will be held the 'c' and
+		the state of one additional bit (on the left hand side)
+		will be returned
+
+		for example:
+		let this is 001010000
+		after Rcl2(3, 1) there'll be 010000111 and Rcl2 returns 1
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::Rcl2(uint bits, uint c)
+	{
+		TTMATH_ASSERT( bits>0 && bits<TTMATH_BITS_PER_UINT )
+
+		uint move = TTMATH_BITS_PER_UINT - bits;
+		uint i, new_c;
+
+		if( c != 0 )
+			c = TTMATH_UINT_MAX_VALUE >> move;
+
+		for(i=0 ; i<value_size ; ++i)
+		{
+			new_c    = table[i] >> move;
+			table[i] = (table[i] << bits) | c;
+			c        = new_c;
+		}
+
+		TTMATH_LOG("UInt::Rcl2")
+
+	return (c & 1);
+	}
+
+
+
+
+	/*!
+		this method moves all bits into the right hand side
+		C -> this -> return value
+
+		the highest *bits* will be held the 'c' and
+		the state of one additional bit (on the right hand side)
+		will be returned
+
+		for example:
+		let this is 000000010
+		after Rcr2(2, 1) there'll be 110000000 and Rcr2 returns 1
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::Rcr2(uint bits, uint c)
+	{
+		TTMATH_ASSERT( bits>0 && bits<TTMATH_BITS_PER_UINT )
+
+		uint move = TTMATH_BITS_PER_UINT - bits;
+		sint i; // signed
+		uint new_c;
+
+		if( c != 0 )
+			c = TTMATH_UINT_MAX_VALUE << move;
+
+		for(i=value_size-1 ; i>=0 ; --i)
+		{
+			new_c    = table[i] << move;
+			table[i] = (table[i] >> bits) | c;
+			c        = new_c;
+		}
+
+		TTMATH_LOG("UInt::Rcr2")
+
+	return (c & TTMATH_UINT_HIGHEST_BIT) ? 1 : 0;
+	}
+
+
+
+
+	/*
+		this method returns the number of the highest set bit in x
+		if the 'x' is zero this method returns '-1'
+
+	*/
+	template<uint value_size>
+	sint UInt<value_size>::FindLeadingBitInWord(uint x)
+	{
+		if( x == 0 )
+			return -1;
+
+		uint bit = TTMATH_BITS_PER_UINT - 1;
+		
+		while( (x & TTMATH_UINT_HIGHEST_BIT) == 0 )
+		{
+			x = x << 1;
+			--bit;
+		}
+
+	return bit;
+	}
+
+
+
+
+
+	/*!
+		this method sets a special bit in the 'value'
+		and returns the last state of the bit (zero or one)
+
+		bit is from <0,63>
+
+		e.g.
+		 uint x = 100;
+		 uint bit = SetBitInWord(x, 3);
+		 now: x = 108 and bit = 0
+	*/
+	template<uint value_size>
+	uint UInt<value_size>::SetBitInWord(uint & value, uint bit)
+	{
+		TTMATH_ASSERT( bit < TTMATH_BITS_PER_UINT )
+
+		uint mask = 1;
+
+		if( bit > 1 )
+			mask = mask << bit;
+
+		uint last = value & mask;
+		value     = value | mask;
+
+	return (last != 0) ? 1 : 0;
+	}
+
+
+
+
+
+
+	/*!
+	 *
+	 * Multiplication
+	 *
+	 *
+	*/
+
+
+	/*!
+		multiplication: result_high:result_low = a * b
+		result_high - higher word of the result
+		result_low  - lower word of the result
+	
+		this methos never returns a carry
+		this method is used in the second version of the multiplication algorithms
+	*/
+	template<uint value_size>
+	void UInt<value_size>::MulTwoWords(uint a, uint b, uint * result_high, uint * result_low)
+	{
+	#ifdef TTMATH_PLATFORM32
+
+		/*
+			on 32bit platforms we have defined 'unsigned long long int' type known as 'ulint' in ttmath namespace
+			this type has 64 bits, then we're using only one multiplication: 32bit * 32bit = 64bit
+		*/
+
+		union uint_
+		{
+			struct
+			{
+				uint low;  // 32 bits
+				uint high; // 32 bits
+			} u_;
+
+			ulint u;       // 64 bits
+		} res;
+
+		res.u = ulint(a) * ulint(b);     // multiply two 32bit words, the result has 64 bits
+
+		*result_high = res.u_.high;
+		*result_low  = res.u_.low;
+
+	#else
+
+		/*
+			64 bits platforms
+
+			we don't have a native type which has 128 bits
+			then we're splitting 'a' and 'b' to 4 parts (high and low halves)
+			and using 4 multiplications (with additions and carry correctness)
+		*/
+
+		uint_ a_;
+		uint_ b_;
+		uint_ res_high1, res_high2;
+		uint_ res_low1,  res_low2;
+		
+		a_.u = a;
+		b_.u = b;
+
+		/*
+			the multiplication is as follows (schoolbook algorithm with O(n^2) ):
+
+                                                   32 bits         32 bits
+
+                                             +--------------------------------+
+                                             |   a_.u_.high   |   a_.u_.low   |
+                                             +--------------------------------+
+                                             |   b_.u_.high   |   b_.u_.low   |
+            +--------------------------------+--------------------------------+
+            |           res_high1.u          |           res_low1.u           |
+            +--------------------------------+--------------------------------+
+            |           res_high2.u          |           res_low2.u           |
+            +--------------------------------+--------------------------------+
+
+                          64 bits                          64 bits
+		*/
+
+
+		uint_ temp;
+
+		res_low1.u        = uint(b_.u_.low) * uint(a_.u_.low);
+
+		temp.u            = uint(res_low1.u_.high) + uint(b_.u_.low) * uint(a_.u_.high);
+		res_low1.u_.high  = temp.u_.low;
+		res_high1.u_.low  = temp.u_.high;
+		res_high1.u_.high = 0;
+
+		res_low2.u_.low   = 0;
+		temp.u            = uint(b_.u_.high) * uint(a_.u_.low);
+		res_low2.u_.high  = temp.u_.low;
+
+		res_high2.u       = uint(b_.u_.high) * uint(a_.u_.high) + uint(temp.u_.high);
+
+		uint c = AddTwoWords(res_low1.u, res_low2.u, 0, &res_low2.u);
+		AddTwoWords(res_high1.u, res_high2.u, c, &res_high2.u);                 // there is no carry from here
+
+		*result_high = res_high2.u;
+		*result_low  = res_low2.u;
+
+	#endif
+	}
+
+
+
+
+	/*!
+	 *
+	 * Division
+	 *
+	 *
+	*/
+	
+
+	/*!
+		this method calculates 64bits word a:b / 32bits c (a higher, b lower word)
+		r = a:b / c and rest - remainder
+		
+		*
+		* WARNING:
+		* the c has to be suitably large for the result being keeped in one word,
+		* if c is equal zero there'll be a hardware interruption (0)
+		* and probably the end of your program
+		*
+	*/
+	template<uint value_size>
+	void UInt<value_size>::DivTwoWords(uint a, uint b, uint c, uint * r, uint * rest)
+	{
+	// (a < c ) for the result to be one word
+	TTMATH_ASSERT( c != 0 && a < c )
+
+	#ifdef TTMATH_PLATFORM32
+
+		union
+		{
+			struct
+			{
+				uint low;  // 32 bits
+				uint high; // 32 bits
+			} u_;
+
+			ulint u;       // 64 bits
+		} ab;
+
+		ab.u_.high = a;
+		ab.u_.low  = b;
+
+		*r    = uint(ab.u / c);
+		*rest = uint(ab.u % c);
+
+	#else
+
+		uint_ c_;
+		c_.u = c;
+
+
+		if( a == 0 )
+		{
+			*r    = b / c;
+			*rest = b % c;
+		}
+		else
+		if( c_.u_.high == 0 )
+		{
+			// higher half of 'c' is zero
+			// then higher half of 'a' is zero too (look at the asserts at the beginning - 'a' is smaller than 'c')
+			uint_ a_, b_, res_, temp1, temp2;
+
+			a_.u = a;
+			b_.u = b;
+
+			temp1.u_.high = a_.u_.low;
+			temp1.u_.low  = b_.u_.high;
+
+			res_.u_.high  = temp1.u / c;
+			temp2.u_.high = temp1.u % c;
+			temp2.u_.low  = b_.u_.low;
+			
+			res_.u_.low  = temp2.u / c;
+			*rest        = temp2.u % c;
+
+			*r = res_.u;
+		}
+		else
+		{
+			return DivTwoWords2(a, b, c,  r,  rest);
+		}
+
+	#endif
+	}
+
+
+#ifdef TTMATH_PLATFORM64
+
+
+	/*!
+		this method is available only on 64bit platforms
+		
+		the same algorithm like the third division algorithm in ttmathuint.h
+		but now with the radix=2^32
+	*/
+	template<uint value_size>
+	void UInt<value_size>::DivTwoWords2(uint a, uint b, uint c, uint * r, uint * rest)
+	{
+		// a is not zero
+		// c_.u_.high is not zero
+
+		uint_ a_, b_, c_, u_, q_;
+		unsigned int u3; // 32 bit
+
+		a_.u  = a;
+		b_.u  = b;
+		c_.u  = c;
+
+		// normalizing
+		uint d = DivTwoWordsNormalize(a_, b_, c_);
+
+		// loop from j=1 to j=0
+		//   the first step (for j=2) is skipped because our result is only in one word,
+		//   (first 'q' were 0 and nothing would be changed)
+		u_.u_.high = a_.u_.high;
+		u_.u_.low  = a_.u_.low;
+		u3         = b_.u_.high;
+		q_.u_.high = DivTwoWordsCalculate(u_, u3, c_);
+		MultiplySubtract(u_, u3, q_.u_.high, c_);
+		
+		u_.u_.high = u_.u_.low;
+		u_.u_.low  = u3;
+		u3         = b_.u_.low;
+		q_.u_.low  = DivTwoWordsCalculate(u_, u3, c_);
+		MultiplySubtract(u_, u3, q_.u_.low, c_);
+
+		*r = q_.u;
+
+		// unnormalizing for the remainder
+		u_.u_.high = u_.u_.low;
+		u_.u_.low  = u3;
+		*rest = DivTwoWordsUnnormalize(u_.u, d);
+	}
+
+
+
+	
+	template<uint value_size>
+	uint UInt<value_size>::DivTwoWordsNormalize(uint_ & a_, uint_ & b_, uint_ & c_)
+	{
+		uint d = 0;
+
+		for( ; (c_.u & TTMATH_UINT_HIGHEST_BIT) == 0 ; ++d )
+		{
+			c_.u = c_.u << 1;
+			
+			uint bc = b_.u & TTMATH_UINT_HIGHEST_BIT; // carry from 'b'
+
+			b_.u = b_.u << 1;
+			a_.u = a_.u << 1; // carry bits from 'a' are simply skipped 
+
+			if( bc )
+				a_.u = a_.u | 1;
+		}
+
+	return d;
+	}
+
+
+	template<uint value_size>
+	uint UInt<value_size>::DivTwoWordsUnnormalize(uint u, uint d)
+	{
+		if( d == 0 )
+			return u;
+
+		u = u >> d;
+
+	return u;
+	}
+
+
+	template<uint value_size>
+	unsigned int UInt<value_size>::DivTwoWordsCalculate(uint_ u_, unsigned int u3, uint_ v_)
+	{
+	bool next_test;
+	uint_ qp_, rp_, temp_;
+
+		qp_.u = u_.u / uint(v_.u_.high);
+		rp_.u = u_.u % uint(v_.u_.high);
+
+		TTMATH_ASSERT( qp_.u_.high==0 || qp_.u_.high==1 )
+
+		do
+		{
+			bool decrease = false;
+
+			if( qp_.u_.high == 1 )
+				decrease = true;
+			else
+			{
+				temp_.u_.high = rp_.u_.low;
+				temp_.u_.low  = u3;
+
+				if( qp_.u * uint(v_.u_.low) > temp_.u )
+					decrease = true;
+			}
+			
+			next_test = false;
+
+			if( decrease )
+			{
+				--qp_.u;
+				rp_.u += v_.u_.high;
+
+				if( rp_.u_.high == 0 ) 
+					next_test = true;
+			}
+		}
+		while( next_test );
+
+	return qp_.u_.low;
+	}
+
+
+	template<uint value_size>
+	void UInt<value_size>::MultiplySubtract(uint_ & u_, unsigned int & u3, unsigned int & q, uint_ v_)
+	{
+	uint_ temp_;
+		
+		uint res_high;
+		uint res_low;
+
+		MulTwoWords(v_.u, q, &res_high, &res_low);
+
+		uint_ sub_res_high_;
+		uint_ sub_res_low_;
+
+		temp_.u_.high = u_.u_.low;
+		temp_.u_.low  = u3;
+
+		uint c = SubTwoWords(temp_.u, res_low, 0, &sub_res_low_.u);
+			
+		temp_.u_.high = 0;
+		temp_.u_.low  = u_.u_.high;
+		c = SubTwoWords(temp_.u, res_high, c, &sub_res_high_.u);
+
+		if( c )
+		{
+			--q;
+
+			c = AddTwoWords(sub_res_low_.u, v_.u, 0, &sub_res_low_.u);
+			AddTwoWords(sub_res_high_.u, 0, c, &sub_res_high_.u);
+		}
+
+		u_.u_.high = sub_res_high_.u_.low;
+		u_.u_.low  = sub_res_low_.u_.high;
+		u3         = sub_res_low_.u_.low;
+	}
+
+#endif // #ifdef TTMATH_PLATFORM64
+
+
+
+} //namespace
+
+
+#endif //ifdef TTMATH_NOASM
+#endif
+
+
+
+

ttmath/ttmathuint_x86_64.h

@@ -36,10 +36,19 @@
  */
 
 
+#ifndef headerfilettmathuint_x86_64
+#define headerfilettmathuint_x86_64
+
+
+#ifndef TTMATH_NOASM
+#ifdef TTMATH_PLATFORM64
+
 
 /*!
-	\file ttmathuint.h
-    \brief template class UInt<uint> for 64bit processors
+	\file ttmathuint_x86_64.h
+    \brief template class UInt<uint> with assembler code for 64bit x86_64 processors
+
+	this file is included at the end of ttmathuint.h
 */
 
 
@@ -52,145 +61,6 @@ namespace ttmath
 	*
 	*/
 
-#ifdef TTMATH_PLATFORM64
-
-
-
-	/*!
-		in 64bit platforms we must define additional operators and contructors
-		in order to allow a user initializing the objects in this way:
-			UInt<...> type = 20;
-		or
-			UInt<...> type; 
-			type = 30;
-
-		decimal constants such as 20, 30 etc. are integer literal of type int,
-		if the value is greater it can even be long int,
-		0 is an octal integer of type int
-		(ISO 14882 p2.13.1 Integer literals)
-	*/
-
-	/*!
-		this operator converts the unsigned int type to this class
-
-		***this operator is created only on a 64bit platform***
-		it takes one argument of 32bit
-	*/
-	template<uint value_size>
-	UInt<value_size> & UInt<value_size>::operator=(unsigned int i)
-	{
-		FromUInt(uint(i));
-
-	return *this;
-	}
-
-
-	/*!
-		a constructor for converting the unsigned int to this class
-
-		***this constructor is created only on a 64bit platform***
-		it takes one argument of 32bit
-	*/
-	template<uint value_size>
-	UInt<value_size>::UInt(unsigned int i)
-	{
-		FromUInt(uint(i));
-	}
-
-
-	/*!
-		an operator for converting the signed int to this class
-
-		***this constructor is created only on a 64bit platform***
-		it takes one argument of 32bit
-
-		look at the description of UInt::operator=(sint)
-	*/
-	template<uint value_size>
-	UInt<value_size> & UInt<value_size>::operator=(signed int i)
-	{
-		FromUInt(uint(i));
-
-	return *this;
-	}
-
-
-	/*!
-		a constructor for converting the signed int to this class
-
-		***this constructor is created only on a 64bit platform***
-		it takes one argument of 32bit
-
-		look at the description of UInt::operator=(sint)
-	*/
-	template<uint value_size>
-	UInt<value_size>::UInt(signed int i)
-	{
-		FromUInt(uint(i));
-	}
-
-
-
-	/*!
-		this method copies the value stored in an another table
-		(warning: first values in temp_table are the highest words -- it's different
-		from our table)
-
-		***this method is created only on a 64bit platform***
-
-		we copy as many words as it is possible
-		
-		if temp_table_len is bigger than value_size we'll try to round 
-		the lowest word from table depending on the last not used bit in temp_table
-		(this rounding isn't a perfect rounding -- look at the description below)
-
-		and if temp_table_len is smaller than value_size we'll clear the rest words
-		in the table
-
-		warning: we're using 'temp_table' as a pointer at 32bit words
-	*/
-	template<uint value_size>
-	void UInt<value_size>::SetFromTable(const unsigned int * temp_table, uint temp_table_len)
-	{
-		uint temp_table_index = 0;
-		sint i; // 'i' with a sign
-
-		for(i=value_size-1 ; i>=0 && temp_table_index<temp_table_len; --i, ++temp_table_index)
-		{
-			table[i] = uint(temp_table[ temp_table_index ]) << 32;
-
-			++temp_table_index;
-
-			if( temp_table_index<temp_table_len )
-				table[i] |= temp_table[ temp_table_index ];
-		}
-
-
-		// rounding mantissa
-		if( temp_table_index < temp_table_len )
-		{
-			if( (temp_table[temp_table_index] & TTMATH_UINT_HIGHEST_BIT) != 0 )
-			{
-				/*
-					very simply rounding
-					if the bit from not used last word from temp_table is set to one
-					we're rouding the lowest word in the table
-
-					in fact there should be a normal addition but
-					we don't use Add() or AddTwoInts() because these methods 
-					can set a carry and then there'll be a small problem
-					for optimization
-				*/
-				if( table[0] != TTMATH_UINT_MAX_VALUE )
-					++table[0];
-			}
-		}
-
-		// cleaning the rest of the mantissa
-		for( ; i >= 0 ; --i)
-			table[i] = 0;
-	}
-
 
 
 	/*!
@@ -205,10 +75,9 @@ namespace ttmath
 	template<uint value_size>
 	uint UInt<value_size>::Add(const UInt<value_size> & ss2, uint c)
 	{
-	register uint b = value_size;
-	register uint * p1 = table;
-	register uint * p2 = const_cast<uint*>(ss2.table);
-
+	uint b = value_size;
+	uint * p1 = table;
+	const uint * p2 = ss2.table;
 
 		// we don't have to use TTMATH_REFERENCE_ASSERT here
 		// this algorithm doesn't require it
@@ -218,17 +87,15 @@ namespace ttmath
 		#endif
 
 		#ifdef __GNUC__
+		uint dummy, dummy2;
+
 			/*
 				this part should be compiled with gcc
 			*/
 			__asm__ __volatile__(
 	
-				"push %%rcx						\n"
-			
-				"xorq %%rax, %%rax				\n"
-				"movq %%rax, %%rdx				\n"
-				"subq %%rdi, %%rax				\n"
-
+				"xorq %%rdx, %%rdx				\n"
+				"negq %%rax						\n"     // CF=1 if rax!=0 , CF=0 if rax==0
 
 			"1:									\n"
 				"movq (%%rsi,%%rdx,8), %%rax	\n"
@@ -238,18 +105,15 @@ namespace ttmath
 				"decq %%rcx						\n"
 			"jnz 1b								\n"
 
-				"setc %%al						\n"
-				"movzx %%al,%%rdx				\n"
+				"adcq %%rcx, %%rcx				\n"
 
-				"pop %%rcx						\n"
-
-				: "=d" (c)
-				: "D" (c), "c" (b), "b" (p1), "S" (p2)
-				: "%rax", "cc", "memory" );
+				: "=c" (c), "=a" (dummy), "=d" (dummy2)
+				: "0" (b),  "1" (c), "b" (p1), "S" (p2)
+				: "cc", "memory" );
 
 		#endif
 
-
+		TTMATH_LOG("UInt::Add")
 	
 	return c;
 	}
@@ -274,14 +138,14 @@ namespace ttmath
 			table[1] = 30 + 2;
 			table[2] = 5;
 
-		of course if there was a carry from table[3] it would be returned
+		of course if there was a carry from table[2] it would be returned
 	*/
 	template<uint value_size>
 	uint UInt<value_size>::AddInt(uint value, uint index)
 	{
-	register uint b = value_size;
-	register uint * p1 = table;
-	register uint c;
+	uint b = value_size;
+	uint * p1 = table;
+	uint c;
 
 		TTMATH_ASSERT( index < value_size )
 
@@ -290,12 +154,10 @@ namespace ttmath
 		#endif
 
 		#ifdef __GNUC__
+		uint dummy, dummy2;
 	
 			__asm__ __volatile__(
 
-				"push %%rax						\n"
-				"push %%rcx						\n"
-
 				"subq %%rdx, %%rcx 				\n"
 
 			"1:									\n"
@@ -311,15 +173,13 @@ namespace ttmath
 				"setc %%al						\n"
 				"movzx %%al, %%rdx				\n"
 
-				"pop %%rcx						\n"
-				"pop %%rax						\n"
-
-				: "=d" (c)
-				: "a" (value), "c" (b), "0" (index), "b" (p1)
+				: "=d" (c),    "=a" (dummy), "=c" (dummy2)
+				: "0" (index), "1" (value),  "2" (b), "b" (p1)
 				: "cc", "memory" );
 
 		#endif
 
+		TTMATH_LOG("UInt::AddInt")
 	
 	return c;
 	}
@@ -361,9 +221,9 @@ namespace ttmath
 	template<uint value_size>
 	uint UInt<value_size>::AddTwoInts(uint x2, uint x1, uint index)
 	{
-	register uint b = value_size;
-	register uint * p1 = table;
-	register uint c;
+	uint b = value_size;
+	uint * p1 = table;
+	uint c;
 
 		TTMATH_ASSERT( index < value_size - 1 )
 
@@ -372,10 +232,9 @@ namespace ttmath
 		#endif
 
 		#ifdef __GNUC__
-			__asm__ __volatile__(
+		uint dummy, dummy2;
 
-				"push %%rcx						\n"
-				"push %%rdx						\n"
+			__asm__ __volatile__(
 			
 				"subq %%rdx, %%rcx 				\n"
 				
@@ -396,15 +255,95 @@ namespace ttmath
 				"setc %%al						\n"
 				"movzx %%al, %%rax				\n"
 
-				"pop %%rdx						\n"
-				"pop %%rcx						\n"
-
-				: "=a" (c)
-				: "c" (b), "d" (index), "b" (p1), "S" (x1), "0" (x2)
+				: "=a" (c), "=c" (dummy), "=d" (dummy2)
+				: "0" (x2), "1" (b),      "2" (index), "b" (p1), "S" (x1)
 				: "cc", "memory" );
 
 		#endif
 
+		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;
 	}
@@ -426,9 +365,10 @@ namespace ttmath
 	template<uint value_size>
 	uint UInt<value_size>::Sub(const UInt<value_size> & ss2, uint c)
 	{
-	register uint b = value_size;
-	register uint * p1 = table;
-	register uint * p2 = const_cast<uint*>(ss2.table);
+	uint b = value_size;
+	uint * p1 = table;
+	const uint * p2 = ss2.table;
+	
 
 		// we don't have to use TTMATH_REFERENCE_ASSERT here
 		// this algorithm doesn't require it
@@ -438,14 +378,12 @@ namespace ttmath
 		#endif
 
 		#ifdef __GNUC__
+		uint dummy, dummy2;
+
 			__asm__  __volatile__(
 	
-				"push %%rcx						\n"
-			
-				"xorq %%rax, %%rax				\n"
-				"movq %%rax, %%rdx				\n"
-				"subq %%rdi, %%rax				\n"
-
+				"xorq %%rdx, %%rdx				\n"
+				"negq %%rax						\n"     // CF=1 if rax!=0 , CF=0 if rax==0
 
 			"1:									\n"
 				"movq (%%rsi,%%rdx,8), %%rax	\n"
@@ -455,18 +393,15 @@ namespace ttmath
 				"decq %%rcx						\n"
 			"jnz 1b								\n"
 
-				"setc %%al						\n"
-				"movzx %%al,%%rdx				\n"
-
-				"pop %%rcx						\n"
-
-				: "=d" (c)
-				: "D" (c), "c" (b), "b" (p1), "S" (p2)
-				: "%rax", "cc", "memory" );
+				"adcq %%rcx, %%rcx				\n"
 
+				: "=c" (c), "=a" (dummy), "=d" (dummy2)
+				: "0" (b),  "1" (c), "b" (p1), "S" (p2)
+				: "cc", "memory" );
 
 		#endif
 
+		TTMATH_LOG("UInt::Sub")
 
 	return c;
 	}
@@ -489,14 +424,15 @@ namespace ttmath
 			table[1] = 30 - 2;
 			table[2] = 5;
 
-		of course if there was a carry from table[3] it would be returned
+		of course if there was a carry from table[2] it would be returned
 	*/
 	template<uint value_size>
 	uint UInt<value_size>::SubInt(uint value, uint index)
 	{
-	register uint b = value_size;
-	register uint * p1 = table;
-	register uint c;
+	uint b = value_size;
+	uint * p1 = table;
+	uint c;
+	uint dummy, dummy2;
 
 		TTMATH_ASSERT( index < value_size )
 
@@ -507,9 +443,6 @@ namespace ttmath
 		#ifdef __GNUC__
 			__asm__ __volatile__(
 			
-				"push %%rax						\n"
-				"push %%rcx						\n"
-
 				"subq %%rdx, %%rcx 				\n"
 
 			"1:									\n"
@@ -525,20 +458,106 @@ namespace ttmath
 				"setc %%al						\n"
 				"movzx %%al, %%rdx				\n"
 
-				"pop %%rcx						\n"
-				"pop %%rax						\n"
-
-				: "=d" (c)
-				: "a" (value), "c" (b), "0" (index), "b" (p1)
+				: "=d" (c),    "=a" (dummy), "=c" (dummy2)
+				: "0" (index), "1" (value),  "2" (b), "b" (p1)
 				: "cc", "memory" );
 
 		#endif
 
+		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
@@ -556,8 +575,8 @@ namespace ttmath
 	template<uint value_size>
 	uint UInt<value_size>::Rcl2_one(uint c)
 	{
-	register sint b = value_size;
-	register uint * p1 = table;
+	sint b = value_size;
+	uint * p1 = table;
 	
 
 		#ifndef __GNUC__
@@ -565,13 +584,12 @@ namespace ttmath
 		#endif
 
 		#ifdef __GNUC__
+		uint dummy, dummy2;
+
 		__asm__  __volatile__(
 		
-			"push %%rdx					\n"
-			"push %%rcx					\n"
-
 			"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"
@@ -580,18 +598,15 @@ namespace ttmath
 			"decq %%rcx					\n"
 		"jnz 1b							\n"
 
-			"setc %%al					\n"
-			"movzx %%al, %%rax			\n"
+			"adcq %%rcx, %%rcx			\n"
 
-			"pop %%rcx					\n"
-			"pop %%rdx					\n"
-
-			: "=a" (c)
-			: "0" (c), "c" (b), "b" (p1)
+			: "=c" (c), "=a" (dummy), "=d" (dummy2)
+			: "0" (b),  "1" (c), "b" (p1)
 			: "cc", "memory" );
 	
 		#endif
 
+		TTMATH_LOG("UInt::Rcl2_one")
 
 	return c;
 	}
@@ -614,8 +629,8 @@ namespace ttmath
 	template<uint value_size>
 	uint UInt<value_size>::Rcr2_one(uint c)
 	{
-	register sint b = value_size;
-	register uint * p1 = table;
+	sint b = value_size;
+	uint * p1 = table;
 	
 
 		#ifndef __GNUC__
@@ -623,11 +638,11 @@ namespace ttmath
 		#endif
 
 		#ifdef __GNUC__
+		uint dummy;
+
 		__asm__  __volatile__(
 
-			"push %%rcx						\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"
 			"rcrq $1, -8(%%rbx, %%rcx, 8)	\n"
@@ -635,17 +650,15 @@ namespace ttmath
 			"decq %%rcx						\n"
 		"jnz 1b								\n"
 
-			"setc %%al						\n"
-			"movzx %%al, %%rax				\n"
+			"adcq %%rcx, %%rcx				\n"
 
-			"pop %%rcx						\n"
-
-			: "=a" (c)
-			: "0" (c), "c" (b), "b" (p1)
+			: "=c" (c), "=a" (dummy)
+			: "0" (b),  "1" (c), "b" (p1)
 			: "cc", "memory" );
 
 		#endif
 
+		TTMATH_LOG("UInt::Rcr2_one")
 
 	return c;
 	}
@@ -671,40 +684,37 @@ namespace ttmath
 	{
 	TTMATH_ASSERT( bits>0 && bits<TTMATH_BITS_PER_UINT )
 
-	register sint b = value_size;
-	register uint * p1 = table;
-	register uint mask;
+	uint b = value_size;
+	uint * p1 = table;
+
 
 		#ifndef __GNUC__
 			#error "another compiler than GCC is currently not supported in 64bit mode"
 		#endif
 
 		#ifdef __GNUC__
-		__asm__  __volatile__(
+		uint dummy, dummy2, dummy3;
 
-			"push %%rdx						\n"
-			"push %%rsi						\n"
-			"push %%rdi						\n"
+		__asm__  __volatile__(
 		
 			"movq %%rcx, %%rsi				\n"
 			"movq $64, %%rcx				\n"
 			"subq %%rsi, %%rcx				\n"
 			"movq $-1, %%rdx				\n"
 			"shrq %%cl, %%rdx				\n"
-			"movq %%rdx, %[amask] 			\n"
+			"movq %%rdx, %%r8 				\n"
 			"movq %%rsi, %%rcx				\n"
 
 			"xorq %%rdx, %%rdx				\n"
 			"movq %%rdx, %%rsi				\n"
-
 			"orq %%rax, %%rax				\n"
-			"cmovnz %[amask], %%rsi			\n"
+			"cmovnz %%r8, %%rsi				\n"
 
 		"1:									\n"
 			"rolq %%cl, (%%rbx,%%rdx,8)		\n"
 
 			"movq (%%rbx,%%rdx,8), %%rax	\n"
-			"andq %[amask], %%rax			\n"
+			"andq %%r8, %%rax				\n"
 			"xorq %%rax, (%%rbx,%%rdx,8)	\n"
 			"orq  %%rsi, (%%rbx,%%rdx,8)	\n"
 			"movq %%rax, %%rsi				\n"
@@ -715,16 +725,13 @@ namespace ttmath
 			
 			"and $1, %%rax					\n"
 
-			"pop %%rdi						\n"
-			"pop %%rsi						\n"
-			"pop %%rdx						\n"
-
-			: "=a" (c)
-			: "0" (c), "D" (b), "b" (p1), "c" (bits), [amask] "m" (mask)
-			: "cc", "memory" );
+			: "=a" (c), "=D" (dummy), "=S" (dummy2), "=d" (dummy3)
+			: "0" (c),  "1" (b), "b" (p1), "c" (bits)
+			: "%r8", "cc", "memory" );
 
 		#endif
 
+		TTMATH_LOG("UInt::Rcl2")
 
 	return c;
 	}
@@ -749,9 +756,9 @@ namespace ttmath
 	{
 	TTMATH_ASSERT( bits>0 && bits<TTMATH_BITS_PER_UINT )
 
-	register sint b = value_size;
-	register uint * p1 = table;
-	register uint mask;
+	sint b = value_size;
+	uint * p1 = table;
+	uint dummy, dummy2, dummy3;
 
 		#ifndef __GNUC__
 			#error "another compiler than GCC is currently not supported in 64bit mode"
@@ -761,31 +768,26 @@ namespace ttmath
 		#ifdef __GNUC__
 			__asm__  __volatile__(
 
-			"push %%rdx						\n"
-			"push %%rsi						\n"
-			"push %%rdi						\n"
-			
 			"movq %%rcx, %%rsi				\n"
 			"movq $64, %%rcx				\n"
 			"subq %%rsi, %%rcx				\n"
 			"movq $-1, %%rdx				\n"
 			"shlq %%cl, %%rdx				\n"
-			"movq %%rdx, %[amask]			\n"
+			"movq %%rdx, %%R8				\n"
 			"movq %%rsi, %%rcx				\n"
 
 			"xorq %%rdx, %%rdx				\n"
 			"movq %%rdx, %%rsi				\n"
 			"addq %%rdi, %%rdx				\n"
 			"decq %%rdx						\n"
-
 			"orq %%rax, %%rax				\n"
-			"cmovnz %[amask], %%rsi			\n"
+			"cmovnz %%R8, %%rsi				\n"
 
 		"1:									\n"
 			"rorq %%cl, (%%rbx,%%rdx,8)		\n"
 
 			"movq (%%rbx,%%rdx,8), %%rax	\n"
-			"andq %[amask], %%rax			\n"
+			"andq %%R8, %%rax				\n"
 			"xorq %%rax, (%%rbx,%%rdx,8)	\n"
 			"orq  %%rsi, (%%rbx,%%rdx,8)	\n"
 			"movq %%rax, %%rsi				\n"
@@ -797,23 +799,20 @@ namespace ttmath
 			"rolq $1, %%rax					\n"
 			"andq $1, %%rax					\n"
 
-			"pop %%rdi						\n"
-			"pop %%rsi						\n"
-			"pop %%rdx						\n"
-
-			: "=a" (c)
-			: "0" (c), "D" (b), "b" (p1), "c" (bits), [amask] "m" (mask)
-			: "cc", "memory" );
+			: "=a" (c), "=D" (dummy), "=S" (dummy2), "=d" (dummy3)
+			: "0" (c), "1" (b), "b" (p1), "c" (bits)
+			: "%r8", "cc", "memory" );
 
 		#endif
 
+		TTMATH_LOG("UInt::Rcr2")
 
 	return c;
 	}
 
 
 	/*
-		this method returns the number of the highest set bit in one 32-bit word
+		this method returns the number of the highest set bit in one 64-bit word
 		if the 'x' is zero this method returns '-1'
 
 		***this method is created only on a 64bit platform***
@@ -821,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
@@ -873,10 +874,10 @@ namespace ttmath
 		#endif
 
 		#ifdef __GNUC__
-			__asm__  __volatile__(
+
+			__asm__ (
 
 			"btsq %%rbx, %%rax		\n"
-
 			"setc %%bl				\n"
 			"movzx %%bl, %%rbx		\n"
 			
@@ -901,18 +902,17 @@ namespace ttmath
 
 
 	/*!
-		multiplication: result2:result1 = a * b
-		result2 - higher word
-		result1 - lower word of the result
+		multiplication: result_high:result_low = a * b
+		result_high - higher word of the result
+		result_low  - lower word of the result
 	
 		this methos never returns a carry
+		this method is used in the second version of the multiplication algorithms
 
 		***this method is created only on a 64bit platform***
-
-		it is an auxiliary method for version two of the multiplication algorithm
 	*/
 	template<uint value_size>
-	void UInt<value_size>::MulTwoWords(uint a, uint b, uint * result2, uint * result1)
+	void UInt<value_size>::MulTwoWords(uint a, uint b, uint * result_high, uint * result_low)
 	{
 	/*
 		we must use these temporary variables in order to inform the compilator
@@ -921,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"
@@ -930,7 +930,7 @@ namespace ttmath
 
 		#ifdef __GNUC__
 
-		__asm__ __volatile__(
+		__asm__ (
 		
 			"mulq %%rdx			\n"
 
@@ -941,8 +941,8 @@ namespace ttmath
 		#endif
 
 
-		*result1 = result1_;
-		*result2 = result2_;
+		*result_low  = result1_;
+		*result_high = result2_;
 	}
 
 
@@ -972,20 +972,22 @@ 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
 		*/
 
+		TTMATH_ASSERT( c != 0 )
+
 		#ifndef __GNUC__
 			#error "another compiler than GCC is currently not supported in 64bit mode"
 		#endif
 
 		#ifdef __GNUC__
 		
-			__asm__ __volatile__(
+			__asm__ (
 
 			"divq %%rcx				\n"
 
@@ -1000,6 +1002,12 @@ namespace ttmath
 		*rest = rest_;
 	}
 
-#endif
 
 } //namespace
+
+
+#endif //ifdef TTMATH_PLATFORM64
+#endif //ifndef TTMATH_NOASM
+#endif
+
+